Abstract
Background
Higher plants and algae are able to fix atmospheric carbon dioxide through photosynthesis and store this fixed carbon in large quantities as starch, which can be hydrolyzed into sugars serving as feedstock for fermentation to biofuels and precursors. Rational engineering of carbon flow in plant cells requires a greater understanding of how starch breakdown fluxes respond to variations in enzyme concentrations, kinetic parameters, and metabolite concentrations. We have therefore developed and simulated a detailed kinetic ordinary differential equation model of the degradation pathways for starch synthesized in plants and green algae, which to our knowledge is the most complete such model reported to date.
Results
Simulation with 9 internal metabolites and 8 external metabolites, the concentrations of the latter fixed at reasonable biochemical values, leads to a single reference solution showing βamylase activity to be the ratelimiting step in carbon flow from starch degradation. Additionally, the response coefficients for stromal glucose to the glucose transporter k_{cat }and K_{M }are substantial, whereas those for cytosolic glucose are not, consistent with a kinetic bottleneck due to transport. Response coefficient norms show stromal maltopentaose and cytosolic glucosylated arabinogalactan to be the most and least globally sensitive metabolites, respectively, and βamylase k_{cat }and K_{M }for starch to be the kinetic parameters with the largest aggregate effect on metabolite concentrations as a whole. The latter kinetic parameters, together with those for glucose transport, have the greatest effect on stromal glucose, which is a precursor for biofuel synthetic pathways. Exploration of the steadystate solution space with respect to concentrations of 6 external metabolites and 8 dynamic metabolite concentrations show that stromal metabolism is strongly coupled to starch levels, and that transport between compartments serves to lower coupling between metabolic subsystems in different compartments.
Conclusions
We find that in the reference steady state, starch cleavage is the most significant determinant of carbon flux, with turnover of oligosaccharides playing a secondary role. Independence of stationary point with respect to initial dynamic variable values confirms a unique stationary point in the phase space of dynamically varying concentrations of the model network. Stromal maltooligosaccharide metabolism was highly coupled to the available starch concentration. From the most highly converged trajectories, distances between unique fixed points of phase spaces show that cytosolic maltose levels depend on the total concentrations of arabinogalactan and glucose present in the cytosol. In addition, cellular compartmentalization serves to dampen much, but not all, of the effects of one subnetwork on another, such that kinetic modeling of single compartments would likely capture most dynamics that are fast on the timescale of the transport reactions.
Background
Insolation is the dominating contributor to a sustainable terrestrial energy balance, whether directly captured or transformed into secondary sources such as wind or biomass. Plants have evolved to use this resource to provide themselves with the lowpotential carbon necessary for growth by splitting water and fixing carbon dioxide, a known greenhouse gas. During the light photosynthetic reactions, more carbon can be fixed than can be productively marshalled for growth, and cells store this excess carbon in compact polymers such as starch. Chloroplastic starch is stored in the form of granules [1] that consist of both linear and branched polymers of glucose; the process of phase transfer between the granule and the aqueous chloroplast stroma is not known in great detail, although phosphorylation by glucan water dikinase [2,3] and phosphoglucan water dikinase [4,5] may be involved. Amylopectin, the major component of starch, is moderately branched, comprises the majority of starch mass, and is responsible for the crystallinity of starch granules. Essentially unbranched amylose, on the other hand, is amorphous and constitutes up to 30% by weight of starch, depending on culture status [6]. The backbone of both polymers arises from α1,4 glycosidic bonds; the α1,6 branches of amylopectin occur every 24 to 30 glucose units.
Much remains unknown about the biochemical pathway for starch degradation in plants and algae. Smith, et al. have proposed a pathway of starch degradation in Arabidopsis thaliana leaves, whereby starch is released from the granule in a soluble form, then debranched to yield soluble linear glucans in the chloroplast stroma [1,7]. Two mutually alternative degradation pathways can then cleave the linear glucans. In the first, chloroplastic glucan phosphorylase catalyzes the phosphorolytic release of glucose1phosphate [8,9], which is cleaved to triose phosphate and the latter antiported in exchange with cytosolic inorganic phosphate [10]. In the second, βamylase hydrolyzes linear glucans to maltose and maltotriose. Recent results show this second pathway to be more usual in the Arabidopsis thaliana chloroplast [1,11,12]. βamylase releases maltose from the nonreducing ends of linear glucan chains at each catalytic turnover [1], but cannot act on chains of less than four glucosyl units, leading to maltotriose as a byproduct of βamylolytic degradation. Although generally functioning as a predominantly hydrolytic enzyme in vivo, βamylase from sweet potato has been shown to catalyze the condensation of maltose to maltotetraose in vitro [13].
Once liberated, maltose and maltotriose can enter chloroplastic and cytosolic carbon pathways. Strong experimental evidence suggests that maltose is exported from the chloroplast stroma to the cytosol by the MEX1 transporter [14]. Cytosolic transglucosidase DPE2 [1517] can split the transported sugar, glucosylating a soluble endogenous acceptor [1] and freeing glucose. A possible candidate for this acceptor is a soluble arabinogalactan [18,19] that serves as a glucosylation substrate of cytosolic glucan phosphorylase in vitro with glucose1phosphate as the donor [1,18]. DPE2 and reversible glucan phosphorylase acting together may therefore result in maltosederived glucose1phosphate. The maltotriose product of chloroplastic βamylase may be acted upon by a disproportionating enzyme (α1,4 glucanotransferase, DPE1) [20] catalyzing the disproportionation of two maltotriose molecules to glucose and maltopentaose, that can in turn be cleaved by βamylase to produce maltotriose to reenter the disproportionation reaction and maltose to be transported out of the stroma. At the catabolic end of starch degradation, cytosolic glucose is phosphorylated at C6 by hexokinase [1,21] for entry into general cellular metabolism.
There currently exists no mathematical model of starch degradation pathways that includes the details discussed in the previous paragraphs. We therefore report the development of a detailed ordinary differential equation (ODE) model that includes most of the biochemical reactions discussed above, and detailed kinetic mechanisms captured from the scientific literature, presumed by direct comparisons, or postulated within the range of characterized mechanisms and parameter values. This approach of hypothesizing unknown values differs from flux balance [22,23] or energy balance [24,25] approaches, where extrema in carefully crafted (i.e., the setting of lower and upper bounds, and construction of the objective function(s)) flux spaces are evaluated. Almost all biochemical reactions are catalyzed by enzymes that can saturate, respond nonlinearly to changes in metabolite concentrations, and comprise components of a reaction network capable of dynamic evolution outside of the steadystate assumption. Although insightful results have been obtained from several studies [2628] on specific metabolic pathways incorporating known enzyme kinetics, the most promising features of the current modeling approach are a greater understanding of potential nonlinear network dynamics, and the possibility of characterizing the highdimensional space of metabolic responses with respect to enzyme concentrations and parameters using modern highperformance computing.
The current model focuses on steadystate solubilized starch catabolism. A starch degradation model previously postulated [1] did not include the effective competition for βamylase of maltose condensation to maltotetraose [13]; the effects of this alternative sink reaction on βamylase turnover have been included in the present model by inclusion of an inhibition term as previously formulated by Shiraishi and coworkers [29]. The biochemical reactions in the current starch degradation model have been schematically represented in Figure 1. Enzyme kinetic parameters are taken from reported values, calculated from Haldane relations [30], or assigned reasonable values within the relatively limited range of known values for the particular parameter in question.
Figure 1. Biochemical and transport reactions in current starch degradation model. The carbon flux in the starch degradation pathway from solubilised starch in the chloroplast stroma to hexose phosphates in the cytosol has been schematically represented in this figure. The pathway consists of seven enzyme catalyzed reactions and two transport reactions. The different types of species in this figure have been colorcoded. The metabolite species are shown in black; enzymes are in red and transporters in orange; and, irreversible and reversible reactions are shown as blue and purple arrows, respectively. The dashed arrow showing maltotriose formation from linear glucans reflects that this reaction is substoichiometric, occurring only once per oddlength chain.
Results
Model content
The model detailed above contains 17 metabolites, 6 enzymes, 2 transporter proteins, and 3 inhibitors that participate in 9 reactions characterized by 63 enzyme kinetic and binding parameters. Eight metabolites are at the boundary of the system and therefore act as systemic parameters and are referred to as "external metabolites"; the remaining nine are free, and called "internal metabolites". The model encompasses the chloroplast stroma, cytosol, and chloroplast intermembrane space containing two transporter proteins linking the stromal and chloroplastic metabolite pools. The intermembrane space impacts simulations only by defining the volume affecting these transporters' concentrations.
The initial concentration of the intermediate linear linkage group Starch_{db}_CS was set to zero in all calculations. The pH value of the cytosol is assumed to be 7, so that the proton concentration in the cytosol is fixed at 0.1 μM. All other internal and external metabolites take on concentrations in the molar to nanomolar range. Although extremes of this range are of questionable physiological significance, two conciliatory considerations apply. First, the dynamical system is dictated primarily by the controlling equations, and secondarily by the particular point in concentration space that the system occupies. Thus, the system will evolve toward a steady state in as robust a manner as the underlying phase space permits, dissipating or accumulating excess mass from the external metabolite baths as necessary. Second, we are explicitly interested in the fixed point(s) of the dynamical space arising from the reaction network topology and the structure of the kinetic equations, such that concentrations beyond biologically relevant bounds are desirable to characterize the possible behaviors of the system.
Six external metabolites were held at fixed concentrations to reflect their coupling to a homeostatic cellular reservoir: chloroplastic starch glucosyl residues, and the cytosolic pools of ATP, ADP, phosphate, glucose1phosphate, and glucose6phosphate. Starch and polymerized glucosyl units were also modeled as external metabolites to reflect a large starch reservoir, such as would be relevant to the transition between photosynthetic starch accumulation and biofuelproducing fermentative metabolism. Three species (cytosolic reduced glutathione, cytosolic glucose1,6bisphosphate, and cytosolic 2,3bisphosphoglycerate) act as hexokinase inhibitors but do not otherwise participate in any reaction; they are therefore classified as parameters. It should be noted that from a mathematical and dynamical perspective, external metabolites are also model parameters. Our chosen distinction between "external metabolites" and parameters (the three metabolites above as well as kinetic and binding constants) is based on participation or nonparticipation as a reactant or product in a modeled reaction, and so is less operational than chemically ontological. Enzyme and transporter protein concentrations are held constant to reflect a particular metabolic state.
States explored
Three primary models are explored in detail. The "reference" model or state is described in detail in the Methods section, and is comprised of a single best estimate of concentrations, kinetic parameter values, and protein concentrations. Two other models derived from this reference system were also analyzed to explore the robustness of the reference system to perturbations of kinetic parameters or enzyme concentrations. The first decreased and 10fold, and is named the "parameterperturbed" model or state. The second increased the reference concentrations of βamylase and MEX twoand 10fold, respectively, and is named the "enzymeperturbed" model or state. In addition to these three models, a space of models differing in either initial concentrations of internal metabolites or fixed concentrations of external metabolites was generated by 2way sampling of 14 concentrations (8 internal and 6 external), thus yielding a body of 2^{14 }= 16,384 individual simulations that is analyzed and discussed separately.
Response coefficients with respect to kinetic and binding parameters
Response coefficients quantify the sensitivity of steadystate variables to variation in model parameters, and so pertain to a particular steady state [31]. Here we focus on the response of steadystate metabolite concentrations only, and include as model parameters both kinetic constants and enzyme concentrations. For the reference steady state, the concentration response coefficients with respect to variations in kinetic and binding parameters in the model are presented as a heat map representation in Figure 2. The elements x_{im }of the matrix in Figure 2 are the response coefficients for each internal metabolite species i at steady state, with respect to each parameter m, excepting intermediate species Starch_{db}_CS (S_{db}) for reasons discussed below.
Figure 2. Metabolite response coefficients with respect to kinetic and binding parameters for the reference model. The matrix of response coefficients with respect to the enzyme kinetic and binding parameters in the starch degradation model at the reference steady state is shown as a heat map representation. Only coefficients with at least one magnitude ≥ 0.05 are included.
The βamylase k_{cat }() has the most response coefficients above zero, with the starch Michaelis constant for the same reaction () showing a similar pattern below zero. This inverse relationship is expected due to the definitions of these kinetic parameters. This responsiveness extends through the plastidic starch degradation products maltose and maltotriose, weakens at cytosolic maltose, and is essentially zero for cytosolic glucose and arabinogalactan.
To evaluate the sensitivity of the model's response to the particular kinetic parameters used, a perturbed model with the βamylase starch and MEX decreased by 10fold ("parameterperturbed" model) was simulated to steady state, and yielded the response coefficients shown in Figure 3. These parameters were chosen based on βamylase's role in carbon flux limitation, and MEX's role in coupling the chloroplastic and cytosolic compartments. None of the significant response coefficients in Figure 3 differ in sign from their counterparts in Figure 2. The ratios of parameterperturbed to reference coefficients are mapped in Figure 4. Response coefficients less than a cutoff value (10^{3}) were considered as zero, and the ratio of two such coefficients replaced by 1.0 to avoid numerical noise from dividing small numbers. Figure 4 shows that the response coefficients of cytosolic maltose with respect to βamylase and , and five DPE2 kinetic parameters (k_{cat}, K_{M }for maltose and arabinogalactan, and K_{i }for arabinogalactan and glucosylated AG) are reduced in magnitude relative to those for the reference steady state. However, maltose responsiveness to the DPE2 and hexokinase equilibrium constants increases, due to more negative values upon perturbation. Cytosolic glucose and arabinogalactan also show decreased response with respect to and . Chloroplastic maltopentaose response coefficients exhibit smaller reductions versus , and . Thus, most changes to concentration responsiveness are decreased upon the coupled decrease in βamylase turnover and increase in maltose export from the stroma, primarily relative to βamylase starch degradation kinetics and cytosolic maltose cleavage, with DPE2 and hexokinase K_{eq }being exceptions. This pattern primarily reflects the decreased starch cleavage rate, leading to roughly 100fold lower stromal maltose (240 vs. 2.3 μM), but approximately the same cytosolic maltose (110 vs. 130 mM) concentrations, lower flux through the DPE2 reaction, and a greater sensitivity to the direction of the DPE2 and hexokinase reactions. No ratio is greater than 1.1 or less than 0.08.
Figure 3. Metabolite response coefficients with respect to kinetic and binding parameters for the parameterperturbed model. The response coefficients in Figure 4 are recalculated from the simulation results obtained on reducing both the βamylase turnover number for starch hydrolysis () and the maltose exporter (MEX) Michaelis constant for maltose transport () by a factor of 10 compared to the corresponding values for Figure 2. The change in these two kinetic parameters yields a new steady state and therefore a different set of values for the response coefficients.
Figure 4. Ratios of parameterperturbed and reference response coefficients. The ratios of each response coefficient in Figure 3 to the corresponding coefficient in Figure 2 are shown. The ratio is set to unity if both response coefficients are less than 10^{3}.
Row (ρ) and column norms (κ) of the parametric response coefficient matrix
To evaluate the global sensitivity of specific metabolite concentrations S_{i }with respect to kinetic and binding parameters p_{m }in Figure 2, the Euclidean norm of each metabolite's response coefficient vectors were calculated as the rootmeansquare of response coefficients along row i,
The converse quantity, namely the overall sensitivity of steadystate metabolite concentrations to variation in individual model parameter values p_{m}, was calculated as the Euclidean norm along column m,
For the set of species considered in Figure 2, we show the bar plot of the corresponding ρ_{i }in Figure 5. The calculated κ_{m }for each of the kinetic and binding parameters in Figure 2 are shown in Figure 6, which includes only the twenty largest κ values. The steadystate concentration of stromal maltopentaose is the most generally sensitive metabolite in the present model, while cytosolic glucosylated arabinogalactan shows the least global sensitivity (not shown). Figure 6 indicates that and have the largest aggregate effects on metabolite concentrations at this reference steady state, consistent with these parameters corresponding to the most positive and negative response coefficients in Figure 2.
Figure 5. Metabolite sensitivities of individual metabolites versus kinetic and binding parameters for reference and parameterperturbed states. The collective sensitivity ρ_{i }of individual metabolites to all kinetic and binding parameters are shown for the reference state, and a model with 10fold reduced βamylase and MEX values. The quantities ρ_{i }are defined as the row norms of individual response coefficients, as discussed in the text.
Figure 6. Collective sensitivities with respect to individual kinetic and binding parameters for reference and parameterperturbed states. Collective sensitivity κ_{m }of all steady state metabolite concentrations with respect to individual kinetic or binding parameters are shown for the reference state, and a model with 10fold reduced values of βamylase and MEX . The quantities κ_{m }are defined as the column norms of individual response coefficients, as discussed in the text.
Despite maltotriose and maltose being sensitive to more kinetic and binding parameters than maltopentaose, maltopentaose has the greatest ρ due to the large magnitudes of the maltopentaose parametric response coefficients for reactions involving the conversion of starch and maltopentaose to maltose and maltotriose. A potential reason for this strong coupling is a positive feedback loop: degradation of starch and maltopentaose by βamylase yields maltose and maltotriose, two molecules of which disproportionate to yield maltopentaose. The least sensitive metabolite, glucosylated arabinogalactan, is simply a nontransient intermediate for transfer of a glucosyl unit from maltose to phosphate via a soluble arabinogalactan that is regenerated at the end of the transfer. Modification of flux into this reaction can increase or decrease the AG/GlcAG ratio, but the concentration response of either is bounded by the total fixed arabinogalactan concentration available. The response coefficient of greatest magnitude for glucosylated arabinogalactan corresponds to the equilibrium constant for glucosyl transfer, emphasizing the locality of this metabolite's response.
Figure 5 also displays the row norms for the parameterperturbed steady state. Chloroplastic maltopentaose ρ is reduced relative to that for the reference steady state, due to the lower values of the response coefficients of maltopentaose with respect to three βamylase enzyme kinetic parameters (Figure 4). In spite of the significant reduction in the response coefficients of cytosolic glucose and arabinogalactan versus a number of kinetic and binding parameters for the parameterperturbed state, the ρ values of these two metabolites are comparable between states. This appears to be so because the response coefficients of cytosolic glucose and arabinogalactan are very small in both cases, so large relative changes have little effect on the norm. The ρ value for cytosolic maltose is also not significantly altered upon the perturbation considered, due to compensatory changes in magnitudes of R^{maltose}.
Upon perturbation, the response coefficient column norms with respect to βamylase k_{cat }and K_{M }and several DPE2 kinetic parameters are reduced in magnitude, whereas the column norms with respect to hexokinase and DPE2 equilibrium constants become larger (Figure 6). These trends follow directly from the observations in Figure 4 which have been discussed in detail above.
Comparison of response coefficients for chloroplastic glucose with respect to kinetic and binding parameters
Because glucose is a precursor for fermentative pathways yielding important potential biofuels, including ethanol and hydrogen, its associated response coefficients are of special interest. The glucose steadystate concentration in the stroma is effectively sensitive to only four kinetic parameters, two having direct proportionality (i.e., parameter increase yields metabolite increase), and two inverse proportionality, with all four magnitudes close to 1. is positive and is negative, consistent with expectations regarding k_{cat }and K_{M}. A similar but inverse relationship is evident in the k_{cat }and K_{M }parameters of the glucose transporter, where turnover implies movement of glucose out of the plastid; so, decelerating this transporter's action will increase steadystate chloroplastic glucose concentration. DPE1, maltose transport, and debranching kinetics had little effect on stromal glucose levels.
Rearrangement and integration of the definition for ρ_{i }above allows us to discover the parametric dependence on S_{i }as
In other words, S_{i }varies as the ρ_{i}^{th }power of p_{m}. Because the response coefficient of chloroplastic glucose with respect to is positive and very close to unity, if is varied while keeping all other parameters fixed the steadystate concentration of chloroplastic glucose should increase linearly with the parameter value. This prediction is validated by the plot in Figure 7, where the steadystate concentration of glucose in the chloroplast stroma is shown to vary linearly with βamylase .
Figure 7. Variation in simulated steady state concentration of chloroplast stromal glucose with the βamylase turnover number for starch hydrolysis. The dependence of the simulated chloroplastic glucose steady state concentration on βamylase turnover number for starch hydrolysis () is seen to be linear over a range from ~0.03 to 0.3 s^{1}.
Response coefficients with respect to enzyme and transporter concentrations
The response coefficients of steadystate metabolite concentrations with respect to enzyme and transporter concentrations for the reference steady state in this study are shown as a heat map in Figure 8. Chloroplastic maltopentaose and cytosolic arabinogalactan, glucosylated arabinogalactan, and glucose show responses near zero, implying lack of sensitivity to any enzyme concentration. Because [E] enters linearly into the kinetic equations used, this situation can arise only if (a) enzymes acting on these species are saturated but present in high concentration (i.e., have low turnover numbers), such that small changes in enzyme concentration have a proportionally small effect, or (b) turnover is fast enough to keep reactions near equilibrium, so slowing down interconversion slightly by lowering the enzyme concentration does not result in an observable change. Maltopentaose participates in two reactions in the model, both of which are upstream of most other reactions. The response coefficients of maltopentaose with respect to enzymes catalyzing these other reactions are small, as expected for net flux "downstream" and rate limitation "upstream". If downstream enzymes were saturated, a sensitivity of upstream metabolite concentrations to these enzymatic levels would arise. This relationship can be seen in the DPE1catalyzed disproportionation that yields maltopentaose from maltotriose, and has a higher turnover number (50 s^{1}) than the downstream degradation of maltopentaose by βamylase (~0.1 s^{1}). As expected, maltopentaose levels have a stronger dependence on the kinetics of the ratelimiting degradation step, and so the maltopentaose steadystate concentration is more sensitive to βamylase variation than to that of DPE1.
Figure 8. Metabolite response coefficients with respect to enzyme and transporter levels for the reference model. Metabolites are along the vertical axis, and enzyme and transporter concentrations along the horizontal axis. Only coefficients with at least one instance of magnitude ≥ 0.05 are shown.
Both arabinogalactan and glucosylated arabinogalactan participate in two reactions that are characterized by high enzymatic turnover in both the forward and reverse directions. Hence, small quantities of the enzymes (DPE2 and cytosolic glucan phosphorylase) catalyzing those two reactions are sufficient to maintain the reactions near equilibrium, resulting in the negligible sensitivities of arabinogalactan and glucosylated arabinogalactan steadystate concentrations to these enzyme concentrations. A similar situation may be assigned to cytosolic glucose, which arrives from the chloroplast stroma via a plastidic glucose transporter with a high transport turnover number. Glucose is also formed directly in the cytosol from maltose and arabinogalactan by DPE2, with high turnover numbers in both the forward and reverse directions. Finally, cytosolic glucose is phosphorylated to glucose6phosphate by hexokinase, which has a large forward turnover number. In all three cases, the enzymes and transporters involved are subsaturated, so that reaction kinetics are relatively insensitive to small changes in enzyme levels.
At the other extreme of sensitivity, the chloroplastic βamylase exhibits the highest number of large positive response values, with the product maltose having the largest positive sensitivity. This suggests that hydrolysis of linear starch fragments to maltose is limiting the starch degradation flux near the combination of enzyme concentrations and parameters associated with the reference steady state. The βamylaseassociated response coefficients for cytosolic glucose, arabinogalactan, and glucosylated arabinogalactan have their second largest values for any enzyme, despite these metabolites being the most distant variable metabolites from chloroplastic maltose in the reaction network. To the extent that these numerical results represent cellular behavior, βamylase is a natural target for increasing in vivo activity.
To examine how representative the reference model and steady state are relative to enzyme levels, in Figure 9 are shown the response coefficients in Figure 8 for a model obtained by increasing the βamylase concentration twofold and the maltose exporter concentration tenfold relative to the reference model (the "enzymeperturbed" model). The perturbation significantly affects only the response coefficients of cytosolic maltose with respect to βamylase and DPE2 concentrations, changing them 2.14fold and 2.17fold, respectively. Higher fluxes of maltose production and transport into the cytosol imply increased dependence of the model system on any reaction that consumes cytosolic maltose to achieve a steady state maltose concentration, so increased DPE2 responsiveness is expected. Given the predictability of the qualitative response and the relatively small change in response coefficients, we consider the reference steady state to be a reasonable representation of the metabolism being explored.
Figure 9. Metabolite response coefficients with respect to increased enzyme and transporter levels for the enzymeperturbed model.
Row norms of response coefficient matrix with respect to enzyme and transporter concentrations and response coefficients for chloroplastic glucose
The calculated ρ_{i }values for the metabolite species in Figure 8 indicate that enzyme and transporter concentrations generally have the greatest impact on the reference steadystate concentrations of chloroplastic maltose and glucose (Figure 10). The two response coefficients of largest magnitude for chloroplastic glucose with respect to βamylase and the plastidic glucose transporter both had values of 1.0, with the next largest coefficients effectively zero (~10^{18}). To increase the standing concentration of chloroplastic glucose, one would need either to overexpress the chloroplastic βamylase enzyme, or suppress glucose transport. Simulations with four different concentrations of the βamylase enzyme show a linear variation on steadystate stromal glucose concentration (Figure 11), consistent with a chloroplastic glucose response coefficient with respect to βamylase close to unity. For context, experiments suggest that the level of βamylase activity in leaves of the hba1 Arabidopsis thaliana mutant is about six times the corresponding activity level in wild type Arabidopsis leaves [32], and that the combined effects of prolonged exposure to light and treatment with sucrose solution can enhance the βamylase activity in wildtype Arabidopsis leaves by approximately a factor of four [33].
Figure 10. Collective responses of individual metabolites with respect to all catalyst levels for reference and enzymeperturbed states. The collective sensitivity ρ_{i}, defined as the row norms of individual response coefficients, of individual metabolites i to all enzyme and transporter concentrations for the reference and enzymeperturbed states are shown.
Figure 11. Dependence of simulated steady state chloroplastic glucose concentration on βamylase concentration. The simulated steady state chloroplastic glucose concentration increases linearly with βamylase concentration in the range 0.038 to 0.152 μM.
For the enzymeperturbed model, row norms are plotted alongside the reference state values in Figure 10. The only ρ value significantly different from that of the reference steady state is the value for cytosolic maltose, arising primarily from the significant increase (2.17) of as explained earlier. The contribution to ρ from the 2.14fold change in between the two states is minor because this response coefficient is of small magnitude.
Comparison of κ values from response coefficient matrix with respect to enzyme and transporter concentrations
The calculated κ values for each of the enzymes or transporters in Figure 8 are shown in Figure 12 for the reference and enzymeperturbed states. For the former, variation of chloroplastic βamylase concentration has the maximum aggregate effect on the steadystate concentrations of the small metabolome considered in Figure 8. The plastidic maltose exporter has the second largest effect on the steadystate concentrations of the metabolites, with the major contribution originating from the large negative . For the enzymeperturbed steady state, the κ_{DPE2 }value is the only column norm that differs significantly upon this perturbation. In contrast to the increase in seen, the analogous increase in the κ_{βamylase }value is insignificant, because the response coefficient of maltose vs. βamylase is quite small in magnitude such that even a significant relative change has little impact on κ.
Figure 12. Collective responses of metabolome with respect to individual catalyst levels for reference state. The collective sensitivity κ_{m }of all steady state metabolite concentrations with respect to individual enzyme and transporter concentrations are shown for the reference and enzymeperturbed states. The quantities κ_{m }are defined as the column norms of individual response coefficients, as discussed in the text.
Enumeration of initial metabolic states, convergence thresholds, and clustering analyses
To determine if the starch degradation model constructed was capable of hosting multiple steady states, and to examine the variability of the model behavior with respect to variation of metabolic pool sizes external to the dynamic core, multiple simulations were run from different initial internal, and fixed external, metabolite concentrations (Tables 1 and 2). The concentration of starch is not an independent variable and is related to the concentration of starch glucosyl units. Out of the 17 metabolites, we therefore sampled the initial concentrations of 8 internal metabolites and the fixed concentrations of 6 external metabolites. Sampling n initial or fixed concentrations of each of these 14 species yields n^{14 }combinations of initial conditions; thus, even for our modest sampling number of n = 2, corresponding to initial or fixed concentrations of 0.1 and 1000 μM, 16,384 temporal integrations were needed. For each of 2^{6 }= 64 system parameter combinations, there are thus 2^{8 }= 256 different initial value vectors to test for multistationarity. For reasons of parallel load balancing, we elected a fixed integration time of 10^{7 }virtual seconds as a stopping criterion, rather than test for convergence directly. The final metabolite and time derivative data was then postprocessed. To analyze the results, we chose to take a global view of all the points generated, as well as to focus on the most highly converged points.
Table 1. Concentrations of external metabolites, enzymes, transporters and inhibitors
Table 2. Initial concentrations of internal metabolites
We have therefore clustered the temporal integration end point data at different convergence thresholds, with the following intent. Because each end point represents a trajectory potentially at a different "stage" of evolution, and starting from a different initial point within the metabolite concentration space, we anticipate that selection based on increasingly tight convergences will progressively select for classes of initial points in the metabolite space that lie closer to fixed points of the different systems created by sampling fixed values of the external metabolites. Nevertheless, the long fixed virtual time of the simulations should permit some structure to be observed in the total dataset, and most importantly to classify different clusters based on their associated initial metabolite vectors. If bivalued samples of the six external metabolites completely determine the phase spaces and each such space has one steadystate fixed point, then one expects the endpoint data to fall into 2^{6 }= 64 different clusters, each associated with a different sample of the external metabolite concentrations. Deviations from this null result and multivariate correlations are potentially interesting features of the formulated model system. Centroids with concentration coordinates intermediate between two fixed sample values simply highlight membership within a cluster of trajectories starting from different values of external metabolites, suggesting that despite these different values, the trajectory end points are still near each other in metabolite coordinate space. Correlations among centroid coordinates can elucidate deterministic relationships between metabolites—for example, complete correlation between variables implies that some are dependent on others and hold little explanatory power. Weaker correlations should reflect control relationships among metabolites—an anticorrelation between freely variable × and Y implies an increase of [X] puts negative pressure on [Y].
In Figure 13 is the cumulative distribution of trajectories at different convergence thresholds. All 16,384 trajectories were found to have converged to less than 10^{2 }(μM s^{1}) μM. Notable drops in trajectory number can be seen between 10^{6 }and 10^{8 }(μM s^{1}) μM, and between 10^{14 }(μM s^{1}) μM and the tightest cutoff explored here, 10^{16 }(μM s^{1}) μM. The classes of trajectories found at the more tightly converged ends of these two breakpoints might be expected to differ qualitatively from more loosely converged trajectories and reflect important features of the system. To seek these features, three classes of end points were mapped—all points, those trajectories evolving at less than 10^{8 }(μM s^{1}) μM^{1}, and those at less than 10^{16 }(μM s^{1}) μM^{1}. In each case, correlation coefficients between metabolites over the cluster centroids were calculated, in order to highlight covariation among the species in the model. Data are provided as Additional File 1. Examination of all 16,384 end points shows that stromal maltose, maltotriose, and maltopentaose levels are dictated by the size of the fixed starch pool, with cytosolic maltose levels correlating predominantly with this pool as well. The latter correlates negatively but weakly with the arabinogalactan pool, ATP, and phosphate, and positively with glucosylated arabinogalactan and glucose1phosphate. Cytosolic glucose correlation coefficients, though less than 1.0, show a similar correlation pattern similar those of stromal glucose. The magnitude of 0.3 suggests that most of the cytosolic glucose coupling is to species other than starch and maltose oligomers in the stroma; the uniform magnitude of correlation coefficients to the latter suggest that transport between the stroma and the cytosol may decrease the metabolic coupling between the starch repository and cytosolic glucose levels. Arabinogalactan is anticorrelated with cytosolic maltose, GlcAG, and glucose1phosphate, and correlated positively with cytosolic glucose and free phosphate. GlcAG correlates negatively with cytosolic glucose, AG, and free phosphate, but positively and weakly with cytosolic maltose, ADP, and glucose1phosphate. This pattern can be reasoned from the pathway diagram in Figure 1, with the possible exception of the negative correlation between AG and GlcAG. This latter relationship likely arises due to the conserved pool size of total arabinogalactan—if both varied independently, more GlcAG would push the reversible DPE2 reaction toward AG; since the total pool size is constant, however, more GlcAG must mean less AG. Negative correlations between ATP and cytosolic glucose and maltose are consistent with increases of the latter putting pressure on ATP supply in the hexokinase reaction. Phosphate shows mainly positive correlation with AG and glucose1phosphate, and negative with GlcAG, arising from the glucan phosphorylase and DPE2 reactions together with the fixed pool of AG + GlcAG. Glucose1phosphate couplings reflect these same relationships. The glucose6phosphate centroid coordinates show no strong coupling to any other metabolite in this total data set.
Figure 13. Cumulative distribution of trajectory number as a function of convergence threshold. Each integration starting from a unique point in the metabolite space was stopped at 10^{7 }simulated seconds, and the convergence evaluated as described in the text.
Additional file 1. Clustering analysis. A single Excel 2008 file containing statistics about kmeans clusters for 10^{2 }convergence cutoff with 63 clusters, and 10^{8 }convergence cutoff with 45 clusters. Correlation coefficients are also contained as separate sheets within the file, colorcoded by magnitude of coefficient: gray (1.0), orange (0.50.99), green (0.20.49), and white (< 0.2).
Format: XLSX Size: 42KB Download file
Applying a more stringent threshold of 10^{8 }(μM s^{1}) μM^{1 }and an initial seed of 64 cluster centroids results in 45 final clusters encompassing all 4,771 end points. Comparison of cluster centroid sizes and coordinates to those of the 10^{2 }(μM s^{1}) μM^{1 }dataset show that trajectories with 1000 μM starch are underrepresented at this threshold—each starch concentration represented 8,192 points in the 10^{2 }(μM s^{1}) μM^{1 }dataset, since this included all 16,384 data points, whereas there are only 768 end points with 1000 μM starch versus 4,003 end points with 0.1 μM starch at 10^{8 }(μM s^{1}) μM^{1}. Reaction rates are proportional to substrate concentrations and no starchbased inhibition was present, so the low representation of highconcentration starch trajectories in the end point dataset suggests that these trajectories may have started further from their systems' fixed points in general. On the other hand, samples with lower ATP and orthophosphate concetrations are also selected against; for ATP, there are 3,773 and 998 end points with 1 mM and 0.1 μM fixed concentrations, respectively, and for phosphate 3,334 and 1,437 end points. In this case, rate limitation in the lowerconcentration samples may broadly explain the patterns seen. The remaining external metabolites show roughly equal numbers of clusters at each of the two sampled concentrations (see Additional File 1).
Correlation analysis of this smaller, more tightly converged subset of end points shows patterns among stromal starch and maltooligomers similar to those seen for the total dataset. Here, however, arabinogalactan is strongly coupled to the starchdependent stromal metabolites due to glucose transport and a full correlation of 1000 μM starch glucosyl units with 1.11 mM AG (the lower starch amount of 0.1 does not correlate neatly with any one AG concentration, thus lowering the magnitude of the correlation coefficient). Within this smaller subset of sampled points, the coupling strength between cytosolic and stromal metabolites is larger, and the uniformity of coupling becomes more apparent. This uniformity is likely due in part to the irreversible formulation of the transport reactions, which mimics a cellular situation with high demand for maltose and glucose.
The 213 points representing trajectories with evolution rates less than 10^{16 }(μM s^{1}) μM^{1 }were found to cluster into 8 distinct clusters using a sequential clustering algorithm. Mapping the Euclidean distances between them gives rise to an obviously apparent symmetry, seen in Figure 14. Although there are 29 numerically distinct pairwise distances, they visually cluster into only 6 (zero selfdistance, 2 edge lengths, and 3 diagonals, represented as raw Euclidean distances, rather than logarithms, due to higher visible resolution). The presence of 8 (2^{3}) centroids implies three determining variables; to assess which three metabolites are relevant, the centroid coordinates for all 8 states are examined in Table 3. What is sought is the least number of compounds, or strongly coupled compound groups, that can explain the symmetric array of states. Immediately, species with invariant concentrations may be omitted, as they have no explanatory power; of the remainder, the external cytosolic ADP and glucose 6phosphate pools clearly define a symmetric quartet of coordinate pairs (0.1, 0.1), (0.1, 1000), (1000, 0.1), and (1000, 1000). Compounds with concentrations that correlate with these two compounds can then be eliminated as determinants of the symmetry, as they are not independent variables (e.g., the glucose concentration in row 13 correlates completely with the values of ADP and glucose 6phosphate). The third symmetry determinant appears to be the metabolically coupled pair of arabinogalactan and the glucosylated form, which have only one independent degree of freedom (i.e., their sum is conserved since we neglect biosynthetic or catabolic pathways as mentioned above).
Figure 14. Euclidean distances between steady states obtained with tight convergence threshold. Each steady state is conceptualized as a point in the positive orthant of a multidimensional space, with coordinates equal to steadystate concentrations. The numerical threshold for defining a steady state here was 10^{16 }(μM s^{1}) μM^{1}. Although the distances between different points shown are numerically unique, only 6 are visually distinct due to the compression of small differences by the overall range.
Table 3. Steady state concentration vectors for the 8 distinct steady states obtained for tight convergence threshold.
The symmetry observed in this geometric interpretation of steadystate metabolite vectors is identical to that found in the distance matrix of corner positions for a rectangular parallelepiped with two sides equal. This situation arises from two concentrations sampled per compound (0.1 and 1000 μM), and only 4 compounds defining the total range, two of which comprise a dynamically coupled pair (and so provide one effective degree of freedom), giving 2^{3 }= 8 points. ADP and glucose 6phosphate are external metabolites, and so are fixed at values defining their respective points; the glucosylated arabinogalactanarabinogalactan pair are internal metabolites. The four states on the corners of the square faces differ in the values of fixed concentrations; the two square faces differ in the third independent determinant's final value. The variability of cytosolic maltose, which slightly perturbs the otherwise perfect symmetry and accounts for the numerical uniqueness of the 29 pairwise distances, illustrates modest sensitivity (9.753 to 11.944 μM) with respect to the much larger ranges of ADP, glucose 6phosphate, and AG/GlcAG (0.1 to 1000, 0.1 to 1000, and 177.527 to 1677.080 μM). The AG pool size explains the difference in steadystate maltose concentrations between evenand oddnumbered states—the differences between pairs of states (e.g., S1/S2 vs. S3/S4 in Table 3) correlate with cytosolic glucose levels. This variation is expected from the reversible cytosolic DPE2 reaction.
Although arabinogalactan and glucosylated arabinogalactan are internal metabolites and therefore free to evolve as the simulation does, their sum is conserved and dictated by their initial concentrations. It was of interest to examine how their sampled initial concentrations correlated with the eight distinct solutions found. In Table 4, one can see that if the initial concentrations of both arabinogalactan and glucosylated arabinogalactan are taken to be 0.1 μM, no trajectories that converge to 10^{16 }(μM s^{1}) μM^{1 }are found. However, if one of the pair has an initial concentration of 0.1 μM and the other at 1000 μM, half of the 8 distinct steady states are obtained; the remaining four are found when both initial concentrations are taken as 1000 μM. It is thus the total amount of arabinogalactan (AG + GlcAG) in the system that acts as a third determinant of the symmetry seen in Figure 14—0.2 μM is not represented, 1000.1 μM correlates with the oddnumbered steady states, and 2000 μM with the evennumbered states.
Table 4. Trajectory and steady state number versus initial concentrations of arabinogalactan and its glucosylation product.
Discussion
Several observations arise about starch biochemistry and general metabolic simulation from the modeling and simulation of soluble starch degradation kinetics described above. First, reversible βamylase action on starch should be incorporated to account for lower starch hydrolysis rates upon maltose accumulation, such as that seen in the A. thaliana mex1 mutant lacking the maltose transporter (MEX) that accumulates abnormally high levels of maltose and has reduced rates of starch degradation [1,34]. Maltose has also been reported to inhibit some βamylases at high concentration [35]. Alternatively, the reduction in starch degradation rate in the mex1 mutant might arise from multioligosaccharides inhibiting an enzyme involved in the attack on the starch granule, possibly by competing with granular starch for a starchbinding domain required for attack on the granule [1,36]. This is supported by experimental observations that the Arabidopsis dpe1 mutant, lacking the chloroplastic disproportionating enzyme required for maltotriose metabolism, also exhibits reduced starch degradation rate [1,20]. This mode of inhibition is outside the scope of the model, but is effectively captured by the reversibility of βamylase kinetics acting on soluble starch.
The model herein is based on the starch degradation pathway postulated by Smith, et al. [1], which suggests that starch granules are solubilised to yield soluble branched glucans that are then degraded by debranching and βamylase enzymes. The mechanism by which the solubilization occurs is not well understood and is likely to involve two dikinases—glucan water dikinase (GWD) [2,3] and phosphoglucan water dikinase (PWD) [4,5]. An alternative pathway might involve the direct attack on the starch granule by βamylase [1]. Although βamylase cannot hydrolyze linkages beyond branch points, it could act in tandem with a debranching enzyme to degrade starch granules gradually and directly to maltose and maltotriose. In such a case, the actions of the two dikinases GWD and PWD would determine the extent to which βamylase can attack chains at the granule surface, since the distribution of the phosphate groups added to amylopectin by these enzymes would reduce the degree of crystalline packing of chains inside the starch granules [37]. Such a pathway would result in formation of solutionphase maltooligosaccharides directly from insoluble starch granules without the intermediacy of soluble branched and linear glucans. Although earlier studies indicate that βamylase is incapable of degrading native starch granules [38,39], a chloroplastic βamylase from potato leaves was recently shown [40] to release maltooligosaccharides from potato tuber starch granules, which lends credibility to the alternative starch degradation pathway. If βamylase catalysis of insoluble starch cleavage is possible, the current model can be interpreted as capturing this in an effective way; however, the βamylase kinetics would need to be reexamined and likely reformulated to accurately capture all the subtlety of this more physically complex process.
A second related conclusion is that flow from soluble photosynthetically fixed carbon stores into metabolic pathways of interest in biofuel production is likely primarily limited by the cleavage of the linear polymer to oligomers, and not by subsequent reactions or debranching. This conclusion certainly depends on expression levels and characteristics of enzymes in particular cases, but is supported by our best estimate of biochemically relevant conditions herein, as well as the correlations among simulation endpoints representing a reasonably wide range of conditions (100 nM to 1 mM of 8 free initial and 6 fixed concentrations). Experimental investigations have suggested that the solubilisation of starch granules, rather than the hydrolysis of solubilised starch, might constitute the overall limiting step in starch degradation at low temperature [39,41,42]. Although potentially capturing such direct starch cleavage qualitatively as noted above, the model herein focuses on biochemical processes after starch solubilization. The βamylase rate limitation identified is thus relative to subsequent solution processes only. This conclusion is supported by experimental observations [1,43] that βamylase activity is strongly correlated to a decrease in starch during fruit ripening in banana plants, and that knockout mutants in A. thaliana lacking one of the chloroplastic βamylases show lowered rates of starch degradation.
A third observation is that transport reactions can serve as a kinetic bottleneck, as seen in the strong negative response coefficients of stromal glucose with respect to the glucose transporter k_{cat }(Figure 2) and concentration (Figure 8), and strong positive response coefficients versus K_{M }(Figure 2). Stromal glucose response stands in contrast to that of cytosolic glucose, which is insensitive to any parameters other than the hexokinase equilibrium constant. The latter behavior arises as a consequence of rapid equilibration of cytosolic glucose with the glucose 6phosphate pool mediated by hexokinase, which mimics rapid glucose flux into downstream carbon sinks. The relative sluggishness of transport can thereby dampen the sensitivity of metabolite concentrations in one compartment from the effects on reactions in another organelle. Such a conclusion is intuitive from topological considerations (two networks are connected by few edges, so requiring perturbations to propagate linearly through at least one reaction step), as well as kinetic ones (to the extent that transport is limiting, fast dynamics on one side of the reaction will not be visible to the other). One can surely devise exceptional cases, and a quantitative elucidation of this statement requires further exploration, but the results here lend support to the validity of a "divideandconquer" approach to cellular dynamical studies, with explicit consideration of single membranebound compartments or phases coupled by an effective variation of transport flux at the network boundaries, rather than explicitly by largescale, multicompartment dynamics. It should be noted that the irreversible formulation of transport we have employed naturally limits the degree to which this model can be generalized—situations without strong downstream cytosolic glucose demand would not be well represented, nor special cases in which dynamic glucose or maltose transients in the cytosol occur, since these effects could not be communicated through the transporters to the stromal metabolite pools.
A fourth conclusion may be summarized as a critique of what appears to be an implicit assumption that a single steady state resembling a handful of observations is necessarily the most important. To the degree that enzyme kinetics measured in vitro reflects turnover response to metabolite concentrations in vivo, kinetic models similar to that presented here are large dynamical systems linear in flux, but nonlinear in concentration. Such dynamical systems potentially possess great complexity, not only with respect to bifurcations as parameters vary, but also with respect to the phase space of concentrations under a single assumed set of parameters. A system infinitely "robust" with respect to temporally varying metabolite concentrations would indeed evolve toward a single unique steady state starting from anywhere in the relevant phase space dictated by catalyst concentrations and kinetic properties. The system explored herein behaves robustly with respect to the initial concentrations of internal metabolites, and with respect to perturbation of kinetic parameters and enzyme concentrations. Changing the latter for βamylase and MEX by up to 10fold resulted in quantitatively similar steadystate response coefficients. From a biological control perspective, this robustness is desired—variation of enzyme levels arising from genetic regulation or posttranslation modification should not give rise to catastrophic divergence of the cellular state. This robustness is also favorable for evolutionary or in vitro metabolic engineering, in that changing the nature of the phase space by mutagenic kinetic parameter variation will not lead to lethal cellular phenotypes. Nevertheless, given the complexity of biochemical regulation and chemical dynamics, we expect some surprises as biochemical models grow in fidelity, highperformance simulation and advanced analytical tools become available, and crossdisciplinary fertilization occurs between biochemists and mathematicians with interests in system dynamics.
Simulation of broadly sampled model concentrations and subsequent trajectory endpoint analysis showed correlation patterns consistent with the topology of the model metabolic network. This approach was found to be useful in identifying certain relationships, such as tight coupling of stromal metabolism, relative convergence values between classes of simulations (high and low starch, low and high ATP and orthophosphate), and sensitivity of internal metabolite concentrations (steadystate maltose concentration dependence on arabinogalactan and glucose levels in the cytosol). Although computational studies of largescale kinetic systems of the complexity found in cells is still in its infancy, the dual contemporary interests in parallel and dataintensive computing is opening the door to discovering unforeseen behaviors and patterns in biochemical networks. The work here has only touched on the full content of even this small metabolic model—simply exploring the appropriate cellular context by varying bath metabolite pool sizes, and potential multiplicity of fixed points by varying initial internal metabolite concentrations, generates a substantial quantity of data requiring significant analysis effort. Nevertheless, bringing metabolic engineering on par with traditional engineering disciplines will require thorough quantitative understanding of both system dynamics, and the effects of parametric variability beyond values that nature provides. This transformation will be facilitated by further development and adaptation of analytical and visualization methods for biochemical systems analysis.
Conclusions
Construction and characterization of a kinetically detailed model of starch metabolism shows that βamylase activity is the limiting factor in saccharide production under conditions of high glucose demand, using best estimates for kinetic parameters and enzyme levels. Sensitivity analyses and sampling of internal and external metabolite concentrations and clustering analysis of fixedtime simulation endpoints showed that soluble starch levels are the main determinant of debranched starch and maltooligomer levels in the chloroplast stroma, but that transport reactions partially decouple the cytosolic chemical subsystem directing carbon flow to downstream sinks. The most tightly converged end points illustrate a role for the metabolically coupled arabinogalactan/glucosylated arabinogalactan pair and cytosolic glucose levels in determining steady state maltose levels. No evidence for multistationarity was found. The model and explorations described highlight areas in starch metabolism for deeper study and experimental testing, as well as potential opportunities for methodological advancement.
Methods
Model Formulation: Nomenclature of Metabolites, Enzymes, Transporters, and Parameters
Compound abbreviations are defined in Tables 1 and 2. The intracellular compartment in which a compound or enzyme resides is denoted by a twoletter suffix preceded by the underscore sign, with "_CY" denoting cytosolic and "_CS" a chloroplast stromal species. Where ionization is possible, pools of ionizable species containing all the biologically occurring ionized and unionized forms are appended with a "tot" subscript. For example, the pool of phosphate in the cytosol is represented as Pitot_CY. For polymeric starch, when a single residue of the polymer—the starch glucosyl unit—constitutes a separate model entity, the residue identity is represented as GlcStarch_CS. The aggregate pool of linear linkage groups released from solubilized starch by the action of the debranching enzyme is assigned the abbreviation Starch_{db}_CS. The enzymes are usually represented using their Enzyme Commission (EC) numbers, such that the enzyme names consist of the prefix "ec_" followed by the EC number with the dots substituted with underscores. Thus, the chloroplastic form of the disproportionating enzyme 1 DPE1 with EC number 2.4.1.25 is therefore represented as ec_2_4_1_25_CS.
Kinetic parameters are generally referred to with appropriate formatting in the text, e.g., βamylase for the turnover number of βamylase acting on starch as a substrate. For figure labels, an alternative nomenclature was used for simplicity in formatting. Thus, βamylase is referred to as "betaAmylase_Gn_kcat". Differing parameters for alternative substrates are denoted by superscripted parentheticals, so the k_{cat }parameters for starch (Gn) and maltopentaose (G5) degradation by βamylase are denoted by and , respectively. The names of the enzyme kinetic parameters and equilibrium constants for the cytosolic disproportionation catalyzed by DPE2 all start with the locus tag for this enzyme in Arabidopsis thaliana, AT2G40840, rather than the KEGG [44] reaction identifier.
To treat transport reactions between compartments of different volumes consistently, all reaction rate equations carry the appropriate volumetric factors in the model (see Additional File 2). Thus, each rateofchange is calculated as a mass rateofchange rather than a concentration rateofchange, consistent with the SBML standard [45]. However, we have excluded these factors in the Tables to be more consistent with standard biochemical nomenclature.
Additional file 2. SBML Model File. A single Systems Biology Markup Language file representing the reference model in this study. This model is automatically transformed to C++ by the translation utility of the HighPerformance Systems Biology Toolkit and compiled into highperformance executable programs for sampling and optimization tasks, as well as simple forward integration. This model may also be found in the BioModels database under accession number MODEL1106030000.
Format: XML Size: 71KB Download file
Model Formulation: Biochemical Processes
Degradation of starch to maltose and maltotriose
Degradation of solubilized starch to maltose and maltotriose in the chloroplast stroma is modeled by a modified version of the kinetic expression proposed by Shiraishi, et al., to describe maltose production from soluble starch by the concerted use of βamylase and debranching enzymes [29]. The total mass concentration of starch is divided into two parts as seen in Figure 10 of Ref. [29], namely that hydrolyzable solely by βamylase S_{β0 }and that requiring both βamylase and a debranching enzyme for hydrolysis, S_{df}. Using the representation scheme in this figure, the total mass concentration hydrolyzed solely by βamylase is
and the total mass concentration hydrolyzed by the combined action of βamylase and debranching enzyme is given by
f_{β }= 0.582 is defined as the mass fraction of starch that can be degraded by the action of βamylase alone.
The rate equations of maltose and maltotriose formation resulting from starch degradation are given in Table 5. The mass concentration of linear linkage groups released from starch produced by the action of a debranching enzyme is defined as S_{db}, the time evolution of which was modeled after [29] and is also detailed in Table 5. The initial concentration of such groups is set to 0 at t = 0. Fractions f_{M }= 0.87 and f_{G3 }= 0.13 represent starch ultimately convertible to maltose and maltotriose, respectively.
Table 5. Kinetic equations and parametric assignments for βamylase.
Disproportionation reactions
The disproportionation of two maltotriose molecules to glucose and maltopentaose is catalyzed by disproportionation enzyme 1 (DPE1). The relevant equations and parameters for this biochemical reaction adapted to a single substrate pool are shown in Table 6. The cytosolic transglucosidase (DPE2) reaction splits maltose into two glucosyl units, one of which ends up as free glucose and the other transferred to a heteroglycan (arabinoglycan) acceptor to yield glucosylated arabinogalactan. Both reactions were modeled with PingPong BiBi kinetics [30], because many enzymatic disproportionation reactions have been reported in the literature to proceed by this mechanism [4648]. Kinetic equations and parameters are given in Tables 6 and 7.
Table 6. Kinetic equations for maltotriose disproportionation to glucose and maltopentaose by DPE1.
Table 7. Kinetic equations for maltose disproportionation to glucose and glucosylated arabinogalactan by DPE2.
Degradation of maltopentaose by βamylase
The degradation of maltopentaose formed by the chloroplastic disproportionation reaction is modeled as
where G_{5 }denotes the maltopentaose mass concentration, E_{β }represents βamylase mass concentration, and and represent the turnover number of the βamylase enzyme for maltopentaose degradation and the Michaelis constant for maltopentaose, respectively [29].
Transport reactions
The present model includes two catalyzed transport processes from chloroplast stroma to cytosol. Both are modeled as irreversible MichaelisMenten processes,
with S being either chloroplastic maltose or glucose, and V_{max }factorable as k_{cat }× [transporter]. The value of is calculated as 11.9 μM s^{1 }using the maximal maltose consumption rate for a given residual maltose concentration and the mass concentration of cells in anaerobic sugarlimited yeast (Saccharomyces cerevisiae CBS 8066) chemostat cultures [49]. Based on a value for of 519 μmole (mg of chlorophyll hr)^{1 }[21] and the equivalence between 1 mg of chlorophyll to 30 μL of chloroplast stroma from [50], we obtain = 4806 μM s^{1}_{. }Parameters are given in Table 8.
Table 8. Kinetic factors for maltose and glucose transport between chloroplast stroma and cytosol.
Release of Glc1P from glucosylated heteroglycan and hexokinase The cytosolic glucan phosphorylase (glycosylated heteroglycan → glucose1phosphate) was modeled as a rapidequilibrium random BiBi mechanism [30,51]. The velocity expression for this reaction is given in Table 9. Cytosolic hexokinase forming glucose6phosphate from glucose and ATP is modeled similarly using the same velocity equation and enzyme kinetic and binding parameters as in the model of human erythrocyte hexokinase developed by Mulquiney and Kuchel [51].
Table 9. Kinetic equation and factors for cytosolic glucan phosphorylase.
The rateofchange expressions for individual metabolites in terms of reaction fluxes are contained in Additional File 3.
Additional file 3. Summary graphics and enzyme kinetic models. Integration timecourse for the reference model, closeup of debranched stromal starch evolution, a graphical summary of the response coefficients and norms for the 8 states detailed in Table 8 and rateofchange equations in terms of reaction fluxes.
Format: DOCX Size: 100KB Download file
Simulation Framework and Methodology: Model instantiation and simulation
The starch degradation model was expressed in Systems Biology Markup Language (SBML) [45] and input into our metabolic simulation software framework, the HighPerformance Systems Biology Toolkit (HiPer SBTK) [52]. The time evolution of metabolite concentrations and fluxes was simulated, and the possible space of rate and binding parameters sampled. Biologically feasible concentrations of internal metabolites and fixed concentrations of external metabolites, enzymes, transporters and inhibitors were imposed, then the ODE system integrated. Sensitivities of steadystate internal metabolite concentrations to enzyme concentrations and kinetic parameter values were then computed. Additionally, the structure of the dynamical space was explored by sampling initial concentrations of internal metabolites given a fixed set of enzyme and inhibitor concentrations and kinetic parameter values, and evaluating the distance between steady states so achieved.
A reference system was defined by starting from a best estimate of biologically relevant metabolite, enzyme, transporter and inhibitor concentrations (Tables 1 and 2). This model was integrated for 10^{7 }virtual seconds to a convergence of 1.6 × 10^{13 }(μM/s) μM^{1}. Response coefficients [31] of the internal metabolite steadystate concentrations S_{i }with respect to the kinetic and binding parameters and fixed enzyme and transporter protein concentrations p_{j }were calculated. In addition, "parameterperturbed" and "enzymeperturbed" systems were modeled by decreasing and 10fold relative to the reference model in the former case, and increasing βamylase concentration twofold and the maltose exporter concentration tenfold in the latter. Integrations of the parameterand enzymeperturbed models converged to 1.5 × 10^{16 }and 5.2 × 10^{12 }(μM/s) μM^{1}, respectively.
Additional File 3 contains the integration timecourse for the reference model, closeup of debranched stromal starch evolution, and a graphical summary of the response coefficients and norms for the 8 states detailed in Table 3.
Analysis of States from Varying External or Initial Internal Metabolite Concentrations
To evaluate both the existence of multiple stationary states within individual dynamic systems, and the effects of external metabolite concentrations on the fixed point(s) of the modeled biochemical network, 6 external and 8 initial internal metabolite concentrations were sampled at two values, and in each of the 2^{14 }cases integrated as described above. Each ODE integration was run for 10^{7 }virtual seconds, and convergence of each trajectory quantified as
where [S_{i}] is the concentration of metabolite i, and C carries units of (μM/s) μM^{1}. Final metabolite vectors at convergence cutoffs of 10^{2 }and 10^{8 }were clustered by the kmeans algorithm (as implemented in the SciPy package). Each row of the data matrix comprised a final point in the logarithm of metabolite concentration space. The variance of each column was normalized to 1.0, then the data clustered with a seed cluster estimate of 64 (2^{6}, for two sampled fixed concentrations of six external metabolites). Cluster centroid distances were checked to be greater than the sum of their respective cluster standard deviations. Establishing unit variance for the most tightly converged metabolite vectors was problematic, and so data in this case were clustered sequentially in the direct space of concentrations, using a distance threshold of 10^{8 }μM. For this latter case, a cluster could be defined by the position of the first point added, rather than the centroid, due to the relatively tight bunching observed.
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
AN developed and encoded the starch degradation model, ensured model selfconsistency, simulated the model behavior, analyzed results from serial and parallel simulations, and drafted the manuscript. ML implemented the parallel parameter sampling routines and carried out parallelized simulations. PG wrote the sequential clustering routine, supervised the software design and implemented the core simulation components. CHC developed the modeling strategy, wrote and ran the kmeans clustering routine and analyzed the results, codrafted and edited the manuscript, and supervised the research. All authors have read and approved the final manuscript.
Authors' information
AN is a research associate in the Computational Science Center at the National Renewable Energy Laboratory (NREL) with research interests including mathematical biology and theoretical immunology. His current explorations include metabolome representation and expression for potential biofuelproducing organisms, enzymatic hydrolysis of biomass particles containing multiple polysaccharide types, flux balance analysis, and modeling regulatory networks using datamining, bioinformatics analysis and correlation of multiple types of highthroughput "omics" data.
ML is a research associate in the Computational Science Center at NREL. His research interests are numerical search, high performance computing (HPC) and uncertainty quantification. In addition to supporting systems biology, Monte is also working on exploring more efficient building designs using parallel computing.
PG is a senior scientist in the Computational Science Center at NREL. His research interests center on the application of mathematics and HPC to renewable energy challenges, particularly optimizations involving expensive simulations. In addition to developing kinetic simulation software for HPC, he studies multiscale device simulation and optimization and highthroughput materials screening, inverse material design by electronic structure data mining and optimization, and computational battery modeling.
CHC is a senior scientist in the Computational Science Center at NREL. His research focuses on systems biochemistry and computational chemistry from molecular to cellular scales. His current projects touch on dynamical properties of chemical and biochemical systems, electron and energy transfer, and homogeneous transition metal catalysis.
Acknowledgements
This work was supported by the U.S. Department of Energy's Office of Science (DOESC) through the Scientific Discovery through Advanced Computing (SciDAC) program, the Office of Biological and Environmental Research, and the Office of Advanced Scientific Computing Research under contract number DEAC3608GO28308.
References

