A Method of Accounting for Enzyme Costs in Flux Balance Analysis Reveals Alternative Pathways and Metabolite Stores in an Illuminated Arabidopsis Leaf1[OPEN]

A diel genome-scale model of leaf metabolism highlights the inherent flexibility in meeting the demands for daytime carbon storage and light-dependent energy balancing. Flux balance analysis of plant metabolism is an established method for predicting metabolic flux phenotypes and for exploring the way in which the plant metabolic network delivers specific outcomes in different cell types, tissues, and temporal phases. A recurring theme is the need to explore the flexibility of the network in meeting its objectives and, in particular, to establish the extent to which alternative pathways can contribute to achieving specific outcomes. Unfortunately, predictions from conventional flux balance analysis minimize the simultaneous operation of alternative pathways, but by introducing flux-weighting factors to allow for the variable intrinsic cost of supporting each flux, it is possible to activate different pathways in individual simulations and, thus, to explore alternative pathways by averaging thousands of simulations. This new method has been applied to a diel genome-scale model of Arabidopsis (Arabidopsis thaliana) leaf metabolism to explore the flexibility of the network in meeting the metabolic requirements of the leaf in the light. This identified alternative flux modes in the Calvin-Benson cycle revealed the potential for alternative transitory carbon stores in leaves and led to predictions about the light-dependent contribution of alternative electron flow pathways and futile cycles in energy rebalancing. Notable features of the analysis include the light-dependent tradeoff between the use of carbohydrates and four-carbon organic acids as transitory storage forms and the way in which multiple pathways for the consumption of ATP and NADPH can contribute to the balancing of the requirements of photosynthetic metabolism with the energy available from photon capture.

