The minimal model of Hahn for the Calvin cycle

TL;DR

This paper analyzes a minimal mathematical model of the Calvin cycle, demonstrating how photorespiration can stabilize the cycle's operation and providing rigorous proofs of the model's steady states and their stability.

Contribution

The paper provides the first rigorous analysis of Hahn's minimal Calvin cycle model, including stability of steady states with and without photorespiration.

Findings

In the no-photorespiration model, exactly one unstable positive steady state exists.

In the photorespiration model, two positive steady states exist, one stable and one unstable.

Photorespiration can stabilize the Calvin cycle, preventing concentrations from diverging.

Abstract

There are many models of the Calvin cycle of photosynthesis in the literature. When investigating the dynamics of these models one strategy is to look at the simplest possible models in order to get the most detailed insights. We investigate a minimal model of the Calvin cycle introduced by Hahn while he was pursuing this strategy. In a variant of the model not including photorespiration it is shown that there exists exactly one positive steady state and that this steady state is unstable. For generic initial data either all concentrations tend to infinity at lates times or all concentrations tend to zero at late times. In a variant including photorespiration it is shown that for suitable values of the parameters of the model there exist two positive steady states, one stable and one unstable. For generic initial data either the solution tends to the stable steady state at late times or all concentrations tend to zero at late times. Thus we obtain rigorous proofs of mathematical statements which together confirm the intuitive idea proposed by Hahn that photorespiration can stabilize the operation of the Calvin cycle. In the case that the concentrations tend to infinity we derive formulae for the leading order asymptotics using the Poincar\'e compactification.

Equations

$$\frac{dx}{dt}=-\alpha x-2\beta x^2+3y^2,$$

$$\frac{dy}{dt}=2\alpha x+3\beta x^2-5y^2-y.$$

$$\left[{\begin{array}[]{cc}-\alpha-4\beta x&6y\\ 2\alpha+6\beta x&-10y-1\\ \end{array}}\right]$$