Smith AM, Zeeman SC, Smith SM: Starch Degradation.
Annu Rev Plant Biol 2005, 56:7398. PubMed Abstract  Publisher Full Text

Mikkelsen R, Baunsgaard L, Blennow A: Functional Characterization of αGlucan, Water Dikinase, the Starch Phosphorylating Enzyme.
Biochem J 2004, 377:525532. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Ritte G, Lloyd JR, Eckermann N, Rottmann A, Kossmann J, Steup M: The Starchrelated R1 Protein Is an AlphaGlucan, Water Dikinase.
Proc Natl Acad Sci USA 2002, 99:71667171. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Baunsgaard L, Lutken H, Mikkelsen R, Glaring MA, Pham TT, Blennow A: A Novel Isoform of Glucan, Water Dikinase Phosphorylates PrePhosphorylated AlphaGlucans and Is Involved in Starch Degradation in Arabidopsis.
Plant J 2005, 41:595605. PubMed Abstract  Publisher Full Text

Kotting O, Pusch K, Tiessen A, Geigenberger P, Steup M, Ritte G: Identification of a Novel Enzyme Required for Starch Metabolism in Arabidopsis Leaves. The Phosphoglucan, Water Dikinase.
Plant Physiol 2005, 137:242252. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Ball SG, Deschamps P: Starch Metabolism. In The Chlamydomonas Sourcebook. Volume 2. Edited by Stern DB, Harris EH. Amsterdam: Academic Press; 2009::140.

