BACKGROUND
the calvin-benson cycle is a central part of the carbon metabolism in oxygenic photosynthesis, involving  <dig> different enzymes that catalyze  <dig> reactions  <cit> .
the cycle is an open system, connected to light photosynthetic reactions, co <dig> fixation and other parts of carbon metabolism , e.g., starch and sucrose synthesis.
it is this complexity that motivates the use of mathematical modeling to unravel the dynamic regulation that underlies experimental observations of the calvin-benson cycle.
there are two common approaches used in modeling the calvin-benson cycle: kinetic modeling, e.g.,  <cit>  and stoichiometric modeling, e.g.,  <cit> .
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 co <dig> 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  <cit> .
even though it is a core part in a variety of photosynthetic models or projects modeling the entire photosynthesis  <cit> , kinetic modeling has mostly focused on the stability of the calvin-benson cycle itself and continues to attract considerable scientific interest  <cit> .
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 .
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  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.
methods
to investigate mathematical and computational models for their suitability and possible errors all models were encoded both in matlab, using the simbiology toolbox  and at the same time in copasi  <cit> .
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 .
RESULTS
available kinetic models for the calvin-benson cycle
there are studies that successfully use a simplified description of the calvin-benson cycle, for instance the photosynthetic oscillations in the chlorophyll fluorescence  <cit> .
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  <cit> , 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.
hahn models
this fundamental model was presented by hahn  <cit>  who also published further related models of the calvin-benson cycle  <cit> .
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  <cit>  was extended by photorespiration in the next model generation  <cit> ; together with the realistic response to the changes in the o <dig> and co <dig> concentrations  <cit> .
this assumption is based on the fact that the equations describing the calvin-benson cycle are the same in both models  <cit> ; however, some values of the kinetic parameters differ significantly.
the most prominent example is the rate of the co <dig> fixation expressed by the rate constant k <dig>  the value of this calculated parameter is in both models identical  <cit> , 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 k <dig> is  <dig>  - fold higher in the later model  <cit> .
hahn mentions that the value of k <dig> was, in order to compensate the neglected photorespiration in the earlier model  <cit> , decreased .
however, it is unlikely that the original and unknown value of k <dig> in the earlier model  <cit>  is  <dig>  - fold lower due to photorespiration because the author himself assumed in the later model  <cit>  the ratio for the carboxylase/oxygenase activity of rubisco  equal to  <dig> .
nevertheless, we have encoded the earlier hahn model  <cit>  into biochemical simulators according its original description.
having only the rate constants, we employed mass action kinetics as the author did  <cit> .
for initial concentrations we used the steady-state values calculated by hahn  <cit> .
in a pilot simulation we have encountered the very low level of atp in the steady-state; atp · - <dig> =  <dig> .
we have therefore performed a stability test in the same condition as author did , i.e., all state variables were set to their calculated steady-state values, except rubp  which was set to zero, simulation ran  <dig> minutes instead of  <dig> s in the original work.
one can draw two major conclusions from the results of the stability test, which is shown in figure  <dig>  at first, the phosphorylation of pga to bpga  is much faster than the rate of co <dig> fixation; moreover the co <dig> fixation is too slow, such that we can even observe a gradual accumulation of rubp.
secondly, the atp level is dropping down in  <dig> s of the simulation which leads to an accumulation of ru5p  and consequent suspension of the cycle.
we note that the stability test performed by the author lasted only  <dig> s.
the reason for the drop in the atp level in the certain time point  was a depletion of both intra- and extracellular phosphates .
finally, we achieved a stable solution by fixing  the pi_ext .
this approach prevents the suspension of the cycle  but the steady-state of the system is still far away from known values  or results presented by author .
pettersson model
the pettersson model  <cit>  reflects a significant effort in the elaboration of the kinetic behavior and control properties of the calvin-benson cycle.
our analysis of this model revealed two problems: an insufficient rate of the atp synthesis  and we did not reach the same results as the authors did .
for instance, the steady state concentration for pga is  <dig>  <dig> - fold lower  in comparison to the calculated value  <dig>  mm reported in the original publication  <cit> , see figure 3a1-a <dig>  finally, if we move to the physiological range for the atp concentration given by atp · - <dig> ratio  <cit> , concentration of pga further decreases .
this shift indicates that the chosen set of kinetic parameters in the pettersson model does not describe the reality in which the steady-state concentration of pga was reported to be in the range  <dig>  mm  <cit>  -  <dig>  mm  <cit> .
poolman and zhu-laisk models
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  <cit> .
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 , the original results  <cit>  are not reproducible.
therefore, if one wants to study the cycle itself or expand the model by adding other part 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  or they are made available as a supplement to a publication.
for the calvin-benson cycle, there are two main options, the poolman model  <cit>  and the zhu model  <cit> .
the poolman model is available as a .xml file  and can be therefore easily imported to any biochemical simulator supporting sbml  level  <dig>  we have employed both copasi  <cit>  and the simbiology toolbox for matlab  to study this model.
the poolman model is structurally very close to the pettersson model  <cit>  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  <cit>  and stability analysis  rather  <cit>  than for kinetic modeling of the particular metabolite.
the initial conditions were exactly the same as in the original work  <cit> .
pilot simulations based on the poolman model showed that the system reached the steady state at  <dig>  s from the start of the simulation.
if the initial concentration of pga  is increased 10-times, the steady state is reached after  <dig>  s. this feature of the model is based on the assumption that metabolites in the system are maintained at a fast equilibrium  <cit> .
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  are equal to  <dig>  <dig>  <dig> l·mmol -1·s- <dig>  <cit>  and the flux through these reactions is  <dig>  mmol·s- <dig> .
this flux is in sharp contrast to the proposed  and measured  maximal rates for these reactions.
the poolman model reaches its equilibrium at  <dig>  s of the simulation also in other cases .
on the other hand, it is known that the real system is limited by rubisco  activation whose induction phase takes  <dig> s  <cit>  -  <dig> s  <cit> .
we note that the activation of rubisco was also considered in a kinetic model with the time constant equals to  <dig> s  <cit> .
the calvin-benson cycle is an open system with many inputs and outputs .
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  concentration to variable concentration .
we note that the original setting for other fixed metabolites, i.e., for inputs  was not changed.
analysis of sinks showed an expected accumulation of the cytosolic gap , dhap  and pga, see figure  <dig>  but the system multiplied its initial metabolites pool 61-fold at  <dig> s of the simulation.
more strikingly, allowing changes in the starch concentration revealed that the starch is not only degraded but reached a high negative concentration, - <dig> -fold of the initial pool of metabolites at  <dig> s of the simulation, as shown in figure  <dig>  the negative flux through the starch degradation pathway increased the rate of sinks accumulation, i.e., the mass production, by 40% .
it was caused by error in the equation describing the starch degradation .
this problem with negative starch concentration can be either bypassed by switching off this reaction, which would set the system for the slow starch synthesis, or solved  by a replacement of pi_ch by starch in the original  <cit>  rate equation:
stpase_vm×pi_chpi_ch+stpase_km× 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  <cit> .
zhu and coworkers presented a model of c <dig> photosynthesis, extending the laisk model  <cit>  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 .
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  <dig> matlab files.
the metabolites in the zhu model reach the equilibrium at  <dig> s of the simulation , which is in agreement with previous findings  <cit>  and in contrast to  <dig>  s based on the poolman model  <cit> .
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 ke <dig> and ki <dig> in appendix c and the nowhere used constant km <dig>  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.
analysis of the zhu model and problems occurring in kinetic modeling of the calvin-benson cycle
the programming environments used for kinetic modeling of the calvin-benson cycle include turbo pascal, used in  <cit> ; scampi, used in  <cit>  or matlab, used in  <cit> .
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  <cit>  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  <cit> .
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  and another one encoded in the matlab file psrate.m in supplement of the paper  <cit> .
• wrong dimension for equations v <dig> and v <dig>
• kinetic parameters were incorrectly taken from kinetic characterization of transaldolase instead of transketolase.
• the value of the equilibrium constant ke <dig> in equation v <dig> is  <dig> times higher than what is known from literature  <cit> .
in order to analyze these problems, we have rewritten the  <dig> zhu model, originally encoded in matlab, for use with copasi  <cit>  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  <dig> zhu model.
we then compared the steady-state values of the metabolites in dependence of atp · - <dig> ratio, reflecting the energy limitation  <cit> .
since all errors are related to the transketolase, we discuss the consequence only for two substrates of transketolase, f6p  and s7p .
at first, we analyzed the key problem of the zhu model - that there are actually two versions of this model due to the different descriptions of reactions catalyzed by transketolase, i.e., equations v <dig> and v <dig>  the first version of the zhu model can be found in appendix a  <cit>   and the other one is in the supplemental matlab files  <cit> .
it is not clear which encoding was actually used for simulations in the main part of the published paper.
if the equations from appendix a are employed, the steady-state concentrations of f6p and s7p are  <dig> -fold and  <dig> -fold higher, respectively, in comparison to the model encoded in matlab, as indicated in figure  <dig> .
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  <cit> , is the consequence of a wrong dimension for equations describing the reactions catalyzed by transketolase  and has its origin in the laisk model  <cit> .
the problem is that the modifiers e4p  and ri5p  were multiplied instead of added in one of the equations  employed in the aforementioned reactions.
even if this particular error has an negligible impact in the calculation itself , the incorrect rate equations do not allow a computation of the simulations in biochemical simulators that test the dimension of equations .
an alternative but not plausible approach was employed in the case of kinetic parameters for rate equation v <dig> in the zhu model.
zhu and co-workers employed different kinetic parameters  <cit>  in comparison to laisk model  <cit> .
the reason why this is not justified is that the new kinetic parameters were taken from transaldolase  <cit>  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 v <dig> from the laisk model  <cit>  were not employed in the zhu model.
strikingly, the original parameters from the laisk model  <cit>  block the cycle and lead to an accumulation of f6p and s7p, see figure  <dig> .
this behavior is caused primarily by big difference between the original  <cit>  and new  <cit>  value of the parameter km <dig>   <dig>  mmol·l- <dig> and  <dig>  mmol·l- <dig>  respectively.
in order to sustain a stable solution, the value of the parameter km <dig> in the rate equation v <dig> must be equal or higher than  <dig>  mmol·l- <dig> .
this problem occurs, however, only for model version from appendix a  <cit> ,.
the matlab model version  <cit>  can employ the original parameter km <dig> =  <dig>  mmol·l- <dig> from laisk model  <cit>  and sustains the stable solution sensitive to the changes of atp level, see figure  <dig> .
finally, and probably the most visible error in the zhu model is an incorrect value for the equilibrium constants ke <dig>  that is,  <dig>  <cit>  instead of  <dig>   <cit>  or  <dig>   <cit> .
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  and the matlab version of zhu model, both with the correct value of ke <dig> =  <dig>   <cit> .
at first, we can compare original versions of zhu model after correction on ke <dig>  the steady state concentrations of f6p and s7p are now 6- and  <dig> -fold higher, respectively, if we compare the original  <cit>  and corrected appendix a model .
this is shown in figure  <dig> in black and green.
in the case of the original  <cit>  and here corrected matlab version of the zhu model , the steady state concentrations of f6p is increased by  <dig> -fold but s7p is  <dig> -fold decreased, respectively; see figure  <dig> - 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  <dig> - green and blue.
calvin-benson cycle as a part of the metabolic network: pga sink implementation
we have revealed several problems with the commonly cited model of the calvin-benson cycle.
having in mind that the zhu model, and others models as well  <cit> , considers also the phosphate translocator, an essential question arises of how the errors or generally any changes in the calvin-benson cycle can influence other parts of carbon metabolism.
the starting point for our analysis was naturally the laisk model  <cit>  where one can find the most complex description of the phosphate translocator.
the approach employed in laisk model for the phosphate translocator between stroma of chloroplast and cytosol, is very complex.
we therefore focused our effort only on the analysis of the pga export, pgastroma ↔ pgacytosol, described by equation vpgaout in the laisk model  <cit> :
vm11**opkmop+pgakmpga+gap+dhapkmt3p+ opckmopc+pgackmpgac+gapc+dhapckmt3pc+ subst − vm11**opkmop+pgakmpga+gap+dhapkmt3p+ opckmopc+pgackmpgac+gapc+dhapckmt3pc+ substsubst= *
the approach employed in the laisk model has its problems as well.
at first, the stromal pga is gradually accumulated without any inhibition effect in the rate of translocator, see figure  <dig>  this fact contradicts what was affirmed before  <cit> .
moreover, the export of pga, based on the equation for vpgaout  <cit> , does not depend, due its complexity and reversibility, on the value of maximum rate of the reaction vm <dig>  even if this parameter is part of the equation .
on the other hand, the original equation vpgaout is sensitive to the changes in metabolite concentrations.
for instance, if the concentration of cytosolic phosphate is increased 100-times, in comparison to stromal phosphate, the stromal pga is decreased  <dig> -times but still without any sign of saturation as shown in figure  <dig>  the problem with pga accumulation in the stroma can be partially solved by adding another reaction: pgastroma → pgacytosol → sink.
this solution stabilizes the stromal pga and defines the expected saturation, see figure  <dig>
besides the approach employed in the laisk model for the phosphate translocator, one can find another and simple one, introduced in the pettersson model  <cit>  and adopted in the next model generation  <cit> .
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  <cit> , the reaction pgacytosol → sink is redundant in this case but was employed anyway  <cit> .
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  <dig> hours in the steady state conditions.
the results of our analysis summarized in the table  <dig> show several interesting outcomes.
at first, it is clear that the content of metabolites in the calvin-benson cycle depends significantly both on the changes  within the cycle itself, as discussed in the section above, but also on the level of regulation of pga sink .
a similar pattern can be observed also in the case of starch synthesis.
more strikingly, the steady state production of pga depends only slightly on how/if it is regulated or if any changes occurred in equilibrium within the calvin-benson cycle.
finally, if we compare the efficiency of the implemented approaches, which can be measured by total mass production, unregulated approach represented by mass action kinetics is less efficient in comparison to two regulated approaches.
this suggests that "export" of metabolites from the calvin-benson cycle is, as expected, regulated also in cyanobacteria, which do not have the phosphate translocator.
however, there are no significant differences between more complex regulation, employed in poolman and zhu models, or here introduced simplified approach for pga sink, see table  <dig>  taken together, the employed approach for pga sink has cardinal influence on the calvin-benson cycle and starch synthesis but not for the export of pga itself.
the results show several approaches in the  <dig> - phosphoglycerate sink implementation and its consequence on the system efficiency , regulation and content of the pool of metabolites as well as the changes within the calvin-benson after  <dig> hours of simulation .
the numbers indicates percentage changes in comparison to corrected zhu model , matlab version.
the response of system to the changes has been tested by using the original zhu model  which has significantly different equilibrium in comparison to the corrected model versions.
meaning of particular abbreviations: c-b pool - metabolic pool of the calvin-benson cycle, a - appendix a model version, m - matlab model version, ma - mass action law kinetics, mm -michaelis-menten kinetics , cb - designed for cyanobacteria.
the bold font indicates the highest value.
minor sinks for adjacent pathways, photorespiration and system efficiency
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  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  <cit> , 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 .
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  <dig>  adjusted for highest efficiency of the mass production  <cit> .
we have taken the model a , particularly the modified subversion for cyanobacteria  and used this model as a template for our further modeling efforts.
photorespiration not directly implemented but considered within pga sink
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 c <dig> and c <dig>  models c were "extended" by pgca  sink  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  <dig> % of the rubp  synthesis flux  <cit>  multiplied by a stoichiometric factor  <dig>  .
both models  were adapted for the same rubp synthesis flux in the steady state  and flux through f6p sinks,  <dig> %  <cit>  of the rubp  synthesis flux in the steady state within the calvin-benson cycle.
the key difference between models c <dig> and c <dig> is that in the case of the model c <dig> , sinks for dhap, e4p and ri5p were not encoded in the model but their basic fluxes were multiplied by a stoichiometric factor q , summed and added to the basic flux of the pga sink which simulates the other sinks within the pga sink .
the basic fluxes for ri5p, dhap and e4p sinks are  <dig> %,  <dig> % and  <dig> % of the rubp synthesis flux , respectively.
on the other hand, model c <dig>  was extended by sinks for dhap, e4p and ri5p with basic fluxes according the above mentioned proportions .
in the case of pga sink, the maximal and applied stable flux  was found to be equal to  <dig> % in the case of c <dig> model  and  <dig> %  instead of proposed theoretical  <dig> % of the rubp synthesis flux  <cit> .
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  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.
the comparison of results based on c models is striking.
having in mind that the only change between the models c <dig> and c <dig> was redistribution of the relatively small flux  and slight difference in the pga sink , the ratios of both the total mass production and sum of sinks show dynamic changes for several hours of the simulation, see figure  <dig>  one can even observe a difference around 45% or even an oscillation in the ratio curve for particular metabolite  which, however, does occur for separate time-series data of e4p .
finally, the presented difference in the mass production and e4p was a response to dynamic changes within the models, changes which are in certain extent inevitable even in the controlled environment.
photorespiration implemented directly as o <dig> fixation
models d, in contrast to models c, implement the photorespiration directly as o <dig> fixation.
at first, we have employed for the description of o <dig> fixation rate a slightly modified equation applied in the zhu model  <cit> , based on previous experimental and theoretical findings  <cit>  :
rubp×vo×o2o2+ko×rubp+kr×
however, our simulations showed that it is not necessary to use two separated complex equations for the carboxylase and oxygenase  activity of rubisco as it was done by zhu and coworkers  <cit> .
it is possible to encode only one equation  and the reaction for co <dig> fixation can be extended by additional product pgca: rubp + co <dig> → a × pga + b × pgca; in our case a =  <dig>  and b =  <dig>  .
the differences between these two approaches are negligible both in and out of the steady state .
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 .
both models  were again adapted for the same rubp synthesis flux in the steady state .
the maximal and applied stable flux for pga sink  was found to be equal to  <dig> % in the case of d <dig> model  and  <dig> %  which shows slightly lower efficiency in comparison to c models.
the results based on d <dig> model  and d <dig> model  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  <dig>  if a model does not consider photorespiration described as oxygenase activity of rubisco, the difference in the export efficiency is reaching 10% .
furthermore, even if the time series data of particular metabolite shows the same qualitative pattern, see an example of e4p oscillation in figure  <dig>  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  if the photorespiration is implemented as oxygenase activity of rubisco.
taken together, the comparisons indicate that it really matters how the photorespiration is incorporated into the model.
therefore, the majority of models, above all those models employing michaelis-menten kinetics  which was also employed in our models, are not suitable for simulating the dynamic responses  focused on analyzing the changes in metabolite content of the calvin-benson cycle  or on mass production .
on the other hand, the zhu model can be employed for such analysis or modeling the photorespiration if the aforementioned corrections of the calvin-benson cycle  are applied.
finally, one can also conclude that the minor sinks, even with small relative fluxes, should be considered in the models of the calvin-benson cycle for analyses out of the steady state.
the scheme of the model d <dig>  which we recommend for next modeling efforts, is shown in figure  <dig>
CONCLUSIONS
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 , 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  computational methods employed in the original works  <cit> .
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  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 co <dig> 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.
competing interests
the authors declare that they have no competing interests.
authors' contributions
ideas and concepts were jointly discussed among all authors.
the manuscript was written by jj and ow contributed to the writing of the final manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
corrected model of the calvin-benson cycle and starch synthesis based on the zhu model_appendix a version.
model a, sbml l2v4; model does not consider the photorespiratory pathway
click here for file
additional file 2
corrected model of the calvin-benson cycle and starch synthesis based on the zhu model_matlab version.
model b, sbml l2v4; model does not consider the photorespiratory pathway
click here for file
additional file 3
model_based on the corrected zhu model_ appendix a version.
model c <dig>  sbml l2v4; modified model a for the purpose of comparison with model c <dig> -analysis of minor sinks and mass production efficiency; model does not consider the photorespiratory pathway but pgca sink is considered within the pga sink
click here for file
additional file 4
minor sinks model_based on the corrected zhu model_ appendix a version.
model c <dig>  sbml l2v4; model does not consider the photorespiratory pathway but pgca sink is considered within pga sink, minor sinks for dhap, e4p and ri5p are implemented
click here for file
additional file 5
model_based on the corrected zhu model_ appendix a version.
model d <dig>  sbml l2v4; modified model a for the purpose of comparison with model d <dig> -analysis of sinks and mass production efficiency; model does consider the first step of the photorespiratory pathway
click here for file
additional file 6
minor sinks model_based on the corrected zhu model_ appendix a version.
model d <dig>  sbml l2v4; model does consider the first step of the photorespiratory pathway , the minor sinks for dhap, e4p and ri5p are implemented
click here for file
acknowledgements
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  <dig>
this work presented here has been supported by the german research foundation  as part of promics research group  <dig>

 the calvin-benson cycle is a central part of the carbon metabolism in oxygenic photosynthesis, involving different enzymes that catalyze reactions. the cycle is an open system, connected to light photosynthetic reactions, co fixation and other parts of carbon metabolism, e.g., starch and sucrose synthesis. it is this complexity that motivates the use of mathematical modeling to unravel the dynamic regulation that underlies experimental observations of the calvin-benson cycle. there are two common approaches used in modeling the calvin-benson cycle: kinetic modeling, e.g., and stoichiometric modeling, e.g.,. 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 co 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. even though it is a core part in a variety of photosynthetic models or projects modeling the entire photosynthesis, kinetic modeling has mostly focused on the stability of the calvin-benson cycle itself and continues to attract considerable scientific interest. 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. 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 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. methods to investigate mathematical and computational models for their suitability and possible errors all models were encoded both in matlab, using the simbiology toolbox 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.  available kinetic models for the calvin-benson cycle 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, 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. hahn models this fundamental model was presented by hahn who also published further related models of the calvin-benson cycle. 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 o and co concentrations. this assumption is based on the fact that the equations describing the calvin-benson cycle are the same in both models; however, some values of the kinetic parameters differ significantly. the most prominent example is the rate of the co fixation expressed by the rate constant k the value of this calculated parameter is in both models identical, 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 k is - fold higher in the later model. hahn mentions that the value of k was, in order to compensate the neglected photorespiration in the earlier model, decreased. however, it is unlikely that the original and unknown value of k in the earlier model is - fold lower due to photorespiration because the author himself assumed in the later model the ratio for the carboxylase/oxygenase activity of rubisco equal to. 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. in a pilot simulation we have encountered the very low level of atp in the steady-state; atp · - =. we have therefore performed a stability test in the same condition as author did, i.e., all state variables were set to their calculated steady-state values, except rubp which was set to zero, simulation ran minutes instead of s in the original work. one can draw two major conclusions from the results of the stability test, which is shown in figure at first, the phosphorylation of pga to bpga is much faster than the rate of co fixation; moreover the co fixation is too slow, such that we can even observe a gradual accumulation of rubp. secondly, the atp level is dropping down in s of the simulation which leads to an accumulation of ru5p and consequent suspension of the cycle. we note that the stability test performed by the author lasted only s. the reason for the drop in the atp level in the certain time point was a depletion of both intra- and extracellular phosphates. finally, we achieved a stable solution by fixing the pi_ext. this approach prevents the suspension of the cycle but the steady-state of the system is still far away from known values or results presented by author. pettersson model the pettersson model reflects a significant effort in the elaboration of the kinetic behavior and control properties of the calvin-benson cycle. our analysis of this model revealed two problems: an insufficient rate of the atp synthesis and we did not reach the same results as the authors did. for instance, the steady state concentration for pga is - fold lower in comparison to the calculated value mm reported in the original publication, see figure 3a1-a finally, if we move to the physiological range for the atp concentration given by atp · - ratio, concentration of pga further decreases. this shift indicates that the chosen set of kinetic parameters in the pettersson model does not describe the reality in which the steady-state concentration of pga was reported to be in the range mm - mm. poolman and zhu-laisk models 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. 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, the original results are not reproducible. therefore, if one wants to study the cycle itself or expand the model by adding other part 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 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 and can be therefore easily imported to any biochemical simulator supporting sbml level we have employed both copasi and the simbiology toolbox for matlab 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 rather than for kinetic modeling of the particular metabolite. 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 s from the start of the simulation. if the initial concentration of pga is increased 10-times, the steady state is reached after 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 are equal to l·mmol -1·s- and the flux through these reactions is mmol·s-. this flux is in sharp contrast to the proposed and measured maximal rates for these reactions. the poolman model reaches its equilibrium at s of the simulation also in other cases. on the other hand, it is known that the real system is limited by rubisco activation whose induction phase takes s - s. we note that the activation of rubisco was also considered in a kinetic model with the time constant equals to s. the calvin-benson cycle is an open system with many inputs and outputs. 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 concentration to variable concentration. we note that the original setting for other fixed metabolites, i.e., for inputs was not changed. analysis of sinks showed an expected accumulation of the cytosolic gap, dhap and pga, see figure but the system multiplied its initial metabolites pool 61-fold at s of the simulation. more strikingly, allowing changes in the starch concentration revealed that the starch is not only degraded but reached a high negative concentration, - -fold of the initial pool of metabolites at s of the simulation, as shown in figure the negative flux through the starch degradation pathway increased the rate of sinks accumulation, i.e., the mass production, by 40%. it was caused by error in the equation describing the starch degradation. this problem with negative starch concentration can be either bypassed by switching off this reaction, which would set the system for the slow starch synthesis, or solved by a replacement of pi_ch by starch in the original rate equation: stpase_vm×pi_chpi_ch+stpase_km× 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 c 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. 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 matlab files. the metabolites in the zhu model reach the equilibrium at s of the simulation, which is in agreement with previous findings and in contrast to 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 ke and ki in appendix c and the nowhere used constant km 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. analysis of the zhu model and problems occurring in kinetic modeling of the calvin-benson cycle 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 and another one encoded in the matlab file psrate.m in supplement of the paper. • wrong dimension for equations v and v • kinetic parameters were incorrectly taken from kinetic characterization of transaldolase instead of transketolase. • the value of the equilibrium constant ke in equation v is times higher than what is known from literature. in order to analyze these problems, we have rewritten the 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 zhu model. we then compared the steady-state values of the metabolites in dependence of atp · - ratio, reflecting the energy limitation. since all errors are related to the transketolase, we discuss the consequence only for two substrates of transketolase, f6p and s7p. at first, we analyzed the key problem of the zhu model - that there are actually two versions of this model due to the different descriptions of reactions catalyzed by transketolase, i.e., equations v and v the first version of the zhu model can be found in appendix a and the other one is in the supplemental matlab files. it is not clear which encoding was actually used for simulations in the main part of the published paper. if the equations from appendix a are employed, the steady-state concentrations of f6p and s7p are -fold and -fold higher, respectively, in comparison to the model encoded in matlab, as indicated in figure. 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 and has its origin in the laisk model. the problem is that the modifiers e4p and ri5p were multiplied instead of added in one of the equations employed in the aforementioned reactions. even if this particular error has an negligible impact in the calculation itself, the incorrect rate equations do not allow a computation of the simulations in biochemical simulators that test the dimension of equations. an alternative but not plausible approach was employed in the case of kinetic parameters for rate equation v 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 v 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. this behavior is caused primarily by big difference between the original and new value of the parameter km mmol·l- and mmol·l- respectively. in order to sustain a stable solution, the value of the parameter km in the rate equation v must be equal or higher than mmol·l-. this problem occurs, however, only for model version from appendix a,. the matlab model version can employ the original parameter km = mmol·l- from laisk model and sustains the stable solution sensitive to the changes of atp level, see figure. finally, and probably the most visible error in the zhu model is an incorrect value for the equilibrium constants ke that is, instead of or. 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 and the matlab version of zhu model, both with the correct value of ke =. at first, we can compare original versions of zhu model after correction on ke the steady state concentrations of f6p and s7p are now 6- and -fold higher, respectively, if we compare the original and corrected appendix a model. this is shown in figure in black and green. in the case of the original and here corrected matlab version of the zhu model, the steady state concentrations of f6p is increased by -fold but s7p is -fold decreased, respectively; see figure - 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 - green and blue. calvin-benson cycle as a part of the metabolic network: pga sink implementation we have revealed several problems with the commonly cited model of the calvin-benson cycle. having in mind that the zhu model, and others models as well, considers also the phosphate translocator, an essential question arises of how the errors or generally any changes in the calvin-benson cycle can influence other parts of carbon metabolism. the starting point for our analysis was naturally the laisk model where one can find the most complex description of the phosphate translocator. the approach employed in laisk model for the phosphate translocator between stroma of chloroplast and cytosol, is very complex. we therefore focused our effort only on the analysis of the pga export, pgastroma ↔ pgacytosol, described by equation vpgaout in the laisk model: vm11**opkmop+pgakmpga+gap+dhapkmt3p+ opckmopc+pgackmpgac+gapc+dhapckmt3pc+ subst − vm11**opkmop+pgakmpga+gap+dhapkmt3p+ opckmopc+pgackmpgac+gapc+dhapckmt3pc+ substsubst= * the approach employed in the laisk model has its problems as well. at first, the stromal pga is gradually accumulated without any inhibition effect in the rate of translocator, see figure this fact contradicts what was affirmed before. moreover, the export of pga, based on the equation for vpgaout, does not depend, due its complexity and reversibility, on the value of maximum rate of the reaction vm even if this parameter is part of the equation. on the other hand, the original equation vpgaout is sensitive to the changes in metabolite concentrations. for instance, if the concentration of cytosolic phosphate is increased 100-times, in comparison to stromal phosphate, the stromal pga is decreased -times but still without any sign of saturation as shown in figure the problem with pga accumulation in the stroma can be partially solved by adding another reaction: pgastroma → pgacytosol → sink. this solution stabilizes the stromal pga and defines the expected saturation, see figure 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. 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 hours in the steady state conditions. the results of our analysis summarized in the table show several interesting outcomes. at first, it is clear that the content of metabolites in the calvin-benson cycle depends significantly both on the changes within the cycle itself, as discussed in the section above, but also on the level of regulation of pga sink. a similar pattern can be observed also in the case of starch synthesis. more strikingly, the steady state production of pga depends only slightly on how/if it is regulated or if any changes occurred in equilibrium within the calvin-benson cycle. finally, if we compare the efficiency of the implemented approaches, which can be measured by total mass production, unregulated approach represented by mass action kinetics is less efficient in comparison to two regulated approaches. this suggests that "export" of metabolites from the calvin-benson cycle is, as expected, regulated also in cyanobacteria, which do not have the phosphate translocator. however, there are no significant differences between more complex regulation, employed in poolman and zhu models, or here introduced simplified approach for pga sink, see table taken together, the employed approach for pga sink has cardinal influence on the calvin-benson cycle and starch synthesis but not for the export of pga itself. the results show several approaches in the - phosphoglycerate sink implementation and its consequence on the system efficiency, regulation and content of the pool of metabolites as well as the changes within the calvin-benson after hours of simulation. the numbers indicates percentage changes in comparison to corrected zhu model, matlab version. the response of system to the changes has been tested by using the original zhu model which has significantly different equilibrium in comparison to the corrected model versions. meaning of particular abbreviations: c-b pool - metabolic pool of the calvin-benson cycle, a - appendix a model version, m - matlab model version, ma - mass action law kinetics, mm -michaelis-menten kinetics, cb - designed for cyanobacteria. the bold font indicates the highest value. minor sinks for adjacent pathways, photorespiration and system efficiency 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 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, 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. 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 adjusted for highest efficiency of the mass production. we have taken the model a, particularly the modified subversion for cyanobacteria and used this model as a template for our further modeling efforts. photorespiration not directly implemented but considered within pga sink 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 c and c models c were "extended" by pgca sink 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 % of the rubp synthesis flux multiplied by a stoichiometric factor. both models were adapted for the same rubp synthesis flux in the steady state and flux through f6p sinks, % of the rubp synthesis flux in the steady state within the calvin-benson cycle. the key difference between models c and c is that in the case of the model c, sinks for dhap, e4p and ri5p were not encoded in the model but their basic fluxes were multiplied by a stoichiometric factor q, summed and added to the basic flux of the pga sink which simulates the other sinks within the pga sink. the basic fluxes for ri5p, dhap and e4p sinks are %, % and % of the rubp synthesis flux, respectively. on the other hand, model c was extended by sinks for dhap, e4p and ri5p with basic fluxes according the above mentioned proportions. in the case of pga sink, the maximal and applied stable flux was found to be equal to % in the case of c model and % instead of proposed theoretical % 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 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. the comparison of results based on c models is striking. having in mind that the only change between the models c and c was redistribution of the relatively small flux and slight difference in the pga sink, the ratios of both the total mass production and sum of sinks show dynamic changes for several hours of the simulation, see figure one can even observe a difference around 45% or even an oscillation in the ratio curve for particular metabolite which, however, does occur for separate time-series data of e4p. finally, the presented difference in the mass production and e4p was a response to dynamic changes within the models, changes which are in certain extent inevitable even in the controlled environment. photorespiration implemented directly as o fixation models d, in contrast to models c, implement the photorespiration directly as o fixation. at first, we have employed for the description of o fixation rate a slightly modified equation applied in the zhu model, based on previous experimental and theoretical findings: rubp×vo×o2o2+ko×rubp+kr× however, our simulations showed that it is not necessary to use two separated complex equations for the carboxylase and oxygenase activity of rubisco as it was done by zhu and coworkers. it is possible to encode only one equation and the reaction for co fixation can be extended by additional product pgca: rubp + co → a × pga + b × pgca; in our case a = and b =. the differences between these two approaches are negligible both in and out of the steady state. 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. both models were again adapted for the same rubp synthesis flux in the steady state. the maximal and applied stable flux for pga sink was found to be equal to % in the case of d model and % which shows slightly lower efficiency in comparison to c models. the results based on d model and d model 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 if a model does not consider photorespiration described as oxygenase activity of rubisco, the difference in the export efficiency is reaching 10%. furthermore, even if the time series data of particular metabolite shows the same qualitative pattern, see an example of e4p oscillation in figure 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 if the photorespiration is implemented as oxygenase activity of rubisco. taken together, the comparisons indicate that it really matters how the photorespiration is incorporated into the model. therefore, the majority of models, above all those models employing michaelis-menten kinetics which was also employed in our models, are not suitable for simulating the dynamic responses focused on analyzing the changes in metabolite content of the calvin-benson cycle or on mass production. on the other hand, the zhu model can be employed for such analysis or modeling the photorespiration if the aforementioned corrections of the calvin-benson cycle are applied. finally, one can also conclude that the minor sinks, even with small relative fluxes, should be considered in the models of the calvin-benson cycle for analyses out of the steady state. the scheme of the model d which we recommend for next modeling efforts, is shown in figure  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, 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 computational methods employed in the original works. 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 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 co 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. competing interests the authors declare that they have no competing interests. authors' contributions ideas and concepts were jointly discussed among all authors. the manuscript was written by jj and ow contributed to the writing of the final manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 corrected model of the calvin-benson cycle and starch synthesis based on the zhu model_appendix a version. model a, sbml l2v4; model does not consider the photorespiratory pathway click here for file additional file 2 corrected model of the calvin-benson cycle and starch synthesis based on the zhu model_matlab version. model b, sbml l2v4; model does not consider the photorespiratory pathway click here for file additional file 3 model_based on the corrected zhu model_ appendix a version. model c sbml l2v4; modified model a for the purpose of comparison with model c -analysis of minor sinks and mass production efficiency; model does not consider the photorespiratory pathway but pgca sink is considered within the pga sink click here for file additional file 4 minor sinks model_based on the corrected zhu model_ appendix a version. model c sbml l2v4; model does not consider the photorespiratory pathway but pgca sink is considered within pga sink, minor sinks for dhap, e4p and ri5p are implemented click here for file additional file 5 model_based on the corrected zhu model_ appendix a version. model d sbml l2v4; modified model a for the purpose of comparison with model d -analysis of sinks and mass production efficiency; model does consider the first step of the photorespiratory pathway click here for file additional file 6 minor sinks model_based on the corrected zhu model_ appendix a version. model d sbml l2v4; model does consider the first step of the photorespiratory pathway, the minor sinks for dhap, e4p and ri5p are implemented click here for file acknowledgements 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 this work presented here has been supported by the german research foundation as part of promics research group 
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. 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. 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. 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. modeling the calvin-benson cycle has a history in the field of theoretical biology. 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. we further show that photorespiration or at least its first step has to be implemented in the model if this model is aimed for analyses out of the steady state.
 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 has to be implemented in the model if this model is aimed for analyses out of the steady state. 