Modeling the Calvin-Benson cycle
© Jablonsky et al; licensee BioMed Central Ltd. 2011
Received: 7 September 2011
Accepted: 3 November 2011
Published: 3 November 2011
Skip to main content
© Jablonsky et al; licensee BioMed Central Ltd. 2011
Received: 7 September 2011
Accepted: 3 November 2011
Published: 3 November 2011
Modeling the Calvin-Benson cycle has a history in the field of theoretical biology. Anyone who intends to model this system will look at existing models to adapt, refine and improve them. With the goal to study the regulation of carbon metabolism, we investigated a broad range of relevant models for their suitability to provide the basis for further modeling efforts. Beyond a critical analysis of existing models, we furthermore investigated the question how adjacent metabolic pathways, for instance photorespiration, can be integrated in such models.
Our analysis reveals serious problems with a range of models that are publicly available and widely used. The problems include the irreproducibility of the published results or significant differences between the equations in the published description of the model and model itself in the supplementary material. In addition to and based on the discussion of existing models, we furthermore analyzed approaches in PGA sink implementation and confirmed a weak relationship between the level of its regulation and efficiency of PGA export, in contrast to significant changes in the content of metabolic pool within the Calvin-Benson cycle.
In our study we show that the existing models that have been investigated are not suitable for reuse without substantial modifications. We furthermore show that the minor adjacent pathways of the carbon metabolism, neglected in all kinetic models of Calvin-Benson cycle, cannot be substituted without consequences in the mass production dynamics. We further show that photorespiration or at least its first step (O2 fixation) has to be implemented in the model if this model is aimed for analyses out of the steady state.
There are two common approaches used in modeling the Calvin-Benson cycle: kinetic modeling, e.g., [2, 3] and stoichiometric modeling, e.g., [4–6]. Kinetic modeling requires to obtain/have available kinetic properties of the enzymes involved. These are mostly known if one assumes conservative kinetic parameters among the species with several different ways of CO2 fixation. On the other hand, stoichiometric modeling allows an analysis of the entire carbon metabolism within the given constraints and without the need of kinetic properties. It would seem plausible to combine both approaches and overlap fluxes from the kinetic model of the Calvin-Benson cycle with the stoichiometric analysis of the metabolic network [7, 8].
Even though it is a core part in a variety of photosynthetic models or projects modeling the entire photosynthesis [9, 10], kinetic modeling has mostly focused on the stability of the Calvin-Benson cycle itself and continues to attract considerable scientific interest [11, 12]. To anyone entering this field, it is clear that the Calvin-Benson cycle is the best-studied plant metabolic system. Since modeling the Calvin-Benson cycle has a long history, the question arises how one should proceed with a further analysis of this system. The standard scientific approach is to build on existing knowledge and models. A natural progression would thus to adapt existing models, to refine and expand them for adjacent metabolic pathways.
What our present work reveals is that, without significant modifications and corrections of the well-established models, it is not possible to re-use existing models. Our analysis shows that one of the existing and widely used models produces an "accumulation" of starch in high negative concentration with consequences to the behavior of the whole system (apart from being biologically implausible). In another case an interchanged parameter leads to 100-times faster equilibrium in one particular reaction, which significantly shifted the steady states of all metabolites. We furthermore found errors, which have propagated from one generation (publication) to another. Our study highlights difficulties with the reproducibility of model-based results, provides solutions and suggestions on how to proceed with extensions to existing models of the Calvin-Benson cycle.
In order to support the evolution of models for the Calvin-Benson cycle and to improve the reproducibility of model-driven results, we discuss ways to improve and correct existing models. We also open a discussion about how to connect the Calvin-Benson cycle with other parts of the carbon metabolism in the case of cyanobacteria and in relation to the system's efficiency, regulation, mass production and implementation of photorespiration.
To investigate mathematical and computational models for their suitability and possible errors all models were encoded both in Matlab, using the SimBiology toolbox (MathWorks Inc., Natick, Massachusetts, USA) and at the same time in Copasi . For transparency and to allow readers to use corrected models of the Calvin-Benson cycle, we make all files and data required for the simulation publicly available in additional files (.xml).
There are studies that successfully use a simplified description of the Calvin-Benson cycle, for instance the photosynthetic oscillations in the chlorophyll fluorescence . However, if one analyses the regulation of carbon metabolism, such study cannot be done without a detailed kinetic model of the Calvin-Benson cycle. Since kinetic modeling of the Calvin-Benson cycle has a long tradition [2, 3, 15–22], the standard approach would be to adapt, refine and expand an existing model. Focusing on the Calvin-Benson cycle and carbon metabolism, we have investigated several well-established models of this system, which are discussed in further detail below.
This fundamental model was presented by Hahn  who also published further related models of the Calvin-Benson cycle [16, 17]. The kinetic parameters in all models were chosen on the basis of a realistic photosynthetic rate in order to give reasonable steady-state values for metabolites. One would assume that the earlier model  was extended by photorespiration in the next model generation ; together with the realistic response to the changes in the O2 and CO2 concentrations . This assumption is based on the fact that the equations describing the Calvin-Benson cycle are the same in both models [15, 17]; however, some values of the kinetic parameters differ significantly.
The most prominent example is the rate of the CO2 fixation expressed by the rate constant k1. The value of this calculated parameter is in both models identical [15, 17], i.e. on the basis of estimated gross photosynthetic rate per unit area of leaf tissue . Since the is the same in both models, it is not clear why value of k1 is 57.3 - fold higher in the later model . Hahn mentions that the value of k1 was, in order to compensate the neglected photorespiration in the earlier model , decreased (without any detail regarding the original value of k1). However, it is unlikely that the original and unknown value of k1 in the earlier model  is 57.3 - fold lower due to photorespiration because the author himself assumed in the later model  the ratio for the carboxylase/oxygenase activity of RuBisCO (Ribulose-1,5-bisphosphate carboxylase oxygenase) equal to 3.7.
Nevertheless, we have encoded the earlier Hahn model  into biochemical simulators according its original description. Having only the rate constants, we employed mass action kinetics as the author did . For initial concentrations we used the steady-state values calculated by Hahn .
The reason for the drop in the ATP level in the certain time point (Figure 2) was a depletion of both intra- and extracellular phosphates (data not shown). Finally, we achieved a stable solution by fixing (having constant concentration) the Pi_ext (extracellular phosphate). This approach prevents the suspension of the cycle (Figure 2) but the steady-state of the system is still far away from known values (e.g., the aforementioned range for PGA) or results presented by author (Figure 2 in ).
Although it is not a problem to reconstruct the established fundamental models from their description, a problem is that our simulations based on the Hahn and Pettersson models demonstrate totally different behaviors in comparison to what was presented by their authors [2, 15]. It is not clear what the problem is because the computational methodology used for running the simulations has not been described in the publications. But if one employs well established computational methods and tools (e.g., Matlab), the original results [2, 15] are not reproducible. Therefore, if one wants to study the cycle itself or expand the model by adding other part(s) of the carbon metabolism, the reasonable cause of action would be to look for available models for which one knows all details about the computational methods. Such models can be nowadays found in model databases (e.g., http://www.ebi.ac.uk/biomodels-main) or they are made available as a supplement to a publication. For the Calvin-Benson cycle, there are two main options, the Poolman model  and the Zhu model .
The Poolman model is available as a .xml file (http://www.ebi.ac.uk/biomodels-main/) and can be therefore easily imported to any biochemical simulator supporting SBML (system biology markup language) level 2. We have employed both Copasi  and the SimBiology toolbox for Matlab (MathWorks) to study this model. The Poolman model is structurally very close to the Pettersson model  but with an altered kinetic parameters. One can speculate that the author encountered the same problems with Pettersson model as we have and as a result, employed different set of kinetic parameters. This model was used as an example to illustrate possible applications in the metabolic modeling  and stability analysis (e.g., number of steady states) rather  than for kinetic modeling of the particular metabolite(s). The initial conditions were exactly the same as in the original work .
Pilot simulations based on the Poolman model showed that the system reached the steady state at 0.6 s from the start of the simulation. If the initial concentration of PGA (3-phosphoglycerate) is increased 10-times, the steady state is reached after 1.4 s. This feature of the model is based on the assumption that metabolites in the system are maintained at a fast equilibrium . Furthermore, to ensure the assumed fast equilibrium, rate constants in the Poolman model were set to extremely high values. For example, the rate constants for reactions catalyzed by transketolase (equations E7 and E10, ) are equal to 500 000 000 l·mmol -1·s-1  and the flux through these reactions is 39.7 mmol·s-1 (our calculation based on the unmodified Poolman model ). This flux is in sharp contrast to the proposed (0.1 mmol·s-1 ) and measured (0.36 mmol·s-1 ) maximal rates for these reactions.
The Poolman model reaches its equilibrium at 0.6 s of the simulation also in other cases (after dark → light switch or changes in CO2 level). On the other hand, it is known that the real system is limited by RuBisCO (ribulose-1,5-bisphosphate carboxylase/oxygenase) activation whose induction phase takes 120 s  - 600 s . We note that the activation of RuBisCO was also considered in a kinetic model with the time constant equals to 250 s .
The Calvin-Benson cycle is an open system with many inputs and outputs (Figure 1). It is therefore important to analyze its output behavior as well. To that end, we have changed the settings for sinks and starch in the Poolman model from the fixed (constant) concentration to variable concentration (responding according the reactions and differential equations). We note that the original setting for other fixed metabolites, i.e., for inputs (CO2, cytosolic Pi, NADPH, NADP+ and H+) was not changed.
The Poolman model was designed for other purposes (stability analysis) but it is clear that this model cannot be used for the purpose of comparison with experimental time-series data without changing both the equations and kinetic parameters.
The remaining option is the Zhu model . Zhu and coworkers presented a model of C3 photosynthesis, extending the Laisk model  by addition of the photorespiratory pathway. Since we do not focus on the complete photorespiratory pathway in this study, we did not analyze the Laisk model separately, except for the difference in the case of phosphate translocator (see the section Calvin-Benson cycle as a part of the metabolic network: PGA sink implementation).
The Zhu model is a natural starting point and suitable template for studying photorespiration, even if it is not easily accessible and understandable because it is encoded in 34 Matlab files. The metabolites in the Zhu model reach the equilibrium at 150 s of the simulation (data not shown), which is in agreement with previous findings  and in contrast to 0.6 s based on the Poolman model . However, several discrepancies and serious errors in the Zhu model, particularly in the part describing the Calvin-Benson cycle, can be found. The discrepancies in the Zhu model include interchanged values of kE13 and kI135 in Appendix C and the nowhere used constant KM8, make the description of the model slightly confusing. The following section is dedicated to the analysis of the errors occurring in the Zhu model, as well as necessary modifications for this model.
The programming environments used for kinetic modeling of the Calvin-Benson cycle include Turbo Pascal, used in ; Scampi, used in  or Matlab, used in . Command-line programming is more prone to errors and also makes it more likely for errors to propagate into new models. For example, in the Laisk model  one error propagated into other models due the missing dimensional analysis in Pascal. This problem can also be found in the last generation of the model containing the Calvin-Benson cycle, in the Zhu model . The following problems occur in the Zhu model:
Two very different descriptions for reactions v7: S7P + GAP = Ri5P + Xu5P and v10: F6P + GAP = E4P + Xu5P), one in the paper (Appendix A) and another one encoded in the MATLAB file PSRate.m in supplement of the paper .
Wrong dimension for equations v7 and v10.
Kinetic parameters were incorrectly taken from kinetic characterization of transaldolase instead of transketolase.
In order to analyze these problems, we have rewritten the 2007 Zhu model, originally encoded in Matlab, for use with Copasi  and also rewrote it for the SimBiology. Since all errors occur in the description of the Calvin-Benson cycle, the photorespiratory pathway was not considered in the rewritten Zhu model. Finally, we improved the Zhu model in a way that all metabolites have their own concentrations, e.g., we consider separated DHAP and GAP instead of T3P pool, in comparison to original 2007 Zhu model. We then compared the steady-state values of the metabolites in dependence of ATP · (ADP + ATP)-1 ratio, reflecting the energy limitation . Since all errors are related to the transketolase, we discuss the consequence only for two substrates of transketolase, F6P (fructose 6-phosphate) and S7P (sedoheptulose-7-phosphate).
The case with two model versions gets more complicated because both of them suffer from other problems. The first one, related only to the version in Appendix A , is the consequence of a wrong dimension for equations describing the reactions catalyzed by transketolase (v7 and v10) and has its origin in the Laisk model . The problem is that the modifiers E4P (erythrose 4-phosphate) and Ri5P (ribose 5-phosphate) were multiplied instead of added in one of the equations (TEMP1, ) employed in the aforementioned reactions. Even if this particular error has an negligible impact in the calculation itself (< ± 1%, data not shown), the incorrect rate equations do not allow a computation of the simulations in biochemical simulators that test the dimension of equations (e.g., SimBiology toolbox).
An alternative but not plausible approach was employed in the case of kinetic parameters for rate equation v10 in the Zhu model. Zhu and co-workers employed different kinetic parameters  in comparison to Laisk model . The reason why this is not justified is that the new kinetic parameters were taken from transaldolase  instead of transketolase. The fact that the source of kinetic parameters is the non-photosynthetic organism Dictyostelium discoideum might be justified. However, we cannot overlook the fact that the transaldolase catalyses the same reaction as transketolase but in opposite direction; it is localized in the pentose phosphate pathway.
One question that arises is why the original kinetic parameters for the rate equation v10 from the Laisk model  were not employed in the Zhu model. Strikingly, the original parameters from the Laisk model  block the cycle and lead to an accumulation of F6P and S7P, see Figure 5 (gray). This behavior is caused primarily by big difference between the original  and new  value of the parameter kM103, 0.015 mmol·l-1 and 0.46 mmol·l-1, respectively. In order to sustain a stable solution, the value of the parameter kM103 in the rate equation v10 must be equal or higher than 0.06 mmol·l-1 (data not shown). This problem occurs, however, only for model version from Appendix A ,. The Matlab model version  can employ the original parameter kM103 = 0.015 mmol·l-1 from Laisk model  and sustains the stable solution sensitive to the changes of ATP level, see Figure 5 (violet).
Finally, and probably the most visible error in the Zhu model is an incorrect value for the equilibrium constants kE7, that is, 10  instead of 0.076  or 0.1 . The consequence is a theoretical 100-fold increase in the forward rate of the reaction F6P + GAP = E4P + Xu5P, which is able to shift the equilibrium of entire model and reduce the starch synthesis. In order to analyze the possible impact, we performed the simulations based on both Appendix A (with kM103 = 0.06 mmol·l-1) and the Matlab version of Zhu model, both with the correct value of kE7 = 0.076 . At first, we can compare original versions of Zhu model after correction on kE7. The steady state concentrations of F6P and S7P are now 6- and 4.8-fold higher, respectively, if we compare the original  and corrected Appendix A model (additional file 1). This is shown in Figure 5 in black and green. In the case of the original  and here corrected Matlab version of the Zhu model (additional file 2), the steady state concentrations of F6P is increased by 6.2-fold but S7P is 9.9-fold decreased, respectively; see Figure 5 - red and blue. The quantitative differences in the behavior between these two versions of the Zhu model suggest that one should really speak about two different models rather than two versions of one model. This conclusion is supported by final comparison of both corrected models - the difference in steady state concentration of S7P between both Zhu models is 121-fold as shown in Figure 5 - green and blue.
Besides the approach employed in the Laisk model for the phosphate translocator, one can find another and simple one, introduced in the Pettersson model  and adopted in the next model generation [3, 21]. In this approach, the phosphate translocator is encoded in the Michaelis-Menten kinetics and regulated, besides the PGA, by concentration of internal/external phosphate, DHAP and GAP. We note that, in comparison to Laisk model , the reaction PGAcytosol → Sink is redundant in this case but was employed anyway . Since this approach can be used only for chloroplast, we also tested a simplified version for cyanobacteria described with the simplest Michaelis-Menten kinetics, regulated only by concentration of PGA and estimated kM. To have an unregulated implementation of the PGA sink to compare with, we employed also irreversible kinetics using the mass action law without modifiers. Since we were interested in how the changes within the Calvin-Benson can influence associated pathways of carbon metabolism and the efficiency of the cycle itself, we have focused our analysis on the content of metabolites within the Calvin-Benson cycle as well as on the accumulation of key products, i.e., starch and PGA sink after 3 hours in the steady state conditions.
Analysis of Zhu model, versions from Matlab and Appendix A, in dependence of employed implementation of PGA sink in the model
A model, MA
A model, MM
A model, MMcb
M model, MA
M model, MM
M model, MMcb
The Calvin-Benson cycle in cyanobacteria is directly linked to other parts of the carbon metabolism. The question therefore arises if all adjacent metabolic pathways have to be considered in the kinetic model. An associated problem to solve is if it is the number of sinks itself (Figure 1) or just the total flux out of the Calvin-Benson cycle which has an impact on the efficiency of the mass production. Finally, even if photorespiration is considered only in the minority of models [3, 17], its flux out of the system is considered, usually in the phosphate translocator reactions, otherwise the Calvin-Benson cycle cannot reach the steady state due to infinite grow of PGA concentration (data not shown). We therefore aimed our focus also in this direction.
The starting point for our analysis was the reconstruction of the metabolic network for cyanobacteria Synechocystis sp. PCC 6803, adjusted for highest efficiency of the mass production . We have taken the model A (corrected Zhu model - version from Appendix A; additional file 1), particularly the modified subversion for cyanobacteria (sinks only for PGA and F6P) and used this model as a template for our further modeling efforts.
As previously mentioned, photorespiration does not have to be considered in the model but its flux out of the system cannot be neglected. Otherwise, the Calvin-Benson cycle cannot reach the steady state due to an accumulation of PGA. In order to analyze the consequences of this approach, together with analysis of the minor sinks, we have developed models C1 and C2. Models C were "extended" by PGCA (phosphoglycolate) sink (photorespiration) in a way which was employed in the majority of models of the Calvin-Benson cycle i.e., indirectly within the phosphate translocator. In our case, the PGCA sink was considered as a part of the PGA sink with the basic flux 4.2% of the RuBP (ribulose 1,5-biphosphate) synthesis flux  multiplied by a stoichiometric factor 1.66 (RuBP has 5 carbons, PGA has 3 carbons).
Both models (C1 and C2) were adapted for the same RuBP synthesis flux in the steady state (Δ flux = 5 × 10-6) and flux through F6P sinks, 2.4%  of the RuBP (ribulose 1,5-biphosphate) synthesis flux in the steady state within the Calvin-Benson cycle. The key difference between models C1 and C2 is that in the case of the model C1 (additional file 3), sinks for DHAP, E4P and Ri5P were not encoded in the model but their basic fluxes were multiplied by a stoichiometric factor q (e.g., for Ri5P, q = 1.66), summed and added to the basic flux of the PGA sink which simulates the other sinks within the PGA sink (the same approach was used for PGCA). The basic fluxes for Ri5P, DHAP and E4P sinks are 0.7%, 0.32% and 0.31% of the RuBP synthesis flux (Henning Knoop, personal communication), respectively. On the other hand, model C2 (additional file 4) was extended by sinks for DHAP, E4P and Ri5P with basic fluxes according the above mentioned proportions (Henning Knoop, personal communication).
In the case of PGA sink, the maximal and applied stable flux (value after 3 hours) was found to be equal to 18.2% in the case of C1 model (after subtraction of other sinks, see above) and 19.3% (C2 model) instead of proposed theoretical 21.7% of the RuBP synthesis flux . We note that the same flux through the PGA might be achieved without disbalancing the RuBP synthesis flux in both C models by using different set of parameters for each model but the parameters may be out of the physiological range. The initial concentrations for both C models (as well as for D models) were taken from the original Zhu model which has significantly different equilibrium in comparison to the corrected version and the system was set out of the steady state.
However, our simulations showed that it is not necessary to use two separated complex equations for the carboxylase and oxygenase (see the equation above) activity of RuBisCO as it was done by Zhu and coworkers . It is possible to encode only one equation (the one for the carboxylase activity) and the reaction for CO2 fixation can be extended by additional product PGCA: RuBP + CO2 → A × PGA + B × PGCA; in our case A = 1.916 and B = 0.084 (i.e., 4.2%, ). The differences between these two approaches are negligible both in and out of the steady state (data not shown). Moreover, our approach is easier for adjusting the level of photorespiration according our expectations and tested hypotheses and it was therefore employed in D models (see below).
Both models (D1 and D2) were again adapted for the same RuBP synthesis flux in the steady state (Δ flux = 3 × 10-4). The maximal and applied stable flux for PGA sink (value after 3 hours) was found to be equal to 17.1% in the case of D1 model (after subtraction of other sinks, see above) and 18.2% (D2 model) which shows slightly lower efficiency in comparison to C models.
The results based on D1 model (additional file 5) and D2 model (additional file 6) show totally different picture of what is going on during the response to perturbation until the steady state is reached, above all in the range of seconds - thousands of seconds, in comparison to C models, see Figure 7. If a model does not consider photorespiration described as oxygenase activity of RuBisCO, the difference in the export efficiency is reaching 10% (Figure 7, difference between gray lines). Furthermore, even if the time series data of particular metabolite shows the same qualitative pattern, see an example of E4P oscillation in Figure 7, the quantitative behavior demonstrates the differences in tenths of percentage. It is apparent that the content of metabolites within the Calvin-Benson cycle is less sensitive to changes induced by minor sinks implementation (Figure 7) if the photorespiration is implemented as oxygenase activity of RuBisCO.
Due to its long history of modelling, the most studied plant metabolic system, the Calvin-Benson cycle, provides a portfolio of kinetic models that can be used as a basis for further studies. What our present work reveals is that it is not possible to adapt existing and well-established models without significant modifications or corrections. We have shown that providing the basis for a look back at this history is important for further development in this field and essential to avoid the propagation of errors into new models generation.
In the present work we analyze two newer models, which are readily available and usable (the Poolman model and the Zhu model), as well as two fundamental models, referred to as the Hahn model and Pettersson model. Our study reveals that it is not possible to reproduce the simulations based on these fundamental models, probably due different (and unknown) computational methods employed in the original works [2, 15]. Furthermore, in the case of two newer models, we suggest corrections for these models and strongly recommend the use of two modeling tools, in our case Copasi and Matlab, as a standard approach. This approach significantly reduces the risk of systematic errors caused, e.g. by a lack of dimensional analysis.
We also discuss possible approaches for how the 3-phosphoglycerate sink should be implemented in the model. We analyze it with respect to system efficiency, regulation and content of the pool of metabolites in the inner cycle. We further show that photorespiration or at least its first step (O2 fixation) has to be implemented in the model if this model is aimed for analyses out of the steady state - we suggest a simplified modification of the model, particularly of chemical equation for CO2 fixation without changing the rate equation or kinetic parameters. Finally, the minor adjacent pathways of the carbon metabolism have been analyzed and we show that they cannot be neglected and should be considered in the model.
Our study shows that modeling photosynthesis, in combination with a kinetic model and FBA together, constrains the model and we get more accurate results and better predictions. Further combinations of methods, for instance with non-steady state metabolic control analysis in the case of detailed focus on dynamic response, can provide even more information about the properties of the model and above all about the biological system. This approach is a promising avenue for further research.
We would like to thank to Henning Knoop from Humboldt University in Berlin for providing data for the minor sinks fluxes out of the Calvin-Benson cycle on the basis of FBA for metabolic network of Synechocystic sp. PCC 6803.
This work presented here has been supported by the German Research Foundation (DFG) as part of PROMICS research group 1186.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.