Nakamura Y: Some Properties of Starch Debranching Enzymes and Their Possible Role in Amylopectin Biosynthesis.
Plant Sci 1996, 121:118. Publisher Full Text

Lin TP, Spilatro SR, Preiss J: Subcellular Localization and Characterization of Amylases in Arabidopsis Leaf.
Plant Physiol 1988, 86:251259. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zeeman SC, Northrop F, Smith AM, Rees T: A Starchaccumulating Mutant of Arabidopsis thaliana Deficient in a Chloroplastic Starchhydrolysing Enzyme.
Plant J 1998, 15:357365. PubMed Abstract  Publisher Full Text

Hausler RE, Schlieben NH, Schulz B, Flugge UI: Compensation of Decreased Triose Phosphate/Phosphate Translocator Activity by Accelerated Starch Turnover and Glucose Transport in Transgenic Tobacco.
Planta 1998, 204:366376. PubMed Abstract  Publisher Full Text

Zeeman SC, Thorneycroft D, Schupp N, Chapple A, Weck M, Dunstan H, Haldimann P, Bechtold N, Smith AM, Smith SM: Plastidial alphaGlucan Phosphorylase Is Not Required for Starch Degradation in Arabidopsis Leaves but Has a Role in the Tolerance of Abiotic Stress.
Plant Physiol 2004, 135:849858. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lao NT, Schoneveld O, Mould RM, Hibberd JM, Gray JC, Kavanagh TA: An Arabidopsis Gene Encoding a Chloroplasttargeted BetaAmylase.
Plant J 1999, 20:519527. PubMed Abstract  Publisher Full Text

