Mass action stoichiometric simulation models: incorporating kinetics and regulation into stoichiometric models

Abstract

The ability to characterize biological dynamics is important for understanding the integrated molecular processes that underlie normal and abnormal cellular states. The availability of metabolomic data, in addition to new developments in the formal description of dynamic states of networks, has enabled a new data integration approach for building large-scale kinetic networks. We show that dynamic network models can be constructed in a scalable manner using metabolomic data mapped onto stoichiometric models, resulting in mass action stoichiometric simulation (MASS) models. Enzymes and their various functional states are represented explicitly as compounds, or nodes in a stoichiometric network, within this formalism. Analyses and simulations of MASS models explicitly show that regulatory enzymes can control dynamic states of networks in part by binding numerous metabolites at multiple sites. Thus, network functional states are reflected in the fractional states of a regulatory enzyme, such as the fraction of the total enzyme concentration that is in a catalytically active versus inactive state. The feasible construction of MASS models represents a practical means to increase the size, scope, and predictive capabilities of dynamic network models in cell and molecular biology.

Figure 1

Data-driven kinetic model construction through the integration of disparate data types. Top panel: Conceptual integration of various data. Stoichiometric network models in combination with equilibrium constants for the reactions in the network and metabolomic data can be used to solve the steady-state mass balance equations for kinetic rate constants. These constants, using bilinear kinetics, can then be used to define the dual Jacobian matrices. Bottom panel: The steps involved in constructing kinetic models, and challenges/quality control (QC) issues associated with each step. The main QC issue associated with each step is highlighted on the right side of the equations for an m×n network with stoichiometric matrix rank r. Failure of the QC step results in the need for iterative revisions of the model, which may be accomplished by, e.g., adjusting flux/metabolite measurements or equilibrium constant approximations based on experimental errors.

Figure 2

Map of the MASS model plus mechanisms for five enzymes. Each enzyme reaction scheme was added and tested before incorporation into the full network. The bottom-right corner of the figure shows the stoichiometric matrix for the mass action model, S_rbc. To construct the regulated model, a single column (HK in the figure) becomes expanded to an additional stoichiometric matrix, S_hk, corresponding to the reaction mechanism for HK. Thus, for each regulated enzyme reaction, a single column in the stoichiometric matrix became multiple columns, and corresponding rows were added to the matrix as well to account for the new reaction intermediates. This modularity is illustrated for HK in the bottom right-hand side of the figure. S_rbc represents the stoichiometric matrix for the mass action network (in the center of the figure). The first column (corresponding to HK) becomes expanded into a new substoichiometric matrix, S_hk. S_hk represents a “subnetwork”, reflected in the illustrated reaction scheme. This process was carried out for all of the enzymes in the network, and they were then all integrated with S_rbc to produce a final stoichiometric matrix for the regulated red cell. The reaction scheme for HK is as described by Mulquiney and Kuchel (29). hkE represents the unbound, free form of the enzyme. The individual steps and interactions account for the catalytic steps in PFK in addition to the allosteric interactions. The reaction scheme for the Rapoport-Leubering shunt is carried out by the same enzyme (27). The reaction schemes for G6PDH are as described by Mulquiney and Kuchel (29). g6pdER represents the active enzyme. The reaction schemes for AK are as described by Hawkins and Bagnara (36). akE represents the enzyme in the catalytic state, and akET represents the tense or inactive form of the enzyme. See the Supporting Material for full metabolite abbreviations.

Figure 3

Tiled phase-plane diagram for all of the fluxes involving the enzyme PFK for the cell in response to a transient, instantaneous decrease in ATP by 60%. The lower-left triangle depicts the trajectories of the corresponding variables. The upper-right triangle shows the correlation coefficient over the simulated time course for the two corresponding entries (when the correlation coefficient is >0.85, the slope of the line appears below it). Note that not all of the steps are correlated with one another. Indeed, the net flux through the enzyme is actually a subnetwork in and of itself. The catalytic step is PFK3. PFK14–PFK17 reflect ATP binding to the tense form of the enzyme, and PFK18–PFK21 represent Mg binding to the tense form of the enzyme. Correlations among some of the fluxes are observed, such as the sequential binding of multiple ATP molecules to the enzyme. Other fluxes (e.g., PFK15 and PFK18) move dynamically independently from the other PFK fluxes.

Figure 4

Network responses to pulsed energy loads. The energy charge (ATP + 0.5^∗ADP/(ATP + ADP + AMP)) (43) is shown as a function of time during a rectangular pulse load on ATP lasting from t = 5 to t = 6 h (magnitude 4 mM/h). Network responses to pulsed redox loads, and the changes in glutathione and NADPH/(NADPH + NADP) as a function of time when the cell is exposed to a first-order load on NADPH (100^∗NADPH) from t = 5 to t = 10 h are also shown.

Figure 5

Enzyme-state plots for HK and AK under energy (left panel) and redox (right panel) perturbations from the simulations in Fig. 4. The majority of the enzyme is in the inactivated state; however, in response to either load, the amount in the inactivated state decreases and the amount in the activated state increases to increase the flux through glycolysis or the pentose pathway.