A generalized stoichiometric model of C3, C2, C2+C4, and C4 photosynthetic metabolism

A model based solely on mass-balance constraints refines our understanding of the trade-offs, energy requirements, leaf-level fluxes, and plasticity mechanisms in different biochemical types of assimilation.


Introduction
Interest in engineering a biochemical carbon-concentrating mechanism (CCM, abbreviations listed in Table 1) to suppress photorespiration in crops is increasing (Furbank et al., 2015;Long et al., 2015). The metabolic activities of a CCM are shared between mesophyll (M) and bundle sheath (BS) cells. Structurally, the leaf parenchyma is organized in concentric layers of cells, with an outer mesophyll encircling one or two layers of BS cells. In some species, the BS cells are isolated from the surroundings by a gas-tight suberized cell wall (Lundgren et al., 2014). Biochemically, the compartmentalization of glycine decarboxylase (GDC) activity in the BS allows plants to take advantage of photorespiratory CO 2 release (Sage et al., 2012;Mallmann et al., 2014), giving rise to a mechanism that delivers CO 2 around the Rubisco in the BS -the so-called C 2 shuttle (Schulze et al., 2013;Keerberg et al., 2014). The 'C 4 cycle' is a further adaptation involving an energy-dependent carboxylation-decarboxylation cycle. CO 2 is initially fixed into four-carbon (C 4 ) organic (amino) acids in the M by phosphoenolpyruvate (PEP) carboxylase (PEPC). These then diffuse to the BS where they are decarboxylated. C 4 plants have traditionally been grouped into three subtypes (Table 2) according to the main decarboxylating acid (PGA). Further, pyruvate phosphate dikinase (PPDK) has traditionally thought to be confined to the M; however, PPDK was shown to be present and active in the BS (Aoyagi and Nakamoto, 1985;Majeran et al., 2010), although the role of PPDK in the BS remains elusive.
Detailed quantification of ATP and NADPH supply and demand in the BS and M is critical for understanding the physiology and regulation of photosynthesis. In terms of supply, the partitioning of ATP and NADPH production varies between cell types depending on the light available locally in the BS or M (Bellasio and Griffiths, 2014c). The dependence of ATP generation upon anatomical traits in the evolutionary continuum from C 3 to C 4 has recently been studied (Bellasio and Lundgren, 2016), and will not be addressed here. In this work, I shall concentrate on ATP and NADPH demand.
Quantifying ATP and NADPH demand in the M and BS requires detailed mechanistic understanding of assimilatory biochemistry. Mathematical modelling offers valid support for integrating knowledge at the systems level (Morandini, 2013;Singh et al., 2014). Classical photosynthetic models have allowed the simulation of leaf-level assimilation in C 3 , C 2 , C 2 +C 4 , and C 4 plants using a mechanistic description based on Rubisco PEPC kinetics (von Caemmerer, 1989(von Caemmerer, , 2000(von Caemmerer, , 2013. These models are based on several simplifications, which limit their applicability. First, they do not account for spatial segregation of biochemical processes, offering only limited support when the separate requirements of the BS and M are being studied. Second, classical model(s) do not distinguish between biochemical subtypes, making it difficult to evaluate the particular requirements of each subtype. Finally, the models were primarily developed to predict leaf-level CO 2 exchange, while the stoichiometry, energetics, and fluxes between the BS and M, which represent a critical bottleneck for C 4 photosynthesis (Pick et al., 2011), are not treated sufficiently, and a dedicated model is consequently needed.
The aim of this work was to develop a stoichiometric model of assimilation (hereafter SMA) which (i) generalizes all pathways of assimilation (C 3 , C 2 , C 2 +C 4 , and C 4 including C 4 photosynthetic subtypes); (ii) is based only on stoichiometry and therefore does not rely on kinetic measurements; and (iii) includes all main photosynthetic reactions, but is user-friendly for non-specialists. The theoretical underpinnings of the SMA are described in detail and a range of simulations to exemplify the model rationale and behaviour are provided. Numerous topics are covered, including (i) an analysis of C 3 and C 2 photosynthesis and all subtypes of C 2 +C 4 and C 4 (NADP-ME, NAD-ME, and PEPCK, in different combinations) in operational conditions; (ii) the energetics involved in manipulating photorespiration in a C 3 plant; (iii) the consequences of progressively upregulating a C 2 shuttle in a background of C 3 photosynthesis; and (iv) the consequences of progressively upregulating a C 4 cycle in a background of C 2 photosynthesis. Results quantify ATP and NADPH demand, which link dark and light reactions; refine our understanding of the stoichiometry of photosynthesis and the trade-offs inherent to different photosynthetic subtypes; and represent a useful framework for the integration of existing biochemical models.

SMA development
The SMA was developed on the basis of a stoichiometric model of NADP-ME C 4 photosynthesis (Bellasio and Griffiths, 2014c;McQualter et al., 2016) to augment all the pathways of carbon assimilation in a single tool. The SMA calculates key reaction rates, and ATP and NADPH requirements, in the M and BS as well as fluxes between the BS and M when the following parameters are known: the locality of GDC and Rubisco, leaf-level Rubisco rates of carboxylation and oxygenation (V O , V C ), and PEP carboxylation rate (V P ). The ATP and NADPH requirements are SMA outputs and they are not related to light reactions at this stage. Reactions are typically grouped by the biochemical function of the pathways, of which only the entry point is calculated. The SMA is based only on well-established reaction stoichiometry. The SMA accounts for the interactions between C 2 and C 4 cycles, including the fluxes associated with amino group rebalancing, and for the NADPH and ATP demand of assimilatory processes. The theory underpinning the SMA can be followed in Fig. 1. Owing to space limitations, the full description of the SMA is reported in Supplementary file 1.

Parameterization
The SMA has 12 input quantities (Table 3): three define Rubisco activity and assimilation (net assimilation, A; respiration in the light, R LIGHT ; and Rubisco rate of oxygenation relative to carboxylation, r O/C ), two define the activity of the CCM (PEPC engagement, as V P ; and PEPCK activity relative to V P , r PEPCK ), and seven partition key processes between the BS and M (f C , for Rubisco carboxylation rate; f O , for Rubisco oxygenation rate; f GDC , for glycine decarboxylase; f RLIGHT , for respiration in the light; f PR , for PGA reduction; f CS , for carbohydrate synthesis; f PPDK , for pyruvate phosphate dikinase). In addition, for NAD-ME subtypes it is possible to constrain malate dehydrogenase (MDH) activity to zero. Inputs can be constrained in different ways depending on the research questions. When the goal is to refine the analysis of a particular metabolic state of the leaf, input quantities may represent realistic biochemistry, otherwise inputs can be freely manipulated to simulate bioengineering or explore hypothetical scenarios. A, R LIGHT , and r O C / can be measured with accuracy (Bellasio et al., 2014;Bellasio et al., 2016a, b). Quantification of V P is less straightforward, and can be achieved through an in vitro assay of PEPC activity (Pfeffer and Peisker, 1998). This is complicated by PEPC-sensitive regulation and feedback inhibition, and the consequent necessity of reproducing physiological metabolite and ion concentrations in the reaction mixture. In alternative A, r O C / and V P can be predicted with biochemical models. Two types of formulations exist, either based on enzyme kinetics (and referred to as enzyme-limited), or based on the rate of ATP and NADPH made available by light reactions (and referred to as light-limited). Integrating enzyme-limited formulations in the SMA is straightforward, because the A, r O C / , and V P output by the biochemical model can be directly input into the SMA (see references below). Integrating light-limited formulations is more complicated and will be addressed in a dedicated paper. Physiological values are available for r PEPCK , gained through extensive biochemical work (Kanai and Edwards, 1999;Koteyeva et al., 2015). The traits underpinning f O , f C , f GDC , and f RLIGHT may require evolution or long acclimation periods to change (Sage et al., 2012;Christin and Osborne, 2014) and are therefore considered constant during gas exchange experiments. In the SMA, f O , f C , and f GDC represent the fraction of enzyme activity in the BS, rather than the physical distribution of the enzyme, but when enzyme compartmentalization is complete (     Aoyagi and Nakamoto (1985), Majeran et al. (2005), and Keerberg et al. (2014)] and are generally known for model species (Edwards and Ku, 1987). Cases of intermediate GDC distribution are rare (Sage et al., 2014), and, even in these cases, the enzyme distribution may be confidently taken as f GDC because substrate concentrations and regulation may be similar in the M and BS. For Rubisco, when compartmentalization is incomplete (as in C 2 +C 4 species), predicting f C and f O requires that the increased CO 2 concentration in the BS be modelled, for instance using the validated models for C 4 , C 3 , C 2 , and C 2 +C 4 photosynthetic subtypes (Farquhar et al., 1980;von Caemmerer, 1989von Caemmerer, , 2000von Caemmerer, , 2013von Caemmerer and Furbank, 1999). f RLIGHT is often assumed to be 0.5 in C 4 plants and 0.2 in C 2 +C 4 plants (von Caemmerer, 1989(von Caemmerer, , 2000, and may be determined from the relative BS/M mitochondrial abundance, or simply from the relative BS/M volume (Edwards and Ku, 1987). Circumstantial evidence suggests that f CS , f PR , r PEPCK , and f PPDK may rapidly change and this has likely contributed to the difficulty of distinguishing them experimentally. Exact empirical parameterization may therefore be of limited interest, and it may be more informative to define physiological maxima and minima (Bellasio and Griffiths, 2014c). Within these limits, the adjustment of f CS , f PR , r PEPCK and f PPDK may help the plant maximize assimilation under transient environmental inputs, such as changes in light quality, which may unbalance the partitioning of ATP and NADPH supply (Bellasio and Griffiths, 2014c). The effect of varying these inputs will be simulated in the subsection below dedicated to C 4 photosynthesis to show how these mechanisms can re-balance ATP and NADPH demand.

SMA simulations
The following simulations were selected to illustrate the capabilities, rationale, and behaviour of the SMA, while at the same time making some considerations of interest for the theme of this special issue. Simulations are grouped for photosynthetic types (C 3 , proto-Kranz and C 2 , C 2 +C 4 , and C 4 ). For each type, two sets of simulations are presented: static scenarios, where the SMA is calculated for one set of inputs, representing realistic operational conditions (Supplementary Table S1), and dynamic simulations where inputs are varied (Table 3).

Simulation 1. C 3 photosynthesis Simulation 1.1. Operational conditions
Outputs for C 3 photosynthesis in operational conditions are shown in Supplementary Fig. S1. The C 4 CCM is not operational. Carbon fixation and the reductive pentose phosphate (RPP) cycle operate in the M, while no carbon fixation is occurring in the BS. Although photorespiration and GDC activity are high, no glycine (GLY) is exported in the BS, and there is no net metabolite flux at the M-BS interface. All ATP and NADPH demand is in the M.

Simulation 1.2. Dynamic simulations
A dynamic scenario for a C 3 photosynthetic type was simulated by varying r O/C between 0 and 1 while keeping other quantities at 0 ( Table 3). The resultant ATP and NADPH demand as a function of r O/C is plotted in Supplementary Fig.  S2A while the ratio between ATP and NADPH demand is plotted in Supplementary Fig. S2B.

Simulation 2. C 2 and proto-Kranz Simulation 2.1. Operational conditions
Outputs for C 2 photosynthesis are shown in Supplementary  Fig. S3. To operate a C 2 shuttle, GDC is absent in the M and the photorespiratory GLY produced in the M is decarboxylated in the BS and fixed by a small fraction of Rubisco located in the BS. Serine (SER) diffuses back to the M, exporting half of the GLY amino groups. The excess ammonia produced by GDC in the BS is fixed by glutamine synthetase/glutamine oxoglutarate aminotransferase (GS/GOGAT), transaminated to alanine (ALA), and diffuses back to the M so as to re-balance the amino groups. For the C 2 shuttle to operate, an import flux of pyruvate (PYR) equimolar to R is required in the M. The C 4 CCM is not operational. In this example, PGA is reduced mainly in the M and a substantial BS↔M triose phosphate exchange occurs. Because the CCM is not operational, there is no net export of NADPH from M to BS and the NADPH demand in BS of 1.1 μmol m −2 s −1 must be met by linear electron flow (LEF); however, cyclic electron flow (CEF) may be preponderant as the ratio of ATP BS to NADPH BS was ~8.

Simulation 2.2. Dynamic simulations
f GDC was varied from 0 to 1 in an idealized C 3 plant, thus simulating the transition to C 2 photosynthesis through intermediate states of GDC compartmentalization, which are generally referred to as proto-Kranz (Sage et al., 2014). At each f GDC level, f C was fitted so that CO 2 leakage was zero (Supplementary Equation S34). The resultant f C as a function of f GDC is plotted in Fig. 2A. Values of f C above the curve will result in a net influx of CO 2 into the BS driven by Rubisco fixation, and in a CO 2 concentration in the BS lower than that in the M. Values of f C below the curve will result in a net CO 2 efflux out of the BS, and in a CO 2 -concentrating effect of the C 2 shuttle. Moving Rubisco and the photosynthetic carbon oxygenation (PCO) cycle to the BS results in changing the locality of ATP and NADPH demand ( Fig. 2A) and requires a substantial traffic of metabolites (Fig. 2B). The flux of PGA out of the BS and the opposite flux of dihydroxyacetone phosphate (DHAP; which is lower than that of PGA by a value equal to carbohydrate synthesis, CS) results from setting f PR = 0 in this simulation. This constraint also determines an excess of reducing power in the BS because the NADH produced by GDC in the BS is not used by PR. In these conditions the SMA predicts MAL to diffuse from the BS to M to shuttle the excess reducing power in the BS. MAL is oxidized to oxaloacetate (OAA) in the M and transaminated to ASP, which diffuses back to the BS and is transaminated back to OAA to supply MDH in the BS. The pair PYR/ALA balances the amino groups resulting from OAA/ASP transamination, and the amino groups resulting from the flux of GLY and SER are directly dependent on the operation of the C 2 shuttle. Alternative scenarios may involve NADH resulting from GDC activity sustaining a minimal level of PR in the BS (shown in Supplementary  Fig. S3).
Simulation 3. C 2 +C 4 C 2 +C 4 photosynthesis has complete GDC compartmentalization to the BS and intermediate states of CCM engagement (Table 2), corresponding to moderate V P and incomplete Rubisco compartmentalization. The CCM subtype is defined by the engagement of MDH in the M and PEPCK in the BS.

Simulation 3.1. Operational conditions
An example of NADP-ME C 2 +C 4 is shown in Supplementary  Fig. S4A. The C 4 CCM is operational and CCM activity is sufficient to exceed Rubisco carbon fixation in the BS, and there is a net CO 2 efflux from BS (L > 0). PR is mainly located in the M, and triose phosphate trafficking is higher than for the C 2 photosynthetic type. Because a large part of Rubisco is mainly located with GDC in the BS, and photorespiration is reduced by the activity of the CCM, the effectiveness of the C 2 shuttle (as R) is reduced relative to C 2 photosynthesis. PR in the BS consumes all NADPH available through MDH M , therefore T ≈ 0; however, there is no residual NADPH deficit, and LEF is not required in the BS. In a NAD-ME C 2 +C 4 photosynthetic subtype ( Supplementary Fig. S4B), the C 4 CCM results in the same export of CO 2 to the BS as in the NADP-ME C 2 +C 4 photosynthetic subtype; however, MDH activity in the M is zero, all OAA is transaminated, and the CCM does not export reducing power to the BS. PR and glycolate recycling in the BS consumes NADPH at a rate of 2.32 μmol m −2 s −1 , which must be generated through LEF.
In a PEPCK C 2 +C 4 photosynthetic subtype ( Supplementary  Fig. S4C), the C 4 CCM results in the same export of CO 2 to the BS as in the NADP-ME and NAD-ME C 2 +C 4 photosynthetic subtypes. However, although MDH is present in the M, PEPCK activity in the BS diverts OAA produced by PEPC in the M, making OAA unavailable for MDH activity in the M, which is consequently zero. Furthermore, because MDH activity is zero, the CCM does not export reducing power to the BS, and the NADPH must be supplied through LEF in the BS.

Simulation 3.2. Dynamic simulations
Intermediate states of 'C 4 ness', represented by intermediate degrees of f C and V P , were explored to simulate the transition from C 2 photosynthesis to C 4 photosynthesis. f C was incrementally varied, and, at each f C level, V P was iteratively fitted so that the CO 2 leak rate (L) remained zero (Supplementary Equation S34). The resultant V P as a function of f C is plotted in Fig. 3A. Values of V P above the curve will result in an effective CCM. Values of V P below the curve are insufficient to sustain Rubisco fixation, which will drive a net influx of CO 2 into the BS. For f C < 0.2, Rubisco fixation in the BS is supplied by the C 2 shuttle, meaning the predicted V P is zero. Moving Rubisco to the BS causes the locality of ATP and NADPH demand to change (Fig. 3A). Reducing power demand in the BS is low and is only used by the PCO cycle by setting f PR at zero. This requires PGA to diffuse out of the BS and in an opposite flux of DHAP, which is lower than that of PGA by a value equal to CS. In addition, ALA and ASP are used by the CCM to bypass MDH in the M. The concave trend of ALA and ASP fluxes reflects a decreasing recruitment for the C 2 shuttle at low f C and an increasing recruitment for the CCM at high f C . GLY and SER are recruited only by the C 2 shuttle and decrease to zero with f C . MAL and PYR are used by the C 4 cycle and their fluxes increase linearly with f C .

Simulation 4. C 4 photosynthetic subtypes
Here, Rubisco and GDC are completely compartmentalized to the BS and V P exceeds V C .  Supplementary Fig. S5A. The activity of the CCM is strong enough to exceed Rubisco carbon fixation in the BS and a positive CO 2 leakage out of the BS. PR is mainly located in the M and metabolite traffic is up to five-times that of gross assimilation (GA). The net effect of the C 2 shuttle (as R), which depends on Rubisco oxygenating activity in the M, is zero. The CCM supplies all NADPH BS (6.54 μmol m −2 s −1 ), no LEF is required in the BS, and there is excess OAA that is not reacted with by MDH and is subsequently transaminated (T).  Table 3.
Simulation 4.1.2. Supplementary Fig. S5B shows a typical NADP-ME subtype with engagement of PEPCK. PEPCK activity regenerates part of the PEP required by PEPC, driving a positive PEP flux out of the BS, which reduces the activity of PPDK in the M. PEPCK consumes half the ATP of PPDK, resulting in a 4% lower ATP/GA than for the NADP-ME subtype. However, PEPCK activity is generally low, all NADPH BS is supplied by the CCM, and no LEF is required in the BS. Simulation 4.1.3. Supplementary Fig. S5C shows a typical NAD-ME C 4 photosynthetic subtype with no engagement of PEPCK. The CCM exports CO 2 at the same rate as the NADP-ME subtype; however, MDH activity in M is zero, all OAA is transaminated, and the NADPH demand in the BS has to be generated in the BS through LEF. Simulation 4.1.4. Supplementary Fig. S5D shows a background NADP-ME metabolism using PEPCK as a sole decarboxylase (r PEPCK = 1). The CCM exports CO 2 at the same rate as the other subtypes. Because MDH is present in the M, the CCM would export reducing power to the BS if r PEPCK < 1, but, here, PEPCK activity in the BS diverts OAA produced by PEPC in the M, making OAA unavailable for MDH activity in the M, which is consequently zero. Furthermore, because MDH activity is zero, the CCM does not export reducing power to the BS, and the NADPH must be supplied through LEF in the BS. Simulation 4.1.5. Supplementary Fig. S5E shows SMA output for a background NAP-ME metabolism using PEPCK as a sole decarboxylase (r PEPCK = 1). The CCM exports CO 2 at the same rate as the other subtypes. Here, MDH in M is not operational, and the CCM does not have the potential to export reducing power to the BS.

Simulation 4.2. Dynamic simulations (C 4 plasticity mechanisms)
For the next four simulations, r O/C was set at a typical C 4 value of 0.05 (Bellasio et al., 2014); V P was set at 10.85 estimated after von Caemmerer (2000); and gradual transitions of f CS , f PR , r PEPCK , and f PPDK were simulated by calculating the SMA for 21 discrete values of f CS , f PR , r PEPCK , and f PPDK between 0 and 1. Other model inputs are listed in Table 3. An increase in f PR caused an increase in the ATP and NADPH demand in the BS, which occurred in both NADP-ME and NAD-ME subtypes (Fig. 4A, B). The NADPH demand to be supplied by LEF in the BS followed different trends. In the NADP-ME subtype, with an f PR of up to 0.5, the NADPH demand for PR was met by the CCM through the MAL shuttle, and the resultant NADPH demand through LEF was zero. Additional levels of PR (f PR > 0.5) required the engagement of LEF in the BS. In the NAD-ME subtype, the demand for LEF in the BS increased linearly for f PR > 0. Simulation 4.2.2. Manipulating PEPCK activity. The engagement of PEPCK was manipulated through the input parameter r PEPCK . f PR was set at 0.25 to highlight differences between NADP-ME and NAD-ME subtypes, while f PPDK and f CS were kept at zero (Table 3). An increase in r PEPCK increased the ATP demand in the BS in both the NADP-ME and NAD-ME subtypes (Fig. 4C, D), but the partitioning of NADPH demand was unaffected. The NADPH demand to be supplied by LEF in the BS, however, followed different trends. In the NADP-ME subtype, with an r PEPCK of up to 0.6, the NADPH demand through LEF was zero, and increased linearly up to 4.5 μmol m −2 s −1 for r PEPCK > 0.6. In the NAD-ME subtype, the NADPH demand to be supplied by LEF was constant at 4.5 μmol m −2 s −1 , regardless of r PEPCK . Simulation 4.2.3. Partitioning CS. In the simulation, the BS fraction of CS was manipulated to increase through the input parameter f CS , while r PEPCK , f PPDK , and f PR were kept at zero. Increasing f CS determined a marginal increase in ATP BS /ATP M from 0.26 to 0.29 in both the NADP-ME and NAD-ME subtypes, while the NADPH demand was unaffected (Fig. 4E, F). Simulation 4.2.4. Effect of PPDK engagement in the BS.The effect of PPDK engagement in the BS was studied by manipulating the input parameter f PPDK , while r PEPCK , f CS , and f PR were kept at zero. The ATP demand in the BS increased substantially, while the NADPH demand was unaffected in both the NADP-ME and NAD-ME subtypes (Fig. 4G, H).  Table 3.

Discussion
Flux-balance analysis models are constraint-based models in which steady state fluxes are predicted in a metabolic network by applying mass-balance constraints based on reaction stoichiometry (Sweetlove and Ratcliffe, 2011). Lately the complexity of models has grown to embrace a suite of photosynthetic processes (Laisk et al., 2009;Wang et al., 2014a;Wang et al., 2014b) and reconstruct genome-wide metabolism (Dal'Molin et al., 2010;Saha et al., 2011). Given this complexity, these models are only available for a few well-studied species and modifying them requires considerable coding effort, meaning they are not ideal for studying bio-manipulation or testing hypotheses. The metabolic model developed by Bellasio and Griffiths (2014c) facilitated straightforward modifications and changes to metabolic pathways. For instance, the model computed the effects of partitioning biochemical work between the BS and M on the locality of ATP and NADPH demand in photosynthesizing C 4 leaves, and was modified to support the interpretation of biochemical, gas exchang,e and transcriptomic data in engineered sugarcane (McQualter et al., 2016). Based on these theoretical underpinnings, the SMA was developed to account for the interactions between C 2 and C 4 cycles, including the fluxes associated with amino group rebalancing, recently described in Mallmann et al. (2014), and their effect on NADPH and ATP availability. The SMA integrates assimilatory metabolism and energetics as well as calculating key reaction rates, metabolite traffic fluxes, and ATP and NADPH requirements in the  Table 3.
M and BS when the locality of GDC (as f GDC ) and Rubisco (as f O , and f C ) activity, leaf-level Rubisco rates of carboxylation and oxygenation (as r O/C ), and PEP engagement (as V P ) are known. While previously published models are tailored for a particular species or a discrete photosynthetic type that is fixed a priori, the SMA allows rates to vary continuously between boundaries defined by what is biochemically feasible. The SMA is designed for integration with existing biochemical models. However, these models rely on assumptions that are specific to a particular photosynthetic type and, as such, they were not included in the SMA at this stage. The SMA is based upon logical constraints and well-established reaction stoichiometry, and is, consequently, of general use. Compared with other recent flux-balance models [e.g. Dal'Molin et al. (2010)], the SMA has several distinctive features. It is focused on assimilation; it generalizes all pathways of carbon assimilation in a single tool, except CAM, which requires an explicit temporal dynamic (Owen and Griffiths, 2013); it allows straightforward modification by the user; and it explicitly accounts for NADPH and ATP demand, allowing the study of the effect of biochemical regulation on energetics.
The following discussion refers to the simulations described in simulations 4.2.1-4.2.4 (C 4 plasticity mechanisms, Fig. 4), which are of particular interest for this special issue. Possible mechanisms that C 4 plants can exploit to stabilize their energy and redox state were recently reviewed (Stitt and Zhu, 2014). Reversible reactions linking PGA and triose phosphates, the reversibility of the PPDK reaction, the inter-conversion of PEP and PGA, and the transient build-up and utilization of PGA and triose phosphate require explicit temporal dynamics to be modelled and cannot be addressed by the SMA. Nevertheless, the futile cycles involved probably lower the biochemical efficiency of assimilation and are likely to be downregulated under steady state conditions. The other processes reviewed by Stitt and Zhu (2014) are shown in the simulation results. Specifically, the rate of ASP/MAL decarboxylation was mechanistically linked to reducing power requirement in the BS through Supplementary Equation S21, and can be followed in Supplementary file 2 (workbook cell O14); CO 2 leakage and leakiness are calculated by balancing all CO 2 fluxes in and out of the BS (Supplementary Equation S34) and can be followed in Supplementary file 2 (workbook cells O19 and O6, respectively); the exchange of triose phosphate and PGA and the effect of switching decarboxylase (ME versus PEPCK) are covered in simulations 4.2.1 and 4.2.2. In addition to the list of Stitt and Zhu (2014), here I investigated a possible role of CS by manipulating f PR (simulation 4.2.3), and I suggest a possible role for PPDK in the BS, which was investigated by manipulating f PPDK (simulation 4.2.4). The simulations were carried out as a single factor analysis in which only one input was varied at a time (Table 3). This approach has the benefit of allowing a comparison of the effectiveness of different plasticity mechanisms. Combining effects, however, is also straightforward: when multiple inputs are manipulated, the total effect is the sum of single effects. The four simulations are discussed separately.

Partitioning PR and regulation of NADPH demand
When BS engagement in PR was manipulated to increase, ATP and NADPH demand in the BS also increased. In particular, in NADP-ME subtypes above a f PR of 0.5, there was a substantial demand for NADPH through LEF in the BS. While ATP is produced by BS photophosphorylation in both NADP-ME and NAD-ME plants, BS chloroplasts in NADP-ME plants are generally not able to oxidize water and produce substantial amounts of NADPH. In NADP-ME plants, high levels of f PR are therefore unlikely, consistent with the observation that, in maize, PGA is mainly reduced in the M . Confining O 2 production to M chloroplasts, away from Rubisco carboxylating sites, contributes to photorespiration suppression (Majeran and van Wijk, 2009), but imposes functional trade-offs. The inability to sustain high levels of PR in the BS requires PGA to diffuse to M cells, and DHAP to return to BS cells in order to complete the RPP cycle (Fig. 1). This imposes substantial traffic across the M-BS interface that is five-times that of GA ( Supplementary Fig. S5A-E). For this traffic, high BS conductance would be beneficial, but imposes high levels of CO 2 retrodiffusion, and limits the effectiveness of the CCM (Bellasio and Griffiths, 2014a). Because of this trade-off, BS conductance is thought to be tightly regulated (Kromdijk et al., 2014), and was observed to scale with assimilation to optimize the efficiency of the CCM (Ubierna et al., 2013;Bellasio and Griffiths, 2014b;von Caemmerer and Furbank, 2016). The available evidence, however, is indirect and further investigation is required.

Manipulating PEPCK activity
While ATP is directly used in PEPCK reactions, PEPCK activity does not directly require NADPH, as shown by the plot of NADPH demand in Fig. 4C. The effect on NADPH demand through LEF is therefore indirect and depends on the subtype (NADP-ME or NAD-ME). NADP-ME subtypes have the potential to export reducing power to the BS through MDH activity in the M: when PEPCK activity increases, increasing levels of OAA are required, which are not available for MDH in the M and cannot be used to shuttle reducing power from the M to the BS through MAL. In other words, the NADPH potentially available through the CCM decreases. Ultimately, when PEPCK is fully engaged, there is no surplus of OAA available for exporting NAPDH to the BS and all NADPH demand in the BS must be met by LEF. In NAD-ME subtypes, which do not have the capacity to export reducing power, all NADPH demand has to be met by LEF. When PEPCK is the only decarboxylating enzyme, irrespective of the presence of MDH in the M (that is, in both NAD-ME and NADP-ME subtypes), no reducing power can be exported to the BS (Fig. 2C, D). Wang et al. (2014a) noted considerable variability in PEPCK engagement. For instance, in the NADP-ME subtypes there is a gradient from Sorghum bicolor with virtually no PEPCK engagement, through low engagement in Flaveria species, intermediate engagement in maize and sugarcane (Saccharum species), to virtually complete engagement in the atypical C 4 Alloteropsis semialata (Gutierrez et al., 1974;Ueno and Sentoku, 2006;Wang et al., 2014a;Lundgren et al., 2016;Dunning et al., in review). Similarly, in NAD-ME subtypes the gradient spans Panicum species with virtually no PEPCK engagement, Cleome C 4 species with intermediate PEPCK activity, and Spartina maritima where PEPCK engagement is virtually complete (Gowik and Westhoff, 2011;Bellasio and Griffiths, 2014c;Lundgren et al., 2016). Regardless of how this variability is classified (Furbank, 2011;Wang et al., 2014a;Koteyeva et al., 2015), the SMA allows rates to vary continuously between boundaries set by what is biochemically realistic.
Given that PEPCK regenerates PEP with half the ATP required by PPDK, moderate levels of PEPCK engagement could potentially increase the biochemical efficiency of assimilation (compare ATP TOT /GA in Supplementary  Fig. S5A with Fig. S5D). The rapid regulation of PPDK  could contribute to the flexibility and efficiency of the CCM, supported by the observation that redundant decarboxylating pathways appeared multiple times at late evolutionary stages . At high levels of PEPCK engagement ( Supplementary Fig. S5D, E), the ATP/GA predicted is substantially lower than for the other photosynthetic subtypes. However, biochemical efficiency may be reduced by the need to hydrolyse part of the newly synthetized PEP to drive the PEPCK reaction, which is close to thermodynamic equilibrium and may be too slow to support physiological decarboxylation rates (Huber and Edwards, 1975). Clarification is still required, and there may be interspecific variability (Smith and Woolhouse, 1983); however, circumstantial evidence gained in comparative experiments have shown lower quantum efficiency of PEPCK plants [Furbank (2011) and references therein].
PEPCK is required for the CCM to work, and when PEPCK is the only decarboxylating enzyme it is likely to be modulated solely by the requirements of the CCM. As a consequence, in PEPCK plants, PEPCK activity cannot contribute to finetuning ATP demand, and hence the biochemical plasticity of the CCM is lower. The additional ATP used by PEPCK can therefore be considered part of the minimum ATP demand in the BS, ATP BSMIN = RuP phospBS + PEPCK + V OBS + R, corresponding to a minimum ATP demand partitioning ratio of ATP ATP  (Carnal et al., 1993). This process is less energy efficient than photophosphorylation (Buckley and Adams, 2011;Kramer and Evans, 2011), and may contribute to the generation of ATP BSMIN only when illumination of the BS chloroplast is insufficient.

Partitioning CS
The main photosynthetic products in C 4 species are starch and sucrose. Both can be synthesized in the M and BS (Majeran and van Wijk, 2009;Friso et al., 2010), but sucrose is preferentially synthesized in the M, while starch disproportionately accumulates in the BS (Furbank et al., 1985;Lunn and Furbank, 1997). It has been observed that accumulation of sucrose in the leaf does not directly influence the partitioning between sucrose and starch (Lunn and Hatch, 1995), apparently discounting a role for starch synthesis as a sink for carbon overspill when sucrose synthesis is inhibited by sucrose build-up (Lunn and Hatch, 1995). Experiments conducted with Panicum species examined the effect of illumination on the distribution of sucrose-phosphate synthase activity between the M and BS, showing contrasting activation patterns between subtypes (Ohsugi and Huber, 1987).

Effect of PPDK engagement in the BS and fast regulation of ATP demand
PPDK was shown to be present and active in the BS (Aoyagi and Nakamoto, 1985;Friso et al., 2010), although its elusive role is generally not included in textbook descriptions of C 4 photosynthesis. The engagement of PPDK increases ATP demand in the BS in the same way as PR, as both processes require two ATP per catalytic turnover. However, PPDK does not require NADPH, whereas PR does. I propose that the tight regulation of PR and PPDK activities (and, in some cases, PEPCK, see above) contribute to fine-tuning ATP and NADPH demand in the BS in response to illumination of the BS chloroplast. Light availability in the BS is determined by anatomy and changes dynamically according to light intensity and quality. The locality of ATP production is therefore largely independent of metabolic control [Bellasio and Lundgren (2016) and references therein]. Local ATP imbalances cannot be rebalanced by ATP diffusion because maintaining ATP concentration and a high ATP to ADP ratio in each cell compartment is critical. The only way to balance supply and demand in each compartment is to regulate the locality of ATP demand (Evans et al., 2007;Bellasio and Griffiths, 2014c). For instance, blue-rich diffuse sky radiation is strongly absorbed in the superficial M and results in preferential ATP production in the M (Evans et al., 2007). In these conditions BS activity can be downregulated, and countered by a proportional increase in M activity, to maintain assimilation and biochemical efficiency. However, operating the RPP, the C 2 , and PCO cycles in the BS impose a limit, and a threshold of ATP demand in the BS (ATP BSMIN = RuP phospBS + V OBS + R), corresponding to an ATP partitioning ratio of ATP ATP BSMIN M ≈ 0 3 . . If ATP supply in the BS is lower than the minimum demand, C 4 photosynthesis ceases, leading to stunted phenotypes (Bellasio and Griffiths, 2014c;McQualter et al., 2016). More penetrating light (e.g. red light or direct sunlight) may drive higher ATP synthesis in the BS. In this case, PR and PEPCK may be upregulated, but are limited by the amount of NADPH available in the BS through the MAL shuttle or through LEF, as discussed above. PPDK may be engaged to take advantage of additional ATP, which may be made available under transient exposure to even more penetrating light qualities (e.g. green-enriched canopy light).

Conclusion
A SMA has been developed as a modelling framework based solely on mass-balance constraints, which generalizes all pathways of assimilation (except CAM) in a single tool. A range of examples detailed the energetics and metabolite fluxes involved in the gradual activation of the C 2 shuttle and the CCM along a spectrum of photosynthetic subtypes from C 3 to C 4 photosynthesis. This knowledge is important for basic and applied research, to support advanced breeding techniques, or to study natural variability of biochemical traits. For instance, by providing quantitative estimates for fluxes and energy requirements in the BS and M, it is possible to set realistic targets for bioengineering projects. Future work will integrate biochemical models to allow the SMA to respond directly to environmental variables.

Availability
SMA is coded in an Excel workbook that is freely available to download. Macros are avoided.

Supplementary data
Supplementary data are available at JXB online. File 1. SMA development and equations S1-S34; simulation of static scenarios under physiological operational conditions: Table S1 and Figures S1, S3-S5E; simulation 1.2, a dynamic scenario for C 3 photosynthesis: Figure S2.
File 2. Excel workbook coding the SMA.