Hehre EJ, Okada G, Genghof DS: Configurational Specificity: Unappreciated Key to Understanding Enzymic Reversions and de Novo Glycosidic Bond Synthesis. I. Reversal of Hydrolysis by Alpha, Betaand Glucoamylases with Donors of Correct Anomeric Form.
Arch Biochem Biophys 1969, 135:7489. PubMed Abstract

Niittyla T, ComparotMoss S, Lue WL, Messerli G, Trevisan M, Seymour MDJ, Gatehouse JA, Villadsen D, Smith SM, Chen JC, et al.: Similar Protein Phosphatases Control Starch Metabolism in Plants and Glycogen Metabolism in Mammals.
J Biol Chem 2006, 281:1181511818. PubMed Abstract  Publisher Full Text

Chia T, Thorneycroft D, Chapple A, Messerli G, Chen J, Zeeman SC, Smith SM, Smith AM: A Cytosolic Glucosyltransferase Is Required for Conversion of Starch to Sucrose in Arabidopsis Leaves at Night.
Plant J 2004, 37:853863. PubMed Abstract  Publisher Full Text

Boos W, Shuman H: Maltose/Maltodextrin System of Escherichia coli: Transport, Metabolism, and Regulation.
Microbiol Mol Biol Rev 1998, 62:204229. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lloyd JR, Blennow A, Burhenne K, Kossmann J: Repression of a Novel Isoform of Disproportionating Enzyme (stDPE2) in Potato Leads to Inhibition of Starch Degradation in Leaves but Not Tubers Stored at Low Temperature.
Plant Physiol 2004, 134:13471354. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Fettke J, Eckermann N, Poeste S, Pauly M, Steup M: The Glycan Substrate of the Cytosolic (Pho 2) Phosphorylase Isozyme from Pisum sativum L.: Identification, Linkage Analysis and Subcellular Localization.
Plant J 2004, 39:933946. PubMed Abstract  Publisher Full Text

