METHODS AND SYSTEMS FOR GENOME-SCALE KINETIC MODELING
Embodiments of the present invention generally relate to the construction, analysis, and characterization of dynamical states of biological networks at the cellular level. Methods are provided for analyzing the dynamical states by constructing matrices using high-throughput data types, such as fluxomic, metabolomic, and proteomic data. Some embodiments relate to an individual, while others relate to a plurality of individuals.
Latest The Regents of the University of California Patents:
- METHOD TO IMPROVE PERFORMANCES OF TUNNEL JUNCTIONS GROWN BY METAL ORGANIC CHEMICAL VAPOR DEPOSITION
- MULTI-PHASE HYBRID POWER CONVERTER ARCHITECTURE WITH LARGE CONVERSION RATIOS
- TARGETED GENE DEMETHYLATION IN PLANTS
- Integrated Programmable Strongly Coupled Three-Ring Resonator Photonic Molecule with Ultralow-Power Piezoelectric Control
- PROGRAMMABLE SMART-SENSING SHELF LINER APPARATUS FOR AUTONOMOUS RETAIL
The present application is a continuation of U.S. application Ser. No. 12/918,317, filed on Aug. 18, 2010, which is a national stage application under 35 U.S.C. §371 of PCT Application No. PCT/US2009/034587, filed on Feb. 19, 2009, which published in English as WO 2009/105591 A3 on Aug. 27, 2009, and which claims benefit of U.S. Provisional Application No. 61/029,891, entitled Methods and Systems for Genome-scale Kinetic Modeling, filed Feb. 19, 2008, the entire contents of which applications and publication are herein incorporated by reference in their entirety.
STATEMENT OF GOVERNMENT SUPPORTThis invention was made with government support under grant number BNGBP47 awarded by NIH. The government has certain rights in the invention.
BACKGROUND OF THE INVENTION1. Field of the Invention
Embodiments of the present invention generally relate to the construction, analysis, and characterization of dynamical states of biological networks at the cellular level. Methods are provided for analyzing the dynamical states by constructing matrices using high-throughput data types, such as fluxomic, metabolomic, and proteomic data. Some embodiments relate to an individual, while others relate to a plurality of individuals.
2. Description of the Related Art
There have been significant advances in the development of statistical analysis of genome-scale networks (Slonim, 2002), that have been propelled by the availability of genome-scale high-throughput datasets and the successes of constraint-based modeling approaches (Kummel et al., 2006; Palsson, 2006; Pharkya et al., 2004; Price et al., 2004; Reed et al., 2006). The foundation of such genome-scale analysis is built upon the descriptions of the biochemical transformations in a network in a self-consistent and chemically accurate matrix format. Much progress has been made with the genome-scale network reconstruction process, and a growing number of genome-scale metabolic reconstructions are now available (Feist et al., 2007; Jamshidi and Palsson, 2007; Oh et al., 2007; Reed et al., 2006; Resendis-Antonio et al., 2007).
Reconstructions of genome-scale biochemical reaction networks (Reed et al., 2006) have been analyzed by topological (Barabasi and Oltvai, 2004) and constraints-based (Price et al., 2004) methods, but dynamic models, at this scale, have been unachievable to date. It turns out that the matrix representation of the biochemical transformations in a network are not only a requisite for dynamic models but also a major determinant in their properties, and thus, it is important to have well curated reconstructions available. The growing availability of metabolomic and fluxomic data sets (Breitling et al., 2007; Goodacre et al., 2004; Hollywood et al., 2006; Sauer, 2006) and methods to estimate the thermodynamic properties (Henry et al., 2007; Henry et al., 2006; Mavrovouniotis, 1991) of biochemical reactions has opened up the possibility to formulate large-scale kinetic models.
The growing amount of disparate data sets (including proteomic, fluxomic, and metabolomic) has produced the unmet need to integrate and analyze this data in a functional context that is tailored to specific individuals. The availability of annotated genomes enabled the reconstruction of genome-scale networks, and now the availability of high-throughput metabolomic and fluxomic data along with thermodynamic information opens the possibility to build genome-scale kinetic models. This invention describes a framework for building and analyzing such models, thus satisfying those needs.
SUMMARY OF THE INVENTIONThe structure of the workflow that leads to large-scale dynamic models is described in this invention along with the description of the framework which enables the integration of data to construct genome-scale kinetic models.
A description of the types of data matrices which are needed to construct large-scale kinetic models of biological processes is provided.
A description of the way to construct the data matrices needed for the formulation of kinetic models is provided. A description of important properties of these data matrices which enable approximations of the values in the data matrices when partial or incomplete data is available is also provided.
A method for the analysis of kinetic models and for providing a method for interpreting and predicting physiological and pathophysiological states of biological reaction networks is provided.
A method is provided for developing data-driven dynamic models of biological networks, the method comprising providing a data structure relating a plurality of biological network reactants to a plurality of biological network reactions, wherein each of the biological network reactions comprises one or more biological network reactants identified as a substrate of the reaction, one or more biological network reactants identified as a product of the biological network reaction and a stoichiometric coefficient relating the substrate and the product; providing a constraint set for a plurality of biological network reactions; providing a constraint set for the plurality of concentrations associated with the plurality of biological network reactions; providing a constraint set for a plurality of kinetic constants associated with the plurality of biological network reactions; providing a constraint set for a plurality of equilibrium constants associated with the plurality of biological network reactions; integrating the sets of constraints, thereby forming the stoichiometric, gradient, and Jacobian matrices; providing an objective function; and predicting a physiological function related to a biological network by determining at least one flux distribution that minimizes or maximizes the objective function when the constraint set is applied to the data structure.
A method is provided wherein the constraints for the plurality of reactions, the constraints for the plurality of concentrations, the constraints for the plurality of kinetic constants, and the constraints for the plurality of equilibrium constants are uniquely defined.
A method is provided wherein providing a constraint set for a plurality of concentrations comprises measuring metabolite levels at steady state.
A method is provided, wherein providing a constraint set for a plurality of kinetic constants comprises solving a steady state mass conservation relationship comprising a stoichiometric matrix and flux vector.
A method is provided further comprising defining a unique set of matrices using the objective function, wherein the constraints within at least one of the sets selected from the reactions constraint set, the concentration constraint set, the kinetic constant constraint set, and the equilibrium constraint set are not unique across the set.
A method is provided further comprising defining a unique set of gradient and Jacobian matrices using order of magnitude approximations, wherein the constraints within at least one of the sets selected from the reactions constraint set, the concentration constraint set, the kinetic constant constraint set, and the equilibrium constraint set are not unique across the set.
A method is provided further comprising determining a set of gradient and Jacobian matrices corresponding to the plurality of constraint conditions, wherein the constraints within at least one of the sets selected from the reactions constraint set, the concentration constraint set, the kinetic constant constraint set, and the equilibrium constraint set are not unique across the set.
A method is provided wherein constraints within at least one of the concentration constraint set and the kinetic constant constraint set are not unique across the set and are estimated by order of magnitude approximations.
A method is provided wherein the stoichiometric, gradient, and Jacobian matrices refer to a biological network process.
A method is provided wherein the biological network process comprises at least one of metabolism, signal transduction or signaling, transcription, and translation.
A method is provided wherein the sets of constraints correspond to an individual.
A method is provided wherein the sets of constraints correspond to a group or plurality of individuals.
A method is provided wherein the model is used to determine the effects of changes from aerobic to anaerobic conditions.
A method is provided wherein the biological network reactions are obtained from a metabolic reaction database that includes the substrates, products, and stoichiometry of a plurality of biological reactions.
A method is provided for the analysis of disease states due to genetic and or environmental influences, the method comprising providing a data structure relating a plurality of reactants to a plurality of reactions, wherein each of the reactions comprises one or more reactants identified as a substrate of the reaction, one or more reactants identified as a product of the reaction and a stoichiometric coefficient relating the substrate and the product; providing a constraint set for the plurality of reactions; providing a constraint set for a plurality of concentrations associated with the plurality of reactions; providing a constraint set for a plurality of kinetic constants associated with the plurality of reactions; providing a constraint set for a plurality of equilibrium constants associated with the plurality of reactions; integrating the sets of constraints, thereby forming the stoichiometric, gradient, and Jacobian matrices; providing an objective function; modifying at least one of the reaction constraint set, the concentration constraint set, the constraint set for a plurality of kinetic constants and the set for a plurality of equilibrium constants based at least partly on altered environmental conditions; and predicting a physiological function related to the disease state by determining at least one flux distribution that minimizes or maximizes the objective function when at least one constraint set is applied to the data structure.
A method is provided wherein the model is used to determine the dynamics and consequences of genetic defects selected from the group consisting of deficiencies in metabolic enzymes and metabolic transporters.
A method is provided wherein the genetic defects comprise the functioning of a compound selected from the group consisting of phosphofructokinase, phosphoglycerate kinase, phosphoglycerate mutase, lactate dehydrogenase adenosine deaminase, ABC transporters, the SLC class of transporters, and the cytochrome P450 class of enzymes.
A method is provided wherein providing a constraint set for a plurality of concentrations comprises measuring metabolite levels at steady state.
A method is provided wherein providing a constraint set for a plurality of kinetic constants comprises solving a steady state mass conservation relationship comprising a stoichiometric matrix and flux vector.
A method is provided wherein the modified sets of constraints correspond to an individual.
A method is provided wherein the modified sets of constraints correspond to a group or plurality of individuals.
A method is provided wherein a therapeutic regime is applied to regulate a physiological function based on at least one of the modified sets based at least partly on genetic variations.
A method is provided wherein a therapeutic regime is applied to regulate a physiological function based on at least one of the modified sets based at least partly on altered environmental conditions.
A method is provided for the analysis of therapies for drug treatments, the method comprising providing a data structure relating a plurality of reactants to a plurality of reactions, wherein each of the reactions comprises one or more reactants identified as a substrate of the reaction, one or more reactants identified as a product of the reaction and a stoichiometric coefficient relating the substrate and the product; providing a constraint set for the plurality of reactions; providing a constraint set for a plurality of concentrations associated with the plurality of reactions; providing a constraint set for a plurality of kinetic constants associated with the plurality of reactions; providing a constraint set for a plurality of equilibrium constants associated with the plurality of reactions; integrating the sets of constraints, thereby forming the stoichiometric, gradient, and Jacobian matrices; providing an objective function; providing a experimental constraint set for the plurality of reactions, concentrations, or kinetic constants for an altered set of constraints due to the treatment with a drug or chemical agent; and predicting a physiological function related to a drug treatment by determining at least one flux distribution that minimizes or maximizes the objective function when a constraint set is applied to the data structure.
A method is provided wherein a set of fluxes and concentrations are measured in an individual before and after giving a dose of a therapeutic agent to the individual and constructing kinetic network models according to the measured data specific to the individual.
A method is provided wherein providing a constraint set for a plurality of concentrations comprises measuring metabolite levels at steady state.
A method is provided wherein providing a constraint set for a plurality of kinetic constants comprises solving a steady state mass conservation relationship comprising a stoichiometric matrix and flux vector.
A method is provided wherein a physiological function related to a drug treatment is instead predicted by determining at least one concentration that maximizes or minimizes the objective function when the experimental constraint set is applied to the data structure.
A method is provided wherein a drug or chemical agent is applied to regulate a physiological function based on the at least one experimental constraint set for an altered state of constraints due to the treatment with a drug or chemical agent.
A method is provided wherein a drug or chemical agent is applied to regulate a physiological function based on the at least one experimental constraint set which maximizes or minimizes the objective function when the experimental constraint set is applied to the data structure.
A method is provided wherein the sets of constraints correspond to an individual.
A method is provided wherein the sets of constraints correspond to a group or plurality of individuals.
A method is provided for accounting for small metabolite regulation, the method comprising: providing a data structure relating a plurality of reactants to a plurality of reactions, wherein each of the reactions comprises one or more reactants identified as a substrate of the reaction, one or more reactants identified as a product of the reaction and a stoichiometric coefficient relating the substrate and the product, in which the reactions may reflect regulatory processes; providing a constraint set for the plurality of reactions; providing a constraint set for a plurality of concentrations associated with the plurality of reactions; providing a constraint set for a plurality of kinetic constants associated with the plurality of reactions; providing a constraint set for a plurality of equilibrium constants associated with the plurality of reactions; integrating the sets of constraints, thereby forming the stoichiometric, gradient, and Jacobian matrices; providing an objective function; and predicting a physiological function related to the small metabolite regulation by determining at least one flux distribution that minimizes or maximizes the objective function when the constraint set is applied to the data structure.
A method is provided wherein providing a constraint set for a plurality of concentrations comprises measuring metabolite levels at steady state.
A method is provided wherein a physiological function related to the small metabolite regulation is instead predicted by determining at least one concentration that minimizes or maximizes the objective function when the constraint set is applied to the data structure.
A method is provided wherein reactants and products include proteins that interact with a plurality of other components.
A method is provided wherein reactants or products activate a particular enzyme or plurality of enzymes.
A method is provided wherein reactants or products inhibit a particular enzyme or plurality of enzymes.
A method is provided for determining one or more parameters to characterize different functional states, the method comprising: analyzing a network function; determining time scales of interest; and identifying parameters that describe the network function.
A method is provided wherein the parameters comprise at least one of a reaction, a concentration and a rate constant.
A method is provided further comprising characterizing the dynamic behavior of the network using the time scale determination; and identifying reactions or concentrations that determine the network function.
A computer readable medium is provided, the medium containing software that, when executed, causes the computer to perform the acts of: providing a data structure relating a plurality of reactants to a plurality of reactions, wherein each of the reactions comprises one or more reactants identified as a substrate of the reaction, one or more reactants identified as a product of the reaction and a stoichiometric coefficient relating the substrate and the product; providing a constraint set for the plurality of reactions; providing a constraint set for a plurality of concentrations associated with the plurality of reactions; providing a constraint set for a plurality of kinetic constants associated with the plurality of reactions; providing a constraint set for a plurality of equilibrium constants associated with the plurality of reactions; integrating the sets of constraints, thereby forming the stoichiometric, gradient, and Jacobian matrices; providing an objective function; and predicting a physiological function by determining at least one flux distribution that minimizes or maximizes the objective function when the constraint set is applied to the data structure.
A computer implemented system is provided comprising: a computing environment; a storage in data communication with the computing environment and configured to store data structure relating a plurality of reactants to a plurality of reactions; and a software program operating on the computing environment and configured to: provide the data structure relating the plurality of reactants to the plurality of reactions, wherein each of the reactions comprises one or more reactants identified as a substrate of the reaction, one or more reactants identified as a product of the reaction and a stoichiometric coefficient relating the substrate and the product; provide a constraint set for the plurality of reactions; provide a constraint set for a plurality of concentrations associated with the plurality of reactions; provide a constraint set for a plurality of kinetic constants associated with the plurality of reactions; provide a constraint set for a plurality of equilibrium constants associated with the plurality of reactions; integrate the sets of constraints, thereby forming the stoichiometric, gradient, and Jacobian matrices; provide an objective function; and predict a physiological function by determining at least one flux distribution that minimizes or maximizes the objective function when the constraint set is applied to the data structure.
A system is provided comprising means for providing a data structure relating a plurality of reactants to a plurality of reactions, wherein each of the reactions comprises one or more reactants identified as a substrate of the reaction, one or more reactants identified as a product of the reaction and a stoichiometric coefficient relating the substrate and the product; means for providing a constraint set for the plurality of reactions; means for providing a constraint set for a plurality of concentrations associated with the plurality of reactions; means for providing a constraint set for a plurality of kinetic constants associated with the plurality of reactions; means for providing a constraint set for a plurality of equilibrium constants associated with the plurality of reactions; means for integrating the sets of constraints, thereby forming the stoichiometric, gradient, and Jacobian matrices; means for providing an objective function; and means for predicting a physiological function by determining at least one flux distribution that minimizes or maximizes the objective function when the constraint set is applied to the data structure.
Methods for constructing and analyzing dynamic in silico models of biological networks are provided herein. A dynamic model can be used to simulate different aspects of cellular behavior of human cells under different environmental and genetic conditions, thereby providing valuable information for a range of industrial and research applications. Models and methods disclosed herein can also be extended to simulate the activity of multiple interacting cells, including organs, physiological systems, and whole body metabolism. Developing kinetic models of biological networks can be used to inform and guide the research process, potentially leading to the identification of new enzymes, medicines or metabolites of commercial importance.
As used herein, the term “concentration” refers to a numerical value with physical units of mass*volume−1, such as molar and millimolar. Such quantities include but are not limited to molecules in a biochemical network, such as glucose, pyruvate, lactate, enzymes that carry out biochemical transformations or transport reactions such as hexokinase, adenosine deaminase, and sodium-potassium ATPase pump, or the number of cells in a volume of tissue.
As used herein, the term “initial conditions” refers to the set of initial variable values that need to be specified for carrying out dynamic simulations in a system of dynamic differential equations.
As used herein, the term “dynamic mass balance” refers to a linear first order differential equation with respect to time which accounts for the summed reactions producing or consuming each component in a specified biological network.
As used herein the term “flux” or “reaction flux” refers to a single flux in a flux distribution. Individual fluxes can be represented as the difference between a forward and reverse flux, with units mass, or number of moles, or number of molecules per unit volume per unit time.
As used herein, the term “flux distribution” refers to a directional, quantitative list of values corresponding to the set of reactions in a network, representing the mass flow per unit time for each reaction in the network being analyzed.
As used herein, the term “concentration distribution” refers to a non-negative, quantitative list of values corresponding to the set of concentrations of variables, including but not limited to metabolites and enzymes, in a network, representing the mass, or number of moles, or number of molecules per unit volume for each variable.
As used herein, the term “biological network” refers to assembled reactions reflecting biological processes. Examples of biological networks include but are not limited to metabolic networks, signaling networks, and regulatory networks.
As used herein, the term “stoichiometric matrix” or “S” refers to a matrix with the stoichiometric coefficients for reactions represented by the columns and the substrates and products in the rows. The stoichiometric matrix may be written for any biological network. Examples include but are not limited to metabolic networks and signaling networks.
As used herein, the term “gradient matrix” or “G” refers to a matrix whose row entries are the first order partial derivatives of reaction fluxes with respect to the metabolites (see equation 2). The gradient matrix is specific for a particular biological network. The number of rows in the gradient matrix is equal to the number of reactions in the network. The number of columns in the gradient matrix is equal to the number of metabolites or compounds in the network.
As used herein, the term “Jacobian”, “Jacobi matrix”, “Jacobian matrix”, “J”, or “Jx” refers to a matrix of the first order partial derivatives of the linearized dynamic mass balances, in which the variables are the metabolites or compounds. The Jacobian matrix is derived when a system of differential equations is linearized around a reference state such that only the first order terms remain. The Jacobian is a square matrix and has the same number of rows and columns as there are metabolites or compounds in the network. “J” can be written as the matrix product of S and G (see equations 1 and 3).
As used herein, the term “flux Jacobian” or “Jv” refers to a matrix of the first order partial derivatives of the linearized dynamic mass balances, in which the variables are the fluxes. The flux Jacobian is a square matrix and has the same number of rows and columns as there are reactions or fluxes in the network. The flux Jacobian can be written as the matrix product of G and S (see equation 5).
As used herein, the term “kappa”, “kappa matrix”, or κ refers to a diagonal matrix composed of kinetic constants. This matrix can be calculated in a variety of ways. For example the rows of G can be scaled or normalized. The factorization may be carried out but not limited to the following manner, by factoring out the 1-norm, 2-norm, or sum of the rate constants with non-zero entries in the row.
As used herein, the term “gamma”, “gamma matrix”, or Γ refers to resulting matrix after κ has been factored out of G.
As used herein, the term “transpose”, transpose of a matrix, orT refers to a standard mathematical operation in which the elements of a matrix are reflected about the diagonal entries of the matrix, such that following the transformation the entries in the rows of the matrix now appear in the columns of the matrix and the entries in the columns of the matrix appear in the rows of the matrix.
As used herein, the term “time constant” refers to a numerical value with physical units of time, such as days, hours, minutes, and microseconds. Examples of time constants include but are not limited to the negative reciprocal of eigenvalues or the reciprocal entries of the kappa matrix, which may correspond to a particular reaction or set of reactions.
As used herein, the term “kinetic constant” or reaction rate constant refers to a time constant characteristic of a particular reaction or set of reactions. The physical units are [concentration]n+1*time−1, where n is the kinetic order of the corresponding reaction. For example, a zero order, first order, and second order kinetic constant has units [concentration]*time−1, time−1, and concentration−1*time−1 respectively.
As used herein, the term “steady state” refers to any of a number of conditions for which the flux distribution and concentration distribution do not change over time. Equilibrium is a special case which also satisfies the steady state. In general, unless otherwise specified, references to the steady state imply a non-equilibrium steady state.
As used herein, the term “concentration deviation variable”, “deviation variable”, or x′ is vector equal to the concentration variable vector minus a reference state. This reference state may be one of any number of concentration distributions including but not limited to steady state concentrations.
As used herein, the term “time scale decomposition” refers to the established practice of characterizing network dynamics by decomposing the Jacobian matrix into its constituent eigenvectors and eigenvalues, from which dynamically independent modes and time constants can be derived.
As used herein, the term “modal matrix” or “modes” refers to the inverse of the eigenvector matrix of J.
The reconstruction of a genome-scale reaction network requires the identification of all its chemical components and the chemical transformations that they participate in. This process primarily relies on annotated genomes and detailed bibliomic assessment. See Reed et al., 2006, which is hereby incorporated by reference in its entirety. The result of this process is the stoichiometric matrix, S, that is used in the dynamic mass balances
dx/dt=S·v(x)x=x0 at time t=0 (1)
that are the basis of all kinetic models. Here d(.)/dt denotes the time derivative, x is the vector of the concentrations of the compounds in the network and v(x) is the vector of the reaction rates. All biochemical transformations are fundamentally uni- or bi-molecular. Such reactions can be represented by mass action kinetics, or generalizations thereof. See Segel, 1975, which is hereby incorporated by reference in its entirety. The net reaction rate for every elementary reaction in a network can be represented by the difference between a forward and reverse flux (e.g. see
The characterization of the dynamic states of networks can be studied through numerical simulation or through using mathematical analysis. A simulation is context-dependent and represents a case study. Mathematical methods for the analysis of model characteristics typically rely on studying the properties of the transformation between the concentrations and fluxes. The analysis of such fundamental properties normally relies on the linearization of the governing equations at a defined condition. The linearization of the dynamic mass balance equations comes down to the linearization of the reaction rate vector to form the gradient matrix
G=[gij]=[∂vi/∂xj] (2)
in which ∂·/∂xj is the partial derivative operator. The Jacobian matrix at a reference state Xref can then be described:
dx′/dt=S·G·x′=J·x′ (3)
where x′=x−xref and J is the well-known Jacobian matrix. Analysis of the characteristics of Jacobian matrix is standard procedure in mathematical analysis of system dynamics (e.g. Strogatz, 1994, which is hereby incorporated by reference in its entirety). The application of these methods to biochemical networks has been carried out for decades. See Heinrich et al., 1977; Heinrich et al., 1991; Heinrich and Sonntag, 1982; Palsson et al., 1987, all of which are hereby incorporated by reference in their entireties. In recent years there has been renewed interest and recently further developments in the dynamic analysis of the properties of J have appeared. See Bruggeman et al., 2006; Famili et al., 2005; Grimbs et al., 2007; Kauffman et al., 2002; Steuer et al., 2006; Teusink et al., 2000, all of which are hereby incorporated by reference in their entireties.
The integration of genomic, bibliomic, fluxomic, and metabolomic data with thermodynamic information into dynamic models of metabolism is illustrated in
The availability of genomic and bibliomic data can be used to define S (
Table 1A compares properties of S and G, highlighting important features. For example, the formation of pools and reaction co-sets is determined by S. See Heinrich and Schuster, 1996; Jamshidi and Palsson, 2006; Papin et al., 2004, all of which are incorporated by reference in their entireties. However, time scale separation is in the realm of G. This decomposition factors the various underlying data needed for the formulation of genome-scale kinetic models. Furthermore, it illustrates the various underlying disciplinary interests that need to be considered and integrated to successfully achieve the analysis of dynamic states of genome-scale networks.
As mentioned supra, simulation is context-dependent and the respective values will depend on the domain being modeled.
For example, in the bioinformatic context, S is primarily derived from an annotated genomic sequence fortified with any direct bibliomic data on the organisms' gene products. The construction of G will rely on fluxomic and metabolmic data, in addition to direct kinetic characterization of individual reactions and assessment of thermodynamic properties.
In the physico-chemical context, S represents chemistry (i.e. stoichiometry of reactions), while G represents kinetics and thermodynamics. The chemical information is relatively easy to obtain, the thermodynamics information harder but possible (see Henry et al., 2007, which is incorporated by reference in its entirety), and the kinetic information is the hardest to acquire. The former two represent hard physico-chemical constraints, while the third represents biologically manipulable numbers; i.e., reaction rates are accelerated by enzymes.
In the biological and genetic considerations context, the matrix S is reconstructed based on the content of a genome and is a property of a species. It has thus been used productively for the analysis of distal causation (Fong and Palsson, 2004; Harrison et al., 2007; Ibarra et al., 2002; Pal et al., 2006, which are all incorporated by reference in their entireties). Distal (or ultimate) causation results from (genomic) changes that occur from generation to generation, in contrast to proximal (or proximate) causation which occurs against a fixed genetic background (i.e. an individual) (Mayr, 1961, which is incorporated by reference in its entirety). G is genetically derived and can represent the variations amongst individuals in a biopopulation. It is important in studying proximal causation and the differences in phenotypes of individuals in a biopopulation. For example, many disease states in higher order organisms result from disruptions or deficiencies in enzyme kinetics (Jamshidi et al., 2002, which is incorporated by reference in its entirety). These changes are reflected in G since this contains the information about kinetics, consequently the analysis of disease states, inter-individual differences, and transitions from a healthy to disease state in a particular individual will generally focus on G.
Generally speaking, the equations (S) will represent the pertinent physical laws and the matrices (G) will represent data. Whereas S is a ‘perfect’ matrix comprised of integers (i.e., digital), G is an analog matrix whose entries are real numbers, and values of G may only be known within an order of magnitude. From a numerical and mathematical standpoint, S is a well conditioned matrix comprised of integers (−2,−1,0,1,2), whereas G is an ill-conditioned matrix of real numbers that can differ up to 10 orders of magnitude in their numerical values. This property leads to time scale separation. Both matrices are sparse; that is most of their elements are zero.
The elements of G, in the absence of comprehensive fluxomic, metabolomic, or thermodynamic information, may be approximated to the appropriate order of magnitude in order to enable the construction of a kinetic model for analysis.
The wide differences in the numerical values of the entries of G leads to its factorization (G=κ·Γ) by scaling it by the length of its rows to yield a factorization of J into three matrices:
J=S·κ·Γ (4)
where κ is a diagonal matrix of the norms of the rows in G (FIG. 1.B.1). We emphasize that the ith column of S contains the stoichiometric coefficients of the ith reaction in the network and the ith row in G contains the forward and reverse rate constants of that same reaction. Thus, the reciprocal of the diagonal entries, 1/κii, correspond to characteristic time constants of the corresponding reactions. Their numerical values will differ significantly.
The factorization of the Jacobian in equation (4) shows that the study of the dynamic properties of biochemical networks can be formally decomposed into chemistry (represented by S), kinetics (represented by κ), and thermodynamic driving forces (represented by Γ). The effects of each can thus be formally determined. Chemistry and thermodynamics are physico-chemical properties, whereas the kinetic constants are biologically set through a natural selection process. The particular numerical values (chosen through the selection process) lead to the formation of biologically meaningful dynamic properties of the network. These biological features of the network can be assessed through time scale decomposition.
The glycolytic pathway (
It follows from Equation (2) that if sij=0 then gji is also zero; that is if a compound does not participate in a reaction, it has no kinetic effect on it. Conversely, if sij≠0 then gji is also not zero. Further inspection reveals the property that S is structurally similar to −GT as illustrated in FIG. 1.B.2. Thus, the non-zero entries in S have corresponding non-zero elements in −GT, but with a different numerical value. This fundamental feature shows that the topology of the network as reflected in S has a dominant effect on its dynamic features; providing another example of the biological principle that structure has a dominant effect on function.
The duality relationship illustrates that fluxes or concentrations can be used as the independent variables: A flux deviation variable, v′ can be defined such that v′=G·x′, from which it follows that
dv′/dt=G·S·v′=Jv·v′ (5)
This transformation illustrates the switch from concentrations to fluxes as the independent variables. Fluxes can tie together the multiple parts of a network to form its overall functions. Furthermore, the ability to relate the fluxes and concentrations independently of a specific rate law formulation, if the elements of G can be approximated, has significant implications for the construction and analysis of kinetic networks.
The two Jacobian matrices, G·S and S·G, are similar and share eigenvalues. Equation (5) is the ‘dynamic’ flux balance equation since the variables in it are the fluxes, v′. One can thus analyze network properties either in terms of concentrations or fluxes as illustrated within FIG. 1.B.3. The fluxes are ‘network’ variables, as they tie all the components together, while the concentrations are ‘component’ variables.
The properties of the Jacobian matrix that determine the characteristics of the network dynamics are its eigen-properties. The eigenvalues are network-based time constants (in contrast to the reaction based time constants in κ). Formally, the standard eigen-analysis is performed by the diagonalization of the Jacobian matrix as:
J=M·Λ·M−1 (6)
where M is the matrix of eigenvectors, Λ is a diagonal matrix of the eigenvalues, and M−1 is the matrix of eigenrows. If the decomposition of equation (6) is introduced into equation (3) we obtain differential equations for the modes (m=M−1·x′)
dm/dt=Λ·m, or dmi/dt=λl·mi, or mi(t)=mi,0exp(λlt) (7)
or a set of completely decoupled dynamic variables; that is, each mi moves on it own time scale defined by λi independent of all the other mj. The eigenrows give lumped or aggregate variables that move independently on each time scale since m is a set of dynamically independent variables. The eigenvectors form a natural coordinate system to describe the dynamics motion of the modes. We note that this decomposition can be applied to J or Jv (
The distribution of the numerical values of the eigenvalues is the basis for time scale separation. Time scale separation forms the basis for decomposition of biochemical reaction networks in time and the interpretation of the biochemical events that take place on the various time constants. Time scale separation has been analyzed in the context of biological networks and shown to lend insight and enable the simplification of theses networks (Heinrich and Sonntag, 1982; Palsson et al., 1987; Reich and Selkov, 1981, which are all incorporated by reference in their entireties).
Methods described herein enable the explicit construction of kinetic models using the incorporation of different data types. The properties, completeness, and accuracy of the data can be explicitly traced to dynamic properties. This decomposition will not only tie the models directly back to the data, but also explicitly gives us the parts of the model that are under biological control and subject to change with adaptation or evolution. Measurement uncertainties are primarily in κ and are subject to evolutionary changes. These ‘biological’ design parameters will likely need to be dealt with through the use of methods that bracket the range of values within uncertainty limitations.
Embodiments show how data can be used for the full delineation of the chemical equations that make up a network and ideally their representation as net combinations of elementary reactions (i.e. vnet=v+−v−). In this format the structure of the gradient matrix can be determined, but integration of multiple networks is also possible and the explicit analysis of the effects of regulatory molecules is also enabled. Furthermore, the underlying bilinear kinetic nature of network dynamics is explicitly recognized by embodiments disclosed herein, as all chemical reactions are combinations of bilinear interactions, including regulatory processes (Segel, 1975).
The duality between flux distributions and concentration distributions enabled by the transformations of the data matrices, G·S and S·G allows representation of network kinetics in terms of concentration variables or flux variables. The systems biology paradigm of ‘components→networks→computational models→physiological states’ naturally leads to the use of fluxes as variables to characterize the functional states of a network. Fluxomic data ties components in a network together and is interpreted through network models, whereas concentration data is component data. Fluxes have been widely used for steady state analysis and can now be used to study dynamic states as well.
Time scale decomposition has been known and studied for several decades (Heinrich et al., 1977; Heinrich and Sonntag, 1982; Palsson et al., 1987; Reich and Selkov, 1981, which are all hereby incorporated by reference in their entireties). Such studies have primarily been performed for small-scale models today, but their conceptual foundation has been established. At larger scales, new issues will arise. These are likely to include, data sensitivity, course-graining and modularization of physiological functions in time. New methods to study the bases vectors in M and M−1 that directly relate them to biochemistry and physiological functions need to be established. The promise of the elucidation of (dynamic) structure-(physiological) function relationships (Palsson et al., 1987, which is hereby incorporated by reference in its entirety) may now grow to large-scale networks.
As an example, models derived by methods disclosed herein can be used to determine the effects of changes from aerobic to anaerobic conditions, such as occurs in skeletal muscles during exercise or in tumors, or to determine the effect of various dietary changes. These modes can also be used to determine the dynamics and consequences of genetic defects, such as deficiencies in metabolic enzymes such as phosphofructokinase, phosphoglycerate kinase, phosphoglycerate mutase, lactate dehydrogenase and adenosine deaminase. As another example, the dynamics and consequences of defects in metabolic transporters, such as ABC transporters, the SLC class of transporters, as well as the cytochrome P450 class of enzymes, may be determined.
The kinetic network models can also be used to choose appropriate targets for drug design in individuals or groups of individuals. Such targets include genes, proteins or reactants, which when modulated positively or negatively in a simulation produce a desired therapeutic result. Models and methods disclosed herein can also be used to predict the effects of a therapeutic agent or dietary supplement on a cellular function of interest. Likewise, models and methods can be used to predict both desirable and undesirable side effects of the therapeutic agent on an interrelated cellular function in the target cell, as well as the desirable and undesirable effects that may occur in other cell types. Thus, models and methods disclosed herein can make the drug development process more rapid and cost effective than is currently possible.
The kinetic models can be used to develop individualized drug therapies and doses. Measuring a set of fluxes and concentrations in an individual before and after giving a dose of a therapeutic agent to the individual and constructing kinetic network models according to the measured data specific to the individual can be constructed. Analysis of the dynamic and functional properties of the individual specific network can be used to determine whether the therapeutic agent is effective in the desired outcome and the appropriate dose ranges if applicable.
As used herein, the term reaction is intended to mean a conversion that consumes a substrate or forms a product that occurs in a biological network. The term can include a conversion that occurs due to the activity of one or more enzymes that are genetically encoded by the human genome. The term can also include a conversion that occurs spontaneously in a cell, such as in human or mammal. Conversions included in the term include, for example, changes in chemical composition such as those due to nucleophilic or electrophilic addition, nucleophilic or electrophilic substitution, elimination, isomerization, deamination, phosphorylation, methylation, glycolysation, reduction, oxidation or changes in location such as those that occur due to a transport reaction that moves one or more reactants within the same compartment or from one cellular compartment to another. In the case of a transport reaction, the substrate and product of the reaction can be chemically the same and the substrate and product can be differentiated according to location in a particular cellular compartment. Thus, a reaction that transports a chemically unchanged reactant from a first compartment to a second compartment has as its substrate the reactant in the first compartment and as its product the reactant in the second compartment. It will be understood that when used in reference to an in silico model or data structure, a reaction is intended to be a representation of a chemical conversion that consumes a substrate or produces a product.
As used herein, the term reaction is intended to mean a chemical that is a substrate or a product of a reaction that occurs in a biological network. The term can include substrates or products of reactions performed by one or more enzymes encoded by human gene(s), reactions occurring in cells that are performed by one or more non-genetically encoded macromolecule, protein or enzyme, or reactions that occur spontaneously in a cell. Metabolites are understood to be reactants within the meaning of the term. It will be understood that when used in reference to an in silico model or data structure, a reactant is intended to be a representation of a chemical that is a substrate or a product of a reaction that occurs in a cell.
As used herein the term “substrate” is intended to mean a reactant that can be converted to one or more products by a reaction. The term can include, for example, a reactant that is to be chemically changed due to nucleophilic or electrophilic addition, nucleophilic or electrophilic substitution, elimination, isomerization, deamination, phosphorylation, methylation, reduction, oxidation or that is to change location such as by being transported across a membrane or to a different compartment.
As used herein, the term “product” is intended to mean a reactant that results from a reaction with one or more substrates. The term can include, for example, a reactant that has been chemically changed due to nucleophilic or electrophilic addition, nucleophilic or electrophilic substitution, elimination, isomerization, deamination, phosphorylation, methylation, reduction or oxidation or that has changed location such as by being transported across a membrane or to a different compartment.
As used herein, the term “stoichiometric coefficient” is intended to mean a numerical constant correlating the number of one or more reactants and the number of one or more products in a chemical reaction. Typically, the numbers are integers as they denote the number of molecules of each reactant in an elementally balanced chemical equation that describes the corresponding conversion. However, in some cases the numbers can take on non-integer values, for example, when used in a lumped reaction or to reflect empirical data.
As used herein, the term “plurality,” when used in reference to reactions or reactants is intended to mean at least 2 reactions or reactants. The term can include any number of reactions or reactants in the range from 2 to the number of naturally occurring reactants or reactions for a particular cell. Thus, the term can include, for example, at least 10, 20, 30, 50, 100, 150, 200, 300, 400, 500, 600 or more reactions or reactants. The number of reactions or reactants can be expressed as a portion of the total number of naturally occurring reactions for a particular cell such as at least 20%, 30%, 50%, 60%, 75%, 90%, 95% or 98% of the total number of naturally occurring reactions that occur in the particular cell.
As used herein, the term “data structure” is intended to mean a physical or logical relationship among data elements, designed to support specific data manipulation functions. The term can include, for example, a list of data elements that can be added combined or otherwise manipulated such as a list of representations for reactions from which reactants can be related in a matrix or network. The term can also include a matrix that correlates data elements from two or more lists of information such as a matrix that correlates reactants to reactions. Information included in the term can represent, for example, a substrate or product of a chemical reaction, a chemical reaction relating one or more substrates to one or more products, a constraint placed on a reaction, a stoichiometric coefficient, a rate constant, a gradient matrix coefficient, a gamma matrix coefficient, or a Jacobian matrix coefficient.
As used herein, the term “constraint”, reaction constraint, or flux constraint is intended to mean an upper or lower boundary for a reaction. A boundary can specify a minimum or maximum flow of mass, electrons or energy through a reaction. A boundary can further specify directionality of a reaction. A boundary can be a constant value such as zero, infinity, or a numerical value such as an integer and non-integer. Alternatively, a boundary can be a variable boundary value as set forth below.
As used herein, the term “concentration constraint” is intended to mean an upper and lower boundary for the concentration of a metabolite or chemical species in a biological network. The boundary can specify the amount of mass or number of moles, or number of molecules per unit volume. This may be within a tissue, a sample of biological fluid, a cell, a subcompartment of a cell, or another confined quantity. The boundary is a non-negative real number including zero, integers, and non-integers. Alternatively, a boundary for a concentration constraint can be a variable boundary as set forth below, with the stipulation that negative values are not permitted.
As used herein, the term “kinetic constraint” is intended to mean an upper and lower boundary for the kinetic constant for a reaction in a biological network. The boundary can specify the time constant characteristic of a particular reaction or set of reactions. The physical units are [concentration]−n+1*time−1, where n is the kinetic order of the corresponding reaction. For example, a zero order, first order, and second order kinetic constant has units [concentration]*time−1, time−1, and concentration−1*time−1 respectively. The boundary is a non-negative real number including zero, integers, and non-integers. Alternatively, a boundary for a concentration constraint can be a variable boundary as set forth below, with the stipulation that negative values are not permitted.
As used herein, the term “equilibrium constant constraint” is intended to mean an upper and lower boundary for the equilibrium constant for a reaction in a biological network. The boundary can specify the equilibrium constant of a particular reaction or set of reactions. The boundary is a non-negative real number including zero, integers, and non-integers. Alternatively, a boundary for a concentration constraint can be a variable boundary as set forth below, with the stipulation that negative values are not permitted.
As defined herein, the term “uniquely defined constraint” is intended to mean a constraint whose lower bound is equal in value to the upper bound.
As used herein, the term “variable,” when used in reference to a constraint is intended to mean capable of assuming any of a set of values in response to being acted upon by a constraint function. The term “function,” when used in the context of a constraint, is intended to be consistent with the meaning of the term as it is understood in the computer and mathematical arts. A function can be binary such that changes correspond to a reaction being off or on. Alternatively, continuous functions can be used such that changes in boundary values correspond to increases or decreases in activity. Such increases or decreases can also be binned or effectively digitized by a function capable of converting sets of values to discreet integer values. A function included in the term can correlate a boundary value with the presence, absence or amount of a biochemical reaction network participant such as a reactant, reaction, enzyme or gene. A function included in the term can correlate a boundary value with an outcome of at least one reaction in a reaction network that includes the reaction that is constrained by the boundary limit. A function included in the term can also correlate a boundary value with an environmental condition such as time, pH, temperature or redox potential.
As used herein, the term “activity,” when used in reference to a reaction, is intended to mean the amount of product produced by the reaction, the amount of substrate consumed by the reaction or the rate at which a product is produced or a substrate is consumed. The amount of product produced by the reaction, the amount of substrate consumed by the reaction or the rate at which a product is produced or a substrate is consumed can also be referred to as the flux for the reaction.
As used herein, the term “constraint set” refers to the collection of flux, concentration, and equilibrium constraints for a biological network. Constraint sets for kinetic constants may be found by solving a steady state mass conservation relationship comprising a stoichiometric matrix and flux vector.
As used herein, the term “allosteric regulation” refers to the regulation of a protein or enzyme by a molecule that binds to a site other than the primary active site, thereby altering the activity of the enzyme and the corresponding reaction that it regulates.
As used herein, the term “activate” or activation refers to an effect a compound has on another compound, serving to alter the constraints in a positive manner, such as increasing the activity of a reaction. This includes but is not limited to allosteric and non-allosteric regulation of enzymes.
As used herein, the term “inhibit” or inhibition refers to an effect a compound has on another compound, serving to alter the constraints in a negative manner, such as decreasing the activity of a reaction. This includes but is not limited to allosteric and non-allosteric regulation of enzymes.
As used herein, the term “growth” refers to the production of a weighted sum of metabolites identified as biomass components.
As used herein, the term “energy production” refers to the production of metabolites that store energy in their chemical bonds, particularly high energy phosphate bonds such as ATP and GTP.
As used herein, the term “redox equivalent” refers to a metabolite that can alter the oxidation state of other metabolites, generally through the transfer of electrons and or protons.
As used herein, the term “biomass” refers to a set of metabolites that are required for cellular functions including replication, maintenance, and/or survival. As such, the term “biomass precursor” refers to a metabolite that can be converted to or is a member of the class of biomass metabolites.
As used herein, the term “bioactive small molecule” refers to a metabolite that can inhibit or activate another biological compound, including metabolites, proteins, and nucleic acid polymers.
As used herein, the term “cofactor” refers to a metabolite or chemical compound that is required for a catalytic activity of an enzyme to be carried out.
As used herein, the term “optimization problem” refers to a problem whose solution is obtained through the calculation of a minimum (or maximum) value of the objective function.
As used herein, the term “objective function” refers to a physiological function that can be described in terms of the data structure of a network, generally in terms of one or more weighted sum or product of reactions.
Methods and models disclosed herein can be applied to any H. sapiens cell type at any stage of differentiation, including, for example, embryonic stem cells, hematopoietic stem cells, differentiated hematopoietic cells, skeletal muscle cells, cardiac muscle cells, smooth muscle cells, skin cells, nerve cells, kidney cells, pulmonary cells, liver cells, adipocytes and endocrine cells (e.g., beta islet cells of the pancreas, mammary gland cells, adrenal cells, and other specialized hormone secreting cells).
Methods and models disclosed herein can be applied to normal cells or pathological cells. Normal cells that exhibit a variety of physiological activities of interest, including homeostasis, proliferation, differentiation, apoptosis, contraction and motility, can be modeled. Pathological cells can also be modeled, including cells that reflect genetic or developmental abnormalities, nutritional deficiencies, environmental assaults, infection (such as by bacteria, viral, protozoan or fungal agents), neoplasia, aging, altered immune or endocrine function, tissue damage, or any combination of these factors. The pathological cells can be representative of any type of mammalian, primate or non-primate pathology, including, for example, various metabolic disorders of carbohydrate, lipid or protein metabolism, obesity, diabetes, cardiovascular disease, fibrosis, various cancers, kidney failure, immune pathologies, neurodegenerative diseases, and various monogenetic diseases described found in humans and documented in the Online Mendelian Inheritance in Man database (Center for Medical Genetics, Johns Hopkins University (Baltimore, Md.) and National Center for Biotechnology Information, National Library of Medicine (Bethesda, Md.).
Methods and models disclosed herein can also be applied to cells undergoing therapeutic perturbations, such as cells treated with drugs that target participants in a reaction network, cells treated with gene-based therapeutics that increase or decrease expression of an encoded protein, and cells treated with radiation. As used herein, the term “drug” refers to a compound of any molecular nature with a known or proposed therapeutic function, including, for example, small molecule compounds, peptides and other macromolecules, peptidomimetics and antibodies, any of which can optionally be tagged with cytostatic, targeting or detectable moieties. The term “gene-based therapeutic” refers to nucleic acid therapeutics, including, for example, expressible genes with normal or altered protein activity, antisense compounds, ribozymes, DNAzymes, RNA interference compounds (RNAi) and the like. The therapeutics can target any reaction network participant, in any cellular location, including participants in extracellular, cell surface, cytoplasmic, mitochondrial and nuclear locations. Experimental data that are gathered on the response of cells to therapeutic treatment, such as alterations in gene or protein expression profiles, can be used to tailor a network for a pathological state of a particular cell type.
Methods and models disclosed herein can be applied to H. sapiens cells as they exist in any form, such as in primary cell isolates or in established cell lines, or in the whole body, in intact organs or in tissue explants. Accordingly, the methods and models can take into account intercellular communications and/or inter-organ communications, the effect of adhesion to a substrate or neighboring cells (such as a stem cell interacting with mesenchymal cells or a cancer cell interacting with its tissue microenvironment, or beta-islet cells without normal stroma), and other interactions relevant to multicellular systems.
The reactants to be used in a reaction network data structure of the invention can be obtained from or stored in a compound database. As used herein, the term “compound database” is intended to mean a computer readable medium containing a plurality of molecules that includes substrates and products of biological reactions. The plurality of molecules can include molecules found in multiple organisms, thereby constituting a universal compound database. Alternatively, the plurality of molecules can be limited to those that occur in a particular organism, thereby constituting an organism-specific compound database. Each reactant in a compound database can be identified according to the chemical species and the cellular compartment in which it is present. Thus, for example, a distinction can be made between glucose in the extracellular compartment versus glucose in the cytosol. Additionally each of the reactants can be specified as a metabolite of a primary or secondary metabolic pathway. Although identification of a reactant as a metabolite of a primary or secondary metabolic pathway does not indicate any chemical distinction between the reactants in a reaction, such a designation can assist in visual representations of large networks of reactions.
As used herein, the term “compartment” is intended to mean a subdivided region containing at least one reactant, such that the reactant is separated from at least one other reactant in a second region. A subdivided region included in the term can be correlated with a subdivided region of a cell. Thus, a subdivided region included in the term can be, for example, the intracellular space of a cell; the extracellular space around a cell; the periplasmic space; the interior space of an organelle such as a mitochondrium, endoplasmic reticulum, Golgi apparatus, vacuole or nucleus; or any subcellular space that is separated from another by a membrane or other physical barrier. Subdivided regions can also be made in order to create virtual boundaries in a reaction network that are not correlated with physical barriers. Virtual boundaries can be made for the purpose of segmenting the reactions in a network into different compartments or substructures.
As used herein, the term “substructure” is intended to mean a portion of the information in a data structure that is separated from other information in the data structure such that the portion of information can be separately manipulated or analyzed. The term can include portions subdivided according to a biological function including, for example, information relevant to a particular metabolic pathway such as an internal flux pathway, exchange flux pathway, central metabolic pathway, peripheral metabolic pathway, or secondary metabolic pathway. The term can include portions subdivided according to computational or mathematical principles that allow for a particular type of analysis or manipulation of the data structure.
The reactions included in a reaction network data structure can be obtained from a metabolic reaction database that includes the substrates, products, and stoichiometry of a plurality of biological reactions. The reactants in a reaction network data structure can be designated as either substrates or products of a particular reaction, each with a stoichiometric coefficient assigned to it to describe the chemical conversion taking place in the reaction. Each reaction is also described as occurring in either a reversible or irreversible direction. Reversible reactions can either be represented as one reaction that operates in both the forward and reverse direction or be decomposed into two irreversible reactions, one corresponding to the forward reaction and the other corresponding to the backward reaction.
Reactions included in a reaction network data structure can include intra-system or exchange reactions. Intra-system reactions are the chemically and electrically balanced interconversions of chemical species and transport processes, which serve to replenish or drain the relative amounts of certain metabolites. These intra-system reactions can be classified as either being transformations or translocations. A transformation is a reaction that contains distinct sets of compounds as substrates and products, while a translocation contains reactants located in different compartments. Thus, a reaction that simply transports a metabolite from the extracellular environment to the cytosol, without changing its chemical composition is solely classified as a translocation, while a reaction such as the phosphotransferase system (PTS) which takes extracellular glucose and converts it into cytosolic glucose-6-phosphate is a translocation and a transformation.
Exchange reactions are those which constitute sources and sinks, allowing the passage of metabolites into and out of a compartment or across a hypothetical system boundary. These reactions are included in a model for simulation purposes and represent the metabolic demands placed on a cell. While they may be chemically balanced in certain cases, they are typically not balanced and can often have only a single substrate or product. As a matter of convention the exchange reactions are further classified into demand exchange and input/output exchange reactions.
Input/output exchange reactions are used to allow extracellular reactants to enter or exit the reaction network represented by a model of the invention. For each of the extracellular metabolites a corresponding input/output exchange reaction can be created. These reactions can either be irreversible or reversible with the metabolite indicated as a substrate with a stoichiometric coefficient of one and no products produced by the reaction. This particular convention is adopted to allow the reaction to take on a positive flux value (activity level) when the metabolite is being produced or removed from the reaction network and a negative flux value when the metabolite is being consumed or introduced into the reaction network. These reactions will be further constrained during the course of a simulation to specify exactly which metabolites are available to the cell and which can be excreted by the cell.
A demand exchange reaction is always specified as an irreversible reaction containing at least one substrate. These reactions are typically formulated to represent the production of an intracellular metabolite by the metabolic network or the aggregate production of many reactants in balanced ratios such as in the representation of a reaction that leads to biomass formation, also referred to as growth.
A demand exchange reaction can be introduced for any metabolite in a model of the invention. Most commonly these reactions are introduced for metabolites that are required to be produced by the cell for the purposes of creating a new cell such as amino acids, nucleotides, phospholipids, and other biomass constituents, or metabolites that are to be produced for alternative purposes. Once these metabolites are identified, a demand exchange reaction that is irreversible and specifies the metabolite as a substrate with a stoichiometric coefficient of unity can be created. With these specifications, if the reaction is active it leads to the net production of the metabolite by the system meeting potential production demands. Examples of processes that can be represented as a demand exchange reaction in a reaction network data structure and analyzed by the methods of the invention include, for example, production or secretion of an individual protein; production or secretion of an individual metabolite such as an amino acid, vitamin, nucleoside, antibiotic or surfactant; production of ATP for extraneous energy requiring processes such as locomotion; or formation of biomass constituents.
In addition to these demand exchange reactions that are placed on individual metabolites, demand exchange reactions that utilize multiple metabolites in defined stoichiometric ratios can be introduced. These reactions are referred to as aggregate demand exchange reactions. An example of an aggregate demand reaction is a reaction used to simulate the concurrent growth demands or production requirements associated with cell growth that are placed on a cell, for example, by simulating the formation of multiple biomass constituents simultaneously at a particular cellular growth rate.
Depending upon the particular environmental conditions being tested and the desired activity, a reaction network data structure can contain smaller numbers of reactions such as at least 200, 150, 100 or 50 reactions. A reaction network data structure having relatively few reactions can provide the advantage of reducing computation time and resources required to perform a simulation. When desired, a reaction network data structure having a particular subset of reactions can be made or used in which reactions that are not relevant to the particular simulation are omitted. Alternatively, larger numbers of reactions can be included in order to increase the accuracy or molecular detail of the methods of the invention or to suit a particular application. Thus, a reaction network data structure can contain at least 300, 350, 400, 450, 500, 550, 600 or more reactions up to the number of reactions that occur in a particular cell or that are desired to simulate the activity of the full set of reactions occurring in the particular human cell.
A reaction network data structure or index of reactions used in the data structure such as that available in a metabolic reaction database, as described herein, can be annotated to include information about a particular reaction. A reaction can be annotated to indicate, for example, assignment of the reaction to a protein, macromolecule or enzyme that performs the reaction, assignment of a gene(s) that codes for the protein, macromolecule or enzyme, the Enzyme Commission (EC) number of the particular metabolic reaction or Gene Ontology (GO) number of the particular metabolic reaction, a subset of reactions to which the reaction belongs, citations to references from which information was obtained, or a level of confidence with which a reaction is believed to occur in a particular human cell. A computer readable medium of the invention can include a gene database containing annotated reactions. Such information can be obtained during the course of building a metabolic reaction database or model of the invention as described below.
Thus, a method is provided for making or expanding a data structure relating a plurality of H. sapiens reactants to a plurality of H. sapiens reactions in a computer readable medium. The method includes the steps of: (a) identifying a plurality of H. sapiens reactions and a plurality of H. sapiens reactants that are substrates and products of the H. sapiens reactions; (b) relating the plurality of H. sapiens reactants to the plurality of H. sapiens reactions in a data structure, wherein each of the H. sapiens reactions includes one or more reactants identified as a substrate of the reaction, one or more reactants identified as a product of the reaction and a stoichiometric coefficient relating the substrate and the product; (c) making a constraint set for the plurality of H. sapiens reactions; (d) providing an objective function; (e) determining at least one flux distribution that minimizes or maximizes the objective function when the constraint set is applied to the data structure, (f) if at least one flux distribution is not predictive of human physiology, then adding a reaction to or deleting a reaction from the data structure and repeating step (e), and (g) if at least one flux distribution is predictive of human physiology, then storing the data structure in a computer readable medium.
The majority of the reactions occurring in human reaction networks are catalyzed by enzymes/proteins, which are created through the transcription and translation of the genes found on the chromosome(s) in the cell. The remaining reactions occur through non-enzymatic processes. Furthermore, a reaction network data structure can contain reactions that add or delete steps to or from a particular reaction pathway. For example, reactions can be added to optimize or improve performance of a H. sapiens model in view of empirically observed activity. Alternatively, reactions can be deleted to remove intermediate steps in a pathway when the intermediate steps are not necessary to model flux through the pathway. For example, if a pathway contains three nonbranched steps, the reactions can be combined or added together to give a net reaction, thereby reducing memory required to store the reaction network data structure and the computational resources required for manipulation of the data structure.
A reaction network data structure as described herein can be used to determine the activity of one or more reactions in a plurality of H. sapiens reactions independent of any knowledge or annotation of the identity of the protein that performs the reaction or the gene encoding the protein. A model that is annotated with gene or protein identities can include reactions for which a protein or encoding gene is not assigned. While a large portion of the reactions in a cellular metabolic network are associated with genes in the organism's genome, there are also a substantial number of reactions included in a model for which there are no known genetic associations. Such reactions can be added to a reaction database based upon other information that is not necessarily related to genetics such as biochemical or cell based measurements or theoretical considerations based on observed biochemical or cellular activity. For example, there are many reactions that are not protein-enabled reactions. Furthermore, the occurrence of a particular reaction in a cell for which no associated proteins or genetics have been currently identified can be indicated during the course of model building by the iterative model building methods as disclosed herein.
The reactions in a reaction network data structure or reaction database can be assigned to subsystems by annotation, if desired. The reactions can be subdivided according to biological criteria, such as according to traditionally identified metabolic pathways (glycolysis, amino acid metabolism and the like) or according to mathematical or computational criteria that facilitate manipulation of a model that incorporates or manipulates the reactions. Methods and criteria for subdividing a reaction database are described in further detail in Schilling et al., J. Theor. Biol. 203:249-283 (2000), which is hereby incorporated by reference in its entirety. The use of subsystems can be advantageous for a number of analysis methods, such as extreme pathway analysis, and can make the management of model content easier. Although assigning reactions to subsystems can be achieved without affecting the use of the entire model for simulation, assigning reactions to subsystems can allow a user to search for reactions in a particular subsystem, which may be useful in performing various types of analyses. Therefore, a reaction network data structure can include any number of desired subsystems including, for example, 2 or more subsystems, 5 or more subsystems, 10 or more subsystems, 25 or more subsystems or 50 or more subsystems.
Flux constraints can be placed on the value of any of the fluxes in the metabolic network using a constraint set. These constraints can be representative of a minimum or maximum allowable flux through a given reaction, possibly resulting from a limited amount of an enzyme present. Additionally, the constraints can determine the direction or reversibility of any of the reactions or transport fluxes in the reaction network data structure.
The ability of a reaction to be actively occurring is dependent on a large number of additional factors beyond just the availability of substrates. These factors, which can be represented as variable constraints in the models and methods of the invention include, for example, the presence of cofactors necessary to stabilize the protein/enzyme, the presence or absence of enzymatic inhibition and activation factors, the active formation of the protein/enzyme through translation of the corresponding mRNA transcript, the transcription of the associated gene(s) or the presence of chemical signals and/or proteins that assist in controlling these processes that ultimately determine whether a chemical reaction is capable of being carried out within an organism. Of particular importance in the regulation of human cell types is the implementation of paracrine and endocrine signaling pathways to control cellular activities. In these cases a cell secretes signaling molecules that may be carried far afield to act on distant targets (endocrine signaling), or act as local mediators (paracrine signaling). Examples of endocrine signaling molecules include hormones such as insulin, while examples of paracrine signaling molecules include neurotransmitters such as acetylcholine. These molecules induce cellular responses through signaling cascades that affect the activity of biochemical reactions in the cell. Regulation can be represented in an in silico H. sapiens model by providing a variable constraint.
The methods described herein can be implemented on any conventional host computer system, such as those based on Intel® or AMID® microprocessors and running Microsoft Windows operating systems. Other systems, such as those using the UNIX or LINUX operating system and based on IBM®, DEC® or Motorola® microprocessors are also contemplated. The systems and methods described herein can also be implemented to run on client-server systems and wide-area networks, such as the Internet.
Software to implement a method or model of the invention can be written in any well-known computer language, such as Java, C, C++, Visual Basic, FORTRAN or COBOL and compiled using any well-known compatible compiler. The software of the invention normally runs from instructions stored in a memory on a host computer system. A memory or computer readable medium can be a hard disk, floppy disc, compact disc, DVD, magneto-optical disc, Random Access Memory, Read Only Memory or Flash Memory. The memory or computer readable medium used in the invention can be contained within a single computer or distributed in a network. A network can be any of a number of conventional network systems known in the art such as a local area network (LAN) or a wide area network (WAN). Client-server environments, database servers and networks that can be used in the invention are well known in the art. For example, the database server can run on an operating system such as UNIX, running a relational database management system, a World Wide Web application and a World Wide Web server. Other types of memories and computer readable media are also contemplated to function within the scope of the invention.
The data matrices constructed by the methods described in this invention can be represented in a markup language format including, for example, Standard Generalized Markup Language (SGML), Hypertext markup language (HTML) or Extensible Markup language (XML). Markup languages can be used to tag the information stored in a database or data structure of the invention, thereby providing convenient annotation and transfer of data between databases and data structures. In particular, an XML format can be useful for structuring the data representation of reactions, reactants and their annotations; for exchanging database contents, for example, over a network or internet; for updating individual elements using the document object model; or for providing differential access to multiple users for different information content of a data base or data structure of the invention. XML programming methods and editors for writing XML code are known in the art as described, for example, in Ray, Learning XML O'Reilly and Associates, Sebastopol, Calif. (2001).
A computer system of the invention can further include a user interface capable of receiving a representation of one or more reactions. A user interface of the invention can also be capable of sending at least one command for modifying the data structure, the constraint set or the commands for applying the constraint set to the data representation, or a combination thereof. The interface can be a graphic user interface having graphical means for making selections such as menus or dialog boxes. The interface can be arranged with layered screens accessible by making selections from a main screen. The user interface can provide access to other databases useful in the invention such as a metabolic reaction database or links to other databases having information relevant to the reactions or reactants in the reaction network data structure or to human physiology. Also, the user interface can display a graphical representation of a reaction network or the results of a simulation using a model of the invention.
Once an initial reaction network data structure and set of constraints has been created, a model disclosed herein can be tested by preliminary simulation. During preliminary simulation, gaps in the network or “dead-ends” in which a metabolite can be produced but not consumed or where a metabolite can be consumed but not produced can be identified. Based on the results of preliminary simulations areas of the metabolic reconstruction that require an additional reaction can be identified. The determination of these gaps can be readily calculated through appropriate queries of the reaction network data structure and need not require the use of simulation strategies, however, simulation would be an alternative approach to locating such gaps.
In the preliminary simulation testing and model content refinement stage an existing model may be subjected to a series of functional tests to determine if it can perform basic requirements such as the ability to produce the required biomass constituents and generate predictions concerning the basic physiological characteristics of the particular organism strain being modeled. The more preliminary testing that is conducted, typically, the higher the quality of the model that will be generated. Typically the majority of the simulations used in this stage of development will be single optimizations. A single optimization can be used to calculate a single flux distribution demonstrating how metabolic resources are routed determined from the solution to one optimization problem. An optimization problem can be solved using linear programming as demonstrated in the Examples below. The result can be viewed as a display of a flux distribution on a reaction map. Temporary reactions can be added to the network to determine if they should be included into the model based on modeling/simulation requirements.
Once a model is sufficiently complete with respect to the content of the reaction network data structure according to the criteria set forth above, the model can be used to simulate activity of one or more reactions in a reaction network. The results of a simulation can be displayed in a variety of formats including, for example, a table, graph, reaction network, flux distribution map or a as a modal matrix.
If no features are observed, or further analysis is required, parameters or steady state conditions are modified in the model 611 so as to simulate a different dynamic, and the simulation 604 observed again until features of interest may be identified.
Alternatively, where features of interest have been identified, physiological materials relevant to them may be prepared 612 and the process completes 613. In the domain of drug therapy, for example, the preparation of physiological materials may comprise the identification and synthesis of compounds which will generate the desired metabolic behavior responsible for generating the features of interest previously modeled. One skilled in the art, however, will readily recognize that the model will lend itself to a wide variety of domains requiring a comprehensive understanding of large and complex metabolic pathways. Accordingly, the information derived from the model may be used in numerous physical applications, from confirming disease diagnosis, to identifying environmental factors impacting individual health.
Thus, methods provided herein can predict physiological function. In one embodiment, a method is provided, the method including providing a data structure relating a plurality of H. sapiens reactants to a plurality of H. sapiens reactions, wherein each of the H. sapiens reactions includes one or more reactants identified as a substrate of the reaction, one or more reactants identified as a product of the reaction and a stoichiometric coefficient relating the substrate and the product; providing a constraint set for the plurality of H. sapiens reactions; providing an objective function, and determining at least one flux distribution that minimizes or maximizes the objective function when the constraint set is applied to the data structure, thereby predicting a H. sapiens physiological function.
As used herein, the term “physiological function,” when used in reference to H. sapiens, is intended to mean an activity of a human cell as a whole. An activity included in the term can be the magnitude or rate of a change from an initial state of a human cell to a final state of the human cell. An activity can be measured qualitatively or quantitatively. An activity included in the term can be, for example, growth, energy production, redox equivalent production, biomass production, development, or consumption of carbon, nitrogen, sulfur, phosphate, hydrogen or oxygen. An activity can also be an output of a particular reaction that is determined or predicted in the context of substantially all of the reactions that affect the particular reaction in a human cell or substantially all of the reactions that occur in a human cell. Examples of a particular reaction included in the term are production of biomass precursors, production of a protein, production of an amino acid, production of a purine, production of a pyrimidine, production of a lipid, production of a fatty acid, production of a cofactor, production of a hormone, production of a bioactive small molecule, or transport of a metabolite. A physiological function can include an emergent property which emerges from the whole but not from the sum of parts where the parts are observed in isolation (see for example, Palsson Nat. Biotech 18:1147-1150 (2000)).
A physiological function of H. sapiens reactions can be determined using dynamic phase plane analysis.
A physiological function of H. sapiens reactions can also be determined using a reaction map to display a flux distribution. A reaction map can be used to view reaction networks at a variety of levels. In the case of a cellular metabolic reaction network, a reaction map can contain the entire reaction complement representing a global perspective. Alternatively, a reaction map can focus on a particular region of metabolism such as a region corresponding to a reaction subsystem described above or even on an individual pathway or reaction.
Methods disclosed herein can be used to determine the activity of a plurality of H. sapiens reactions including, for example, biosynthesis of an amino acid, degradation of an amino acid, biosynthesis of a purine, biosynthesis of a pyrimidine, biosynthesis of a lipid, metabolism of a fatty acid, biosynthesis of a cofactor, production of a hormone, production of a bioactive small molecule, transport of a metabolite and metabolism of an alternative carbon source.
Methods disclosed herein can be used to determine a phenotype of a H. sapiens mutant. The activity of one or more H. sapiens reactions can be determined using the methods described above, wherein the reaction network data structure lacks one or more gene-associated reactions that occur in H. sapiens. Alternatively, methods can be used to determine the activity of one or more H. sapiens reactions when a reaction that does not naturally occur in H. sapiens is added to the reaction network data structure. Deletion of a gene can also be represented in a model of the invention by constraining the flux through the reaction to zero, thereby allowing the reaction to remain within the data structure. Thus, simulations can be made to predict the effects of adding or removing genes to or from H. sapiens. The methods can be particularly useful for determining the effects of adding or deleting a gene that encodes for a gene product that performs a reaction in a peripheral metabolic pathway. In one contemplated embodiment, adenovirus vectors are used for in vivo transfer of genes determined in silico to be required for a desired functioning of the metabolic pathway.
A drug target or target for any other agent that affects human cellular function can be predicted using the methods of the invention. Such predictions can be made by removing a reaction to simulate total inhibition or prevention by a drug or agent. Alternatively, partial inhibition or reduction in the activity a particular reaction can be predicted by performing the methods with altered constraints. The methods can be particularly useful for identifying a target in a peripheral metabolic pathway.
Once a reaction has been identified for which activation or inhibition produces a desired effect on H. sapiens function, an enzyme or macromolecule that performs the reaction in H. sapiens or a gene that expresses the enzyme or macromolecule can be identified as a target for a drug or other agent. A candidate compound for a target identified by the methods of the invention can be isolated or synthesized using known methods. Such methods for isolating or synthesizing compounds can include, for example, rational design based on known properties of the target (see, for example, Decamp et al., Protein Engineering Principles and Practice, Ed. Cleland and Craik, Wiley-Liss, New York, pp. 467-506 (1996), which is hereby incorporated by reference in its entirety), screening the target against combinatorial libraries of compounds (see for example, Houghten et al., Nature, 354, 84-'86 (1991); Dooley et al., Science, 266, 2019-2022 (1994), which describe an iterative approach, or R. Houghten et al. PCT/US91/08694 and U.S. Pat. No. 5,556,762 which describe the positional scanning approach, all four of which are hereby incorporated by reference in their entireties), or a combination of both to obtain focused libraries. Those skilled in the art will know or will be able to routinely determine assay conditions to be used in a screen based on properties of the target or activity assays known in the art.
A candidate drug or agent, whether identified by methods described above or by other methods known in the art, can be validated using an in silico human metabolic model or method disclosed herein. The effect of a candidate drug or agent on human physiological function can be predicted based on the activity for a target in the presence of the candidate drug or agent measured in vitro or in vivo. This activity can be represented in an in silico human metabolic model by adding a reaction to the model, removing a reaction from the model or adjusting a constraint for a reaction in the model to reflect the measured effect of the candidate drug or agent on the activity of the reaction. By running a simulation under these conditions the holistic effect of the candidate drug or agent on human physiological function can be predicted.
Once a candidate has been identified, medical practitioners may readily conceive a plurality of delivery methods that will induce the desired simulated metabolic change. Suitable pharmaceutically acceptable carriers and their formulation are described in standard formulation treatises, e.g. Remington's Pharmaceuticals Sciences by E. W. martin. See also Wang, Y. J. and Hanson, M. A. “Parental formulations of Proteins and Pepties; Stability and Stabilizers,” Journals of Parental Sciences and Technology, Technical Report No. 10, Supp. 42:2 S (1988).
The methods of the invention can be used to determine the effects of one or more environmental components or conditions on an activity of a H. sapiens cell. As set forth above, an exchange reaction can be added to a reaction network data structure corresponding to uptake of an environmental component, release of a component to the environment, or other environmental demand. The effect of the environmental component or condition can be further investigated by running simulations with adjusted values for the metabolic flux vector of the exchange reaction target reaction to reflect a finite maximum or minimum flux value corresponding to the effect of the environmental component or condition. The environmental component can be, for example an alternative carbon source or a metabolite that when added to the environment of a H. sapiens cell can be taken up and metabolized. The environmental component can also be a combination of components present for example in a minimal medium composition. Thus, methods disclosed herein can be used to determine an optimal or minimal medium composition that is capable of supporting a particular activity of H. sapiens.
In some embodiments, a method for determining a set of environmental components to achieve a desired activity for H. sapiens is provided. The method includes the steps of providing a data structure relating a plurality of H. sapiens reactants to a plurality of H. sapiens reactions, wherein each of the H. sapiens reactions includes one or more reactants identified as a substrate of the reaction, one or more reactants identified as a product of the reaction and a stoichiometric coefficient relating the substrate and the product; providing a constraint set for the plurality of H. sapiens reactions; applying the constraint set to the data representation, thereby determining the activity of one or more H. sapiens reactions; determining the activity of one or more Homo sapiens reactions according to steps (a) through (c), wherein the constraint set includes an upper or lower bound on the amount of an environmental component, and repeating steps (a) through (c) with a changed constraint set, wherein the activity determined.
Further contemplated embodiments comprise down-regulating portions of metabolic networks by introducing inhibiting compounds regulating flows found in the model. These compounds may inhibit deleterious feedback systems, such as by suppressing antigen processing or antigen presentation. Other embodiments may up-regulate portions of metabolic networks to account for deletions or omissions in preceding stages of a “normally” functioning network. Such pathways need not be confined to H. sapiens but can be used in foreign entities as well. For example, metabolic pathways in viruses and bacteria can be modeled to determine an exogenous chemical dependency which will permit the virus' live introduction into a host, and subsequent controlled attenuation via reduction of the chemical dependency. Such live introduction facilitates better immune response. Embodiments of the present invention, by modeling the metabolic pathway, permit more controlled attenuation of these vectors, with less risk to the host.
One skilled in the art will readily recognize that these are but a few possible implementations of the invention and that many more exist.
The following examples are intended to illustrate but not limit the invention.
Example 1 A Stoichiometric Model of Glycolysis in Homo sapiensThe stoichiometric matrix as described in this invention was constructed for the glycolytic pathway in Homo sapiens.
A model of glycolysis was constructed by the methods described in this invention. Fluxes and concentrations for a normal human red blood cell are listed in Table 2. Table 2 also includes equilibrium constants for the reactions in the network. Reactions with equilibrium constants much greater than 1 or much less than one are effectively irreversible. The gradient matrix can be factored into the kappa and gamma matrices as described in this invention and illustrated for the glycolytic pathway in Homo sapiens.
By solving the steady state mass balance equations using the stoichiometric matrix in Example 1 and the mass-action rate law formulation and substituting in steady state metabolite concentrations (Table 2), steady state flux concentrations (Table 3), and equilibrium constants (Table 3), the forward rate constants were calculated.
The mass-action rate laws were then written in terms of the equilibrium constants and rate constants, as function of the metabolite concentrations. The gradient matrix was calculated from this.
The structural similarity between the stoichiometric matrix and negative of the transpose of the gradient matrix for glycolysis is immediately apparent by comparing Table 1 and Table 4.
As discussed above, the gradient matrix can be factored into the kappa and gamma matrices as illustrated in Table 5 and Table 6.
The elements in the kappa matrix are spread over approximately ten (log [κmax/κmin]=9.7) orders of magnitude, which conveys the biologically important time scale hierarchy and also conveys the ill-conditioned properties of the gradient matrix.
Example 3 The Dynamic Glycolytic Pathway in Homo sapiensUsing the stoichiometric matrix in Example 1 and the gradient matrix in Example 2, a dynamic model of the glycolytic pathway was constructed by calculating the concentration Jacobian matrix by post-multiplication of the stoichiometric matrix by the gradient matrix (Table 7) and the flux Jacobian matrix by pre-multiplication of the stoichiometric matrix by the gradient matrix (Table 8), as per the methods described above.
The ability to describe the Jacobian matrix in terms of concentration variables or flux variables, and consequently being able to describe the dynamics in terms of fluxes or concentrations, is a result of the duality property described herein.
Example 4 Dynamic Properties of the Glycolytic Pathway in Homo sapiensThe construction of the key dynamic data matrices constructed for the glycolytic pathway in Homo sapiens in Example 1, 2, and 3, now enable analysis of the dynamic properties of the pathway on different time scales.
Pooling of variables (metabolites or fluxes) refers to the formation of aggregate groups of variables, in which the group of variables move together in a concerted manner. Pools that form on very fast time scales reflect the formation of chemical equilibrium pools, whereas pooling that occurs on the slower time scales reflect physiologically relevant interactions. Moving from the very fast time scales to the slower ones (
The construction of the key dynamic data matrices constructed for a full model of metabolites in Homo sapiens.
Currently, the human red cell is the only cell-scale kinetic model available, whose formulation followed a 30 year history of iterative model building (summarized in (Joshi and Palsson, 1989; Joshi and Palsson, 1990), which are both hereby incorporated by reference in their entireties). Analysis of the dynamic structure of this model resulted in the simplification of network dynamics and the description of the cellular functions in terms of physiologically meaningful pools of metabolites. Analysis of the hierarchical dynamics of the human red cell model resulted in a richer and more complex physiological pool formation (Kauffman et al., 2002, which is hereby incorporated by reference) that detailed above for glycolysis alone. An overview of the hierarchical reduction of the network into a functional diagram is schematized in
This procedure characterized the dynamics of a network about a particular steady state. The first step involves specifying a particular steady state flux distribution. Next, the concentrations (generally there will need to be some approximations of these) of the metabolites at this steady state need to be identified. Third, the equilibrium constant approximations need to be determined and the steady state fluxes and concentrations will be incorporated into mass action rate laws. Finally, the forward rate constants can be solved by solving n linear equations for forward rate constants.
For situations in which all of the concentrations are not known or if one wants to trace the effects of particular metabolites in multiple rate constants, rather than n linear equations, one can solve the m equations by solving the steady state mass conservation relationship (S.v˜0). This will result in multiple alternative solutions.
As described previously (5), the dynamics of this model are consistent with the observations from previous studies in which the well-known examples of pool formation between hexose phosphates (HP) and phosphoglycerates (PG). With successive removal of time scales (modes), more and more of the metabolites begin to form aggregate variables and move together in a concerted fashion at fixed ratios. For glycolysis alone, the successive aggregation of chemical moieties (i.e. HP, PG) culminates in, on the slowest time scale, the formation of a physiologically meaningful pool that represents the sum of high-energy phosphate bonds found in the glycolytic intermediates (i.e. their ATP equivalents). The last row of the modal matrix for the flux Jacobian (Supplemental Material) shows that this pool is moved by hexokinase as the input and ATPase as the output, reflecting the catabolism of glucose to generate ATP. The dimensions of the matrices and their corresponding subspaces are shown in Table 11. This is a single pathway model, hence the single dimension if the right null space is expected.
Scaleability of the Approach
Since the approach outlined in
A map of the network is shown in
As expected from the size of the left null space, there are three conserved pools of metabolites, the stoichiometric sum of all of the phosphate containing compounds, the NAD moiety, and the NADP moiety. The amount of glutathione in the network is also conserved, however since the dimeric form is not included as a dynamic variable in the network, it does not appear in the left null space. The size of the right null space increased significantly, by 8 dimensions, reflecting the various alternate pathways that can occur from the Rappoport-Leubering Shunt, the pentose-phosphate pathway, and the connections between the non-oxidative branches of the pentose phosphate pathway and the nucleotide salvage pathways.
The gradient and Jacobian matrices have the characteristics leading towards the expected behavior for the metabolic network, such as hierarchical structure with fast reactions corresponding to isomerization reactions and slower ones reflecting the more physiological behaviors of the network, such as phosphate group transfers. There is good time scale separation (about 20 milliseconds to 5 hours), although the slowest and fastest modes of the traditional red cell model is slightly larger (16, 19). These slower time scales in the kinetic model resulted from the potassium leak channels (20, 21), which were not included in this model. The compositions of the very fast motions are similar to those seen for the red cell (e.g. isomerization between glycolytic intermediates, hexose sugars, pentose sugars, etc). With this implementation, although the diphosphoglycerates move on the slower modes in the network, they are still faster than expected (on the order of 5 minutes as opposed to 10 hours).
Construction of Mechanistic Regulated Model of the Human Erythrocyte
Fundamentally, all interactions between metabolites and enzymes are bilinear, including regulatory reactions (22). The mass action model for the red blood cell was expanded in order in order to include an increased level of detail for the biochemical reactions, accounting for mechanistic regulation in the key regulatory enzymes of the red cell. Many of the association and dissociation constants for the steps involved in enzymatic catalysis for the enzymes in the red cell have been characterized. The five key regulatory enzymes in the red blood cell were described in terms of their mechanisms, these include hexokinase (HK), phosphofructokinase (PFK), diphosphoglyceromutate and diphosphoglycerate phosphatase (DPGM and DPGase carried out by the same enzyme), glucose-6-phosphate dehydrogenase (G6PDH), and adenosine kinase (AK). The mass action model discussed in the previous section was used as a basis for this model. The incorporation of the bilinear interactions reflecting small metabolite regulation of the enzymes resulted in a dramatic increase in the size of the network. Given the complexity and increase in size, the addition of each regulated enzyme was carried out in a systematic, ‘modular’ manner. Thus the parameters for each enzyme were calculated individually. After each individual enzymatic step was completed and test, they were integrated into the whole network model. Since a particular flux state was under consideration, there were no problems with alternative flux distributions in the network.
The mechanism for unordered mechanism for hexokinase along with the association and dissociation constants were used as described in (23).
PFK is a tetramer and as is apparent from the reaction diagram in
The diphosphoglyceromutase and diphosphoglycerol phosphatase reactions are both carried out by the same enzyme and the enzymatic steps were described according to (23). As with HK, the majority of the association and dissociations steps were measured (23). A similar procedure was carried out in order to calculate the concentrations of the enzyme-metabolite intermediates.
The reactions for G6PDH were implemented as described by (23) and solved in an analogous manner as HK. The reaction scheme for AK is outlined in
A summary of the dimensions and sizes of the subspaces for the regulated model appear in Table 13.
Note that this network is more than twice as large as the previous version. Enzyme catalyzed mechanisms were introduced in this model for 5 enzymes and the size of the left null space increased from 3 to 10. The size of the left null space grew by 7 dimensions due to the addition enzymes as reactants; 5 enzymes (note that PFK has an inactive and active form) plus the conservation of total magnesium. The right null space increased by 3 as a result of the non-sequential mechanism of some of the regulated enzymes, such as HK.
Modal matrices for larger networks are much more difficult to interpret and were dependent on involved algorithmic approaches to elucidate the structure (2). As alluded to earlier, this is because as the size of networks increase and the non-integer entries in G which tends to be an illconditioned matrix, interactions become more tightly tied together. However, observing the structure of the modal matrix for the regulated network (see Supplemental Material) a few points can be observed,
The time scale separation has improved from the non-regulated model, with motions ranging from less than microseconds to ten hours.
The slow modes are dominated by the metabolite bound enzymes.
In particular, the reactions that have a larger number of intermediate complexes playa more dominant role, especially at the slower time scales.
The movement of 2,3DPG (enzyme bound forms) have moved to the slowest time scales.
In short, the structure of the modal matrix more closely approximates the features in the full red cell kinetic model. This has encouraging and interesting implications, since the starting point was a basic model of glycolysis, expanded it to a model of red cell metabolism, and then incorporated regulatory interactions in a mechanistic manner. The results from the regulated model of red cell metabolism suggest that one manner in which regulatory enzymes exert their influence over the rest of the network is their ability to bind more than one substrate and to exist in multiple states (active and inactive).
Simulations
Simulations were carried out with the regulated cell scale model in response to instantaneous ATP depletion (
The NADHP/NADP ratio, G6PD enzyme ratio, and flux through the G6PD reaction are all highly correlated under both conditions (energy and redox stresses), suggesting that either of these measurements can be used to infer the states of the others. Correlations among many of the other variables appear to vary under different conditions.
ConclusionsStoichiometric genome-scale models have made significant advancements in recent years (26). The incorporation of dynamic information and the ability to describe the kinetics of these networks is an important, physiologically relevant aspect that needs to be accounted for. The kcone formalism took the initial steps towards developing kinetic models in a data driven fashion through sampling the kinetic parameter space (4). Herein we described an approach which uses stoichiometric models as a scaffold for constructing kinetic models based upon bilinear kinetics, through the incorporation of metabolomic data. This stepwise procedure is scaleable and feasible for networks of increasing size.
In order to further expand the mechanistic detail of stoichiometric models, we demonstrated how complex regulatory schemes, which are described with more fidelity as a scheme of bilinear interactions, can fit within this framework of kinetic model construction. Subsequent analysis of the network suggested that the regulatory influence of enzymes is determined by the number of metabolites that can bind to the enzyme. We have previously shown that stoichiometry is a key aspect of kinetic networks (5); here we further highlight that it is also a key factor in regulatory control as well.
Accounting for regulation allows one to observe the sub-networks of regulation within regulated metabolic networks. These subnetworks can exert influences across the whole network, and this appears at least in part in their ability to bind multiple compounds. As with small metabolites, the ratios of enzymes in different forms (e.g. bound versus unbound) reflect different functional states. For some enzymes these are highly correlated with the flux and/or related small metabolite ratios, as with glucose-6-phosphate dehydrogenase. However for other enzymes, particularly those with a larger number of different binding states, such ratios will not always be correlated with fluxes.
Increasing availability of metabolomic data in conjunction with the availability of Keq approximations and related methods for extrapolation and approximation of these constants (27-29) will enable the application of this approach to more metabolic networks. During the preparation of this work, Tran et al (10) describe an approach with similarities to the scheme described herein in which the solution space is sampled for a glycolysis/pentose phosphate pathway model in E. coli. The addition of sampling would enable to the exploration of alternative solutions and highlights an area worth further exploration in the future. Taken together, the approach described within this work is scale-able and systematic, thus representing a feasible, realistic path towards genome-scale kinetic descriptions of metabolic networks.
Methods
In general for any reaction, k,+
2Ak1−←→k1+B
the resulting net rate law is given by,
v]=k1+A2−k1−B (1),
with the equilibrium constant,
Keq=k1+/k1− (2)
For the regulated reactions in the network, the enzymes are explicit species and combined with each metabolite, as outlined in the figures.
Steps for Model Construction—Modular Model Construction
Substituting flux rate expressions (1) and equilibrium constants (2), into the steady state mass balances,
S·v=0 (3)
for which S is the m×n stoichiometric matrix and v is the flux vector of length n, allows one to solve for the n forward rate constants k. With the rate constants and equilibrium constants in hand in addition to steady state concentrations, the gradient matrix can be defined. The gradient matrix G is given by dv/dx. Decomposition into the dual Jacobian matrices (5) is given by
Jx−S·G=MΛM′; (4)
mx=Mx−1x′ (5)
for the concentration Jacobian, Jx, and
Jv=G·S=MΛM;′ (6)
mv=Mv−1v′ (7)
for the flux Jacobian, Jv. These calculations were carried out for each of the sub-models that were developed in the studies described in the manuscript. The regulated cell model also included magnesium complexing with ATP, ADP, AMP, and 2,3DPG. The matrices S and G are represented for various studies in
Models were constructed and tested in Mathematica (Wolfram Research). Simulations were also all carried out in Mathematica. However, model construction and simulation can be carried out in any number of other commercial software packages, such as Matlab, free software packages such as Scilab and R, or written in programming and scripting languages such as Java, C++, python, perl, etc. Since the protein concentrations were much smaller than the concentrations of the metabolites, simulations were carried out using normalized enzyme concentrations, in order to reduce the stiffness of the system of differential equations.
When carrying out simulations, the enzyme and bound enzyme variables were normalized, by redefining the enzyme rate constant. This allowed the stiffness of the ordinary differential equations to be significantly reduced. The system was then integrated out to >1000 hours to define the steady state.
This example shows how physiologically meaningful dynamic structures form as a result of the particular numerical values in G and how they overlay on the network structure given by S. Not all sets of numerical values give this dynamic structure. It has been shown that genetic variation, as represented by sequence polymorphrism in particular enzymes (pyruvate kinase and glucose-6-phosphate dehydrogenase), can disrupt this dynamic structure and lead to pathological states (Jamshidi et al., 2002).
Example 6 Simulating Genetic Deficiencies in Homo sapiens Through Altered Metabolite MeasurementsConsidering the dynamic model constructed in the above examples for metabolism in the human erythrocyte, one can simulate genetic deficiencies through measurements of metabolites. For example a model for a cell with a deficiency in Phosphofructokinase can be constructed and simulated based on metabolite measurements. Following the algorithm for model construction, a new set of forward rate constants will be calculated, resulting in a different set of reaction rate expressions and gradient matrix. Subsequent simulations will reflect the resulting altered biochemical phenotype due to the genetic mutation(s) resulting in Phosphofructokinase Deficiency.
Example 7 Simulating Genetic Deficiencies in Homo sapiens Through Altered Rate ConstantsConsidering the dynamic model constructed in the above examples for metabolism in the human erythrocyte, one can simulate genetic deficiencies through measurements of kinetic parameters, such as rate constants. For example, a model for a cell with a deficiency in Phosphofructokinase can be constructed and simulated based on measurements of altered kinetic parameters for the individual regulated steps. The mechanistic representation of the regulation of PFK by other metabolites enables the representation of genetic mutations in a direct, detailed manner. Following the algorithm for model construction, the altered rate constants are adjusted and an altered set of concentrations are recalculated through solving the steady state equations and simulating the time course of the dynamics. This will result in a new set of reaction rate expressions and gradient matrix. Subsequent simulations will reflect the resulting altered biochemical phenotype due to the genetic mutation(s) resulting in Phosphofructokinase Deficiency.
Example 8 Simulating Genetic Deficiencies in Homo sapiens Through Altered Lumped Parameter MeasurementsConsidering the dynamic model constructed in the above examples for metabolism in the human erythrocyte, one can simulate genetic deficiencies through measurements of kinetic parameters, such as rate constants. For example, a model for a cell with a deficiency in Phosphofructokinase can be constructed and simulated based on measurements of kinetic parameters such as Vmax. Following the algorithm for model construction, the steady state equations and dynamic simulations are re-implemented. This will result in a new set of reaction rate expressions and gradient matrix. Subsequent simulations will reflect the resulting altered biochemical phenotype due to the genetic mutation(s) resulting in Phosphofructokinase Deficiency.
Example 9 Simulating the Effects of Drugs with Known Targets in Homo sapiensConsidering the dynamic model constructed in the above examples for metabolism in the human erythrocyte, one can simulate the effects of known drugs on metabolism. For example, the effect of a redox depleting drug which is known to act on a particular enzyme within the network, such as glucose-6-phosphate dehydrogenase, can be simulated by reducing or abrogating the rate constant for glucose-6-phosphate dehydrogenase by 10%, 20%, 40%, 60%, 80%, or 100%. Following the algorithm for model construction, the altered rate constants are adjusted and an altered set of concentrations are recalculated through solving the steady state equations and simulating the time course of the dynamics. This will result in a new set of reaction rate expressions and gradient matrix. The resulting model will reflect a metabolic network treated with the drug at the specified dose.
Example 10 Simulating the Effects of Drugs with Unknown Targets in Homo sapiensConsidering the dynamic model constructed in the above examples for metabolism in the human erythrocyte, one can simulate the effects of known drugs on metabolism and identify potential candidates for the identity of the drug target. For example, metabolite concentration measurements can be made experimentally following the administration of a drug with unknown identity or unknown target. Following the algorithm for model construction, a new set of forward rate constants will be calculated, resulting in a different set of reaction rate expressions and gradient matrix. The resulting model will reflect a metabolic network treated with the drug at the specified dose. Analysis of the Jacobian matrices and dynamic simulations can yield insights into the mechanism of action of the drug and its potential candidate targets.
REFERENCESThe following references are incorporated herein in their entirety:
- 1. Scriver, C. 2000. The Metabolic and Molecular Bases of Inherited Disease. McGraw-Hill Professional.
- 2. Jamshidi, N., and B. O. Palsson. 2008. Top-down analysis of temporal hierarchy in biochemical reaction networks. PLoS Comput BioI4:eI000177.
- 3. Bruggeman, F. J., J. de Haan, H. Hardin, J. Bouwman, S. Rossell, K. van Eunen, B. M. Bakker, and H. V. Westerhoff. 2006. Time-dependent hierarchical regulation analysis: deciphering cellular adaptation. Systems biology 153:318-322.
- 4. Famili, I., R. Mahadevan, and B. O. Palsson. 2005. k-Cone analysis: determining all candidate values for kinetic parameters on a network scale. Biophys J 88:1616-1625.
- 5. Jamshidi, N., and B. O. Palsson. 2008. Formulating genome-scale kinetic models in the post-genome era. Molecular systems biology 4: 171.
- 6. Maurya, M. R., and S. Subramaniam. 2007. A kinetic model for calcium dynamics in RAW 264.7 cells: I. Mechanisms, parameters, and subpopulational variability. Biophys J 93:709-728.
- 7. Saucerman, J. J., L. L. Brunton, A. P. Michailova, and A. D. McCulloch. 2003. Modeling beta-adrenergic control of cardiac myocyte contractility in silico. The Journal of biological chemistry 278:47997-48003.
- 8. Smallbone, K., E. Simeonidis, D. S. Broomhead, and D. B. Kell. 2007. Something from nothing bridging the gap between constraint-based and kinetic modelling FEBS J 274:5576-5585.
- 9. Steuer, R., T. Gross, J. Selbig, and B. Blasius. 2006. Structural kinetic modeling of metabolic networks. Proc Nat! Acad Sci USA 103: 11868-11873.
- 10. Tran, L. M., M. L. Rizk, and J. C. Liao. 2008. Ensemble modeling of metabolic networks. Biophys J 95:5606-5617.
- 11. Wishart, D. S., D. Tzur, C. Knox, R. Eisner, A. C. Guo, N. Young, D. Cheng, K. Jewell, D. Arndt, S. Sawhney, C. Fung, L. Nikolai, M. Lewis, M. A. Coutouly, I. Forsythe, P. Tang, S. Shrivastava, K. Jeroncic, P. Stothard, G. Amegbey, D. Block, D. D. Hau, J. Wagner, J. Miniaci, M. Clements, M. Gebremedhin, N. Guo, Y. Zhang, G. E. Duggan, G. D. Macinnis, A. M. Weljie, R. Dowlatabadi, F. Bamforth, D. Clive, R. Greiner, L. Li, T. Marrie, B. D. Sykes, H. J. Vogel, and L. Querengesser. 2007. HMDB: the Human Metabolome Database. Nucleic acids research 35:D521-526.
- 12. Bennett, B. D., J. Yuan, E. H. Kimball, and J. D. Rabinowitz. 2008. Absolute quantitation of intracellular metabolite concentrations by an isotope ratio-based approach. Nat Protoc 3:1299-1311
- 13. Claus, S. P., T. M. Tsang, Y. Wang, O. Cloarec, E. Skordi, F. P. Martin, S. Rezzi, A. Ross, S. Kochhar, E. Holmes, and J. K. Nicholson. 2008. Systemic multicompartmental effects of the gut microbiome on mouse metabolic phenotypes. Molecular systems biology 4:219.
- 14. Ishii, N., K. Nakahigashi, T. Baba, M. Robert, T. Soga, A. Kanai, T. Hirasawa, M. Naba, K. Hirai, A. Hogue, P. Y. Ho, Y. Kakazu, K. Sugawara, S. Igarashi, S. Harada, T. Masuda, N. Sugiyama, T. Togashi, M. Hasegawa, Y. Takai, K. Yugi, K. Arakawa, N. Iwata, Y. Toya, Y. Nakayama, T. Nishioka, K. Shimizu, H. Mori, and M. Tomita. 2007. Multiple high-throughput analyses monitor the response of E. coli to perturbations. Science (New York, N.Y. 316:593-597.
- 15. Smith, C. A., G. O'Maille, E. J. Want, C. Qin, S. A. Trauger, T. R. Brandon, D. E. Custodio, R. Abagyan, and G. Siuzdak. 2005. METLIN: a metabolite mass spectral database. Ther Drug Moni t 27: 747-7 51.
- 16. Joshi, A., and B. O. Palsson. 1990. Metabolic dynamics in the human red cell. Part IV-Data prediction and some model computations. J Theor Bioi 142:69-85.
- 17. Joshi, A., and B. O. Palsson. 1990. Metabolic dynamics in the human red cell. Part III—Metabolic reaction rates. J Theor Bioi 142:41-68.
- 18. Reich, J., and E. Selkov. 1981. Energy metabolism of the cell: a theoretical treatise. Academic Press.
- 19. Jamshidi, N., J. S. Edwards, T. Fahland, G. M. Church, and B. O. Palsson. 2001. Dynamic simulation of the human red blood cell metabolic network. Bioinformatics 17:286-287.
- 20. Jamshidi, N., and B. O. Palsson. 2006. Systems biology of the human red blood cell. Blood cells, molecules & diseases 36:239-247.
- 21. Kauffman, K. J., J. D. Pajerowski, N. Jamshidi, B. O. Palsson, and J. S. Edwards. 2002. Description and analysis of metabolic connectivity and dynamics in the human red blood cell. Biophys J 83:646-662.
- 22. Segel, 1. 1975. Enzyme Kinetics. John Wiley & Sons.
- 23. Mulquiney, P. J., and P. W. Kuchel. 1999. Model of 2,3-bisphosphoglycerate metabolism in the human erythrocyte based on detailed enzyme kinetic equations: equations and parameter refinement. Biochem J 342 Pt 3:581-596.
- 24. Albe, K. R., M. H. Butler, and B. E. Wright. 1990. Cellular concentrations of enzymes and their substrates. J Theor Bioi 143:163-195.
- 25. Hawkins, C. F., and A. S. Bagnara. 1987. Adenosine kinase from human erytluocytes: kinetic studies and characterization of adenosine binding sites. Biochemistry 26: 1982-1987.
- 26. Feist, A. M., M. J. Herrgard, 1. Thiele, J. L. Reed, and B. O. Palsson. 2009. Reconstruction of biochemical networks in microorganisms. Nat Rev Microbiol 7: 129-143.
- 27. Alberty, R. A. 2003. Thermodynamics of biochemical reactions. Wiley-Interscience, Hoboken, N.J.
- 28. Jankowski, M. D., C. S. Henry, L. J. Broadbelt, and V. Hatzimanikatis. 2008. Group contribution method for thermodynamic analysis of complex metabolic networks. Biophys J 95:1487-1499.
- 29. Mavrovouniotis, M. L. 1991. Estimation of standard Gibbs energy changes of biotransformations. The Journal of biological chemistry 266:14440-14445.
- 30. Mahadevan, R., and C. H. Schilling 2003. The effects of alternate optimal solutions in constraint-based genome-scale metabolic models. Metabolic engineering 5:264-276.
- 31. Mulquiney, P. J., W. A. Bubb, and P. W. Kuchel. 1999. Model of 2,3-bisphosphoglycerate metabolism in the human erythrocyte based on detailed enzyme kinetic equations: in vivo kinetic characterization of 2,3-bisphosphoglycerate synthase/phosphatase using 13C and 31P NMR. Biochem J 342 Pt 3:567-580.
Claims
1. A computer-implemented method for developing a dynamic network model of a biological network, the method comprising:
- providing, at a computer system, a network data structure relating a plurality of reactants to a plurality of reactions, the plurality of reactions comprising one or more flux distributions, wherein each of the reactions comprises one or more reactants identified as a substrate of the reaction, one or more reactants identified as a product of the reaction and a stoichiometric coefficient relating the substrate and the product;
- providing, at a computer system, a constraint set for the plurality of reactions in said biological network associated with the plurality of reactions;
- providing, at a computer system, a constraint set for a plurality of constraints, wherein said constraint set comprises at least one physical constraint;
- providing, at a computer system, a constraint set for a plurality of kinetic constants associated with the plurality of reactions;
- providing, at a computer system, a constraint set for a plurality of equilibrium constants associated with the plurality of reactions;
- integrating, at a computer system, the sets of constraints to form stoichiometric, gradient, and Jacobian matrices;
- providing, at a computer system, an objective function; and
- predicting, at a computer system, a physiological function related to said biological network by determining using the computer system at least one flux distribution, from the one or more flux distributions, that minimizes or maximizes the objective function when at least one constraint set from the sets of constraints is applied to the data structure.
2. The method of claim 1, wherein said at least one physical constraint is mass balance.
3. The method of claim 1, wherein said constraint set for said plurality of reactions in said biological network reflects enzymatic interconversions of small metabolites, enzyme reaction schemes, allosteric enzyme interactions, signaling pathways, or macromolecular biochemical reactions.
4. A computer-implemented method for developing a dynamic network model, the method comprising:
- accessing archetype data for a plurality of reactions of a network;
- developing a stoichiometric data structure for the reactions based on the archetype data;
- accessing physiological data reflecting a specific context;
- performing a steady state simulation using the stoichiometric data structure and the physiological data;
- incorporating metabolic data and equilibrium constant data based on the steady state simulation;
- deriving a gradient matrix and dynamic rate equations based on the stoichiometric data structure and the incorporated data; and
- deriving a concentration Jacobian matrix and a flux Jacobian matrix based on the stoichiometric data structure and the gradient matrix so as to characterize dynamic properties of the network.
5. The method of claim 4, wherein an ith column of the stoichiometric data structure contains stoichiometric coefficients of an ith reaction in the network and an ith row in the gradient matrix contains partial derivatives of the ith reaction with respect to all reactants in the network.
6. The method of claim 4, wherein the stoichiometric data structure is derived from genomic, proteomic and metabolomic data.
7. The method of claim 4, wherein the gradient matrix is derived from fluxomic and metabolomic data, kinetic characterization of individual reactions and assessment of thermodynamic properties.
8. The method of claim 4, wherein the stoichiometric data structure is derived from a stoichiometry of reactions, and the gradient matrix is based on kinetic data and thermodynamic information.
9. The method of claim 4, wherein the stoichiometric data structure is based on the content of a genome and is a property of a species, and the gradient matrix is based on kinetic information and is genetically derived.
10. The method of claim 4, wherein a reaction comprises a representation of a chemical conversion that consumes a substrate or produces a product, and a reactant comprises a representation of a chemical that is a substrate or a product of a reaction that occurs in a cell.
11. The method of claim 4, wherein a reaction comprises a conversion that occurs due to the activity of one or more enzymes that are genetically encoded by the human genome, a conversion that occurs spontaneously in a human or mammalian cell, or conversions based on changes in chemical composition such as those due to nucleophilic or electrophilic addition, nucleophilic or electrophilic substitution, elimination, isomerization, deamination, phosphorylation, methylation, glycolysation, reduction, oxidation or changes in location such as those that occur due to a transport reaction that moves one or more reactants within the same compartment or from one cellular compartment to another.
12. The method of claim 4, wherein the biological network reflects a H. sapiens cell type at any stage of differentiation.
13. The method of claim 4, wherein the method is applied to normal cells or pathological cells.
14. The method of claim 4, wherein the concentration Jacobian matrix is calculated by post-multiplication of the stoichiometric data structure by the gradient matrix, and the flux Jacobian matrix is calculated by pre-multiplication of the stoichiometric data structure by the gradient matrix.
15. A computer-implemented method for developing data-driven dynamic models of biological networks, the method comprising:
- providing, at a computer system, a data structure relating a plurality of biological network reactants to a plurality of biological network reactions, the plurality of biological network reactions comprising one or more flux distributions, wherein each of the biological network reactions comprises one or more biological network reactants identified as a substrate of the reaction, one or more biological network reactants identified as a product of the biological network reaction and a stoichiometric coefficient relating the substrate and the product;
- providing, at a computer system, a constraint set for the plurality of biological network reactions;
- providing, at a computer system, a constraint set for a plurality of concentrations associated with the plurality of biological network reactions;
- providing, at a computer system, a constraint set for a plurality of kinetic constants associated with the plurality of biological network reactions;
- providing, at a computer system, a constraint set for a plurality of equilibrium constants associated with the plurality of biological network reactions;
- integrating, at a computer system, the sets of constraints, thereby forming stoichiometric, gradient, and Jacobian matrices;
- providing, at a computer system, an objective function; and
- predicting, at a computer system, a physiological function related to at least one of the biological networks by determining using the computer system at least one flux distribution, from the one or more flux distributions, that minimizes or maximizes the objective function when at least one constraint set from the sets of constraints is applied to the data structure, wherein one or more of said constraint sets reflects a perturbation.
16. The method of claim 15, wherein one or more of said constraint sets has been altered to reflect genetic or environmental perturbations.
17. The method of claim 15, wherein at least one of said kinetic constants or at least one of said concentrations has been altered to reflect a genetic perturbation.
18. The method of claim 15, additionally comprising modifying, at a computer system, at least one of the reaction constraint set, the concentration constraint set, the kinetic constant set and the equilibrium constant set based at least partly on genetic variations.
19. The method of claim 15, additionally comprising modifying, at a computer system, at least one of the reaction constraint set, the concentration constraint set, the constraint set for a plurality of kinetic constants and the set for a plurality of equilibrium constants based at least partly on altered environmental conditions.
20. The method of claim 15, additionally comprising simulating the effects of adding or removing genes on at least one of the plurality of biological networks.
21. The method of claim 15, additionally comprising simulating the effects of one or more environmental perturbations on at least one of the plurality of biological networks.
22. The method of claim 18, additionally comprising determining dynamics and consequences of genetic variations selected from the group consisting of deficiencies in metabolic enzymes and metabolic transporters.
23. The method of claim 22, wherein the genetic variations comprise functioning of a compound selected from the group consisting of phosphofructokinase, phosphoglycerate kinase, phosphoglycerate mutase, lactate dehydrogenase adenosine deaminase, ABC transporters, the SLC class of transporters, and the cytochrome P450 class of enzymes.
24. The method of claim 15, wherein providing a constraint set for a plurality of concentrations comprises measuring metabolite levels at steady state.
25. The method of claim 15, wherein providing a constraint set for a plurality of kinetic constants comprises solving a steady state mass conservation relationship comprising a stoichiometric matrix and flux vector.
26. The method of claim 19, wherein the modified sets of constraints correspond to an individual.
27. The method of claim 19, wherein the modified sets of constraints correspond to a group or plurality of individuals.
28. The method of claim 18, wherein a therapeutic regime is applied to regulate a physiological function based on at least one of the modified sets based at least partly on genetic variations.
29. The method of claim 15, wherein a genetic variation results in alterations of kinetic parameters or concentrations of variables.
30. The method of claim 19, wherein a therapeutic regime is applied to regulate a physiological function based on at least one of the modified sets based at least partly on altered environmental conditions.
31. The method of claim 15, wherein the altered environmental conditions result in changes in measured variables.
32. The method of claim 15, additionally comprising:
- factorizing the gradient matrix into a kappa matrix and a gamma matrix; and
- analyzing the gradient, Jacobian, kappa and gamma matrices to characterize time scales of the network model.
33. The method of claim 32, wherein the kappa matrix is a diagonal matrix of kinetic constants.
34. The method of claim 15, additionally comprising:
- modifying, at a computer system, at least one of metabolite measurement or a kinetic parameter measurement based on a genetic influence, the kinetic parameter comprising a rate constant;
- generating a new set of reaction rate expressions and a new gradient matrix; and
- analyzing a pathological state of the network model based on a simulation of the new set of reaction rate expressions and the new gradient matrix.
35. A computer-implemented method for the analysis of pathological states due to genetic influences, the method comprising:
- providing, at a computer system, a data structure relating a plurality of reactants to a plurality of reactions, the plurality of reactions comprising one or more flux distributions, wherein each of the reactions comprises one or more reactants identified as a substrate of the reaction, one or more reactants identified as a product of the reaction and a stoichiometric coefficient relating the substrate and the product;
- providing, at a computer system, a constraint set for the plurality of reactions;
- providing, at a computer system, a constraint set for a plurality of concentrations associated with the plurality of reactions;
- providing, at a computer system, a constraint set for a plurality of kinetic constants associated with the plurality of reactions;
- providing, at a computer system, a constraint set for a plurality of equilibrium constants associated with the plurality of reactions;
- integrating, at a computer system, the sets of constraints, thereby forming stoichiometric, gradient, and Jacobian matrices;
- providing, at a computer system, an objective function;
- modifying, at a computer system, at least one of the reaction constraint set, the concentration constraint set, the kinetic constant set and the equilibrium constant set based at least partly on genetic variations; and
- predicting, using a computer system, a physiological function related to a pathological state by determining at least one flux distribution, from the one or more flux distributions, that minimizes or maximizes the objective function when at least one constraint set from the sets of constraints is applied to the data structure.
36. The method of claim 35, additionally comprising determining dynamics and consequences of genetic variations selected from the group consisting of deficiencies in metabolic enzymes and metabolic transporters.
37. The method of claim 36, wherein the genetic variations comprise functioning of a compound selected from the group consisting of phosphofructokinase, phosphoglycerate kinase, phosphoglycerate mutase, lactate dehydrogenase adenosine deaminase, ABC transporters, the SLC class of transporters, and the cytochrome P450 class of enzymes.
38. The method of claim 35, wherein providing a constraint set for a plurality of concentrations comprises measuring metabolite levels at steady state.
39. The method of claim 35, wherein providing a constraint set for a plurality of kinetic constants comprises solving a steady state mass conservation relationship comprising a stoichiometric matrix and flux vector.
40. The method of claim 35, wherein the modified sets of constraints correspond to an individual.
41. The method of claim 35, wherein the modified sets of constraints correspond to a group or plurality of individuals.
42. The method of claim 35, wherein a therapeutic regime is applied to regulate a physiological function based on at least one of the modified sets based at least partly on genetic variations.
43. The method of claim 35, wherein the genetic variations result in alterations of kinetic parameters or concentrations of variables.
44. A computer-implemented method for the analysis of pathological states due to environmental influences, the method comprising:
- providing, at a computer system, a data structure relating a plurality of reactants to a plurality of reactions, the plurality of reactions comprising one or more flux distributions, wherein each of the reactions comprises one or more reactants identified as a substrate of the reaction, one or more reactants identified as a product of the reaction and a stoichiometric coefficient relating the substrate and the product;
- providing, at a computer system, a constraint set for the plurality of reactions;
- providing, at a computer system, a constraint set for a plurality of concentrations associated with the plurality of reactions;
- providing, at a computer system, a constraint set for a plurality of kinetic constants associated with the plurality of reactions;
- providing, at a computer system, a constraint set for a plurality of equilibrium constants associated with the plurality of reactions;
- integrating, at a computer system, the sets of constraints, thereby forming stoichiometric, gradient, and Jacobian matrices;
- providing, at a computer system, an objective function;
- modifying, at a computer system, at least one of the reaction constraint set, the concentration constraint set, the constraint set for a plurality of kinetic constants and the set for a plurality of equilibrium constants based at least partly on altered environmental conditions; and
- predicting, using a computer system, a physiological function related to a pathological state by determining at least one flux distribution, from the one or more flux distributions, that minimizes or maximizes the objective function when at least one constraint set from the sets of constraints is applied to the data structure.
45. The method of claim 44, wherein providing a constraint set for a plurality of concentrations comprises measuring metabolite levels at steady state.
46. The method of claim 44, wherein the modified sets of constraints correspond to an individual.
47. The method of claim 44, wherein the modified sets of constraints correspond to a group or plurality of individuals.
48. The method of claim 44, wherein a therapeutic regime is applied to regulate a physiological function based on at least one of the modified sets based at least partly on altered environmental conditions.
49. The method of claim 44, wherein the altered environmental conditions result in changes in measured variables.
50. A computer-implemented method for the analysis of a pathological state of a data-driven model of a biological network due to a genetic influence, the method comprising:
- providing, at a computer system, a data structure relating a plurality of reactants to a plurality of reactions, the plurality of reactions comprising one or more flux distributions, wherein each of the reactions comprises one or more reactants identified as a substrate of the reaction, one or more reactants identified as a product of the reaction and a stoichiometric coefficient relating the substrate and the product;
- providing, at a computer system, a constraint set for the plurality of reactions;
- providing, at a computer system, a constraint set for a plurality of concentrations associated with the plurality of reactions;
- providing, at a computer system, a constraint set for a plurality of kinetic constants associated with the plurality of reactions;
- providing, at a computer system, a constraint set for a plurality of equilibrium constants associated with the plurality of reactions;
- integrating, at a computer system, the sets of constraints, thereby forming stoichiometric, gradient, and Jacobian matrices;
- modifying, at a computer system, at least one of metabolite measurement or a kinetic parameter measurement based on a genetic influence, the kinetic parameter comprising a rate constant;
- generating a new set of reaction rate expressions and a new gradient matrix; and
- analyzing of the pathological state of the network model based on a simulation of the new set of reaction rate expressions and the new gradient matrix.
51. A computer-implemented method for the analysis of an effect of a drug or chemical agent on a physiological function, the method comprising:
- providing, at a computer system, a data structure relating a plurality of reactants to a plurality of reactions, the plurality of reactions comprising one or more flux distributions, wherein each of the reactions comprises one or more reactants identified as a substrate of the reaction, one or more reactants identified as a product of the reaction and a stoichiometric coefficient relating the substrate and the product;
- providing, at a computer system, a constraint set for the plurality of reactions;
- providing, at a computer system, a constraint set for a plurality of concentrations associated with the plurality of reactions;
- providing, at a computer system, a constraint set for a plurality of kinetic constants associated with the plurality of reactions;
- providing, at a computer system, a constraint set for a plurality of equilibrium constants associated with the plurality of reactions;
- integrating, at a computer system, the sets of constraints, thereby forming the stoichiometric, gradient, and Jacobian matrices;
- providing, at a computer system, an objective function;
- adding a reaction to the data structure, removing a reaction from the data structure, or adjusting a constraint for a reaction in at least one of the sets of constraints to reflect a measured effect of a drug or chemical agent on the activity of the reaction; and
- predicting, at a computer system, a physiological function related to the drug or chemical agent by determining at least one flux distribution, from the one or more flux distributions, that minimizes or maximizes the objective function when a constraint set is applied to the data structure.
52. A computer-implemented method for developing a dynamic network model, the method comprising:
- accessing archetype data for a plurality of reactions of a network;
- developing a stoichiometric data structure for the reactions based on the archetype data;
- accessing physiological data reflecting a specific context; performing a steady state simulation using the stoichiometric data structure and the physiological data;
- incorporating metabolic data and equilibrium constant data based on the steady state simulation;
- deriving a gradient matrix and dynamic rate equations based on the stoichiometric data structure and the incorporated data; and
- deriving a concentration Jacobian matrix and a flux Jacobian matrix based on the stoichiometric data structure and the gradient matrix so as to characterize dynamic properties of the network, wherein at least some of the stoichiometric, metabolic or equilibrium constant data is obtained by sampling a solution space.
53. A computer-implemented method for the analysis of an effect of a drug or chemical agent on a physiological function, the method comprising:
- accessing archetype data for a plurality of reactions of a biological network;
- developing a stoichiometric data structure for the reactions based on the archetype data;
- accessing physiological data reflecting a specific context;
- performing a steady state simulation using the stoichiometric data structure and the physiological data;
- incorporating metabolic data and equilibrium constant data based on the steady state simulation;
- adding a reaction to the stoichiometric data structure, removing a reaction from the stoichiometric data structure, or adjusting a constraint for a reaction in at least a set of constraints for reactions in the stoichiometric data structure to reflect a measured effect of a drug or chemical agent on the activity of the reaction;
- deriving a gradient matrix and dynamic rate equations based on the stoichiometric data structure and the incorporated data;
- deriving a concentration Jacobian matrix and a flux Jacobian matrix based on the stoichiometric data structure and the gradient matrix so as to characterize dynamic properties of the network; and
- predicting, at a computer system, a physiological function related to the drug or chemical agent based on a simulation including at least the added reaction, removed reaction, or adjusted constraint.
54. A computer-implemented method for developing a dynamic network model, the method comprising:
- accessing archetype data for a plurality of reactions of a network;
- developing a stoichiometric data structure for the reactions based on the archetype data; performing a steady state simulation using the stoichiometric data structure and physiological data;
- obtaining metabolic data and equilibrium constant data based on the steady state simulation;
- deriving a gradient matrix and dynamic rate equations based on the stoichiometric data structure and the obtained data; and
- deriving a concentration Jacobian matrix and a flux Jacobian matrix based on the stoichiometric data structure and the gradient matrix so as to characterize dynamic properties of the network, wherein unknown concentrations of metabolites are determined by solving a steady state mass conservation relationship S·v=0, where S is a the stoichiometric data structure and v is a flux vector.
Type: Application
Filed: Sep 25, 2012
Publication Date: Apr 11, 2013
Applicant: The Regents of the University of California (Oakland, CA)
Inventor: The Regents of the University of California (Oakland, CA)
Application Number: 13/626,807
International Classification: G06F 19/12 (20060101);