$$G_1=(-,-)=\{(x,y):x>\frac{3}{\alpha }{y}^2,x<\frac{1}{2\alpha }(5y^2+y),$$

$$G_2=(+,-)=\{(x,y):x<\frac{3}{\alpha }{y}^2,x<\frac{1}{2\alpha }(5y^2+y),$$

$$G_3=(-,+)=\{(x,y):x>\frac{3}{\alpha }{y}^2,x>\frac{1}{2\alpha }(5y^2+y),$$

$$G_4=(+,+)=\{(x,y):x<\frac{3}{\alpha }{y}^2,x>\frac{1}{2\alpha }(5y^2+y).$$

$$N_1=\{0,0\}\cup (0,-)\cup S_1\cup (0,+),$$

$$N_2=\{0,0\}\cup (-,0)\cup S_1\cup (+,0).$$

$$\frac{dX}{d\tau }=\alpha X^2-3X{Z}^2,$$

$$\frac{dZ}{d\tau }=2\alpha X+(\alpha -1)XZ-5Z^2-3Z^3.$$

$$\frac{dX}{d\tau }=2u\frac{du}{d\tau }=\alpha u^4-3u^4{v}^2$$

$$\frac{du}{d\tau }=\frac{1}{2}(\alpha u^3-3u^3{v}^2).$$

$$\frac{dZ}{d\tau }=u\frac{dv}{d\tau }+v\frac{du}{d\tau }=2\alpha u^2-5u^2{v}^2+(\alpha -1)u^{3}v-3u^3{v}^3$$

$$\frac{dv}{d\tau }=2\alpha u-5u{v}^2-(\alpha -1)u^{2}v+3u^2{v}^3-\frac{1}{2}(\alpha u^{3}v-3u^3{v}^3).$$

$$\frac{du}{ds}=\frac{1}{2}(\alpha u^2-3u^2{v}^2),$$

$$\frac{dv}{ds}=2\alpha -5v^2+(\alpha -1)uv-3u{v}^3-\frac{1}{2}(\alpha u^{2}v-3u^2{v}^3).$$

$$\frac{dZ}{d\tau }=\frac{dv}{d\tau }=2\alpha u{v}^2-5v^2+(\alpha -1)u{v}^3-3v^3.$$

$$\frac{dX}{d\tau }=v^2\frac{du}{d\tau }+2uv\frac{dv}{d\tau }=\alpha u^2{v}^4-3u{v}^4.$$

$$\frac{du}{d\tau }=\alpha u^2{v}^2-3u{v}^2-4\alpha u^{2}v+10uv-2(\alpha -1)u^2{v}^2+6u{v}^2.$$

$$\frac{du}{ds}=10u+3uv-4\alpha u^2-(\alpha -2)u^{2}v,$$

$$\frac{dv}{ds}=-5v+2\alpha uv+(\alpha -1)u{v}^2-3v^2.$$

$$\frac{dw}{ds}=-10w-4\alpha w^2+3wv+\frac{15}{2\alpha }v-(\alpha -2){\left(w+\frac{5}{2\alpha }\right)}^{2}v,$$

$$\frac{dv}{ds}=2\alpha wv+(\alpha -1)\left(w+\frac{5}{2\alpha }\right)v^2-3v^2.$$

$$\frac{dw}{ds}=-10w+\frac{5(\alpha +10)}{4{\alpha }^2}v-4\alpha w^2-\frac{2(\alpha -5)}{\alpha }wv-(\alpha -2)w^{2}v.$$

$$\frac{dv}{ds}=\left[\frac{\alpha +10}{4\alpha }+\frac{5\alpha -5}{2\alpha }-3\right]v^2+\dots$$

$$=-\frac{1}{4}{v}^2+\dots$$

$$\beta (y^2+2y{)}^2={\alpha }^2(y^2-y)$$

$$y[\beta y^3+4\beta y^2+(4\beta -{\alpha }^2)y+{\alpha }^2]=yp(y)=0.$$

$$\Delta ={\alpha }^2\beta [96{\beta }^2-131{\alpha }^2\beta +4{\alpha }^4].$$

$$q(y)=4\beta y^3+12\beta y^2+(8\beta -2{\alpha }^2)y+{\alpha }^2=0.$$

$$9\gamma (-4{\gamma }^2+131\gamma -96)=0$$

$$\frac{dx}{dt}=-\alpha x-2\beta x^2+3y^5,$$

$$\frac{dy}{dt}=2\alpha x+3\beta x^2-5y^5-y.$$

$$\frac{dX}{d\tau }=\alpha X^5-3X{Z}^5,$$

$$\frac{dZ}{d\tau }=2\alpha X^4+(\alpha -1)X^{4}Z-5Z^5-3Z^6.$$

Full text

1 Introduction

The Calvin cycle is a part of photosynthesis. There are many mathematical models for this biochemical system in the literature. Reviews of these can be found in [1], [2] and [11]. This is an interesting example in which the relations between different mathematical models for the same biological situation can be investigated. A mathematical comparison of a number of these models was carried out in [18]. There it was pointed out that it would be desirable to look more closely at the minimal model of the Calvin cycle introduced by Hahn [9]. In fact Hahn’s paper contains several related systems of ordinary differential equations of dimensions two and three and the aim of the present paper is to obtain an understanding of the dynamics of the two-dimensional models of Hahn which is as complete as possible. There is also a brief discussion of the relation of the two-dimensional models to the three-dimensional one.

The function of the Calvin cycle is to use carbon dioxide to produce sugars. This process is fuelled by ATP and NADPH produced in the light reactions of photosynthesis where the energy contained in light is captured as chemical energy. A comprehensive introduction to the biochemistry of photosynthesis can be found in [10]. The main step in the Calvin cycle, resulting in the production of PGA (phosphoglycerate), is catalysed by the enzyme Rubisco. Interestingly this enzyme has a dual functionality. It can not only catalyse the reaction of carboxylation, which is the primary way in which carbon dioxide is fixed in the Calvin cycle, but also an oxidation reaction. This second reaction competes with the first and reduces the efficiency with which the Calvin cycle produces sugar. The reason for the existence of this apparently wasteful alternative reaction is not clear. One possible explanation, for which the Hahn model is relevant, is that photorespiration stabilizes the system - it creates the possibility of the existence of a stable positive steady state.

The paper is organized as follows. In section 2 the two-dimensional system of Hahn is introduced. In dimensionless form the equations depend on two non-negative parameters $\alpha$ and $\beta$. The case $\beta>0$ corresponds to including photorespiration in the model. The dynamics of the models is first analysed in the case without photorespiration ($\beta=0$). The main result is Theorem 1 which describes the global asymptotic behaviour of general solutions in detail. There exists a unique positive steady state $S_{1}$ which is unstable. For an open set of initial data which is described in detail all concentrations tend to zero as $t\to\infty$. For another open set of initial data all concentrations tend to infinity as $t\to\infty$. The complement of the union of these two sets is the stable manifold of the steady state $S_{1}$. A formula is derived for the leading order asymptotics of the solutions which tend to infinity. In section 3 the case $\beta>0$ is treated. The main result is Theorem 2. In one open set of parameter space, for which an explicit formula is given, all solutions have the property that the concentrations tend to zero as $t\to\infty$. In the interior of the complement of that set more interesting behaviour is observed. There are two positive steady states, one stable and one unstable. For an open set of initial data all concentrations tend to zero as $t\to\infty$. For another open set of initial data the solutions tend to the stable positive steady state as $t\to\infty$.

Many models of the Calvin cycle contain the fifth power of the concentration of the substance GAP (glyceraldehyde phosphate). This is because in the usual coarse-grained descriptions of the Calvin cycle, where many elementary reactions are combined, there is an effective reaction where five molecules of GAP with three carbon atoms each go in and three molecules of a five-carbon sugar come out. Applying mass-action kinetics to this leads to the fifth power. In deriving the model studied in sections 2 and 3 Hahn replaces the fifth power by the second power. His motivation is to make the model analytically more tractable. He assumes implicitly that this change makes no essential difference to the qualitative behaviour of the solutions but gives no justification for this assumption. In section 4 we show that the solutions of the model with the fifth power do indeed behave in a way which is very similar to the behaviour of the model with the second power. The main difference in the analysis is that for the fifth power no explicit formula is obtained for the boundary between the two generic behaviours in parameter space. The results are summarized in Theorems 3 and 4. In [9] the two-dimensional systems are obtained from a three-dimensional one by informal arguments. In section 5 it is shown how the relation between the three-dimensional system and the two-dimensional system with the fifth power can be formalized in a rigorous way using the theory of fast-slow systems. (For an introduction to this theory we refer to [12].) This also gives some limited information about the dynamics of solutions of the three-dimensional system. A full analysis of the three-dimensional system is a much harder problem which is left to future work.

2 The system of Hahn

The system which will be examined in what follows consists of the equations (41)-(42) in [9]:

[TABLE]

In addition to the system with $\alpha>0$ and $\beta>0$, which we will call the full system, we also treat the cases where $\alpha=0$ (no photosynthesis), $\beta=0$ (no photorespiration) or both. Note for future reference that the derivative of the right hand side of this system at the point $(x,y)$ is

[TABLE]

with determinant $\alpha+4\beta x-2\alpha y+4\beta xy$.

Consider first the case $\alpha=\beta=0$. There any solution satisfies the inequality $\frac{dy}{dt}\leq-y$ and thus $y$ decays exponentially at late times. In particular there are no positive steady states. The non-negative steady states are precisely the points on the $x$-axis. Apart from the zero eigenvalue due to the continuum of steady states the other eigenvalue of the linearization about any of these points is $-1$ and this manifold is normally hyperbolic [12]. It follows that given any $x_{0}>0$ there exists a positive solution with $\lim_{t\to\infty}x(t)=x_{0}$. Consider next the case $\alpha=0$, $\beta\neq 0$. Any positive steady state satisfies $y=\sqrt{\frac{2\beta}{3}}x$ by (1). Substituting this into (2) gives $y(-\frac{1}{2}y-1)=0$. Thus there is no positive steady state. The only non-negative steady state is at the origin. In fact $\frac{d}{dt}(3x+2y)=-y^{2}-y$. Thus $3x+2y$ is a strict Lyapunov function on the non-negative orthant away from the origin and it follows that all solutions converge to the origin as $t\to\infty$.

In the case $\alpha\neq 0$, $\beta=0$ we have the inequality $\frac{d}{dt}(5x+3y)\leq\frac{1}{5}\alpha(5x+3y)$, so that all solutions exist globally in the future. Equation (1) shows that for a steady state $x=\frac{3}{\alpha}y^{2}$. Substituting this into (2) gives $y^{2}-y=0$. Thus the steady states are $S_{0}=(0,0)$ and $S_{1}=\left(\frac{3}{\alpha},1\right)$. Now we carry out a nullcline analysis as described in the Appendix. The nullclines are given by $x=\frac{1}{\alpha}y^{2}$ and $x=\frac{1}{2\alpha}(5y^{2}+y)$. These are the graphs of functions of $y$ and it is clear that the complement of the union of the nullclines has four connected components (cf. Figure 1).

These are

[TABLE]

The complement of $S_{1}$ in one of the nullclines has two connected components which can be distinguished by the sign of the time derivative which does not vanish. We write

[TABLE]

Note that if $\dot{x}=0$ at some time then $\ddot{x}=6y\dot{y}$ and that if $\dot{y}=0$ at some time then $\ddot{y}=(2\alpha+6\beta x)\dot{x}$.

Lemma 1 A solution of (1)-(2) belongs to one of the following three cases.

(i) It converges to $S_{0}$ as $t\to\infty$.

(ii) It converges to $S_{1}$ as $t\to\infty$.

(iii) There is a time $t_{1}$ such that it belongs to $G_{4}$ for $t\geq t_{1}$.

Proof Consider a solution which starts at a point on the boundary of $G_{1}$ other than $S_{0}$ or $S_{1}$. If it is on $N_{1}$ then $\dot{y}<0$ and using the equation for $\ddot{x}$ shows that $\dot{x}$ immediately becomes negative. If a solution starts on $N_{2}$ then $\dot{x}<0$ and $\dot{y}$ immediately becomes negative. It follows that any solution which starts in $G_{1}$ remains in $G_{1}$ and any solution which starts at a point on the boundary of $G_{1}$ other than a steady state immediately enters $G_{1}$. Since $y$ is decreasing for solutions in $G_{1}$ any solution which is ever in $G_{1}$ converges to the origin as $t\to\infty$. Consider next a solution which starts in $G_{2}$. Since $y$ is decreasing on $G_{2}$ this solution is bounded. By Lemma A1 of the Appendix it either converges to $S_{0}$ or to $S_{1}$ as $t\to\infty$ or it reaches another point of the boundary of $G_{2}$ after a finite time. In the latter case it reaches a point of the boundary of $G_{1}$ or $G_{4}$ after a finite time. By an analogous argument we can reach a similar conclusion about a solution which starts in $G_{3}$. Either it converges to $S_{0}$ or $S_{1}$ or it reaches another point of the boundary of $G_{1}$ or $G_{4}$ after a finite time. An analysis similar to that done for $G_{1}$ can be carried out for $G_{4}$. A solution which starts in $G_{4}$ must remain there and a solution which starts on the boundary of $G_{4}$ must immediately enter $G_{4}$. Thus any solution which does not belong to case (i) or case (ii) must enter $G_{4}$ after finite time and then it stays there. $\blacksquare$

Lemma 2 The stable manifold of $S_{1}$ intersects both axes. If a solution starts below the stable manifold of $S_{1}$ it converges to $S_{0}$ as $t\to\infty$. If it starts above the stable manifold it eventually lies in $G_{4}$.

Proof Consider the derivative of the right hand side of the system at $S_{1}$. This matrix has trace $-\alpha-10y-1<0$ and determinant $\alpha(1-2y)<0$. Thus it has one positive and one negative eigenvalue. Its stable manifold $V_{s}$ is one-dimensional and lies in $G_{2}\cup G_{3}$. Along this manifold $\frac{\dot{y}}{\dot{x}}$ is negative. It follows that $V_{s}$ is the graph of a function of $x$. As $x$ decreases along the part of $V_{s}$ to the left of $S_{1}$ the derivative of this function remains bounded and so $V_{s}$ intersects the $x$-axis. As $x$ increases along the part of $V_{s}$ to the right of $S_{1}$ the derivative of the function remains bounded away from zero and so $V_{s}$ intersects the $x$-axis. It can be concluded that the complement of $V_{s}$ in the positive orthant has two connected components $H_{1}$ and $H_{2}$, where $H_{1}$ has compact closure. A point is said to lie below the stable manifold if it belongs to $H_{1}$ and above the stable manifold if it belongs to $H_{2}$. A solution which starts in one of these two components remains in it. A solution which starts in $H_{1}$ cannot reach $G_{4}$ and one which starts in $H_{2}$ cannot reach $G_{1}$. Thus Lemma 2 follows from Lemma 1. $\blacksquare$

Lemma 3 A solution of (1)-(2) which is eventually contained in $G_{4}$ has, after a suitable translation of $t$, the asymptotics $x=\frac{\alpha}{5}e^{\frac{\alpha t}{5}}+\ldots$, $y=\frac{\sqrt{2}\alpha}{5}e^{\frac{\alpha t}{10}}+\ldots$ for $t\to\infty$.

Proof If a solution is eventually contained in $G_{4}$ then $r=\sqrt{x^{2}+y^{2}}$ must tend to infinity as $t\to\infty$. For $r$ is an increasing function of $t$ and if it were bounded the solution would have to converge to a steady state. However there are no steady states in $G_{4}$. It then follows from the defining equations for $G_{4}$ that both $x$ and $y$ tend to infinity as $t\to\infty$. We now consider the Poincaré compactification of the system [15]. Usually this compactification is constructed using two charts, covering neighbourhoods of the $x$- and $y$-axes respectively. For a solution which is in $G_{4}$ for $t\geq t_{1}$ we have seen that $y$ tends to infinity for $t\to\infty$ and hence $\frac{x}{y}$ tends to infinity. This means that the solution eventually leaves a neighbourhood of the origin in the chart covering a neighbourhood of the $y$-axis and lies in the chart covering a neighbourhood of the $x$-axis. Moreover it tends to the origin in the latter chart as $t\to\infty$. The chart we are talking about is defined by the coordinates $X=1/x$ and $Z=y/x$. Define a new time coordinate $\tau$ which satisfies $\frac{d\tau}{dt}=x$. The transformed system is

[TABLE]

Both eigenvalues of the linearization at the origin are zero and so we must blow up the origin to get more information. This will be done by means of a quasihomogeneous blow-up following [4]. In the notation used there the exponents calculated using the Newton polygon are $(2,1)$. There are two transformations to be done, corresponding to the two coordinates. The first of these is given by the correspondence $(X,Z)=(u^{2},uv)$. We have

[TABLE]

and hence

[TABLE]

Furthermore

[TABLE]

and hence

[TABLE]

If we now introduce a new time coordinate $s$ satisfying $\frac{ds}{d\tau}=u$ then the system becomes

[TABLE]

When $v=0$ the derivative of $v$ is positive. Thus no solution can have an $\omega$-limit point on the $u$-axis. Hence the solution must eventually be contained in the chart defined by the second transformation, which is given by $(X,Z)=(uv^{2},v)$. In this case

[TABLE]

Furthermore

[TABLE]

and hence

[TABLE]

If we now introduce a new time coordinate $s$ satisfying $\frac{ds}{d\tau}=v$ then the system becomes

[TABLE]

The $u$- and $v$-axes are invariant, the origin is a steady state, which is a hyperbolic saddle, and there is an additional steady state $(5/2\alpha,0)$. The latter steady state has one negative and one zero eigenvalue. If we transform any positive solution then in the blown-up Poincaré compactification it must tend to that point. To get more details we translate the steady state to the origin using a coordinate transformation. Let $w=u-\frac{5}{2\alpha}$ be a new coordinate. Then the equations become

[TABLE]

The first can be rewritten as

[TABLE]

We now apply centre manifold theory (cf. [15], Section 2.7). The centre manifold can be written in the form $w=\frac{\alpha+10}{8\alpha^{2}}v+r(v)$ with a remainder term $r$ which is at least quadratic. Consider the contributions to the right hand side of the evolution equation for $v$ which are quadratic in $v$. We get

[TABLE]

After translating $s$ if necessary we get $v=\frac{4}{s}+\ldots$. Substituting this into the defining equation for $s$ gives $\tau=\frac{1}{8}s^{2}$ and $v=\sqrt{\frac{2}{\tau}}+\ldots$. It follows that $X=\frac{5}{\alpha\tau}+\ldots$ and $Z=\sqrt{\frac{2}{\tau}}+\ldots$. Next we compute the transformation from $\tau$ to $t$. We have $\frac{dt}{d\tau}=\frac{5}{\alpha\tau}+\ldots$ and hence up to a translation of the time coordinate $t=\frac{5}{\alpha}\log\tau+\ldots$ and $\tau=e^{\frac{\alpha t}{5}}+\ldots$. When written in the original variables these relations give $x=\frac{\alpha}{5}e^{\frac{\alpha t}{5}}+\ldots$, $y=\frac{\sqrt{2}\alpha}{5}e^{\frac{\alpha t}{10}}+\ldots$. $\blacksquare$

Theorem 1 A positive solution of (1)-(2) with $\alpha>0$ and $\beta=0$ belongs to one of the following three classes.

(i) It starts below the stable manifold of $S_{1}$ and $x$ and $y$ converge to zero as $t\to\infty$.

(ii) It starts on the stable manifold of $S_{1}$ and converges to $S_{1}$ as $t\to\infty$.

(iii) It starts above the stable manifold of $S_{1}$ and $x$ and $y$ tend to infinity as $t\to\infty$, with the asymptotics given in Lemma 3.

In particular, every bounded solution converges to a steady state as $t\to\infty$.

Proof This is obtained by combining Lemma 1 - Lemma 3. $\blacksquare$

Note that because after transformation to the Poincaré compactification each solution converges to a steady state there exist no periodic solutions. In other words the system does not exhibit sustained oscillations. The eigenvalues of the linearization of the system about $S_{0}$ are real because the axes are invariant manifolds. The steady state $S_{1}$ has been shown to be hyperbolic with the eigenvalues of the linearization about that point being real. Thus damped oscillations decaying to one of the steady states are also ruled out.

3 The case with photorespiration

In this case we consider the full system where $\alpha$ and $\beta$ are both non-zero.

Lemma 4 Corresponding to positive initial data for (1)-(2) with $\alpha>0$ and $\beta>0$ given at $t=t_{0}$ there exists a solution on the interval $(t_{0},\infty)$ and it is bounded.

Proof Taking a suitable linear combination of the equations gives $\frac{d}{dt}(5x+3y)=-\beta x^{2}+\alpha x-3y$. If $x\geq\alpha/\beta$ then the right hand side is negative. If $x\leq\alpha/\beta$ then $\alpha x\leq\alpha^{2}/\beta$. Thus if also $y\geq\alpha^{2}/3\beta$ the right hand side is negative. If a solution satisfies $5x+3y>\beta^{-1}(5\alpha+\alpha^{2})$ at some time then it must be in one of the regions where the time derivative of $5x+3y$ is negative. Thus the value of $5x+3y$ is bounded by the maximum of its initial value and $\beta^{-1}(5\alpha+\alpha^{2})$. It follows that all solutions of this system can be extended to exist globally in the future and are bounded. $\blacksquare$

Consider now steady states of (1)-(2).

Lemma 5

(i) For $\alpha^{2}/\beta<32$ the only non-negative steady state is the origin, which we once again denote by $S_{0}$.

(ii) For $\alpha^{2}/\beta=32$ there is precisely one positive steady state, which we call $S_{1}$.

(iii) For $\alpha^{2}/\beta>32$ there are precisely two positive steady states. For one of these, which we call $S_{1}$, both coordinates are smaller than the corresponding coordinates of the other steady state, which we call $S_{2}$.

Proof Any steady state satisfies $\alpha x=y^{2}+2y$, so that its $y$-coordinate determines its $x$-coordinate. In fact the $y$-coordinate is a monotone function of the $x$-coordinate. At the same time $\beta x^{2}=y^{2}-y$. Squaring the first of these equations and substituting it into the second gives

[TABLE]

and hence

[TABLE]

The positive steady states are in one to one correspondence with the positive roots of the cubic $p(y)$. Since $p(0)>0$ there is at least one negative root and there are at most two positive roots. If $4\beta-\alpha^{2}>0$ then it follows from Descartes’ rule of signs [14] that there are no positive roots. Further information about the number of positive roots can be obtained by looking at the discriminant of the polynomial $p$. It is given by

[TABLE]

There are two positive values of $\alpha^{2}/\beta$ for which $\Delta$ vanishes, namely $\zeta_{\pm}=\frac{131\pm 125}{8}$. We have $\zeta_{-}=\frac{3}{4}$ and $\zeta_{+}=32$. When $\Delta<0$ the polynomial $p$ has only one real root and it must be negative. When $\Delta\geq 0$ all roots are real. Only in this case can there be more than one positive root. For $\alpha^{2}/\beta=32$ there is a root which is at least double. If this root were negative then it would follow that $4\beta-\alpha^{2}>0$, a contradiction. Thus the double root is positive. For $\alpha^{2}/\beta>32$ there are two positive roots. $\blacksquare$

Both of the nullclines of this system are of the form $f(x)=g(y)$ for monotone increasing functions $f$ and $g$. Thus we can write them as graphs of functions of $x$ or of functions of $y$. Due to Lemma 5 we know that the two nullclines intersect at the origin and at no, one or two points in the positive region, depending on the parameters. It can then be concluded that when there are no, one and two positive steady states the complement of the union of the nullclines has three, four and five connected components, respectively. As in the previous section these components can be labelled with the signs of $\dot{x}$ and $\dot{y}$. In case (i) of Lemma 5 there is one component with each of the labels $(-,-)$, $(-,+)$ and $(+,-)$. In case (ii) there are two components with the label $(-,-)$ and one component with each of the labels $(-,+)$ and $(+,-)$. In case (iii) (cf. Figure 2)

there are two components with the label $(-,-)$ and one component with each of the labels $(-,+)$, $(+,-)$ $(+,+)$. The components of the complements of the set of steady states in the nullclines can be labelled as $(-,0)$, $(0,-)$, $(+,0)$ and $(0,+)$.

At a point where the nullclines cross they are tangent if and only if the linearization at the point has a zero eigenvalue. This happens precisely when the determinant of the linearization at that point is zero. This is the only way that a steady state can fail to be hyperbolic since the fact that the trace is negative rules out the possibility of a pair of complex conjugate imaginary eigenvalues. With the information we already have about steady states it can be concluded that the determinant is zero precisely when

[TABLE]

Taking suitable linear combinations of the equations $p(y)=0$ and $q(y)=0$ leads to two quadratic equations for $y$. These can be combined to give a linear equation for $y$ and substituting this back into $r(y)$ leads to the equation

[TABLE]

where $\gamma=\frac{\alpha^{2}}{\beta}$. It follows that all steady states are hyperbolic when $\gamma>32$. It can be concluded that except in case (ii) the steady states are hyperbolic. In case (ii) the linearization at $S_{1}$ has one zero and one negative eigenvalue. In case (iii) it has one positive and one negative eigenvalue.

Theorem 2 Any positive solution of (1)-(2) with $\alpha>0$ and $\beta>0$ converges a steady state as $t\to\infty$. If $\frac{\alpha^{2}}{\beta}<32$ there are no points $S_{1}$ and $S_{2}$ and all solutions converge to $S_{0}$. If $\frac{\alpha^{2}}{\beta}=32$ there is no point $S_{2}$ and points above or on and below the unique centre manifold of $S_{1}$ converge to $S_{1}$ and $S_{0}$ respectively. If $\frac{\alpha^{2}}{\beta}>32$ then points above, on or below the stable manifold of $S_{1}$ converge to $S_{2}$, $S_{1}$ and $S_{0}$ respectively.

Proof Note that the stable manifold of $S_{1}$ is always one-dimensional. By the same arguments as in the case $\beta=0$ it can be shown that this manifold is the graph of a function of $x$ and that its complement is the union of two components $H_{1}$ and $H_{2}$. Components with the sign combination $(-,-)$ have boundaries with the sign combinations $(-,0)$ and $(0,-)$. Using the information on the signs of $\ddot{x}$ and $\ddot{y}$ shows that these components are invariant. For instance, a solution which satisfies $\dot{x}=0$ and $\dot{y}<0$ at some point satisfies $\ddot{x}<0$. Similarly components with the sign combination $(+,+)$ are invariant. It can be shown as in the case $\beta=0$ that any solution which starts in a component with one of sign combinations $(+,-)$ or $(-,+)$ and does not converge to a steady state as $t\to\infty$ must enter one of the components with the sign combination $(-,-)$ or $(+,+)$ after a finite time. Once it enters a component of this type it must stay there and converge to a steady state as $t\to\infty$. Thus every solution converges to a steady state as $t\to\infty$. It is then straightforward to determine which steady state it converges to in different cases. $\blacksquare$

Note that since every solution converges to a steady state the system exhibits no sustained oscillations. With the information we have about the linearization about the fixed points we can argue as in the previous section that there are no damped oscillations close to the point $S_{0}$. Using Lemma A2 of the Appendix we get the corresponding conclusion for $S_{1}$ and $S_{2}$.

Consider what happens if $\beta$ tends to zero while $\alpha$ has a fixed positive value. The polynomial converges. In the limit there is a unique root, which is $S_{1}$. It is a hyperbolic saddle. Thus it is the limit of a steady state of the system in the general case as $\beta\to 0$. The approximating solution must coincide with the point $S_{1}$ in the general system. Let $\gamma=\frac{\alpha^{2}}{\beta}$ and define $z=\gamma^{-1/2}y$ and $q(z)=p(y)$. Then if $\beta$ tends to zero while $\alpha$ has a fixed positive value the polynomial $q$ converges. Again there is a unique positive root in the limit with $z=1$. It is approximated by a root for positive $\beta>0$ and that corresponds to the steady state $S_{2}$. We conclude that as $\beta$ tends to zero the coordinates of $S_{2}$ have the asymptotic behaviour $x=\alpha\beta^{-1}+\ldots$ and $y=\alpha\beta^{-\frac{1}{2}}+\ldots$.

4 The system with the fifth power

In this section we study the system where $y^{2}$ is replaced in (1)-(2) by $y^{5}$. This is

[TABLE]

The aim is to see to what extent the results obtained for (1)-(2) generalize to (32)-(33).

In the case $\alpha=0$ the analysis of (1)-(2) extends without essential changes to (32)-(33) to give the same qualitative results. When $\alpha\neq 0$ and $\beta=0$ the analysis up to and including Lemma 2 extends easily. Of course the explicit formulae in the definitions of the invariant regions $G_{i}$ are modified by replacing $y^{2}$ by $y^{5}$.

Lemma 6 A solution of (32)-(33) which is eventually contained in $G_{4}$ has, after a suitable translation of $t$, the asymptotics $x=(\frac{4\alpha}{5})^{\frac{1}{4}}e^{\frac{\alpha t}{5}}+\ldots$, $y=2^{\frac{1}{20}}\left(\frac{2\alpha}{5}\right)^{\frac{1}{4}}e^{\frac{\alpha t}{25}}+\ldots$ for $t\to\infty$.

Proof That the arguments from the case of (1)-(2) extend easily is also true of the first part of the proof of Lemma 3 which shows that the late time behaviour can be analysed in one of the charts of the Poincaré compactification. In this case the time coordinate must be rescaled in a different way from what was done previously. Let $\frac{d\tau}{dt}=x^{4}$. In the case of (32)-(33) the transformation to this chart gives

[TABLE]

The linearization of the system at the origin is identically zero. To get more information we do a quasi-homogeneous blow-up. The exponents calculated using the Newton polygon are $(5,4)$. Once again, there are two transformations to be done. The first of these is given by the correspondence $(X,Z)=(u^{5},u^{4}v)$. We have

[TABLE]

and hence

[TABLE]

Furthermore

[TABLE]

and hence

[TABLE]

If we now introduce a new time coordinate $s$ satisfying $\frac{ds}{d\tau}=u^{16}$ then the system becomes

[TABLE]

Just as in the case with the quadratic nonlinearity we see that the solution must eventually be contained in the chart defined by the second transformation which is given by $(X,Z)=(uv^{5},v^{4})$. In that case

[TABLE]

Hence

[TABLE]

Furthermore

[TABLE]

and hence

[TABLE]

In terms of the time coordinate $s$ with $\frac{ds}{d\tau}=v^{16}$ we get

[TABLE]

The axes are invariant and the origin is a hyperbolic saddle. There is a steady state at the point $(u_{0},0)$ with $u_{0}=\left(\frac{5}{2\alpha}\right)^{\frac{1}{4}}$. If we transform any solution then in the blown-up Poincaré compactification it must converge to this point. To get more details we translate the steady state to the origin using a coordinate transformation. Let $w=u-u_{0}$. Then the equations become

[TABLE]

Here we have explicitly retained only those terms which are required for the calculation which will now be done. It follows from the definition of the centre manifold that $w=h(v)$ for a function $h$ with $h(v)=O(v^{2})$. The derivative of this relation with respect to time also holds. Hence $\dot{w}=h^{\prime}(v)\dot{v}$. It follows from (49) that $\dot{v}=O(v^{3})$ and so $\dot{w}=O(v^{4})$. It follows from (48) that $w=O(v^{4})$. Hence $\dot{v}=O(v^{5})$ and $\dot{w}=O(v^{6})$. It can be concluded from the evolution equation for $w$ that

[TABLE]

It follows that $\frac{dv}{ds}=-\frac{1}{10}v^{5}+\ldots$. We see that the flow on the centre manifold is towards the steady state. After translating $s$ if necessary we get $v=\left(\frac{5}{2s}\right)^{\frac{1}{4}}+\ldots$. Substituting this into the defining equation for $s$ gives $s=5^{\frac{1}{5}}\left(\frac{5}{2}\right)^{\frac{4}{5}}\tau^{\frac{1}{5}}+\ldots$ and $v=\left(\frac{1}{2\tau}\right)^{\frac{1}{20}}+\ldots$. Substituting for the original variables gives $X=u_{0}\left(\frac{1}{2\tau}\right)^{\frac{1}{4}}+\ldots$ and $Z=\left(\frac{1}{2\tau}\right)^{\frac{1}{5}}+\ldots$. Next we compute the transformation from $\tau$ to $t$. We have $\frac{dt}{d\tau}=\left(\frac{u_{0}^{4}}{2\tau}\right)+\ldots$. Hence $t=\left(\frac{u_{0}^{4}}{2}\right)\log\tau+\ldots$ and $\tau=e^{4\alpha t/5}+\ldots$. Finally we get $x=(\frac{4\alpha}{5})^{\frac{1}{4}}e^{\frac{\alpha t}{5}}$ and $y=2^{\frac{1}{20}}\left(\frac{2\alpha}{5}\right)^{\frac{1}{4}}e^{\frac{\alpha t}{25}}$. $\blacksquare$

Theorem 3 A positive solution of (32)-(33) with $\alpha>0$ and $\beta=0$ belongs to one of the following three classes.

(i) It starts below the stable manifold of $S_{1}$ and $x$ and $y$ converge to zero as $t\to\infty$.

(ii) It starts on the stable manifold of $S_{1}$ and converges to $S_{1}$ as $t\to\infty$.

(iii) It starts above the stable manifold of $S_{1}$ and $x$ and $y$ tend to infinity as $t\to\infty$, with the asymptotics given in Lemma 6.

In particular, every bounded solution converges to a steady state as $t\to\infty$.

Proof The proof is identical to that of Theorem 1 except that Lemma 3 is replaced by Lemma 6. $\blacksquare$

It is interesting to compare the asymptotics in Lemma 6 with those obtained in [19] for a more elaborate model of the Calvin cycle. In Lemma 6 we see that both unknowns have growing exponential asymptotics but that the exponent for GAP is one fifth of that for the other variables. The main system considered in [19] has five unknowns and has solutions for which all unknowns have growing exponential asymptotics. In that case the exponent for GAP is one fifth of that for the other four unknowns. These four unknowns satisfy a system of the form $\frac{d\bar{x}}{dt}=A\bar{x}+R$ where $R$ is considered as a remainder term and the larger exponent is an eigenvalue of $A$. There is a natural analogue of this equation for the system (32)-(33) with $\beta=0$. It is the equation $\frac{d}{dt}(5x+3y)=\frac{\alpha}{5}(5x+y)-\left(\frac{\alpha}{5}+1\right)y.$ Here the last term is to be considered as the remainder. Note that in the asymptotics of Lemma 6 $y$ is much smaller than $x$ at late times so that this treatment as remainder term is reasonable. Since there is ony one unknown growing at the maximal rate in this case the matrix $A$ is replaced by a number and that number is $\alpha/5$. Thus we see that on a heuristic level the exponents in the two cases agree. The statement of Theorem 3 is stronger than the analogous statement in [19] in the following sense. The description of the asymptotic behaviour in [19] is only obtained for some non-empty subset of initial data which is not further characterized while the set of initial data giving rise to this asymptotic behaviour in Theorem 3 is much more explicit.

Consider next the case where the coefficients $\alpha$ and $\beta$ in (32)-(33) are both positive. The solutions are bounded using the same argument as in the proof of Lemma 4. As in the case of (1)-(2) the nullclines are of the form $f(x)=g(y)$ for monotone increasing functions $f$ and $g$. A steady state satisfies the equations

[TABLE]

Substituting the second of these equations into the first gives

[TABLE]

By Descartes’ rule of signs this equation can have at most two positive solutions. Since the derivative of the polynomial $p$ at zero is positive $p(y)>0$ for $y$ slightly larger than zero. Thus if $p(y)<0$ for some $y>0$ the polynomial has two positive roots. Now $p(2)=1296\beta-28\alpha^{2}$. Thus for fixed $\beta$ if $\alpha$ is large enough we have $p(2)<0$ and $p$ has two positive roots. If we define values of $x$ corresponding to these two values of $y$ we obtain two positive steady states of the system (32)-(33). On the other hand, if $\beta>\alpha^{2}$ there are no positive steady states.

We have not succeeded in obtaining information about the hyperbolicity of steady states of this system which is as complete as that which we obtained in the case of a quadratic nonlinearity. It is, however, possible to show that for generic values of the parameter $\gamma$ all steady states are hyperbolic. We can calculate polynomials $p$ and $q$ as in the case with a quadratic nonlinearity but it is not possible to solve explicitly for their common roots $y$. Instead we can proceed as follows. For any non-hyperbolic steady state we obtain equations of the form

[TABLE]

for certain polynomials $p_{1}$ and $q_{1}$ which do not depend on $\gamma$. Hence

[TABLE]

Since the polynomial $s$ is non-constant this equation has only finitely many solutions $y$. For any given solution $y$ there is at most one corresponding value of $\gamma$. Hence for all but finitely many values of $\gamma$ all steady states are hyperbolic.

With the information on steady states just obtained we can prove an analogue of Theorem 2 for the system with the fifth power using the same techniques. The result is

Theorem 4 Any positive solution of (32)-(33) with $\alpha>0$ and $\beta>0$ converges to a steady state as $t\to\infty$. If $\frac{\alpha^{2}}{\beta}<1$ there are no points $S_{1}$ and $S_{2}$ and all solutions converge to $S_{0}$. If $\frac{\alpha^{2}}{\beta}$ is large enough then points above, on or below the stable manifold of $S_{1}$ converge to $S_{2}$, $S_{1}$ and $S_{0}$ respectively.

By scaling the unknowns $x$ and $y$ by the same factor and $t$ by another factor it is possible to transform the more general system

[TABLE]

for general positive constants $A$ and $B$ into the system (32)-(33). Thus the results obtained for (32)-(33) imply analogous results for (56)-(57). This observation will be used in the next section.

5 Derivation from the three-dimensional system

The system (1)-(2) was derived by Hahn from a three-dimensional system but he did not give a mathematical formulation of the relation between the two systems. The three-dimensional system is, in a modified notation,

[TABLE]

We now consider a limit where $k_{3}$ becomes large. This means that the reaction producing triose phosphate from PGA is very fast. Let $k_{3}=\epsilon^{-1}\tilde{k}_{3}$ and introduce a new variable by $w=y+z$. Then the equations above are equivalent to the system

[TABLE]

This is a fast-slow system in standard form with fast variable $z$ and slow variables $x$ and $w$. The critical manifold is given by $z=w$ and the slow system is

[TABLE]

Replacing $w$ by $y$ and setting $k_{1}=\alpha$, $k_{2}=\beta$, $k_{4}=A$ and $k_{5}+k_{6}=B$ gives the system (56)-(57). The critical manifold is normally hyperbolic and the one normal eigenvalue is negative.

We know that for certain values of the parameters the system (32)-(33) has three steady states $S_{0}$, $S_{1}$ and $S_{2}$. Moreover, for generic values of $\gamma$ the steady states $S_{0}$ and $S_{2}$ are hyperbolic sinks while $S_{1}$ is a hyperbolic saddle with a one-dimensional stable manifold. There are heteroclinic orbits connecting $S_{0}$ to $S_{1}$ and $S_{1}$ to $S_{2}$. Putting this together with the fact that the normal eigenvalue is negative shows that for suitable parameters with $\epsilon$ small the three-dimensional system has three steady states $S_{0}$, $S_{1}$ and $S_{2}$ which converge to those with the corresponding names as $\epsilon=0$. Moreover $S_{0}$ and $S_{2}$ are hyperbolic sinks and $S_{1}$ is a hyperbolic saddle with a two-dimensional unstable manifold.

6 Conclusions and outlook

In this paper we have obtained detailed information on minimal models of the Calvin cycle introduced by Hahn in [9]. A rather complete analysis of the two-dimensional models of Hahn was given. The relation of the two-dimensional to the three-dimensional model of Hahn was discussed but a comprehensive analysis of the three-dimensional model, which is likely to be complicated, was postponed to future work. The models in [9] originated by formal simplification of earlier models due to the same author. The first is a model with 19 chemical species defined in [7]. It did not implement a detailed description of photorespiration and a description of this kind was added in the model of [8], with 33 species. In the future it would be desirable to put the understanding of the relations between these different models on a better mathematical footing. This should also allow conclusions about the dynamics of the higher-dimensional systems to be obtained. Note that there are some general references in the literature about the inheritance of dynamical features from reduced systems (see e.g. [3], [5]).

Another interesting task is to relate the models of Hahn to other models of the Calvin cycle in the literature. Possible relations to a model of Grimbs et al. [6] studied in [19] were already mentioned is section 4 and perhaps these could be extended so as to give a wider view of runaway solutions of models for the Calvin cycle, i.e those solutions where all concentrations tend to infinity. One task is to obtain some kind of characterization of models admitting solutions of this type. Another is to obtain formulae for the asymptotics of these solutions in the case that they do occur. This kind of behaviour can be ruled out if the model admits a suitable conservation law. This is, for instance, the case in a model of [16] whose mathematical properties were studied in [13]. In the model of Hahn with photorespiration boundedness of solutions is obtained without there being a conservation law. Another example of this is given by a model studied in Section 6 of [19] where the original model of Grimbs at al. is modified by including the concentration of ATP explicitly.

It also remains to obtain a comprehensive understanding of solutions where some concentrations tend to zero at late times. As discussed in [13] this can be related to the biological phenomenon of overload collapse. This means intuitively that the production of sugar by the cycle cannot meet the demand for export from the chloroplast. In [13] it was shown that there are solutions of the model of [16] which admit this phenomenon while in a modification of the model due to Poolman [17] these solutions are eliminated. The model of [17] does not include photorespiration but does include the mobilization of glucose from starch.

In this paper various aspects of the dynamics of the minimal model of Hahn for the Calvin cycle have been analysed. The reduction to a model with only two variables makes it possible to get a good overview of the dynamics. We believe that this establishes a good starting point for understanding the dynamics of more detailed models of the Calvin cycle in the future and we have indicated some directions in which this could progress.

7 Appendix: nullcline analysis

In this section we discuss how nullclines can be used to obtain information about the global behaviour of solutions of a two-dimensional dynamical system. Consider the following system of ordinary differential equations.

[TABLE]

where the functions $f$ and $g$ are $C^{1}$ and defined on an open subset $U\subset{\mathbb{R}}^{2}$. The nullclines $N_{1}$ and $N_{2}$ are the zero sets of $f$ and $g$ respectively. Let $G=U\setminus(N_{1}\cup N_{2})$. The open set $G$ is a countable union of connected components $G_{i}$. In what follows we will restrict to the case that the following assumption is satisfied.

Assumption 1 The complement of the nullclines has only finitely many connected components.

When the system satisfies Assumption 1 it defines a directed graph as follows. There is one node for each component $G_{i}$ and there is a directed edge from the node corresponding to $G_{i}$ to that corresponding to $G_{j}$ when there is a solution which starts from a point of $G_{i}$ and later enters $G_{j}$ without entering any component of $G$ other than $G_{i}$ and $G_{j}$ at an intermediate time. Let us call this the succession graph. A cycle in a directed graph is a finite sequence of directed edges such that the initial node of each edge is the final node of the previous one and the final node of the last edge is the initial node of the first one. We now restrict to the case that the following assumption is satisfied.

Assumption 2 There exist only finitely many steady states. Whenever a steady state $S_{i}$ is in the closure of a component $G_{j}$ there is a continuous curve joining a point of $G_{j}$ to $S_{i}$ which does not intersect any other $G_{k}$.

Lemma A1 Consider a solution $(x(t),y(t))$ on a time interval $[t_{0},t_{1})$ which lies in $G_{i}$ for some $i$ when $t=t_{0}$ and which lies entirely in $\bar{G}_{i}$. Then $x(t)$ and $y(t)$ are monotone. They are strictly monotone as long as the solution lies in $G_{i}$. Suppose that $t_{1}$ is maximal. If the solution is bounded then it converges to a point $(x^{*},y^{*})$ for $t\to t_{1}$ which is either a point of $N_{1}\cup N_{2}$ or a point of $\bar{G}\setminus G$. If $(x^{*},y^{*})\in G$ then $(x^{*},y^{*})\in N_{1}\cap N_{2}$ if and only if $t_{1}=\infty$.

Proof On a component $G_{i}$ the signs of $\dot{x}$ and $\dot{y}$ are constant and this implies the monotonicity statements. It follows that the limits of $x(t)$ and $y(t)$ as $t\to\infty$ exist, either as real numbers or as $\pm\infty$. If the solution is bounded then these limits are real numbers $x^{*}$ and $y^{*}$. The point $(x^{*},y^{*})$ belongs to the closure of $G$. Suppose now that $t_{1}$ is maximal and that $(x^{*},y^{*})\in G$. We claim that if $t_{1}=\infty$ then $(x^{*},y^{*})$ is a point of $N_{1}\cap N_{2}$, and hence a steady state. Otherwise at least one of $\dot{x}$ or $\dot{y}$ would tend to a non-zero value, say $c$, as $t\to t_{1}$. Suppose w.l.o.g. that $\dot{x}$ has this property and that $c>0$. It follows that if $t_{2}$ is sufficiently large then $\dot{x}(t)\geq\frac{1}{2}ct$ for all $t\geq t_{2}$. Thus $x$ is unbounded, a contradiction. We conclude that if the interval $[t_{0},t_{1})$ is infinite $(x^{*},y^{*})\in N_{1}\cap N_{2}$. Suppose now conversely that $(x^{*},y^{*})\in N_{1}\cap N_{2}$. If $t_{1}$ were finite it would be possible to extend the solution beyond $t=t_{1}$. But then this solution would have to coincide with the time-independent solution $x(t)=x^{*}$, $y(t)=y^{*}$ a contradiction. Thus if $(x^{*},y^{*})\in N_{1}\cap N_{2}$ the interval $[t_{0},t_{1})$ is infinite.

Lemma A2 Let $S_{i}$ be a steady state. Suppose that Assumption 2 is satisfied, that there is more than one component $G_{j}$ having $S_{i}$ as a limit point and that there is no cycle in the succession graph. Then there is no damped oscillation converging to $S_{i}$.

Proof Suppose there is a solution exhibiting a damped oscillation. Due to Assumption 2 it must intersect each $G_{j}$ having $S_{i}$ as a limit point more than once. Hence the succession graph contains a cycle, a contradiction.