Yang Y, Steup M: Polysaccharide Fraction from Higher Plants which Strongly Interacts with the Cytosolic Phosphorylase Isozyme: I. Isolation and Characterization.
Plant Physiol 1990, 94:960969. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Critchley JH, Zeeman SC, Takaha T, Smith AM, Smith SM: A Critical Role for Disproportionating Enzyme in Starch Breakdown Is Revealed by a Knockout Mutation in Arabidopsis.
Plant J 2001, 26:89100. PubMed Abstract  Publisher Full Text

Weber A, Servaites JC, Geiger DR, Kofler H, Hille D, Groner F, Hebbeker U, Flugge UI: Identification, Purification, and Molecular Cloning of a Putative Plastidic Glucose Translocator.
Plant Cell 2000, 12:787802. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Palsson BØ: Systems Biology: Properties of Reconstructed Networks. New York: Cambridge University Press; 2006.

Lee JM, Gianchandani EP, Papin JA: Flux Balance Analysis in the Era of Metabolomics.
Brief Bioinf 2006, 7:140150. Publisher Full Text

Beard DA, Babson E, Curtis E, Qian H: Thermodynamic Constraints for Biochemical Networks.
J Theor Biol 2004, 228:327333. PubMed Abstract  Publisher Full Text

Beard DA, Liang Sd, Qian H: Energy Balance for Analysis of Complex Metabolic Networks.
Biophys J 2002, 83:7986. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Teusink B, Passarge J, Reijenga CA, Esgalhado E, van der Weijden CC, Schepper M, Walsh MC, Bakker BM, van Dam K, Westerhoff HV, Snoep JL: Can Yeast Glycolysis Be Understood in Terms of in Vitro Kinetics of the Constituent Enzymes? Testing Biochemistry.
Eur J Biochem 2000, 267:53135329. PubMed Abstract  Publisher Full Text

