aa r X i v : . [ q - b i o . M N ] S e p A Calvin bestiary
Alan D. RendallInstitut f¨ur MathematikJohannes Gutenberg-Universit¨atStaudingerweg 9D-55099 MainzGermany
Abstract
This paper compares a number of mathematical models for the Calvincycle of photosynthesis and presents theorems on the existence and sta-bility of steady states of these models. Results on five-variable modelsin the literature are surveyed. Next a number of larger models relatedto one introduced by Pettersson and Ryde-Pettersson are discussed. Themathematical nature of this model is clarified, showing that it is naturallydefined as a system of differential-algebraic equations. It is proved thatthere are choices of parameters for which this model admits more thanone positive steady state. This is done by analysing the limit where thestorage of sugars from the cycle as starch is shut down. There is also adiscussion of the minimal models for the cycle due to Hahn.
The Calvin cycle is a part of photosynthesis and there are many mathematicalmodels for this biochemical system in the literature. Reviews of these can befound in [2], [10] and [3]. The aim of this paper is to survey what is known aboutthe dynamics of these models with a focus on what has been proved rigorously.It should be pointed out right away that the rigorous results constitute a smallisland in an ocean of simulations and heuristics. To start with it is necessaryto fix the boundary of the area to be covered. The models treated are alldeterministic, continuous time evolution equations without delays and spatialvariations are neglected. Thus mathematically we are dealing with systems ofordinary differential equations (ODE) or differential-algebraic equations (DAE).The unknowns are concentrations of chemical substances depending on time.Photosynthesis is a process of central importance in biology and, as a con-sequence, in our daily lives. It consists of two major parts. In the first of these(the light reactions) energy is captured from sunlight and molecular oxygen isproduced. In the second (the dark reactions) carbon dioxide from the air isused to make carbohydrates. For reasons to be described later the second part1s also called the Calvin cycle. The models which are the subject of what fol-lows relate to the Calvin cycle and all describe ordinary chemical reactions insolution together with simple sources, sinks and transport processes betweencellular compartments which fit into the same mathematical framework. Thelight reactions involve electrochemistry on a membrane, a type of process whosemodelling will not be considered here. A comprehensive introduction to thebiochemistry of photosynthesis can be found in [9].If we are describing one biological system here, why should there be manymathematical models for it? This is a consequence of some general features ofthe modelling of biochemical systems which will now be listed. The first is thata biochemical system like the Calvin cycle is in reality coupled to many otherchemical processes (the light reactions, sucrose production etc.) and so we haveto make a choice of the set of chemical species whose concentrations are includedas unknowns in the ODE system. The hope is that these concentrations haveonly a small effect on the concentrations of the other species with which thechosen ones interact. The concentrations of these other species are taken tobe constant and we refer to them as external species while the species whoseconcentrations are the unknowns in the ODE system are referred to as internalspecies. A possible justification for this procedure is that if the concentrationof an external species is very high it will remain approximately constant even ifsome amount of the substance concerned is being produced or consumed by someof the other reactions. There is also a choice of which reactions are consideredto be taking place at an appreciable rate. Usually the stoichiometry of thereactions is known but the same cannot be said of the reaction rates. Therefurther assumptions have to be made. Summing up, different mathematicalmodels arise through different choices of the species and reactions included andthe reaction kinetics. Furthermore it may happen that models are replacedby smaller ones using limits involving time scale separation or elimination ofintermediate species in some reactions.After these preliminary considerations we may look at what a standard text-book on cell biology [1] tells us about the Calvin cycle. The essential featuresof this process were worked out by Melvin Calvin and his collaborators (earningCalvin the Nobel prize for chemistry in 1961). Often the situation of carbondioxide and light saturation is considered. Calvin’s experiments were done underthese circumstances and they are often assumed to hold when doing modelling.This means on the one hand that carbon dioxide is so plentiful that it can beconsidered as an external species. On the other hand the substances ATP andNADPH which are supplied by the light reactions are assumed to be plentiful.Thus CO , ATP and NADPH are taken as external species. The same is trueof ADP and NADP which are produced from ATP and NADPH in certain re-actions. Inorganic phosphate P i is often also treated as an external species. Allthe reactions are catalysed by enzymes but these are usually treated as externalspecies. The substances which remain in the description in [1], and which willbe internal species, are ribulose 1,5-bisphosphate (RuBP), 3-phosphoglycerate(PGA), 1,3-bisphosphoglycerate (DPGA), glyceraldehyde 3-phosphate (GAP)and ribulose 5-phosphate (Ru5P). The simplest assumption is that each of these2ubstances reacts to give the next with a final reaction taking us back from Ru5Pto RuBP. Thus we have a cycle, explaining the other part of the name ’Calvincycle’. While most of these are bona fide reactions, that leading from GAP toRu5P is an effective reaction (or, less respectfully, a fudge reaction) resultingfrom collapsing part of a more complicated network. This can easily be rec-ognized by means of the exotic stoichiometry, with five molecules going in andthree coming out. In addition there are two transport processes in which PGAand GAP are exported to the cytosol from the chloroplast where the Calvincycle takes place. Thus some simple models of photosynthesis (to be consideredin more detail in the next section) have five species and seven reactions.In this paper there is no attempt to present a systematic catalogue of models.Instead it is like an accompanied walk through a zoo, where the visitor is takento see the lions and the elephants but also less familiar exhibits such as the giantanteater or the Tasmanian devil. It starts with the simplest and best knownmodels and is led by the consideration of various issues to ones which have beenstudied less.So many related models are considered in what follows that it would becumbersome to have a name for each of them. Some names will be used butin addition the models will be given numbers according to the pattern Modelm.n.k, where roughly speaking this means the variant k of the model n firstintroduced in section m .The paper is organized as follows. Section 2 introduces the models with fivespecies and is mainly a survey of known results concerning them. These provideinformation about the existence and stability of positive steady states and theexistence of solutions where the concentrations tend to zero or to infinity at latetimes. The process of passing from one model to another by making an internalspecies into an external one is carried out in a simple example. In Section 3models with a larger number of variables (about fifteen) are considered which arevariants of one introduced in [14]. Some known results on the (non)-existence ofpositive steady states and the ways in which concentrations can approach zeroat late times are reviewed and extended. It is shown how the model of [14] itselfcan be given a clear mathematical formulation as a system of DAE. It is alsoshown how ATP can be made into an external species in these models. Section4 contains a proof that the model of [14] admits more than one positive steadystate for suitable values of the parameters. This is related to a stoichiometricgenerator for the network and other generators, which might also be helpful inthe search for steady states, are presented. A similar approach can be appliedto the related reduced Poolman model and this is done in Section 5, where itis shown that there are parameter values for which there exist at least threepositive steady states. The last section is concerned with some prospects forfuture progress and briefly discusses some simplified models due to Hahn.3 The five-species models
In this section some models will be considered where the unknowns are theconcentrations of the five substrates introduced in the previous section. Theconcentration of a substance X is denoted by x X . The reactions are thosementioned in the introduction and they are all assumed to be irreversible. Thesimplest type of kinetics which can be assumed is mass action kinetics. For thiswe need to know the exact stoichiometry. At this point there is an ambiguityresulting from the fudge reaction. In that case the biochemistry does not de-termine unique values for the stoichiometric coefficients. It only determines theratio of the number of molecules going into the reaction and the number comingout. One possible choice, which was made in [19], is to use the coefficients 1and 0 .
6. Since this is only an effective reaction there is no strong argument thatthese coefficients should be integers. Nevertheless in [6] the authors preferredto use 5 and 3. This does make a difference to the evolution equations result-ing from the assumption of mass action kinetics. With the stoichiometry of [6]mass action kinetics leads to nonlinear evolution equations (Model 2.1.1) whilethe kinetics of [19] gives linear equations (Model 2.1.2). An alternative to massaction kinetics is Michaelis-Menten kinetics, either with the coefficients of [6](Model 2.2.1) or those of [19] (Model 2.2.2).We now examine the dynamics of these models. It was shown in [6] that ifthe reaction constants satisfy a certain inequality ( k ≤ k ) Model 2.1.1 doesnot possess any positive steady states while if k > k it possesses preciselyone positive steady state for given values of the parameters. These statementsare obtained by explicit calculation. In the latter case it was shown that thesteady state is unstable. This could be seen as a disappointment since it might besupposed for biological reasons that the cycle can exist in a stable configuration.This is not absolutely clear since a good model need not give a good descriptionof the dynamics globally in time but only on a time scale long enough so asto capture the processes which are to be described. This quantitative line ofthought will not be pursued further here. In [17] it was further shown that inthe case k ≤ k all concentrations tend to zero as time tends to infinity. Thiswas done with the help of a Lyapunov function. In [6] and [17] the inequalityrelating k and k was not given any biological interpretation and no explanationwas given for the Lyapunov function which was found by trial and error. Moreinsight on these questions was obtained while studying a more complicatedmodel of the Calvin cycle in [12]. The reaction constant k controls the rate ofexport of PGA from the chloroplast and in reality this is coupled to the importof inorganic phosphate. Thus intuitively k secretly contains a dependence onthe constant concentration of inorganic phosphate in the cytosol. The positivesteady state disappears when there is too much phosphate in the cytosol. Thenthe production of sugars by the cycle cannot keep up with the export. Thisphenomenon has been called overload breakdown [14]. To obtain some insightinto the Lyapunov function it is helpful to consider the total number of carbonatoms in the system. The reactions within the cycle conserve carbon apart fromthe fact that at one point carbon dioxide is imported. The export processes also4o not conserve carbon. Nevertheless the time derivative of the total amountof carbon only has a few contributions. One of them has a positive sign butmodifying the coefficient of x PGA allows a Lyapunov function L to be obtained.There remains the question of what happens to solutions of Model 2.1.1 at latetimes when k > k . It was shown in [17] that there are solutions whereall concentrations tend to zero at late times (using a modification L of L )and solutions for which the concentrations tend to infinity as t → ∞ (runawaysolutions) and their leading order asymptotics were determined. From whathas just been said it can be seen that a number of facts are known about thedynamics of Model 2.1.1 but there remain open questions, for instance whetherperiodic solutions exist. Information about Model 2.1.2 has also been obtainedin [17]. In particular, there are either no positive steady states or a wholecontinuum of steady states, depending on the values of the reaction constants.Both solutions converging to the origin at late times and runaway solutionsoccur.When mass action kinetics is replaced by Michaelis-Menten kinetics there arestill runaway solutions but there are also more interesting steady states. Con-cerning Model 2.2.2 it stated in [19] that there is at most one steady state whichis ‘physiologically feasible’. This last condition includes restrictions on the valuesof the model parameters. In some cases a parameter interval is chosen centred ata value taken from the experimental literature. Some Michaelis constants K mi are set to fixed values and it is not clear to this author where these values comefrom. The paper [19] uses computer-assisted methods which are claimed to provethe assertion about the limitation on steady states. In [5] a purely analyticalproof of the assertion was given under the assumptions on the Michaelis con-stants made in [19]. It was also proved that the assertion depends essentially onthese assumptions. There are examples with κ = ( K m − K m )( K m − K m ) < κ vanish there is a continuum of steady states. It was alsoshown that there exist cases with two isolated steady states where one of themis stable and the other is unstable. This is proved by showing that there isa bifurcation with a one-dimensional centre manifold. In particular we obtainmodels admitting a stable positive steady state although it is unclear whetherthe parameter values required for this are biologically reasonable. Model 2.2.1also permits the existence of two positive steady states, one of which is stableand the other is unstable.Another type of model, considered in [6], is obtained if each of the basicreactions considered up to now is replaced by a Michaelis-Menten scheme witha substrate, an enzyme and a substrate-enzyme complex (Models 2.3.1 and 2.3.2)and these will be discussed in this section although they contain many more thanfive variables. It is possible to pass from these models to Models 2.2.1 and 2.2.2by a Michaelis-Menten reduction which is well-behaved in the sense of geometricsingular perturbation theory (GSPT) - the transverse eigenvalues have negativereal parts. (For background information on GSPT we refer to [11].) This meansthat we can transfer information on the existence and stability of steady statesfrom Models 2.2.1 and 2.2.2 to Models 2.3.1 and 2.3.2 in a straightforward way.5n fact the existence of more than one steady state of Model 2.3.1 was discovereddirectly in [6] with the help of elementary flux modes.In another model introduced in [6] ATP was made into an internal species,thus producing a six-variable model and diffusion of ATP was included. Re-stricting consideration to spatially homogeneous solutions reduces the resultingsystem of reaction-diffusion equations to a system of ODE (Model 2.4.1). Itturns out that Model 2.4.1 can be analysed as in the cases of Models 2.2.1 and2.2.2, giving the existence of two steady states, one stable and one unstable[5]. Interestingly, the solutions of Model 2.4.1 are bounded although this is non-trivial to prove [17]. In all these models ω -limit points where some concentrationvanishes are strongly restricted. In Models 2.1.1, 2.1.2, 2.2.1 and 2.2.2 the onlypoint which can occur is the origin. In Model 2.3.1 the analogue of this is asituation where all substrates are exhausted and the enzymes are completelyin the unbound form. In Model 2.4.1 the corresponding situation is that allconcentrations except that of ATP are zero and the concentration of ATP takeson its maximal value.The process of making a species into an external species can be illustratedby showing how Model 2.1.1 can be obtained as a limit of Model 2.4.1. Theevolution equations of Model 2.4.1 are dx RuBP dt = k x Ru5P x ATP − k x RuBP , (1) dx PGA dt = 2 k x RuBP − k x PGA x ATP − k x PGA , (2) dx DPGA dt = k x PGA x ATP − k x DPGA , (3) dx GAP dt = k x DPGA − k x − k x GAP , (4) dx Ru5P dt = − k x Ru5P x ATP + 3 k x , (5) dx ATP dt = − k x PGA x ATP − k x Ru5P x ATP + k ( c A − x ATP ) (6)where the constant c A is the total concentration of adenosine phosphates. Wewould now like to consider a situation where ATP is in excess. Equivalentlywe can consider a situation where the concentrations of all substances exceptATP are very small. Define x X = η ˜ x X for each substance X except ATP. Define˜ k = η k . Then making these substitutions gives d ˜ x RuBP dt = k ˜ x Ru5P x ATP − k ˜ x RuBP , (7) d ˜ x PGA dt = 2 k ˜ x RuBP − k ˜ x PGA x ATP − k ˜ x PGA , (8) d ˜ x DPGA dt = k ˜ x PGA x ATP − k ˜ x DPGA , (9) d ˜ x GAP dt = k ˜ x DPGA − k ˜ x − k ˜ x GAP , (10)6 ˜ x Ru5P dt = − k ˜ x Ru5P x ATP + 3˜ k ˜ x , (11) d ˜ x ATP dt = − ηk ˜ x PGA x ATP − ηk ˜ x Ru5P x ATP + k ( c A − x ATP ) . (12)Letting η tend to zero gives a system for which x ATP = c A is an invariant mani-fold and the restriction of the system to that manifold reproduces the equationsof Model 2.1.1. This is a regular limit and it follows from the existence of anunstable hyperbolic positive steady state in the limiting system that Model 2.4.1also has an unstable hyperbolic positive steady state. Note that the perturbedsteady state for η small and positive does satisfy x ATP < c A since otherwiseequation (12) would lead to a contradiction. For this system the informationabout a steady state obtained by the perturbation argument is less than what isalready known by analysing the full system directly. The argument has neverthe-less been presented here since analogous arguments may be useful for obtaininginformation about more complicated systems where no alternative is available. The models considered in this section involve more unknowns than those ofthe previous section. The starting point is a model introduced in a paperby Pettersson and Ryde-Pettersson [14] which we refer to for brevity as thePettersson model. The substances included are roughly speaking those whichCalvin found to appear after the dark reactions had run for a few minutes. Inaddition to those in the five-variable models these are dihydroxyacetone phos-phate (DHAP), fructose 1,6-bisphosphate (FBP), fructose 6-phosphate (F6P),erythrose 4-phosphate (E4P), sedoheptulose 7-phosphate (S7P), sedoheptulose1,7-bisphosphate (SBP), xylulose 5-phosphate (X5P) and ribose 5-phosphate(R5P). In addition the process by which sugars can be stored in the chloro-plast as starch is included. Starch itself is treated as an external species. Theintermediates glucose 6-phosphate (G6P) and glucose 1-phosphate (G1P) areincluded as internal species. In contrast to the models of the previous sectioninorganic phosphate in the chloroplast, P i , is modelled dynamically, as is ATP.On the other hand NADPH is still treated as an external species. Some of thereactions which were treated as irreversible in the models of the previous sectionare treated as reversible in the Pettersson model, for instance the reaction inter-converting PGA and DPGA. The decision, which reactions to treat as reversibleand which as irreversible in the Pettersson model is based on experimental data.The only reactions treated as irreversible are those whose substrates are Ru5P,RuBP, FBP and G1P together with the transport processes to the cytosol andto starch.The Pettersson model (Modell 3.1.1) will now be described. In [14] the timederivatives of the relevant concentrations are expressed in terms of the rates v i of the different reactions. The equations are dx RuBP dt = v − v , (13)7 x PGA dt = 2 v − v − v PGA , (14) dx DPGA dt = v − v , (15) dx ATP dt = v − v − v − v st , (16) dx GAP dt = v − v − v − v − v − v GAP , (17) dx DHAP dt = v − v − v − v DHAP , (18) dx FBP dt = v − v , (19) dx F6P dt = v − v − v , (20) dx E4P dt = v − v , (21) dx X5P dt = v + v − v , (22) dx SBP dt = v − v , (23) dx S7P dt = v − v , (24) dx R5P dt = v − v , (25) dx Ru5P dt = v + v − v , (26) dx G6P dt = v − v , (27) dx G1P dt = v − v st , (28) dx P i dt = v + v + v + v PGA + v GAP + v DHAP + 2 v st − v . (29)The total amount of phosphate in the chloroplast is a conserved quantity andmay be used to eliminate the concentration of inorganic phosphate in the chloro-plast from the equations in favour of the other variables. It is assumed that thereversible reactions are much faster than the irreversible ones. This can beimplemented mathematically by introducing a small parameter ǫ and defining˜ v i = ǫv i for the fast reactions. The slow reactions, whose rates are not rescaled,are those with reaction rates v , v , v , v , v , v PGA , v GAP , v DHAP and v st .For each of the slow reactions an explicit phenomenological expression is givenfor the rate. This incorporates the known experimental facts on the activationand inhibition of certain reactions due to the influence of other substances. Noexpressions are given for the rates of the fast reactions. Instead it is assumedthat these reactions can be taken to be in equilibrium, which gives algebraicequations relating the concentrations. It will be shown below how this can be8mplemented mathematically.At this point we interrupt the discussion of the Pettersson model and in-stead take the reaction network underlying the Pettersson model including itsstoichiometry and apply mass action kinetics to get something which was calledthe Pettersson-MA model in [12] (Model 3.2.1). The strategy adopted in [12]was to study the dynamics of Model 3.2.1 so as to try to obtain some insightsfor tackling the Pettersson model later. There is a related model with an addi-tional reaction which liberates G1P from starch (Model 3.2.2). There the aboveevolution equations are modified by adding a contribution v to the evolutionfor x G1P and a contribution − v to the evolution equation for x P i . This willbe called the Poolman-MA model since a modification of the Pettersson modelincluding a mechanism of this type was first introduced by Poolman [15]. Pool-man himself used the same reaction rates as in [14] for the slow reactions andtreated the liberation of G1P as a slow reaction while taking mass action ki-netics for the fast reactions. We call the resulting system of ODE the Poolmanmodel (Model 3.1.2). A hybrid model can be obtained by taking the reactionsincluded in the Pettersson model with the reaction rates as in the Poolmanmodel (Model 3.3.1). This is obtained from the Poolman model by setting oneof the reaction constants k to zero. In these models the total amount ofphosphate is conserved and every unknown contains some phosphate. Thus theconservation law implies that all solutions are bounded and runaway solutionsare ruled out for these models. It was already mentioned that the export of sug-ars from the chloroplast is coupled to the import of inorganic phosphate (whoseconcentration in the cytosol is assumed constant in the model). It is indicatedin [14] that if the external concentration of phosphate is too high then no pos-itive steady state will exist. This is the phenomenon of overload breakdown.Poolman suggested that overload breakdown could be avoided by introducingthe release of G1P from starch. In the case of Model 3.2.1 it was shown in [12]that if k c A ≤ k there exists a Lyapunov function related to the function L of the last section and this proves that under this condition Model 3.2.1 has nopositive steady states. Here the k i are reaction constants and c A is the totalconcentration of adenosine phosphates. For Model 3.2.2 this construction nolonger works. When k c A > k in Model 3.2.1 it is possible to construct ananalogue of the function L of the last section which gives the conclusion thatthere exist no positive steady states where L is less than a certain numberdepending only on the reaction constants. Here we define L = L −
12 ( x DPGA + x GAP + x DHAP ) . (30)and it satisfies d (5 L ) dt = 12 (2 k x DHAP x GAP + k x FBP x GAP + k x E4P x DHAP + k x S7P x GAP − k x GAP − k x DHAP − k x E4P x X5P − k x X5P x R5P − x PGA − k x FBP − k x SBP ) . (31)When L is sufficiently small the positive terms on the right hand side are9ominated by the negative ones. Information can also be obtained for Model3.3.1 by using the function L . In that case L is decreasing provided thequantity k x ATP − v PGA is negative. Now v PGA = V ex x PGA NK PGA where N = 1 + (cid:18) K P ext x P ext (cid:19) (cid:18) x P i K P i + x PGA K PGA + x PGA K GPA + x DHAP K DHAP (cid:19) (32)and the other quantities which have not previously been defined are positiveconstants. Here we treat the total amount of phosphate as a parameter andthen if k is chosen small enough for fixed values of the other parameters inthe kinetics we get a positive lower bound for x − v PGA . Thus under theseconditions L is decreasing. In particular Model 3.3.1 has no positive steadystates when the parameters are restricted in this way.In [12] conditions were derived for ω -limit points of positive solutions ofModels 3.2.1 and 3.2.2. It was pointed out in [12] that many of the argumentsused apply to the original Poolman model since the only property of the reactionrates which is used is under what circumstances they are positive or zero and thisis not changed when the mass action kinetics is replaced by the more complicatedkinetics of the Poolman model. The same argument applies to the hybrid model.It thus follows from the arguments in [12] that Models 3.1.2, 3.2.1, 3.2.2 and 3.3.1have the property that the only substances whose concentrations may fail tovanish at an ω -limit point of a positive solution where at least one concentrationvanishes are G1P, G6P, F6P, E4P, S7P and P i .In [12] information was also obtained on how these points may be approachedby positive solutions of Models 3.2.1 and 3.2.2. This is done by linearizing aboutsteady states where some concentrations are zero and analysing the eigenvaluesof the linearization. In some cases spectral stability could be determined butother cases remain open. In many cases where the spectral analysis was success-ful it turned out that the centre manifold coincides with the center subspace andthe qualitative behaviour on the centre manifold could be analysed. In othercases the centre manifold is nonlinear and its Taylor expansion not yet beencomputed.Consider now again the evolution equations for the Pettersson model ex-pressed in terms of the reaction rates v i . In [14] five linear combinations y i ofconcentrations are identified whose time derivatives only depend on the slowreaction rates. Suppose we now complement these by a suitable set of concen-trations z i , for instance all those except x RuBP , x F6P , x Ru5P , x DHAP and x ATP ,which we denote by s i . Then the concentrations of the internal species are re-lated to the variables y i and z i by an invertible linear transformation. Considerthe equations of the hybrid model, expressed in terms of the variables y i and z i . If we write them in terms of the ˜ v i then all the terms on the right hand sideof the evolution equations for the y i are regular in ǫ while many of those on theright hand side of the equations for the z i contain a factor ǫ − . Multiplyingthese equations with ǫ and letting ǫ tend to zero gives a system of algebraicequations. When expressed in terms of the ˜ v i these equations are linear andthey imply that the ˜ v i vanish for all the fast reactions at ǫ = 0. This fact can be10ead off from the subnetwork obtained by deleting the slow reactions from thefull Pettersson network. Now the ˜ v i can be obtained from v i by replacing thereaction constants k i by ˜ k i = ǫk i . If the ˜ k i are chosen independent of ǫ the al-gebraic equations (20)-(30) of [14] for the concentrations are obtained. Thus inthe limit ǫ → y i and z i . This is what we refer to as the Petterssonmodel. Without further information it is not even clear that local existenceholds for this system. In [14] it is claimed that the equations (20)-(30) of thatpaper can be used to obtain an explicit closed system of evolution equations forthe variables s i but the calculations presented there are not complete. A similarreduction to a DAE can be carried out for the Poolman model and we call theresult the reduced Poolman model (Model 3.3.2).Any of the models considered in this section may be modified so as to makeATP and inorganic phosphate external species. By analogy with what was donein the case of Model 2.4.1 we can fix the concentration of ATP and c A and rescalethe concentrations of the other substances by a factor η . The fact of havingeliminated the concentration of inorganic phosphate by using the conservationof the total amount of phosphate means that x P i then automatically becomesconstant. In the same way as k had to be scaled by a power of η in Model 2.4.1 itis necessary to rescale the reaction constants in the reactions with two substratesin order to get a non-zero limit. For instance, using the notation of [12] weintroduce ˜ k = ηk . The other reaction constants which should be rescaled in asimilar manner are k , k , k , k and k . To maintain the non-trivial effectsof saturation, activation and inhibition in the slow reactions it is also nessaryto rescale certain Michaelis constants. The constants involved are K m , K i , K m , K mst , K m , K i , K i , K i , K i , K m , K i , K i , K ast , K ast , K ast , K PGA , K GAP and K DHAP . Let us call the results of modifying Models3.1.1 and 3.3.1 in this way Models 3.4.1 and 3.4.2 respectively. In both of thesecases the modified model has an invariant manifold x ATP = c A for η = 0 and therestriction of the system to that submanifold reproduces the model with ATP asan external species. As in the discussion of Model 2.4.1 in the last section thismeans that any hyperbolic positive steady state of Model 3.4.1 or 3.4.2 givesrise to a hyperbolic positive steady state of Model 3.1.1 or 3.3.1, respectively.Thus information on steady states can be obtained from information on thecorresponding models with the concentrations of ATP and P i frozen. Five of the equations in the Pettersson model are evolution equations. The righthand sides of these equations are the functions F i in the equations (54)-(58) of[14] and their vanishing is equivalent to the equations (42)-(47) of [14], whichare, with the definition v = v , v = v, (33)11 = v v st , (34) v = v , (35) v = v, (36) v = 3 v + v st − v PGA , (37) v = 3 v ex + 6 v st (38)where v ex = v PGA + v GAP + v DHAP . In this section we concentrate on Model3.4.1, where all these equations hold except that containing v . Suppose that x DHAP is given. Then ˜ v = 0 fixes the value of x GAP . Then ˜ v = 0 and ˜ v = 0fix the values of x DPGA and x PGA . In addition ˜ v = 0 fixes the value of x FBP .With the information we have it is possible to compute v ext in terms of x DHAP .Suppose now that v st is also fixed. Then with this information it is possible toobtain v and hence v and v . The quantities x SBP , x RuBP and x Ru5P are thenuniquely determined. To ensure the existence of these quantities it suffices toassume that the parameters V , V and V are large enough. The equations˜ v = 0 and ˜ v = 0 fix x R5P and x X5P . Next the equation ˜ v = 0 allows x E4P to be determined. Then ˜ v = 0 can be used to determine x F6P and ˜ v = 0 and˜ v = 0 give x G6P and x G1P . This leaves two consistency conditions, namelythe equation for v and the expression for v st in terms of x G1P . Let these bedenoted abstractly by Φ ( x DHAP , v st ) = 0 and Φ ( x DHAP , v st ) = 0 respectively.A general strategy for looking for positive steady states of the evolutionequations defined by a given network is to look at limiting cases where someof the reaction constants are set to zero and the smaller network obtained bydiscarding the reactions concerned. Trying to do this for a larger network with-out some guiding principles may fail because there are too many possibilities.A concept which can be used as a guiding principle is that of elementary fluxmodes. There is a theory of how to produce steady states using these objects [4]but an alternative is to use elementary flux mode to guess which reaction con-stants to set to zero and then proceed directly with the construction of steadystates. This possibility was used in [17] to give an existence proof of steadystates of Model 2.3.1 and it will also be applied in what follows. The elemen-tary flux modes computed are not required for the proofs themselves but theyhelp to put those proofs into context. They helped to find the proofs and thisapproach might also turn out to be useful for analysing other similar models inthe future.Consider the equations satisfied by the fluxes in steady states of a modeldefined by the Pettersson network. These are2 v − v − v PGA = 0 , (39) v = v , (40) v − v − v − v − v − v GAP = 0 , (41) v − v − v − v DHAP , (42) v = v , (43)12 − v − v = 0 , (44) v = v , (45) v = v , (46) v = v , (47) v + v − v = 0 , (48) v = v , (49) v + v − v = 0 , (50) v = v , (51) v = v , (52) v = v st , (53) v − v − v − v st = 0 . (54)The solutions of these linear equations can be parametrized with the help ofstoichiometric generators. The relevant terminology will now be recalled (cf.[4]). The system of ODE arising from a reaction network can be written inthe form ˙ x = N v ( x ), where N is the stoichiometric matrix and v ( x ) are thereaction rates. In this context reversible reactions are treated as two separatereactions. This means that there are two columns of N corresponding to eachreversible reaction. The column corresponding to the forward reaction is minusthe column corresponding to the backward reaction. If we discard one of thetwo columns corresponding to each reversible reaction a truncated matrix ¯ N isobtained. The kernels of the matrices N and ¯ N are related in a simple waywhich will be described below. The set of reaction rates at a steady state isan element of the kernel of N with non-negative components. We can think ofthis as a point in the space of real-valued functions on the set R of reactions.The set of all non-negative elements of the kernel of N is a positive cone andthus consists of all vectors of the form P i a i w i with a i non-negative coefficientsand w i a finite number of vectors which in this context are called elementaryflux modes [18], [4]. Each of these vectors has the property that setting somebut not all of its components to zero gives a vector which is not in the kernelof N . Another important quantity is the incidence matrix. It has one row foreach complex (quantity on the left or right hand side of a reaction) and onecolumn for each reaction. The element for the left hand side of the reaction is −
1, the element for the right hand side is +1 and all other elements are zero.A vector which is in the kernel of the incidence matrix is in the kernel of N .Elementary flux modes which are not in the kernel of the incidence matrix arecalled stoichiometric generators. For a reversible reaction let us make a choiceof which is the forward direction, so as to get a forward reaction r + and abackward reaction r − . Then the vector which has the components +1 at r + , − r − and all other components zero belongs to the kernel of the incidencematrix. Let us call this a trivial mode. It is an elementary flux mode which isnot a stoichiometric generator. The kernel of N is the joint span of the kernelof ¯ N and the trivial modes.In the case of the Pettersson network the following vectors with components13 i are stoichiometric generators[3 6 6 3 1 1 1 1 1 1 1 2 3 0 0 9 0 0 1 0] , (55)[3 6 6 2 1 1 1 1 1 1 1 2 3 0 0 9 0 1 0 0] , (56)[3 5 5 2 1 1 1 1 1 1 1 2 3 0 0 8 1 0 0 0] , (57)[6 12 12 5 3 3 2 2 2 2 2 4 6 1 1 19 0 0 0 1] . (58)Here the components are written in the order[ v v v v v v v v v v v v v v v v v PGA v GAP v DHAP v st ] . (59)We only write out the components corresponding to the forward reactions. Inother words these vectors belong to the kernel of ¯ N . To get the full mode, whichis an element of the kernel of N , it would be necessary to add zeroes for thebackward reactions. It is easily checked that the vectors defined by (55)-(58) aresolutions of the equations (39)-(54). It can also be shown that any solution forwhich the last four components are zero is zero. Furthermore, any non-negativesolution for which at least one of the last four components in non-zero is a linearcombination of the four solutions above and hence all the components exceptthe last four are zero. This verifies the defining property of elementary fluxmodes that a vector obtained by setting some, but not all, of the componentsof the generator to zero is not a solution. These vectors do not belong to thekernel of the incidence matrix and hence are stoichiometric generators. Finally,all solutions are linear combinations with non-negative coefficients of these gen-erators. This follows from the facts that the generators are linearly independentand that the dimension of the kernel of ¯ N is four. Each of the generators isobtained by shutting off all but one of the output reactions.Each stoichiometric generator defines a subnetwork obtained by setting thosereaction rates to zero for which the corresponding component of the generatoris zero. If mass action kinetics are being considered the desired reaction ratescan be set to zero by setting the corresponding reaction constants to zero. Forother kinetics some more thought is necessary. In most of the slow reactionswe can set v i to zero by setting the corresponding coefficient V i to zero. Theexceptions are the transport reactions to the cytosol. There we have threereaction rates but only one coefficient V ex . Here the desired reaction rates, v X can be set to zero by setting K − X formally to zero. In other words, in the limitcertain constants K X tend to infinity. It also turns out to be helpful to setthe quantity K − i to zero in the limit. If we make the same assumptions onthe kinetics as in the Pettersson model we get a system of DAE correspondingto the subnetwork. Call the system of this type obtained from the first of thefour generators listed above Model 4.1.1. Concretely, it is obtained from thePettersson model by setting the reaction constants k , k , k , k , k , k and k to zero and discarding the variables x G1P , x G6P . A limit is consideredwhere these reaction constants are multiplied by a small constant ζ and theconstants K PGA and K GAP and K i are multiplied by ζ − . Biologically thiscorresponds to a situation where PGA and GAP not only fail to be exported14ut even fail to bind to the transporter and thus do not compete with DHAP.In a similar way it is possible to obtain an analogue of the hybrid model for thesubnetwork. Call it Model 4.2.1. We can freeze the concentrations of ATP andP i in Models 4.1.1 and 4.2.1 to get Models 4.1.2 and 4.2.2. Model 4.1.2 is closeto a modelling approach used in [13] although in that paper no complete systemof equations was written out. The steady states of Model 4.1.2 can be studiedby following calculations in [13]. There are two differences between Model 4.1.2and the situation in [13]. One of these corresponds to setting the coefficients K − i to zero in the expression for v in [14] while the other has to do with thefact that v PGA and v GAP are absent from Model 4.1.2.The equations (33)-(38) for the reaction rates in the Pettersson model aremodified in the subnetwork by the removal of v st and v PGA . For ζ = 0 we haveΦ = v st . Following the computations done above we find that for ζ = 0Φ = Ax DHAP B + Cx DHAP − Dx E + F x (60)where A = V ex , B = K DHAP (cid:16) x Pext K Pext (cid:17) x Pi K Pi , C = K DHAP + 1 + x Pext K Pext , D = V k k k k , E = 1 + K − i x P i and F = k k k k . Using the positivity of the unknownwe see that the equation Φ = 0 is equivalent to the quadratic equation( AF − CD ) x − BDx
DHAP + AE = 0 . (61)This equation has two positive solutions precisely when AF − CD >
AE < B D AF − CD ) . Moreover in that case the derivative of Φ is non-zero ateach of those points. Parameters can be chosen such that these inequalities aresatisfied. For instance, starting from arbitrary positive values of the parameters V can be reduced so as to ensure that the first inequality is satisfied. Then k can be increased to arrange that the second one is satisfied. Perturbing ζ awayfrom zero and applying the implicit function theorem we see that there existtwo positive steady states of Model 3.4.1 for suitable choices of the parameters.This implies in turn the existence of two positive steady states for the Petterssonmodel. Summing up, we get the following result Theorem
There are choices of the parameters for which the Pettersson modelhas at least two positive steady states.
The equations for the steady state fluxes in the Poolman model are similarto those in the case of the Pettersson model. If we make use of the conser-vation law for the total amount of phosphate then the only difference is anadditional summand v in the evolution equation for G1P. Note that this re-action rate belongs to a slow reaction. The explicit expression for this rate is v = V x Pi x Pi + K m (cid:16) x G1P Ki (cid:17) (cf. [15], equation (4.4)). The equations (54)-(58) in1514] are replaced by F = v − v , (62) F = v − v − v st + v , (63) F = v + 2 v − v − v st + v , (64) F = 2 v + v st − v ex − v − v − v , (65) F = v + v PGA − v − v − v st . (66)Correspondingly equations (42)-(47) of [14] are replaced by v = v, (67) v = v v st − v , (68) v = v , (69) v = v, (70) v = 3 v + v st − v PGA , (71) v = 3 v ex + 6 v st − v . (72)Let us freeze ATP and P i in the reduced Poolman model and call the resultModel 4.1.1. In the frozen model v is a function of x G1P alone.We now pass to a limit in a similar way to what was done for the Petterssonmodel, setting the reaction rates v PGA , v PGA and v st to zero. This time weallow v to remain non-zero. The calculations are simplified by setting K − i tozero, so that the expression for v reduces to a constant. Then equation (60)is replaced by Φ = Ax DHAP B + Cx DHAP − Dx E + F x − G (73)where G = V x Pi x Pi + K m . If we start from a choice of parameters which gives twopositive steady states in the Pettersson model and perturb G from being zeroto being positive and sufficiently small then the function Φ has three non-degenerate zeroes. For a sufficiently small perturbation of this type does notdestroy the positive zeroes present for G = 0 and does not make them becomedegenerate. On the other hand is makes the value of Φ at the origin negativeand this leads to a new positive zero. This zero is a deformation of the zerowhich lies at the origin for G = 0 and so it is non-degenerate for small parametervalues. It follows by arguments similar to those in the last section that for theparameter values just considered the reduced Poolman model has three positivesteady states.Stoichiometric generators for the Poolman system have been studied in [16]and we now consider the relation of these to the construction just carried out.There are modes generalizing those for the Pettersson model by augmentingthem by a zero entry for v . Strangely these do not seem to fit with Figure2A in [16] where the reaction with flux v is also shown as being shut off. Thatfigure appears to contradict equation (9) of [14] which says that in a steady16tate v = v + v . We can now proceed as in the last section, with v beingset to zero. We use the mode[3 6 6 3 1 1 1 1 1 1 1 2 3 0 0 10 0 0 1 0 0] . (74)That takes us to the same subnetwork as before and the analogous argumentsshow that there are parameter values for which the reduced Poolman model hastwo positive steady states. As has been shown above the possibility of having v non-zero to get results for the reduced Poolman model which go beyondthose obtained for the Pettersson model.To obtain these modes it is necessary to be careful about the differencebetween N and ¯ N . This time we will make a different choice of which reactionsare considered to be in the forward direction. This results in the reaction rates v and v being replaced by their negatives. This is related to the fact that inthis mode material is flowing from starch to the sugars in the cycle. Considerthe mode [3 6 6 4 0 0 1 1 1 1 1 2 3 1 1 9 0 0 1 0 1] . (75)In this case we get a subnetwork with an inflow from starch but no outflow tostarch.Note finally that in this paper we have not obtained any results on multiplestationary solutions for the Pettersson model itself or for the hybrid model. Inorder to do this it would be necessary to obtain information on the transverseeigenvalues, showing that they are all different from zero, at least for some valuesof parameters for which multiple steady states exist for the reduced models. In this paper we have been concerned with a variety of mathematical models forone biological system, the Calvin cycle. The approach has been to understandas much as possible about the relations between the different models and toobtain as much information as possible about the qualitative behaviour of thesolutions of the equations defined by these models. The scope was restricted toresults obtained by purely analytical and rigorous methods without any appealto numerics or reliance on heuristics. These results are piecemeal and shouldbe complemented by a better conceptual understanding of the key mechanismsdetermining the behaviour. To do this it makes sense to look at models whichare as simple as possible.The approach of studying the simplest possible model has been pursuedby Hahn [8]. His three-variable model includes the important phenomenon ofphotorespiration which is not included in the models discussed up to now inthis paper. Because this seems to the author to be an important direction forfuture developments the work of [8] will now be discussed briefly. The unknownsare x RuBP , x PGA and x TP . TP stands for ’triose phophate’ and compared tothe five-variable models it is obtained by lumping together x DPGA and x GAP . x Ru5P has been considered as an intermediate species and discarded. Let us17gnore photorespiration for the moment ( k = 0 in the notation of [8]). We alsoignore the reaction called dark respiration ( k = 0).The reaction from RuBP to PGA is as in Model 2.1.1 and there is a reactiontaking PGA to TP replacing that from PGA to DPGA in Model 2.1.1. There isa sink reaction starting at TP. There is a reaction from TP to RuBP replacingthat from GAP to Ru5P. The evolution equations are dx RuBP dt = − k x RuBP + 3 k x , (76) dx PGA dt = 2 k x RuBP − k x PGA , (77) dx TP dt = k x PGA − k x − k x TP . (78)It is emphasized by Hahn that a key property which is expected from a modelis that it should have a stable positive steady state. He states that k canreasonably be estimated from experimental data but that k , k and k cannot.Thus he adopts the following strategy. It is assumed that a stable steady stateexists in the model and the concentrations of the substances involved which canbe measured under suitable circumstances are assumed to be the values in thesteady state. Assuming in this way certain values for the coordinates of thesteady state the three remaning reaction constants can be calculated. Havingobtained these values we can then ask if for those parameters there are otherpositive steady states.The system above has exactly one positive steady state for any values of theparameters. It satisfies x TP = (cid:16) k x RuBP k (cid:17) , x PGA = k x RuBP k and x RuBP = "(cid:18) k k (cid:19) k k . (79)Note that in contrast to Model 2.1.1 this model always has a positive steadystate for any choice of the parameters. The difference has to do with the factthat in the Hahn model there is no sink for x PGA . It is shown in [8] that thissteady state is unstable.In the case with photorespiration it is shown in [8] how to reduce the problemof finding steady states to that of solving a ninth degree equation for one of theconcentrations. Numerically it was found that for the biologically motivatedvalues of the parameters this model has two steady states. Moreover one ofthese is stable and the other unstable. This is similar to the results which haveproved for other models discussed in the previous sections. In the model of[8] this situation is only possible in the presence of photorespiration. Howeverthe other models suggest that in this respect Hahn’s three-variable model is notrepresentative of what happens in more detailed models. It is desirable to obtainmore analytical results on the Hahn models and their relations to other moredetailed models of the Calvin cycle. Note that in [7] Hahn had also previouslystudied some larger models which have rise to the models of [8] by a process ofsimplification. 18 eferenceseferences