Computational modeling of metabolism is increasingly used to analyze the complexity of plant metabolic networks and to understand system-level properties such as carbon use efficiency (Sweetlove and Ratcliffe, 2011;Nägele and Weckwerth, 2012;de Oliveira Dal'Molin and Nielsen, 2013;Kruger and Ratcliffe, 2015). Flux balance analysis (FBA), which is a method for predicting steadystate flux distributions using a stoichiometric model of the network, is particularly well suited to this task, because it can be applied to large-scale metabolic networks (Lewis et al., 2012). It is also computationally efficient, meaning that models of different cell types (de Oliveira Dal'Molin et al., 2010), different temporal phases (Cheung et al., 2014), and different tissues (Borisjuk et al., 2013;Grafahrend-Belau et al., 2013) can be combined.
FBA can generate accurate predictions of plant metabolic fluxes (Williams et al., 2010;Hay and Schwender, 2011;Cheung et al., 2013), but the analysis is complicated by the presence of alternative pathways that share the same function within the network. For example, mitochondria and chloroplasts have several potential mechanisms for maintaining energetic homeostasis, including alternative electron flow pathways, metabolite shuttles for the transfer of reducing power or ATP, and uncoupling mechanisms (Millar et al., 2011;Taniguchi and Miyake, 2012). More generally, the distributed robustness of metabolic networks means that they have the inherent property of being able to achieve cellular objectives in different ways (Wagner, 2005). However, FBA does not automatically identify these alternative flux distributions, because the immediate output of the analysis is a single flux distribution that satisfies the constraints and objectives applied to the model. This has the effect of masking the potential contribution of alternative pathways, and to avoid this, it is necessary to extend the analysis in a way that will reveal them.
The most commonly used approach to this problem is flux variability analysis (FVA), which defines the permissible flux ranges for each reaction in the optimal flux space (Mahadevan and Schilling, 2003). Another possibility is random sampling of the optimal flux space, using a uniform sampling algorithm that was originally introduced to characterize the entire feasible flux solution space (Price et al., 2004). While both approaches are useful for exploring the capability of the metabolic system in achieving the cellular objectives, they do not give any indication of which alternative optimal flux solutions are more likely to be found in vivo, and they do not generate flux distributions that represent the biological reality in which alternative pathways may be operating simultaneously.
Here, we develop a methodology that permits alternative pathways to be explored efficiently and that allows the consequences of the simultaneous operation of alternative pathways on the rest of the metabolic network to be examined. Our approach emerged from a reconsideration of the use of flux minimization as the objective function. Minimization of the sum of the absolute flux values supported by all the reactions in the network is often used as an objective function in FBA on the principle that cells have evolved to minimize the costs for the synthesis of the enzymes and membrane transporters that support growth and cell maintenance (Holzhütter, 2004). However, no weighting is applied when calculating the sum of fluxes, so there is an implicit assumption that the machinery cost per unit of flux is the same for all reactions. This assumption is invalid in vivo, as enzymes vary in terms of their size, number of subunits, and catalytic capacity. Ideally, each reaction should be weighted by its enzyme machinery costs per unit of flux, but such information is not available for the majority of the reactions in large-/ genome-scale metabolic models. Here, we develop a modeling method, cost-weighted FBA, that avoids the invalid assumption of equal costs and that allows the evaluation of alternative metabolic routes in a complex network. The method was used to demonstrate the flexibility of leaf metabolism in meeting the metabolic requirements of an Arabidopsis (Arabidopsis thaliana) leaf in the light. A genome-scale diel FBA model, in which the light and dark phases of the diel cycle were solved as a single optimization problem, was used for the analysis because this approach currently provides the most realistic constraints-based framework for modeling leaf metabolism (Cheung et al., 2014;de Oliveira Dal'Molin et al., 2015).

Cost-Weighted Flux Minimization
The use of flux-weighting factors to explore alternative pathways was first examined using a small metabolic network ( Fig. 1) with the input flux (W_in) set to carry one unit of flux. Conventional FBA using flux minimization as the objective function arbitrarily chooses one of the three parallel reactions that convert W to X, here R1, and selects R4 for the conversion of X to Z because the reaction involves only a single step (Table  I). Applying FVA reveals the equivalence of R1, R2, and R3, with each taking a flux value of between 0 and 1, emphasizing both the inability of the constraints to produce a unique solution and the potential existence of alternative flux distributions to the one picked by conventional FBA. Geometric FBA (Smallbone and Simeonidis, 2009), which identifies a unique flux solution that corresponds to the center of the optimal flux solution space, predicts an equal division of flux between R1, R2, and R3 and a flux of 1 unit through R4. As an alternative to these established methods, using 1,000 sets of randomly chosen weighting factors and averaging the flux solutions gave predicted fluxes of approximately onethird through R1, R2, and R3, five-sixths through R4, and one-sixth through R5 and R6 (Table I). The costweighting procedure correctly identifies the ambiguity in the prediction of the fluxes between W and X, and it also emphasizes the possibility that the route via R5 and R6 might be important, depending on the relative and unknown costs of supporting the fluxes through R4, R5, and R6. In essence, the new method captures alternative metabolic routes that would be neglected by conventional flux minimization because of their higher unweighted sum of fluxes. Note also that an analytical solution confirms the predicted fluxes for the model analyzed here (Supplemental File S1) but that, in general, it would be impractical to obtain such solutions for all the alternative pathways in a large-scale network.
The flux distributions predicted by cost-weighted FBA can be used to identify reactions that operate together or in parallel by calculating the Pearson correlation coefficient (r) for the set of reaction fluxes. This approach is similar to the concepts of reaction correlation coefficient (Poolman et al., 2007) and flux correlation (Poolman et al., 2009) used in conventional FBA, and it leads to a set of r values for the small metabolic network (Table II). Reactions that operate together have a positive r (e.g. between R5 and R6); reactions that work against each other (e.g. in parallel pathways) have a negative r (e.g. between R1, R2, and R3); and reactions that support fluxes with no necessary correlation have r close to 0 (e.g. between R1 and R4). These reaction correlation coefficients complement the flux range information from FVA by characterizing sets of reactions that operate together or in parallel. While alternative metabolic routes can be easily identified by inspection in a small metabolic network (Fig. 1), this is seldom possible in large-scale models with hundreds or even thousands of reactions.
The solutions for the small network generated by conventional FBA with flux minimization and costweighted FBA can be used to compare the limitations Figure 1. Model for investigating the utility of the cost-weighted flux minimization objective function. The conversion of W to X can be carried out by three single-step reactions, R1, R2, and R3, while the conversion of X to Z can be achieved by a single step, R4, or by two steps, R5 and R6. The input (W_in) and output (Z_out) fluxes are defined, and FBA is used to predict the fluxes through the six internal steps, R1 to R6. of the two approaches. Conventional FBA generates a single solution, but FVA shows the equivalence between R1, R2, and R3, so it is unclear whether the flux between W and X actually uses one, two, or all three routes in some combination. Conventional FBA also ignores the possibility that the flux between X and Y could occur via R5 and R6. In contrast, cost-weighted FBA generates multiple solutions, and each one has the same limitations as a conventional FBA solution (i.e. each solution will be the one with the lowest flux cost given the defined weightings). However, by varying the weightings and generating an averaged solution, it is made clear that there are multiple ways in which the network might function. The averaged solution is not an optimal solution for the conventional FBA problem of flux minimization with uniform weighting (Table I shows that the sum of fluxes for a model with equal weightings is 2 units, whereas the averaged solution gives 2.161 units), but it does emphasize the potential contribution of alternative pathways. These may or may not operate simultaneously, but without further biochemical knowledge, it is not possible to rule out the possibility that all the pathways in the averaged solution contribute to the flux distribution. Ultimately, costweighted FBA explores a larger flux space satisfying the applied constraints and objective function and, thus, highlights the way in which alternative pathways can be harnessed for a particular objective.

Application of Cost-Weighted FBA to a Diel Leaf Model
The flexibility of the metabolic network in a C 3 leaf was investigated by applying flux-weighting factors to a diel Arabidopsis genome-scale metabolic model (Cheung et al., 2014). Flux solutions were obtained for 1,000 sets of randomly chosen weighting factors (Supplemental Table  S1), and most of the reactions in photosynthesis, photorespiration, Suc synthesis, and nitrogen assimilation were found to be well defined with little difference in flux between solutions. However, high flux variability was observed for several reactions in the Calvin-Benson cycle, including the steps catalyzed by sedoheptulose bisphosphatase and Fru-1,6-bisphosphatase. A scatterplot of the predicted reaction fluxes for the two phosphatases in the light showed that the solutions lay in three distinct clusters ( Fig. 2A) and that the values for the Pearson correlation coefficients for the fluxes catalyzed by sedoheptulose bisphosphatase, Fru-1,6bisphosphatase, and transaldolase were greater than 0.999 (Supplemental Table S2). Inspection of the flux solutions (Supplemental Table S1) showed that they corresponded to three distinct flux modes within the network of reductive pentose phosphate metabolism (Fig. 2, B-D). One flux mode, with equal flux through the two phosphatases, corresponds to the conventional Calvin-Benson cycle (Fig. 2C); while the other two flux modes correspond to situations in which one of the phosphatases works in tandem with transaldolase ( Fig. 2, B and D).
Cost-weighted FBA showed a high variability in the accumulation of storage carbon during the day, with a tight tradeoff (r = 20.9992; Supplemental Table S2) between carbohydrate (starch and soluble sugars) and four-carbon organic acids (malate and fumarate; Fig.  3A). This suggests that organic acids can complement carbohydrates for carbon storage during the day, but interestingly, organic acids did not completely replace the carbohydrate store, since there was a minimum requirement for carbohydrate storage at a rate of 0.429 mmol m 22 s 21 . Note that if the model was prevented from storing carbohydrate, it stored organic acids alone, but this required a higher photon input (186.9 versus 182.1 mmol m 22 s 21 ) and a higher sum of fluxes (621.1 versus 598.2 mmol m 22 s 21 ) than the solutions obtained when the model could choose between the different storage forms (data not shown). In fact, the accumulation of four-carbon organic acids was positively correlated (r . 0.93) with net carbon dioxide uptake, Rubisco, and phosphoenolpyruvate carboxylase, indicating that the accumulation of four-carbon organic acids required a higher carbon assimilation rate than the storage of carbohydrates. Table II. Pearson correlation coefficients for the set of fluxes predicted by cost-weighted FBA for the small metabolic network shown in Figure 1 Pearson correlation coefficients between reaction fluxes were calculated from 1,000 flux solutions computed using cost-weighted FBA. Positive values represent reactions that operate together; negative values represent reactions in parallel pathways and/or reactions that act in opposite directions; values close to 0 indicate that the reactions are independent of each other.  Table I. Comparisons of flux predictions for the small metabolic network shown in Figure 1 using different flux minimization modeling methods The FBA column contains the flux solution computed from standard FBA using flux minimization with uniform weighting factors. FVA was performed with flux minimization with uniform weighting factors as the primary objective function. Flux ranges from FVA are represented as (v min ,v max ), where v min and v max are the minimum and maximum feasible flux values compatible with the optimal objective value. Geometric FBA was used to determine the flux solution corresponding to the center of the optimal flux space. The averaged fluxes from the cost-weighted analysis were calculated from 1,000 flux solutions.
The scatterplot (Fig. 3A) presents the results for 10 5 sets of random weighting factors. Generally, 1,000 sets of random weighting factors were sufficient to sample the cost-weighted solution space, but here, a more detailed sampling was undertaken to investigate whether the apparent threshold for a minimum carbohydrate store might be overcome in rare cases. This analysis revealed that a tiny fraction of the solutions (approximately 250 out of 10 5 ) corresponded to a significantly reduced total storage of carbon (carbohydrate + organic acids). Inspection of the network fluxes in these solutions showed that they had reduced synthesis and accumulation of citrate during the night, emphasizing that the availability of carbon skeletons synthesized at night for nitrogen assimilation during the day (Gauthier et al., 2010;Cheung et al., 2014) depends on the relative costs of the alternative pathways for the provision of carbon skeletons that are available in the day and the night.
Although light influx into the model was unconstrained for the simulations shown in Figure 3A, choosing minimization of the sum of fluxes as the objective function tends to produce solutions with a minimal light Figure 2. Cost-weighted analysis of the reductive pentose phosphate pathway in a diel C 3 leaf model. A, Scatterplot for the sedoheptulose 1,7-bisphosphatase and Fru-1,6-bisPase fluxes in the light obtained from 1,000 runs of cost-weighted FBA showing three distinct clusters. B to D, Flux modes in the reductive pentose phosphate pathway corresponding to the three clusters observed in the scatterplot. BPGA, 1,3-Bisphosphoglycerate; DHAP, dihydroxyacetone phosphate; E4P, erythrose 4-phosphate; F6P, Fru-6-P; FBP, Fru-1,6-bisP; GAP, glyceraldehyde 3-phosphate; PGA, 3-phosphoglycerate; P i , inorganic phosphate; R5P, Rib-5-P; Ru5P, ribulose 5-phosphate; RuBP, ribulose 1,5-bisphosphate; S7P, sedoheptulose 7-phosphate; SBP, sedoheptulose 1,7-bisphosphate; X5P, xylulose 5-phosphate. requirement. Given that the accumulation of organic acids requires a higher photon influx, this may explain the threshold for the maximum amount of organic acid that can be accumulated. To test this, a new set of simulations was run in which light influx was constrained to a range of different values ( Fig. 3B; Supplemental Results S1). This showed that increasing the light intensity lowered the minimum requirement for carbohydrate storage, although still not to 0, reduced the clustering around particular combinations of storage forms (Supplemental Results S1), and increased the number of solutions with a significantly reduced total storage of carbon (43.1% of 10 5 solutions in Fig. 3B). These results confirm that increasing the light intensity gives the network more freedom to select less energetically efficient routes, with the bulk of the observed changes occurring in the range up to 226 mmol m 22 s 21 .

Effect of Light Intensity on Leaf Metabolism
The effect of light intensity on daytime leaf metabolism was modeled by varying the photon influx while fixing Suc and amino acid export and the cellular maintenance costs. This corresponds to a carbonlimited scenario where an increase in light intensity does not lead to an increase in net carbon fixation. In keeping with the fixed constraints on productivity and maintenance, net CO 2 assimilation rate in the light remained constant, with only 3.52% deviation from the average value, when the photon influx was varied between 180 and 272 mmol m 2 2 s 2 1 (Figs. 4 and 5; Supplemental Tables S3 and S4). Unweighted (Fig. 4A) and cost-weighted (Fig. 5A) flux minimization gave the same result for the net CO 2 assimilation rate, and at the minimum photon flux required for cell growth and maintenance, the only electron transport processes were photosynthetic LEF (Figs. 4B and 5B) and the mitochondrial ETC (cytochrome pathway; Figs. 4C and 5C). As the light intensity increased, multiple additional fluxes were engaged to match the supply of ATP and reducing power with the invariant demand of the model, and a major difference emerged between the two modeling approaches: conventional flux minimization with no reaction weighting solved the problem by selecting the minimum number of extra fluxes (Fig.  4, C-F), as was also observed when the photon flux was varied in a rice (Oryza sativa) genome-scale model , whereas the use of fluxweighting factors emphasized the potential contribution of alternative fluxes (Fig. 5, C-F). The latter provides a more biologically relevant view of the network, with a number of alternative pathways such as the mitochondrial alternative electron transport pathways that are known to be expressed in Arabidopsis leaves (Baerenfaller et al., 2008) carrying flux, in contrast to conventional FBA, which predicts zero flux for these reactions.
Photosynthetic LEF increased linearly with photon influx (Figs. 4B and 5B), and the use of flux-weighting factors showed how different processes could contribute to the disposal of excess reductant at different light intensities. As the photon influx increased from 180 to . Predicted effects of light intensity on the daytime fluxes of leaf metabolism obtained using flux minimization with uniform reaction weightings. A, Net CO 2 assimilation. B, Photosynthetic linear electron flow (LEF) flux through PSII (blue circles) and PSI (red triangles). C, Fluxes through the mitochondrial electron transport chain (ETC; blue circles), cyclic electron flow (CEF; red triangles), and water-water cycle (WWC; green crosses) normalized to two electrons. D, Fluxes through the chlorophyll cycle (blue circles) and the xanthophyll cycle (red triangles). E, Fluxes through the alternative NADH dehydrogenase (blue circles), the alternative oxidase (red triangles), and the uncoupling protein (green crosses). F, Fluxes through futile cycles involving carbamoyl phosphate (blue circles), formate-tetrahydrofolate (THF; red triangles), chloroplastic phosphofructokinase (green crosses), and peroxisomal palmitate transport (yellow crosses).
195 mmol m 22 s 21 , the model predicted a decrease in flux through the mitochondrial ETC (Fig. 5C) with a corresponding decrease in mitochondrial ATP synthesis. This was compensated by ATP synthesis in the chloroplast, partly from LEF and partly from CEF and the WWC (Fig. 5C). As the photon influx increased from 195 to 215 mmol m 22 s 21 , the mitochondrial ETC flux continued to decrease, and the CEF and WWC passed through maxima as the ATP produced by the LEF flux continued to increase. Cost-weighted FBA also predicted roles for the mitochondrial alternative respiratory pathway, which includes the alternative NADH dehydrogenase and the alternative oxidase (Fig. 5E) as well as the mitochondrial uncoupling protein (Fig. 5E). Above 215 mmol m 22 s 21 , the CEF and WWC fluxes increased again, the alternative respiratory pathway became more important than the mitochondrial ETC, and a role emerged for several ATP-dissipating reactions (Fig. 5F). These futile cycles may not all operate in vivo, but the use of flux-weighting factors reveals the full capability of the metabolic network, in contrast to the output of flux minimization with uniform weighting (Fig. 4F). Note that the relative contribution of the predicted fluxes that continue to increase above 220 mmol m 22 s 21 remained constant, indicating that the optimal flux solution was now independent of the light intensity, whereas at lower intensities, the relative contributions of the predicted fluxes varied as the model phased out the mitochondrial electron transport chain as a source of ATP.

Flux Weighting as a Method for Identifying Alternative Metabolic Pathways
The optimization objective functions that are commonly used in the constraints-based analysis of plant metabolic networks produce solutions that minimize the utilization of alternative pathways and, thus, overlook their potential contribution to network function. The method developed here, cost-weighted FBA, uses randomly chosen flux-weighting factors to overcome this problem, thus permitting the straightforward identification of alternative flux solutions. By running the model thousands of times, each with a different set of flux-weighting values, and then averaging the solutions, it is possible to assess feasible alternative metabolic routes in the network simultaneously. This is an improvement on random sampling of the flux solution space (Price et al., 2004) because the combination of random weighting with flux minimization increases the biological relevance of the flux predictions. Costweighted FBA is also superior to geometric FBA (Smallbone and Simeonidis, 2009), which only identifies alternative pathways with an equal number of reactions (Table I). Moreover, the use of flux-weighting factors allows a more nuanced exploration of alternative pathways than FVA, which only generates a permissible range for each flux satisfying the objective function and constraints. In contrast, cost-weighted FBA accepts that the real costs of each step in the Figure 5. Predicted effects of light intensity on the daytime fluxes of leaf metabolism obtained using cost-weighted flux minimization. A, Net CO 2 assimilation. B, Photosynthetic LEF flux through PSII (blue circles) and PSI (red triangles). C, Fluxes through the mitochondrial ETC (blue circles), CEF (red triangles), and WWC (green crosses) normalized to two electrons. D, Fluxes through the chlorophyll cycle (blue circles) and the xanthophyll cycle (red triangles). E, Fluxes through the alternative NADH dehydrogenase (blue circles), the alternative oxidase (red triangles), and the uncoupling protein (green crosses). F, Fluxes through futile cycles involving carbamoyl phosphate (blue circles), formate-tetrahydrofolate (THF; red triangles), chloroplastic phosphofructokinase (green crosses), and peroxisomal palmitate transport (yellow crosses). The error bars represent SD values obtained in a robustness analysis based on 10 independently calculated averaged solutions, each generated from 1,000 sets of randomly chosen weighting factors (Supplemental Table S4).
network are different, reflecting differences in enzyme size, turnover, and catalytic properties, and predicts flux distributions that show how combinations of alternative pathways can satisfy the constraints on the model. Indeed, it is arguable that cost-weighted FBA gives the most legitimate solution to the flux minimization problem because it is unbiased by the unknown machinery costs of supporting the fluxes.
One consequence of this is to smooth out the abrupt and physiologically unrealistic on-off switches and discontinuities that are commonly observed when FBA is used to model the effects of a variable external parameter such as light intensity ( Fig. 4; Nogales et al., 2012;Poolman et al., 2014). More generally, while the efficiency of network structure, as assessed by the minimization of total flux through the network, has been shown to lead to realistic flux distributions in plant metabolic networks (Williams et al., 2010;Hay and Schwender, 2011;Cheung et al., 2013), costweighted flux minimization highlights the important caveat that greater predictive accuracy might be achievable with a knowledge of the unknown costs of supporting fluxes through different metabolic steps. This is to be expected, because lack of information, whether in the form of kinetic parameters for kinetic models or constraints for constraints-based models, is often encountered when predicting fluxes in complex networks. Here, the focus is on the enzyme machinery costs, but other methods are possible, including a technique known as lumped hybrid cybernetic modeling (Song and Ramkrishna, 2010) in which the flux distribution is determined by choosing elementary modes with optimal yields for ATP and biomass production. This method has been applied to yeast metabolism, but the optimization of ATP and biomass production is not immediately relevant to the metabolism of a mature leaf under high light intensities, and in any case, the method is not easily scalable given the combinatorial explosion of elementary modes that occurs as the size of the model increases.

Alternative Flux Modes in Central Metabolism
The application of randomly chosen weighting factors might be expected to lead to a continuum of alternative solutions. However, in practice, solutions will switch between flux modes within different reaction subnetworks and, hence, will cluster around discrete flux distributions. This is clearly revealed within the reductive pentose phosphate pathway, for which three discrete flux distributions were obtained, each corresponding to a different flux mode (Fig. 2). One of these flux distributions is the classic Calvin-Benson cycle, while the other two are formed by using transaldolase to bypass either Fru-1,6-bisPase or sedoheptulose 1,7bisphosphatase. The physiological significance of the alternative flux modes involving transaldolase in the light may be marginal, since there is some evidence for the light inhibition of transaldolase (Anderson, 1981; but contrast Caillau and Quick, 2005), and transgenic manipulation of sedoheptulose 1,7-bisphosphatase supports a dominant role for the enzyme in CO 2 assimilation (Lefebvre et al., 2005). The potential involvement of transaldolase was revealed previously using elementary modes analysis of the Calvin-Benson cycle in tandem with kinetic modeling (Poolman et al., 2000), but this approach can only be used on relatively small metabolic networks, whereas cost-weighted FBA can reveal much more complex flux modes in largescale networks.
An example of this is the use of alternative stores of daytime-synthesized carbon to support night metabolism in the diel leaf model. Depending on the networkweighting factors, the model either stored carbon exclusively as carbohydrate or as a mixture of carbohydrate and four-carbon organic acids, in a tight tradeoff. Two features of the observed relationship (Fig. 3) reveal the operation of alternative flux modes. First, a subset of solutions stored less carbon overall. This involved an alteration in the source of carbon skeletons for daytime nitrogen assimilation from stored citrate to newly synthesized citrate, necessitating different flux distributions in the day and night. Second, it is apparent that the distribution of solutions in Figure 3 is not a continuum: there was a tendency for solutions to cluster around particular combinations of the storage forms (Supplemental Results S1). These combinations correspond to the operation of different flux modes in the wider network as a consequence of the daytime synthesis and nighttime mobilization of different proportions of carbohydrate and organic acids.
The other aspect of the analysis of alternative carbon stores that is remarkable is that there was no combination of weighting factors that allowed the storage of only organic acids. In modeling terms, this is because such solutions lead to a higher sum of fluxes as a result of the requirement to assimilate more CO 2 when storing organic acids than carbohydrate. This higher CO 2 assimilation rate reflects the high flux through the anaplerotic phosphoenolpyruvate carboxylase reaction for the net synthesis of four-carbon organic acids. This high daytime flux through phosphoenolpyruvate carboxylase would compete with Rubisco for available CO 2 , and this only commonly occurs in C 4 plants where the two enzymes are spatially segregated. Nevertheless, the model demonstrates that it is entirely possible, although metabolically costly, to replace part of the leaf carbohydrate store with four-carbon organic acids; in fact, it is well established that starch and fumarate are used as alternative carbon storage forms in Arabidopsis leaves (Chia et al., 2000;Pracharoenwattana et al., 2010). In this instance, the factor that determines the unusually high accumulation of an organic acid as a carbon storage form in the day is the role played by the acid in the regulation of pH during rapid growth on nitrate, emphasizing the point that there can be other considerations, beyond those captured in the modeling analysis, that can influence the actual flux distribution through a network. In passing, it is also interesting that the proposed optimal synthetic carbon fixation pathway (Bar-Even et al., 2010) is based on organic acid metabolism, and the implementation of this in a photosynthetic context would pose the additional challenge of engineering sufficient capacity for the storage of a non-pH-perturbing carbon form.

Alternative Routes of Electron Flow and Energy Rebalancing in Photosynthetic Metabolism
Leaf metabolism in the light is known to involve multiple routes for balancing the supply and demand for ATP and NADPH (Noctor and Foyer, 1998;Allen, 2003;Kramer and Evans, 2011;Taniguchi and Miyake, 2012). One immediate problem is that the ratio of ATP to NADPH arising from light-driven LEF does not match the requirements for the fixation of carbon dioxide in the Calvin-Benson cycle (Noctor and Foyer, 1998;Allen, 2003); another problem is the need to dispose of excess energy or reducing equivalents under conditions in which an imbalance arises between supply and demand. It is known that alternative pathways exist that alter the proportional synthesis of ATP and NADPH in the chloroplast, including CEF and the WWC (Allen, 2003). It is also apparent that there are many other potential mechanisms for rebalancing ATP and NADPH, such as the operation of the mitochondrial ETC and other pathways that consume ATP and/ or reducing power, such as the chlorophyll cycle and the xanthophyll cycle. The contributions of these pathways are sensitive to environmental and physiological conditions, reflecting the impact of these factors on energy supply and the metabolic demand of the system. However, the relative contribution of the different pathways to the overall energy balance of the leaf is not well understood.
Recently, FBA was used to analyze the contribution of alternative electron flows to photosynthetic metabolism in the cyanobacterium Synechocystis sp. PCC6803, and it was concluded that these pathways were necessary for optimal growth under both light-limiting and carbon-limiting conditions (Nogales et al., 2012). However, the limitations of conventional FBA meant that the alternative pathways were analyzed individually, with only one extra pathway permitted to carry flux in each simulation, and this was not ideal because it is known that more than one pathway can operate simultaneously. This problem can be avoided using fluxweighting factors, so we used this new approach to assess the likely relative importance and simultaneous operation of alternative pathways for energy rebalancing in the photosynthetic metabolism of an Arabidopsis leaf. To do this, the light input into the diel leaf model was varied and the network fluxes were predicted using cost-weighted FBA while keeping the carbon and nitrogen demands for the export of Suc and amino acids to the phloem constant. The biosynthetic and maintenance constraints in the model required the light intensity to exceed a threshold value (Figs. 4 and 5), and at the lowest feasible light intensity, alternative electron flows were required to rebalance the ATP to NADPH ratio generated by LEF. This is in agreement with the conclusion of Nogales et al. (2012) that alternative electron flows are required for the optimal performance of a photosynthetic network. In the Arabidopsis model, the role of these pathways at the lowest light intensity was to convert NAD(P)H to ATP, and the mitochondrial ETC was predicted to be the major route for this purpose (Fig. 5C). The choice of the mitochondrial ETC over other pathways can be explained by the higher ATP yield of the mitochondrial ETC: 0.945 ATP per photon, compared with 0.45 and 0.32 for CEF and the WWC, respectively (Kramer and Evans, 2011). Efficient conversion of reductant to ATP is important when the energy flux (light) into the system is low and the predicted mitochondrial ETC flux is consistent with experimental evidence for the importance of the mitochondrial ETC in sustaining photosynthetic carbon assimilation (Raghavendra and Padmasree, 2003;Nunes-Nesi et al., 2008), implying that mitochondria fulfill an important role in rebalancing reductant and ATP in the light.
As the light intensity increased, flux through the photosynthetic LEF pathway increased proportionally, generating more ATP and NADPH and, thus, changing the rebalancing requirement in relation to the fixed demand. Between 180 and 195 mmol m 22 s 21 , the flux through the mitochondrial ETC decreased in favor of the photosynthetic CEF and, to a lesser extent, fluxes through the alternative NADH dehydrogenases and the alternative oxidase (Fig. 5, C and E). The changing balance between the proton-pumping and non-protonpumping mitochondrial pathways reflects the fact that the mitochondrion is being used to consume excess reductant while needing to generate progressively less ATP. These predictions are consistent with the hypothesis that the mitochondrion can act as a sink for excess reductant (Sweetlove et al., 2006;Miyake, 2010). At a photon influx of approximately 195 mmol m 22 s 21 , flux through photosynthetic CEF and the WWC peaked, reflecting the rising ATP production by photosynthetic LEF and a declining need for ATP synthesis via the action of CEF and the WWC.
The predicted contribution of the non-protonpumping mitochondrial pathways increased up to a light intensity of approximately 220 mmol m 22 s 21 , and there was also a minor contribution from the mitochondrial uncoupling protein that matched the residual (nonzero) contribution from the proton-pumping ETC. The main function of these mitochondrial pathways in this phase is the oxidation of mitochondrial NADH generated by the action of Gly decarboxylase as part of the photorespiratory cycle. In contrast to a simple flux minimization objective function with uniform weightings, cost-weighted flux minimization allows for the possibility that the cytochrome pathway might contribute to the disposal of the reducing equivalents generated by photorespiration, and the resulting small protonpumping flux, which would otherwise lead to the synthesis of ATP, is largely matched by the predicted flux through the uncoupling protein (Fig. 5, C and E). Above 220 mmol m 22 s 21 , the alternative mitochondrial respiratory pathways plateaued and the mitochondrial ETC flux decreased to a low level of approximately 0.2 mmol m 22 s 21 . Gly decarboxylase is the main source of mitochondrial NADH, and the observation that the relative contribution of the pathways for the disposal of mitochondrial NADH stabilizes at moderate light intensities reflects the fact that the rate of photorespiration, and hence Gly decarboxylase, was constant in the model. At higher light intensities, the continued oxidation of excess NADPH generated by LEF was dealt with by chloroplastic pathways, presumably because the transfer of this reductant to mitochondria via metabolite shuttles would require more steps and, thus, a higher sum of fluxes. After the mitochondrial pathways had plateaued, the photosynthetic CEF and WWC began to increase in flux again, with the CEF becoming one of the main routes for NADH rebalancing above a photon influx of 220 mmol m 22 s 21 . This model prediction is consistent with experimental findings that CEF is only important under high light levels (Munekage et al., 2008;Walker et al., 2014). In a CEF mutant of Chlamydomonas reinhardtii, the deficiency in CEF could be compensated by the combined effect of the mitochondrial ETC and WWC (Dang et al., 2014). WWC was shown to be the major alternative electron flow in strains of the photosynthetic dinoflagellate Symbiodinium (Roberty et al., 2014), but there is controversy as to whether it is a major route of electron flow in leaves under illumination (Heber, 2002). The model prediction suggests that WWC would play a minor role in leaves in the light compared with CEF.
At higher light intensities, the model also predicts a major role for the chlorophyll and xanthophyll pigment cycles in dissipating excess reductant, with the chlorophyll cycle having an even greater flux than CEF above light intensities of 240 mmol m 22 s 21 . The major accepted function of the chlorophyll cycle is the synthesis and catabolism of chlorophyll a/b (Tanaka and Tanaka, 2011), and the xanthophyll cycle is well known for its role in photoprotection by dissipating excess light energy as heat (Demmig-Adams and Adams, 1996). The model suggests that both cycles provide efficient mechanisms for maintaining a balanced photosynthetic system through the disposal of excess reducing equivalents.
Another poorly understood area highlighted by the model was the role of futile cycles in dissipating excess ATP. Above a light intensity of 220 mmol m 22 s 21 (where the ATP demand of the modeled metabolic system can be met entirely by the photosynthetic LEF), excess ATP was metabolized by the activation of futile cycles involving the synthesis and degradation of carbamoyl phosphate in the chloroplast by carbamoyl-phosphate synthetase and carbamate kinase and of formate tetrahydrofolate by formate-tetrahydrofolate ligase and formyltetrahydrofolate deformylase in the chloroplast. A palmitate transport cycle across the peroxisomal membrane also makes a minor contribution. Very little is known about the potential operation of these ATP-dissipating mechanisms during photosynthesis, and as such, these are novel predictions that would need to be verified experimentally. Additionally, the model predicts that ATP can be dissipated by futile cycling between Fru-6-P and triose phosphates via the chloroplastic isoforms of phosphofructokinase and Fru-1,6-bisPase. Although the hexose phosphate-triose phosphate futile cycle is well recognized in nonphotosynthetic tissues, where it involves cytosolic isoforms of the enzymes (Kruger and Ratcliffe, 2015), its operation in leaves has not been recognized previously.

Cost-Weighted FBA and Experimentally Defined Flux Weightings
The objective function of flux minimization is a proxy for minimizing the enzyme machinery costs to carry out the metabolic function(s) defined by the model constraints. We have shown that the use of flux-weighting factors in conjunction with the minimization of total flux is an effective way of exploring feasible alternative metabolic pathways in a genome-scale metabolic model. In particular, the application of cost-weighted FBA has identified alternative flux modes in the Calvin-Benson cycle, revealed the potential for alternative transitory carbon stores in leaves, and made predictions about the light-dependent contribution of alternative electron flow pathways and futile cycles in energy rebalancing. In this analysis, the weighting factor for each reaction flux in the objective function corresponds to the enzyme cost per unit of flux. Ideally, a model would be constrained with the actual enzyme cost per unit of flux, but this information is difficult to obtain for most reactions in a genome-scale metabolic model, as it involves not only the size of the enzymes (cost of synthesis) but also their kinetic properties (required amount of the enzyme) and their rate of protein turnover. This has been done for Rubisco in a recent FBA study of Arabidopsis leaf metabolism (Arnold and Nikoloski, 2014), but sufficient data do not yet exist for plants to define the costs of the entire metabolic proteome, as has been achieved for Escherichia coli (O'Brien et al., 2013). An alternative approach is to empirically measure the rate of synthesis of each enzyme in the system by analyzing steady-state protein abundance and protein turnover rates. Recently, a combination of 15 N labeling and proteomics was used to acquire such data for individual proteins (Li et al., 2012;Nelson et al., 2014), and in the future, it may be possible to use such methods to obtain reliable data about protein synthesis rates in absolute terms for the majority of the individual isozymes that make up the metabolic network. Nevertheless, the acquisition of such data are a substantial effort, and although the use of experimentally defined reaction weightings may provide a biologically verified result, it will also tend to lead to model predictions that minimize the use of alternative pathways, because the optimization algorithm will favor the metabolic route(s) with the lowest cost given the defined weightings.
Moreover, the application of experimental weightings in combination with flux minimization will also give spurious results for proteins such as Rubisco, which is the main carboxylating enzyme in leaves despite its high cost. If the goal is to understand and explore alternative pathways, the use of randomly chosen flux weightings is a superior and easier approach.

CONCLUSION
In conclusion, cost-weighted FBA avoids an inherent bias that occurs with unweighted reaction fluxes when the actual costs of enzymes are unknown. This method can be applied to any stoichiometric model of a metabolic network, not just the diel model of leaf metabolism used here, and by exploring a larger feasible flux space, it reveals the potential contribution of alternative pathways in the network. Averaging sets of solutions provides a convenient way to observe the operation of alternative pathways simultaneously and provides a picture that may be more biologically relevant. As well as revealing alternative flux modes in the Calvin-Benson cycle and alternative daytime carbon stores in the Arabidopsis leaf, this approach has highlighted the variable role of alternative pathways for NADPH and ATP rebalancing at different light intensities. The analysis demonstrated the potential importance of several such pathways, including cyclic photosynthetic electron transport, the WWC, mitochondrial uncoupling pathways, and a number of chloroplastic futile cycles that do not carry flux in conventional FBA.

FBA and Cost-Weighted Flux Minimization
FBA was performed by solving the linear programming problem formulated as follows: maximize or minimize where Z is the objective function, c is a vector of weighting factors, v is a vector of metabolic fluxes, N is the stoichiometric matrix, v i is the ith element of v, and min i and max i are the minimum and maximum constraints on v i .
The commonly applied objective function of minimizing the norm (minimization of total flux with equal weighting factor) can be defined as minimize S|c i 3 v i |, where the weighting factor, c i , is identical for all reactions. To allow for the unknown costs of each flux, c i was assigned a random value between 0 and 1 from a uniform distribution for each reaction, and a flux solution was obtained by minimizing S|c i 3 v i |. Multiple flux solutions, typically 1,000 but up to 10 5 when a more detailed analysis was required, were computed by repeatedly solving the optimization problem, each time with a different set of randomly chosen weighting factors. An averaged flux solution was calculated by averaging the multiple flux solutions generated with different sets of c i . Note that the number of flux solutions required depends on the nature of the problem, and this can be determined by testing the robustness of the analysis.
The linear programming problem was solved using the ScrumPy software package (Poolman, 2006) with the Gnu Linear Programming Kit (http://www.gnu.org/software/glpk). FVA and geometric FBA were performed in ScrumPy following the algorithms described by Mahadevan and Schilling (2003) and Smallbone and Simeonidis (2009), respectively.

Quantifying the Relatedness between Reactions by the Pearson Correlation Coefficient
The Pearson correlation coefficient (r) between two reaction fluxes across the multiple flux solutions computed with the flux-weighting factors was used as a measure of relatedness between the two reactions. A value of r = 1 between two irreversible reactions means that the two reactions vary linearly in the same direction, whereas r = 21 indicates that they are competing reactions. For reversible reactions, the sign of r depends on the definition of the reaction directions. A value of r close to 0 indicates that the fluxes are independent of each other. Note that r can only be calculated between reactions carrying variable fluxes across the multiple flux solutions. Reactions with nonvariable flux over the multiple flux solutions are, by definition, correlated with each other.

Metabolic Model and Flux Constraints
A diel Arabidopsis (Arabidopsis thaliana) genome-scale metabolic model was used to model mature C 3 leaf metabolism by setting all fixed carbon to be exported into the phloem and, therefore, no biomass synthesis, as before (Cheung et al., 2014). The model was constrained with values for Suc and amino acid export to the phloem, nitrate uptake, cell maintenance costs, a 3:1 ratio for the carboxylase and oxygenase activities of Rubisco, and other C 3 -specific constraints as described by Cheung et al. (2014). The photon influx into the model was unconstrained unless stated otherwise. For cost-weighted flux minimization, reactions representing CO 2 and oxygen diffusion across membranes, accumulation of storage compounds, biomass synthesis, and generic ATPases and NADPH oxidases for maintenance costs were set to have a weighting of 0, as these reactions were included for modeling purposes and have no directly associated enzyme machinery costs. To relate the model simulations to physiological values, the daytime net carbon assimilation rate in the averaged flux solution from 1,000 sets of flux-weighting factors (Supplemental Table S1) was scaled to the experimental value of 11.8 mmol m 22 leaf area s 21 for Arabidopsis grown at 200 mmol m 22 s 21 (Eckardt et al., 1997).
The response of leaf metabolism to changes in light intensity was explored by defining 101 equally spaced sample points between the minimum possible photon influx (180.14 mmol m 22 s 21 ), calculated by running the model with minimization of photon flux as the objective function, and an arbitrary photon influx (271.98 mmol m 22 s 21 ) above the value where the flux distribution had stabilized. For each sample point, an averaged flux solution was calculated from 1,000 sets of flux-weighting factors. To test the robustness of the analysis, the process was repeated 10 times for 21 of the 101 photon influxes, and SD values were calculated for the means of the averaged flux solutions.

Supplemental Data
The following supplemental materials are available.
Supplemental Table S1. Flux solutions for a diel model of a mature C 3 leaf computed from 1,000 sets of flux-weighting factors.
Supplemental Table S2. Pearson correlation coefficients between selected reactions in a C 3 leaf calculated from the multiple flux solutions obtained by cost-weighted flux minimization.
Supplemental Table S3. Averaged flux solutions obtained by costweighted flux minimization for a mature C 3 leaf with varying photon influx.
Supplemental Table S4. Robustness analysis for the averaged flux solutions obtained by cost-weighted flux minimization for a mature C 3 leaf with varying photon influx.
Supplemental Results S1. Effect of light intensity on the tradeoff between daytime storage of carbohydrate and organic acids in a diel C 3 leaf model.
Supplemental File S1. Analysis of the effect of pathway length on the predicted flux distribution in a small metabolic network.