Beard DA: A Biophysical Model of the Mitochondrial Respiratory System and Oxidative Phosphorylation.
PLoS Comp Biol 2005, 1:e36. Publisher Full Text

Wu F, Yang F, Vinnakota KC, Beard DA: Computer Modeling of Mitochondrial Tricarboxylic Acid Cycle, Oxidative Phosphorylation, Metabolite Transport, and Electrophysiology.
J Biol Chem 2007, 282:2452524537. PubMed Abstract  Publisher Full Text

Shiraishi F, Kawakami K, Yuasa A, Kojima T, Kusunoki K: Kinetic Expression for Maltose Production from Soluble Starch by Simultaneous Use of betaAmylase and Debranching Enzymes.
Biotechnol Bioeng 1987, 30:374380. PubMed Abstract  Publisher Full Text

Segel IH: EnzymeKinetics. New York: John Wiley & Sons; 1976.

Klipp E, Herwig R, Kowald A, Wierling C, Lehrach H: Systems Biology in Practice: Concepts, Implementation and Application. WileyVCH Verlag GmbH & Co; 2005.

Mita S, Hirano H, Nakamura K: Negative Regulation in the Expression of a SugarInducible Gene in Arabidopsis thaliana. A Recessive Mutation Causing Enhanced Expression of a Gene for βAmylase.
Plant Physiol 1997, 114:575582. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Laby RJ, Kim D, Gibson SI: The ram1 Mutant of Arabidopsis Exhibits Severely Decreased βAmylase Activity.
Plant Physiol 2001, 127:17981807. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Niittyla T, Messerli G, Trevisan M, Chen J, Smith AM, Zeeman SC: A Previously Unknown Maltose Transporter Essential for Starch Degradation in Leaves.
Science 2004, 303:8789. PubMed Abstract  Publisher Full Text

Lizotte PA, Henson CA, Duke SH: Purification and Characterization of Pea Epicotyl βAmylase.
Plant Physiol 1990, 92:615621. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Witt W, Sauter JJ: Purification and Properties of the Starch Granuledegrading αAmylase from Potato Tubers.
J Exp Bot 1996, 47:17891795. Publisher Full Text

Blennow A, Nielsen TH, Baunsgaard L, Mikkelsen R, Engelsen SB: Starch Phosphorylation: A New Front Line in Starch Research.
Trends Plant Sci 2002, 7:445450. PubMed Abstract  Publisher Full Text

Beck E, Ziegler P: Biosynthesis and Degradation of Starch in HigherPlants.
Annu Rev Plant Physiol Plant Mol Biol 1989, 40:95117. Publisher Full Text

Dunn G: Model for Starch Breakdown in HigherPlants.
Phytochem 1974, 13:13411346. Publisher Full Text

Scheidig A, Frohlich A, Schulze S, Lloyd JR, Kossmann J: Downregulation of a Chloroplasttargeted betaAmylase Leads to a Starchexcess Phenotype in Leaves.
Plant J 2002, 30:581591. PubMed Abstract  Publisher Full Text

Yano R, Nakamura M, Yoneyama T, Nishida I: Starchrelated alphaGlucan/Water Dikinase Is Involved in the Coldinduced Development of Freezing Tolerance in Arabidopsis.
Plant Physiol 2005, 138:837846. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Colonna P, Leloup V, Buleon A: Limiting Factors of Starch Hydrolysis.
Eur J Clin Nutr 1992, 46:S17S32. PubMed Abstract

do Nascimento JRO, Vieira Junior A, Bassinello PZ, Cordenunsi BR, Mainardi JA, Purgatto E, Lajolo FM: BetaAmylase Expression and Starch Degradation during Banana Ripening.
Postharvest Biol Technol 2006, 40:4147. Publisher Full Text

Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, Katayama T, Kawashima S, Okuda S, Tokimatsu T, Yamanishi Y: KEGG for Linking Genomes to Life and the Environment.
Nucl Acids Res 2008, 36:D480484. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hucka M, Finney A, Sauro HM, Bolouri H, Doyle JC, Kitano H, Arkin AP, Bornstein BJ, Bray D, CornishBowden A, et al.: The Systems Biology Markup Language (SBML): A Medium for Representation and Exchange of Biochemical Network Models.
Bioinf 2003, 19:524531. Publisher Full Text

Nakamura A, Haga K, Yamane K: The Transglycosylation Reaction of Cyclodextrin Glucanotransferase Is Operated by a PingPong Mechanism.
FEBS Lett 1994, 337:6670. PubMed Abstract  Publisher Full Text

Everett RR, Soedjak HS, Butler A: Mechanism of Dioxygen Formation Catalyzed by Vanadium Bromoperoxidase. Steady State Kinetic Analysis and Comparison to the Mechanism of Bromination.
J Biol Chem 1990, 265:1567115679. PubMed Abstract  Publisher Full Text

SauraValls M, Faure R, Ragas S, Piens K, Brumer H, Teeri TT, Cottaz S, Driguez H, Planas A: Kinetic Analysis Using LowMolecular Mass Xyloglucan Oligosaccharides Defines the Catalytic Mechanism of a Populus Xyloglucan Endotransglycosylase.
Biochem J 2006, 395:99106. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Weusthuis RA, Adams H, Scheffers WA, van Dijken JP: Energetics and Kinetics of Maltose Transport in Saccharomyces cerevisiae: A Continuous Culture Study.

Zhu XG, de Sturler E, Long SP: Optimizing the Distribution of Resources between Enzymes of Carbon Metabolism Can Dramatically Increase Photosynthetic Rate: A Numerical Simulation Using an Evolutionary Algorithm.
Plant Physiol 2007, 145:513526. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Mulquiney PJ, Kuchel PW: Model of 2,3Bisphosphoglycerate Metabolism in the Human Erythrocyte Based on Detailed Enzyme Kinetic Equations: Equations and Parameter Refinement.
Biochem J 1999, 342(Pt 3):581596. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Chang CH, Graf P, Alber DM, Kim K, Murray G, Posewitz M, Seibert M: Photons, Photosynthesis, and HighPerformance Computing: Challenges, Progress, and Promise of Modeling Metabolism in Green Algae.

Tewari YB, Goldberg RN, Sato M: Thermodynamics of the Hydrolysis and Cyclization Reactions of α, β, and γCyclodextrin.
Carb Res 1997, 301:1122. Publisher Full Text

Kakefuda G, Duke SH: Characterization of Pea Chloroplast DEnzyme (4alphadGlucanotransferase).
Plant Physiol 1989, 91:136143. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Leemhuis H, Dijkstra BW, Dijkhuizen L: Thermoanaerobacterium thermosulfurigenes Cyclodextrin Glycosyltransferase.
Eur J Biochem 2003, 270:155162. PubMed Abstract  Publisher Full Text

Alberty RA: Thermodynamics of Biochemical Reactions. Hoboken, New Jersey: John Wiley & Sons; 2003.

Gold AM, Johnson RM, Sanchez GR: Kinetic Mechanism of Potato Phosphorylase.
J Biol Chem 1971, 246:34443450. PubMed Abstract  Publisher Full Text

Vinnakota K, Kemp ML, Kushmerick MJ: Dynamics of Muscle Glycogenolysis Modeled with pH Time Course Computation and pHdependent Reaction Equilibria and Enzyme Kinetics.
Biophys J 2006, 91:12641287. PubMed Abstract  Publisher Full Text  PubMed Central Full Text