In the following we describe the reconstruction of an S. bicolor four-cell metabolic model (iTJC1414x4) with the capacity to simulate C₄ photosynthesis as well as to integrate a day phase with a night phase during which storage compounds accumulated during the day are consumed. We based the metabolic reaction network on iEB2140x2 (Bogart and Myers, 2016), a two-cell high confidence representation of Zea mays core metabolism with 635 reactions of which 469 are associated to a total of 2,140 maize genes. To relate the model to S. bicolor, we derived genome wide synteny/orthology associations between maize and sorghum. Predicted protein sequences for amino acid sequences and chromosomal order of 63,480 Z. mays protein encoding genes (B73 v3, Schnable et al., 2009) and of 34,129 S. bicolor protein encoding genes (S. bicolor cultivar BTx623, genome assembly v3.1.1, McCormick et al., 2018) were obtained from Phytozome genomic resource.¹ Z. mays protein sequences were used as queries to identify orthologs in S. bicolor using the SynOrths tool (version 1.0, Cheng et al., 2012). In short, this tool identifies likely pairs of ortholog genes in two related species based on similarities in protein sequences with additional support from homologous flanking genes. For Z. mays genes for which SynOrths did not predict orthologs, protein BLAST (version 2.8.1, Altschul et al., 1997) was used to align amino acid sequences and top hits were assumed to be orthologs if they had at least 90% alignment and 70% identity. Since maize is tetraploid (Swigonova et al., 2004) while S. bicolor is diploid (Price et al., 2005), we expect quite frequently two or more maize genes to map to the same sorghum locus. Indeed, out of over 24,000 S. bicolor genes for which orthology relations were identified, in more than 60% of cases, an S. bicolor gene was orthologous to more than one Z. mays gene (Supplementary Figure 1A).

The Z. mays modeling file (iEB2140x2) was obtained from Bogart and Myers (2016) and transferred to S. bicolor in a spreadsheet format (COBRA toolbox, Microsoft Excel 2016). The model was extended from a two-cell model (x2) connecting a bundle-sheath cell (BSC) and a mesophyll cell (MC) to a diel model (x4). Per sub-model, reactions were expanded or validated using the Kyoto Encyclopedia of Genes and Genomes,² Plant Metabolic Pathway Databases,³ and Phytozome v12.1. The network reactions were revised according to literature, as presented in Supplementary Table 1. Briefly, we moved the reactions for de novo fatty acid biosynthesis (through C₁₈) from the cytosol to the plastid. Acetyl-coenzyme A (CoA) in the plastid can be diverted directly from the CBB cycle by the C₃→C₂ pathway, involving conversion of 3-phosphoglycerate (3-PGA) via 2-PGA, phosphoenolpyruvate and pyruvate to acetyl-CoA (Joyard et al., 2010). TAG was not represented in the initial maize model, so we added its synthesis from phosphatidylcholine (PC) and/or diacylglycerol (DAG), as well as its breakdown into glycerol and free fatty acids. The lipid classes were originally represented by generic molecules (Bogart and Myers, 2016). To simulate metabolic effects resulting from differences in fatty acid chain lengths, we took an approach similar to that taken by Simons et al. (2014) and defined four majorly abundant fatty acid species (C16:0, C18:1, C18:2, and C18:3) to be associated with cytosolic PC, DAG, and TAG molecular species without assigning stereospecific positioning (e.g., sn-1 and sn-2). We added reactions to synthesize these lipid classes, expanded how PC-bound fatty acids are desaturated by cytosol/endoplasmic reticulum-localized desaturase enzymes, and added a PC-dependent acyl editing cycle (Li-Beisson et al., 2013). Next, to allow for testing the effects of futile cycles of lipid synthesis and degradation, we added peroxisomal β-oxidation to recycle cytosolic free fatty acids as well as CoA ligation reactions (Long-Chain Acyl-CoA Synthetase, EC 6.2.1.3) to reuse free fatty acids to synthesize TAG.

For numerical simulations the model was converted to Systems Biology Markup Language (SBML) using a Python⁴ script derived from Hay et al. (2014). Linear programming analyses were performed using the constraint-based reconstruction and analysis (COBRA) toolbox version 3.1 (Heirendt et al., 2019) and GLPK solver⁵ within the MATLAB R2018b environment (The MathWorks, Natick, MA, United States). Mathematical flux ratio constraints were imposed using the COBRA function “addCOBRAConstraints.” The S. bicolor iTJC1414 single-cell model is provided in Supplementary File 1. The full leaf diel iTJC1414x4 model with simulation codes are provided in Supplementary File 2.

Several numerical constraints were added to the model for model simulations to be consistent with the NADP-ME subtype of C₄ photosynthesis (Supplementary Table 2). In plants with the NADP-ME subtype, the malate/pyruvate transfer is the main mechanism to shuttle CO₂ from MCs to BSCs, while a secondary aspartate/malate shuttle mechanism has been reported to be of significance as well in NADP-ME C₄ plants like sorghum (Yin and Struik, 2018). These transport flows through plasmodesmata are generally thought to be driven by strong concentration gradients between the cell types (Wang et al., 2014a; Arrivault et al., 2017), which means that malate, for example, can be assumed to be always moving from MCs to BSCs and never in the opposite direction. Accordingly, among the plasmodesmata transport reactions available in iTJC1414x4, we allowed transport of malate and aspartate from MCs toward BSCs as well as pyruvate and alanine in the opposite direction (Figure 1 and Supplementary Table 2). NADP-ME type species tend to contain little capacity for photoreduction of 3-PGA in the BSCs (Hatch, 1987), which suggests the operation of a 3-PGA/triose phosphate shuttle between BSCs and MCs so that part of the 3-PGA reduction can take place in the MCs. Accordingly, we allowed 3-PGA and glyceraldehyde 3-phosphate to move between the cell types. Due to the presence of various isoforms of ME in the model, cytosolic isoforms and the plastidic NADH-dependent isoform were inactivated so that only the BSC plastidic NADP-ME will be active in decarboxylation of malate in the C₄ cycle. Phosphoenolpyruvate carboxykinase is another decarboxylating enzyme in C₄ metabolism known to be of relevance in species of the NADP-ME subtype. However, the reaction was inactivated in our model because phosphoenolpyruvate carboxykinase activity was reported to be missing in sorghum (Gutierrez et al., 1974; Walker et al., 1997; Döring et al., 2016). Some CO₂ liberated by decarboxylation of malate is expected to leak from BSCs to the MCs due to the concentration gradient between the two cells (von Caemmerer and Furbank, 2003; de Oliveira Dal’Molin et al., 2010; Wang et al., 2014b). Jenkins et al. (1989) estimated this leakage to comprise 13.4% of the inorganic CO₂ liberated by C₄ acid decarboxylation in the BSCs. Since essentially all the CO₂ fixed in BSCs passes through ribulose bisphosphate carboxylase/oxygenase (RubisCO), we applied the 13% leakage estimate on the rate of RubisCO carboxylation. As well, from the expected CO₂ and O₂ availabilities in BSCs, RubisCO carboxylase activity is expected to be 29.2-fold greater than its oxygenase activity in NADP-ME-type C₄ plants (Jenkins et al., 1989). Consequently, the CO₂ leakage and photorespiration reactions were constrained to match these literature estimates (Supplementary Table 2).

CO₂ assimilation by C₄ photosynthesis and allocation of reduced carbon was modeled by considering the interactions of day and night metabolism in a diel Flux Balance Analysis (FBA) model similar to other reported diel plant models (Cheung et al., 2014; de Oliveira Dal’Molin et al., 2018; Shaw and Cheung, 2018). In essence, the S. bicolor C₄ core metabolism model with interconnected MC and BSC sub-models (iTJC1414x2) was duplicated into representations of day and night. In the resulting diel model with four cellular sub-models (iTJC1414x4), the BSC day and night cells share a common pool that represents transitory starch (Figure 1), while MCs are assumed to not store starch. This is according to literature findings that starch is observed mainly in the BS chloroplasts of C₄ species like Z. mays and sorghum (Lunn and Furbank, 1997). Light is only available as an energy source for the day sub-models which contribute to the starch pool while the night models rely on starch degradation to fuel dark metabolism and respiration (Supplementary Table 2). In general, flux states were obtained based on the principle of light-limiting conditions. Given a fixed biosynthetic rate of a photosynthetic product, total light uptake fluxes were minimized. Alternatively, given fixed light uptake fluxes, the biosynthetic rate of a photosynthetic product was maximized. To simulate C₄ photosynthesis, only the MC is allowed to exchange CO₂ with the external environment and RubisCO is only active in the BSC (Supplementary Table 2). Net CO₂ uptake (V_netCO2) is defined as the difference between daily uptake of CO₂ (V_CO2uptake) and CO₂ loss at night (V_CO2export):

By linking the four sub-models (BSC-day, MC-day, BSC-night, and MC-night) through the shared starch pool, a single linear optimization procedure predicts fluxes that are in steady-state over the integrated day-night cycle.

Model analyses in this study assume cells are mature, thus no longer growing, and TAG and sucrose are the primary net carbon sinks. Although we modeled non-growing cells in this study, we also updated the iTJC1414x4 biomass composition constraints using published S. bicolor carbohydrate, lignin, amino acid, and lipid ratios (Supplementary Table 3) and confirmed that the model was functional if maximizing total biomass was used as the objective function. The fatty acid composition of TAG was defined as reported for wild-type sorghum leaves (Vanhercke et al., 2019a). Transgenic studies aimed at increasing TAG accumulation in vegetative tissues have regularly reported TAG compositions that are distinctly different from the wild-type (Xu and Shanklin, 2016; Vanhercke et al., 2019b). As described in section “Results,” changes in TAG fatty acid composition that have been observed in sorghum high-oil lines (Vanhercke et al., 2019a) do not majorly affect biosynthetic energy demands and were therefore not considered here. Similarly, we did not consider here the possibility that increased leaf TAG levels may in part be attributable to re-balancing fatty acids between membrane lipids and TAG. For example, in a recent study on high oil accumulating tobacco lines (Zhou et al., 2020) there was a dramatic decrease in levels of chloroplastic galactolipids in the transgenic line which likely benefited TAG production. Our model does not consider the potentially involved lipid trafficking mechanisms between chloroplast and cytosol/endoplasmic reticulum. To our knowledge these mechanisms are not yet fully elucidated. We verified that sucrose and precursors for TAG biosynthesis can be generated by canonical pathways. In short, the primary CO₂ fixation product in the CBB-cycle in the chloroplast, 3-PGA, can be directly transformed into pyruvate by a series of glycolytic reactions, which can then be transformed into acetyl-CoA by the plastid-localized pyruvate dehydrogenase complex (PDHp) (Bassham and Calvin, 1960; Li-Beisson et al., 2013). Alternatively, PGA can be reduced to triose phosphates (Bassham and Calvin, 1960), which are then transformed into transitory starch or sucrose for export (Raines, 2003).

For modeling carbon allocation scenarios, we consider here that the primary products of photo-assimilation are exported via phloem to other organs or that, due to an envisioned transgenic intervention to engineer TAG accumulation, a part of the assimilated carbon is diverted into oil droplets and stored in the leaf. Apart from changes in leaf TAG levels, the carbon allocation scenarios simulated here do not consider changes in leaf levels of other biomass constituents such as sugars, starch, or amino acids. While for a mature leaf, sucrose can be understood as the predominant exported photo-assimilate, it can be considered that amino acids are exported as well, which is supported by studies of phloem sap composition in various plants (Wilkinson and Douglas, 2003). We therefore defined photo-assimilate based on the composition of sucrose and amino acids found in phloem sap of barley leaves (Winter et al., 1992). The combined assimilate export is comprised of sucrose and 17 amino acids (Supplementary Table 4) with a sucrose to amino acid ratio of 94:6 (by weight). Also, while TAG can be considered to be the predominant product of the envisioned metabolic engineering effort, oleosin proteins are major components of oil droplets that are stored as intracellular particles (Ohlrogge and Chapman, 2011). Oleosins were estimated to comprise 1–5% of the oil droplet weight in maize embryos while the remainder is mostly TAG (Ting et al., 1996). Consequently, in carbon allocation simulations we define the deposition of oil droplets as biosynthesis of TAG and oleosin at a 95:5 ratio (by weight) (Supplementary Table 5). Here, the oleosin has an amino acid composition that matches Sesamum indicum oleosin-L (Tai et al., 2002; Supplementary Table 5). This protein was over-expressed in sorghum leaves in the study of Vanhercke et al. (2019a).

Note that nitrate reduction and ammonia assimilation into amino acids take place in the simulated assimilate export as well as in oleosin synthesis. Nevertheless, since nitrogen plays only a minor role, we will keep here the term “carbon allocation” for the model simulations of assimilate export, and TAG and oleosin deposition. For simulations on TAG accumulation, the objective function used was to maximize the TAG and oleosin accumulation in MC and BSC during the day. This allows TAG and oleosin biosynthesis to take place in day metabolism for both cell types. This procedure appears to be justified because when WRI1, DGAT, and oleosin were expressed under a constitutive promoter in sorghum, oil droplets were found to accumulate in both cell types (Vanhercke et al., 2019a). For simulations of carbon allocation scenarios, different carbon allocation ratios were enforced based on the equation

where P_oil is the percentage of net CO₂ uptake transformed into TAG and oleosin, and V_Aex is the assimilate export flux to the phloem, with the average photo-assimilate being 12.77162 mol carbon per mol. For any chosen carbon allocation ratio P_oil > 0%, the model can be optimized by maximization of the total oil droplet accumulation rate. Because steady-state day-night cycling is integrated into the model, the simulations predict fluxes in units of mol/m²/day, including the total accumulation rate of TAG.

Given a daily oil droplet accumulation rate (V_oilDroplet), one can determine the time period needed (t) in days for a full-grown leaf to accumulate a desired amount of TAG per dry weight:

where C_TAG is the desired TAG content in the leaf expressed as percent dry weight, L is the initial dry weight per leaf area which was set to 60 g/m², as reported for sorghum (Zhao et al., 2005). M_TAG is the molecular weight of TAG (866.811 g/mol) based on the molecule species composition as defined in the model biosynthesis reactions. Note that the reaction equation for the daily oil droplet accumulation rate (V_oilDroplet) is defined so that 1 mol TAG is produced (866.81 g) along with 45.55 g oleosin, so that oleosin is 5% of the total in weight (Supplementary Table 5).

We used our model to explicitly test two types of futile cycles that will be referred to as TAG cycling and FA cycling. TAG cycling is defined as the ratio of fluxes toward TAG synthesis and TAG storage in oil droplets (detailed equations in Supplementary Table 6):

Succinctly, the TAG cycling ratio indicates to what extent TAG is over-produced for the predicted rate of TAG storage. The excess TAG must be hydrolyzed into glycerol and free fatty acids. As will be seen from the model simulations (see section “Results”), under limited energy input (light) and with increasing TAG cycling there will be reduced flux into TAG storage. Like TAG cycling, FA cycling is the ratio of fluxes toward fatty acid synthesis compared to that which is used for TAG storage in oil droplets (detailed equations in Supplementary Table 6):

By defining FA cycling in this way, the over-produced portion of fatty acids cannot be utilized and is therefore subjected to degradation (β-oxidation). Products and energy cofactors generated by the degradation process are recycled.

The default configuration of iTJC1414x4 allows photon flux to freely distribute between the BSC and MC and between the linear and cyclic components of photosynthetic electron transport. To model the light absorption more realistically, we generated sorghum-specific light flux distribution settings from a model for cell-type specific electron transport in C₄ photosynthesis (Yin and Struik, 2018), which integrates multiple photosynthesis-related biophysical and biochemical characteristics (see Supplementary File 3). One predicted parameter, a_BS,M, defines the ratio at which total received light absorptance is partitioned between the BSC and MC (Yin and Struik, 2018). This ratio was applied to iTJC1414x4 as a linear dependency between the two light uptake fluxes:

The model by Yin and Struik (2018) also predicts the fraction of the absorbed light that drives linear electron transport (LET) for each of the two cell types. Here we designate these as f_LET,M and f_LET,BS for the MC and BSC, respectively. To add the f_LET parameters as numerical constraints to iTJC1414x4 we derived the following equation (for details see Supplementary File 3):

with v_PSI and v_PSII being the rates of the reactions representing photosystem I and II in iTJC1414x4, respectively (“PhotosystemImodchloroplast,” “PhotosystemIImodchloroplast”). Applying the energy budget parameters a_BS,M, f_LET,BS and f_LET,M to iTJC1414x4 by Equations 6, 7 will cause the photosynthetic supply of ATP and NADPH to be at fixed proportions in the BSC and MC and at a fixed overall ATP/NADPH supply ratio. If the photosynthetic supply is fixed and the biosynthetic output has a fixed composition, there can be an imbalance between supply and metabolic demands of ATP and NADPH. Such an imbalance causes one of the two energy cofactors to be overproduced, i.e., it will either accumulate or dissipate in a process that is not related to biosynthesis. In the context of iTJC1414x4, no accumulation of metabolites can materialize since FBA is a steady-state modeling approach. Instead, we characterize imbalanced scenarios in iTJC1414x4 by detecting net hydrolysis of ATP that does not drive biosynthesis and by detecting mitochondrial oxidative phosphorylation that leads to more ATP production at the expense of oxidation of photosynthetically generated reducing equivalents. Specifically, ATP surplus is quantified by maximizing an ATP-consuming dummy reaction in a secondary optimization (Supplementary Figure 2). The size of the flux through complex IV of the mitochondrial electron transport chain (“ComplexIVmodmitochondrion”) reveals the amount of photosynthetically produced reducing equivalents that are transferred back to oxygen. Both ATP surplus and NADPH surplus are quantified by a Flux Variability procedure as indicated in Supplementary Figure 2. The ATP surplus and NADPH surplus quantified this way reveals misalignment between photosynthetic supply of ATP and NADPH and biosynthetic demands. Tests showed that if the energy budget parameters are not applied, no ATP or NADPH surplus is detected. Furthermore, as detailed in Supplementary File 3, if sucrose is the biosynthetic product, values for a_BS,M, f_LET,BS and f_LET,M can be derived from the Yin and Struik model and applied to iTJC1414x4, leading to a balanced energy budget, i.e., no ATP or NADPH surplus is detected.

To set the photosynthetic assimilation rate in iTJC1414x4, we searched literature to define physiologically realistic values for daily photosynthetic rates of CO₂ fixation for a sun-lit sorghum leaf under non-stressed field conditions. While direct measurements of total CO₂ uptake have been reported for bioenergy grasses such as miscanthus and switchgrass (Dohleman et al., 2009), we could not find such data for sorghum. Reviewing literature for photosynthetic performance data specifically on S. bicolor, we found measurements of maximal photosynthetic CO₂ fixation rates at mid-day. We therefore first derived an approximation to convert maximal photosynthetic CO₂ fixation rates at mid-day to daily rates. Dohleman et al. (2009) measured leaf CO₂ uptake rates for miscanthus and switchgrass for upper canopy sunlit leaves every 2 h on multiple days during the growing season. By plotting these uptake rates against time of day, the authors proposed a geometric method for estimating the daily photosynthesis rate. Here, we estimate the time course plot can be approximated with a parabolic function that peaks at mid-day and can be integrated to determine the daily rate. Accordingly, we defined the integral of daily CO₂ fixation as 2/3 × (maximal photosynthetic CO₂ fixation) × daylength (Supplementary Figure 3) and used this relation to estimate daily CO₂ uptake rates for sorghum leaves based on mid-day photosynthetic rates. Around mid-day, under non-stressed conditions and at highest photosynthetic photon flux density, CO₂ uptake rates measured for 22 S. bicolor genotypes ranged from 23.5 to 44.8 μmol CO₂/m²/s, with an average of 38.2 μmol CO₂/m²/s (Peng et al., 1991). In another study, under similar conditions, a range of 43.3–56.9 μmol CO₂/m²/s was measured for 26 sorghum genotypes (Balota et al., 2008). Based on the two studies, we chose 40 μmol CO₂/m²/s as a conservative estimate for the maximal photosynthetic rate. Using the parabola approximation with a daylength of 14 h, the CO₂ uptake from the environment into the leaf during the daylight period is 1.344 mol CO₂/m²/day.

In addition to the photosynthetic rates, it is important to consider energetic losses associated with cellular maintenance activities. This refers, for example, to energetic expenditures for maintenance of membrane gradients or continuous turnover of biopolymers. In plant FBA models such energetic costs are often classified as unspecified ATP as well as NADPH consumption (Sweetlove et al., 2013). To mimic energetic burdens of cellular maintenance, each cellular sub-model of iTJC1414x4 includes a generic reaction to hydrolyze ATP (“GenericATPasemod”) and one to oxidize NADPH (“GenericNADPOxidasemod”). Based on the setup of other diel models of C₃ and crassulacean acid metabolism (CAM) leaf metabolism (Cheung et al., 2014; Shaw and Cheung, 2018), we assumed here the maintenance costs to be the same in the light and dark phases, and we adopted a 3:1 ratio between ATP and NADPH drain fluxes (Cheung et al., 2013). To estimate dark respiration during the daytime we used a light response curve of sorghum leaves shown in Figure 4 of Wang et al. (2019). From this figure we determined the maximal photosynthetic rate for adaxial illuminated leaves was 38 μmol CO₂/m²/s, which is close to the maximal photosynthetic rate of 40 μmol CO₂/m²/s we considered here. The light response curve shows a linear increase in photosynthetic response up to a light intensity of about 500 μmol photon/m²/s with a slope of 0.064 mol CO₂/mol photon and the dark respiration intersect of the linear was determined to be at −2.6 μmol CO₂/m²/s. In another study evaluating five sorghum genotypes under field-grown conditions, respiration rates of dark-adapted leaves were on average 1.4 μmol CO₂/m²/s (Li et al., 2021). We therefore consider the dark respiration intersect of 2.6 μmol CO₂/m²/s to represent a rather conservative estimate. Assuming this value is constant over the 14 h daytime, the daily CO₂ export rate is 0.131 mol CO₂/m²/day. The model was adjusted to this dark respiration rate using the generic maintenance reactions (“GenericATPasemod,” “GenericNADPOxidasemod”), as detailed in Supplementary Figure 4.

To give a summarizing overview of the flux solution space in iTJC1414x4, we defined a set of network projections that map selected groups of iTJC1414x4 reactions onto a lumped reaction network as described before (Schwender and Hay, 2012). In short, a bilevel optimization as in Flux Variability Analysis is applied. However, the objective in the secondary optimizations is not to minimize/maximize a single reaction but a linear combination of reactions so that in effect the minimum and maximum bounds for the total flux between lumped metabolite pools is revealed. For example, to aggregate all conversions between oxaloacetate and malate by all malate dehydrogenase enzyme isoforms in the BSC-day sub-model, and taking into consideration the directionality of each isoform, the objective function is to maximize/minimize the following sum: “bsDayMALATEDEHYDROGENASENADPRXNchloroplast – bsDayMalateDHPeroxisome – bsDayMalateDH – bsDayMALATEDEHRXNmitochondrion.” All used summation terms of this type are listed in Supplementary Table 7.

Independent of specific metabolic pathways, any oxygenic photosynthetic process converting CO₂ and water into a carbon-reduced product and oxygen can be defined by its chemical and redox balance. Accordingly, we derived a general calculation scheme to derive the chemical balance of formation and the average oxidation state of carbon and applied it to various carbon-reduced compounds (Supplementary File 3). In addition to the chemical balances, we derived the overall balances of formation from CO₂ and water by plant specific biochemical pathways. For this purpose, we defined a minimal stoichiometric model representing canonical pathways of the CBB cycle of C₃ photosynthesis, the photorespiratory cycle, and sucrose, fatty acid, and TAG synthesis. This one-compartment model can only exchange sucrose, TAG, CO₂, oxygen, and water with the environment. Using the software tool METATOOL (Schuster et al., 1999, 2000), the mass balances of all possible conversions of external metabolites into sucrose or TAG can be determined as Elemental Flux Modes (see Supplementary File 3 for details and Supplementary File 4 for METATOOL input files). Reaction coupling with dummy species “PSATP” and “PSH2” reveals the biosynthetic demands in ATP and reducing equivalents.

The sorghum model iTJC1414x4 is comprised of four sub-models representing leaf MC and BSCs under day and night conditions (Figure 1). Each of the sub-models has over 1,000 reactions, 750 metabolites, and 7 intracellular compartments (Table 1). The day and night parts of the model are connected by a transitory starch pool. The fraction of assimilated carbon that is converted into (transitory) starch is consumed at the same rate in the dark model by biosynthetic or respiratory processes (Figure 1). Note that the simulated flux distributions are to be understood as integrals over the entire day. iTJC1414x4 was constructed using the Z. mays metabolic reconstruction of leaf C₄ metabolism iEB2140x2 (Bogart and Myers, 2016) as a template. We added 380 reactions and modified over 100 reactions per sub-model, thereby correcting mostly for subcellular localization or cofactor specificity of fatty acid biosynthetic reactions. Effort at manual curation was particularly focused at lipid metabolism, which was done mostly in reference to plant-specific reactions as outlined in Li-Beisson et al. (2013). In short, chloroplast localized fatty acid biosynthesis was defined to produce palmitic acid (C16:0, a 16-carbon fatty acid with zero double bonds), stearic acid (C18:0) and oleic acid (C18:1), which are exported to the cytosol as free fatty acids (Li-Beisson et al., 2013). Cytosol/endoplasmic reticulum-localized lipid metabolism was set up as a sub-network specifying the biosynthesis of PC, DAG, and TAG with distinct molecular species that are composed of five acyl chain types, which are palmitate, stearate, oleate, and the polyunsaturated linoleic acid (C18:2) and linolenic acid (C18:3). The linoleate and linolenate species derive by desaturation reactions of acyl chains bound to PC (Li-Beisson et al., 2013). Degradation of TAG and fatty acids is given by TAG lipase reactions, by core activities of β-oxidation of the modeled fatty acid species (Li-Beisson et al., 2013) and by functions of the glyoxylate cycle (Graham, 2008). Besides lipid components, the model is able to synthesize protein, cellulose, hemicellulose, soluble sugars, and lignin (Supplementary Table 3).

Of the 2,133 Z. mays genes associated with reactions in the iEB2140 single-cell model, 85% of the genes had an S. bicolor ortholog and 92% of the reactions had S. bicolor orthologs for most of the associated maize genes (Supplementary Figures 1B,C). Together, these support the core iEB2140x2 model being transferable from maize to sorghum. To integrate information on cell type specific metabolism, we analyzed differential gene expression data for S. bicolor MCs and BSCs, that was generated by Döring et al. (2016) using Illumina sequencing and SuperSage analysis. Of the 1,414 sorghum genes associated with reactions in iTJC1414x4 (Table 1), 52 were found to be differentially expressed if either sequencing method was used. To be consistent with this data, we constrained nine MC and three BSC reactions to be inactive because all of their associated genes were found to be strongly preferentially expressed in the opposing cell (Supplementary Table 8).

Prior to using iTJC1414x4 to model TAG biosynthesis, we assessed energy content and biosynthetic costs of photo-assimilation independent from specific biochemical pathways. For a process of oxygenic photosynthesis that converts CO₂ and water into a reduced product, the physicochemical limits are given by the chemical balance. Table 2 provides a list of reduced theoretical products, including organic acids, carbohydrates, lipids, and fossil fuel derived reduced carbon compounds, sorted by carbon reduction states and physicochemical demands in reducing equivalents. Five of the TAG molecular species listed are representative of the major fatty acids that are commonly accumulated in TAGs in crop plants, in particular palmitate, stearate, oleate, linoleate, and linolenate (Table 2). The amounts of required reducing equivalents and energy contents listed in Table 2 are derived from the elemental composition of each listed compound (for details see Supplementary File 3). Therefore, by averaging the listed TAG species, the properties of any molecular species composed of the mentioned five major fatty acids can be derived. On a per weight basis, the TAG species are on average 2.39-times more energy dense than sucrose (ratio ranging from 2.37 to 2.44; Table 2). However, since this study aims at analyzing partitioning of photo-assimilate into different reduced carbon compounds, it is important to also consider the molar quantity of assimilated CO₂ (i.e., carbon) as a basis for comparing energy content and biosynthetic demands. If considered on a per carbon basis, sorghum wild-type leaf TAG only has 1.29-fold higher energy content than sucrose, while high-oil leaf TAG is 1.31-fold higher and different TAG molecular species are on average 1.3-fold higher in energy content than sucrose (ranging from 1.27 to 1.34; Table 2). As expected, because the biosynthetic process can be understood as the transfer of reducing equivalents onto carbon, the differences in per carbon energy content and per carbon reducing equivalents demands are similar.

In addition, per CO₂ incorporated, sorghum wild-type TAG synthesis requires 1.37-fold the amount of reducing equivalents as sucrose (values for other TAG compositions range from 1.35 to 1.43; Table 2). Since the demands for reducing equivalents listed in Table 2 derive from chemical balances, they represent a chemical minimum of energetic investments. As detailed in Supplementary File 3, the overall biosynthetic demands in reducing equivalents will raise above this minimum if metabolic pathways contain steps where electrons are transferred from pathway intermediates back to oxygen. For synthesis of TAG with sorghum fatty acid composition, the chemical minimum demand is 2.74 reducing equivalents per carbon (Table 2), while biosynthetic pathways demand 2.98 reducing equivalents. As detailed in Supplementary File 3, the difference is accounted for by the transfer of electrons to oxygen at the fatty acid desaturase steps. Relative to the total biosynthetic demands, 8.1% of reducing equivalents are transferred back to oxygen.

To analyze possible photo-assimilatory pathways for sucrose and TAG synthesis, we first used the single-cell model iTJC1414 to simulate three scenarios of sucrose or TAG formation under limiting light conditions (Supplementary File 1). Table 3 shows the simulations in the case of C₃ photosynthesis without photorespiration being active (see Supplementary Table 9 for consideration of photorespiration). In each case the overall demands in CO₂ and H₂O as well as O₂ production are shown in Table 3 and are as expected by the chemical balances (see Supplementary File 3). For each scenario, the optimal flux solution space was further characterized by Flux Variability Analysis with subsequent classification of flux bounds into flux variability types as described before (Hay and Schwender, 2011). For optimal sucrose production, 48 reactions must be active (Table 3, scenario 1, “essential reactions”), while, in between alternative optimal flux states, 25 reactions can have zero or non-zero flux values (“non-essential reactions”). For TAG synthesis, a much larger number of reactions are non-essential than are essential (Table 3, scenario 2). While this might not be unexpected due to the complexity of the TAG biosynthetic network, inspection of the 618 non-essential reactions in scenario 2 revealed that the optimum TAG solution space allows many reactions to be active that should not be required for photo-assimilatory TAG synthesis (Supplementary Table 10). For example, invertase (“RXN1461”) was active in what appears to constitute a cycle of synthesis and degradation of sucrose. The Oxidative Pentose Phosphate Pathway (OPPP) can be active (e.g., “GLU6PDEHYDROGRXNChloroplast”) although in the photoautotrophic context, activity of OPPP and CBB cycle constitute a futile cycle with net ATP consumption (Sharkey and Weise, 2016). Altogether, the optimum solution space seems to feature numerous alternative solutions with ATP consuming futile cycles across the network that dissipate an apparent over-supply in ATP. To confirm this suspicion, we showed that this ATP consuming activity can be aggregated by maximizing a generic ATP consuming reaction, which we added to the network (“ATPsurplus”). If fixed to its maximal possible value under optimality, subsequent Flux Variability Analysis resulted in a much smaller solution space in terms of number of active reactions and futile cycles were no longer detectable (Table 3, scenario 3). In addition to the observations made in Table 3, simulations that include photorespiratory activity led to essentially the same conclusions on the occurrence of ATP surplus for TAG synthesis (Supplementary Table 9). The explanation for an over-supply (“surplus”) in ATP is found in the rigidity of the photosynthetic light-driven production of ATP and NADPH. As we will explore in more detail below, the lowest possible ratio at which ATP and NADPH can be supplied by the photosynthetic apparatus is at ATP/NADPH = 1.5. If the overall metabolic demands ratio for a particular biosynthetic product is below 1.5, stoichiometrically feasible optimal flux solutions will require the dissipation of surplus ATP. This model finding has the equivalence of an ATP/NADPH imbalance in the real process. The photosynthetic apparatus can be considered to deliver ATP and NADPH at fixed proportions, which should precisely match the proportions of the metabolic demands (Kramer and Evans, 2011). Chloroplasts are said to have limited capacity to adapt to variation in the relative demands in ATP and NADPH and a misalignment can stall the overall process (Kramer and Evans, 2011).

To further explain the results on an ATP/NADPH imbalance, we assessed the metabolic demand independent of supply (Figure 2A). We defined a minimal stoichiometric model for CO₂ assimilation, consisting of 44 reactions representing plant canonical pathways of the CBB cycle of C₃ photosynthesis (Sharkey and Weise, 2012), the photorespiratory cycle (Peterhansel et al., 2010), sucrose synthesis (Stitt et al., 1987), the generation of acetyl-CoA directly from 3-PGA (Joyard et al., 2010), and a net equation for chloroplast fatty acid synthesis. To obtain the latter, the net balance from acetyl-CoA and energy cofactors was derived according to generic plant pathways (Li-Beisson et al., 2013) by assembling the stoichiometries of 47 reaction steps. Moreover, the minimal model is a subset of iTJC1414 (see Supplementary File 3 for details). Applying Elementary Flux Modes analysis to the minimal model, the net balance for the conversion of CO₂ and water to oxygen and sucrose or TAG can be obtained along with the energy cofactor requirements. Distinct flux modes were obtained for sucrose and TAG production; in each case one with the photorespiratory pathway being active and one without (Supplementary File 3). Figure 2A summarizes the resulting scenarios for sucrose and TAG synthesis if photorespiration is not active. Flux rates and energy cofactor demands are shown relative to the net assimilation of one mol of CO₂. It is noticeable that, relative to sucrose production, CBB cycle fluxes are expected to be substantially higher if TAG is the product (Figure 2A). This is because, at the PDHp step, one-third of the carbon that is processed from 3-PGA toward fatty acid biosynthesis is released as CO₂ and then re-fixed by RubisCO. Per CO₂ assimilated into TAG, ATP-, and NADPH-demands are higher by 1.3- and 1.49-fold, respectively (Figure 2A). Similar increases are observed if photorespiration is active (see Supplementary File 3).

Considering ATP and NADPH relative to each other, the ATP/NADPH demand ratio for TAG synthesis is lower than that for sucrose (Figure 2A). The basic concept of an energy cofactor budget is outlined by the supply/demand scheme shown in Figure 2B. If the ATP/NADPH supply ratio is higher than the demand ratio for a biosynthetic product, then there remains an ATP surplus which must dissipate in a process other than biosynthesis. If the ATP/NADPH supply ratio is lower than the demand ratio, then there will be a surplus in reducing equivalents, which can be dissipated by activity of oxidative phosphorylation or by another oxidative process (Figure 2B). Based on the stoichiometric setup of light reactions in iTJC1414, the ATP/NADPH ratio supplied by the LET chain is 1.5 (Supplementary File 3). Note that there are some uncertainties about the true ratio of ATP/NADPH output from LET found in plants, but most literature sources determine the ratio as being close to 1.5 (Noctor and Foyer, 1998; Allen, 2003; Kramer and Evans, 2011). Thus, considering LET to be active with only a very small additional ATP contribution by CET, the photosynthetic supply is sufficient to match the demands of 1.54 ATP per NADPH for sucrose synthesis (Figure 2B). Supposing that the ATP and NADPH supply is tailored this way for sucrose synthesis with a rigid supply ratio of 1.54, the demand ratio of 1.34 for TAG (Figure 2A) will cause an ATP surplus. If NADPH supply and demand are aligned, the ATP surplus amounts to 13% of the production [(1.54−1.34)/1.54]. Based on the principle outlined in Figure 2B we modeled and assessed the effects of a rigid energy cofactor supply for the two cell types of C₄ metabolism in iTJC1414x4. This allowed to define a reference scenario for CO₂ assimilation into sucrose, where photosynthetic ATP and NADPH production are in balance with metabolic demands in both cell types. Based on this reference state, the energy balance in TAG biosynthesis scenarios can be evaluated.

In C₄ photosynthesis, different metabolic activities must be coordinated across the two cell types. Before applying sorghum-specific constraints to the distribution of light energy among the two cell types in iTJC1414x4, we first verified the basic configuration where the light uptake fluxes and the relative activities of LET and CET were freely adjustable. If sucrose synthesis was simulated by minimization of total light uptake, it was found that with the overall photosynthetic energy cofactor supply is 3.647 ATP and 2.431 NADPH per mol CO₂ assimilated (Table 4, scenario 1, sucrose). Since sucrose synthesis requires 2 NADPH per CO₂ assimilated (Table 2), there is an oversupply of 0.431 NADPH. Further inspection of the flux scenario revealed that this amount of reducing equivalents is consumed by mitochondrial oxidative phosphorylation, which generates additional ATP, resulting in an adjusted overall cofactor supply of 5.083 ATP and 2 NADPH per mol CO₂ (Table 4, scenario 1, sucrose). This adjusted supply meets the above established general expectation for the basic photosynthetic demands of 3.083 ATP and 2 NADPH per mol CO₂ (Figure 2A), after it is also considered that the operation of the C₄ cycle incurs an additional 2 ATP/CO₂ at the pyruvate-phosphate dikinase (PPDK) step (Figure 1). Potentially, the additional ATP could be generated by CET. However, the basic sucrose synthesis flux scenario predicts CET to be inactive in both cell types (Supplementary Table 7). Flux modes analysis of the isolated ATP and NADPH generating system shows that ATP production via combined LET and oxidative phosphorylation activity has higher quantum use efficiency than ATP production via CET (Supplementary File 3).

To develop a sorghum-specific configuration of light energy distribution, we assumed that a typical mature sorghum leaf must have a photosynthetic energy cofactor supply that is well balanced with the metabolic demands of sucrose synthesis. To define a typical sorghum light distribution pattern with energy balance for sucrose synthesis, we relied on a model of cell-type specific electron transport in C₄ photosynthesis by Yin and Struik (2018). The analytical model predicts the energy budgets for the major sub-types of C₄ photosynthesis based on the cofactor demands of the CBB cycle, taking sub-type specific variations in the C₄ cycles into consideration and by estimating photosynthetic production capacity based on leaf morphology and the distribution of chlorophyll and photosystems I and II between BSCs and MCs. We used the Yin and Struik (2018) model to predict three parameters to be incorporated into our sorghum FBA model: the ratio of light absorptance between the BSC and MC (a_BS,M) and, for the BSC and MC, respectively, the fraction of absorbed light that drives LET (f_LETBS, f_LETM). With sorghum specific parameters and defining sucrose synthesis as the energy consuming metabolic activity, the Yin and Struik (2018) model predicts a balanced C₄ energy budget at a_BS,M = 0.398, f_LETBS = 0.032, and f_LETM = 0.933 (see Supplementary File 3 for details). After the three energy budget parameters were applied to iTJC1414x4 (Equations 6, 7), no ATP surplus or oxidative phosphorylation adjustment was detected for the two cell types (Table 4, scenario 2, sucrose). If, instead of sucrose, TAG is the product in the BSC or the MC, then an ATP surplus of 2.044 or 1.097 ATP per CO₂ fixed can be detected, respectively (Table 4, scenario 2, TAG). This means that using these energy budget parameters, a state was defined in iTJC1414x4 where cell-type specific photosynthetic ATP and NADPH supplies are in balance with cell-type specific demands for sucrose synthesis. Figure 3 illustrates the functioning of C₄ metabolism under this condition. With the 3-PGA/triose phosphate shuttle operating, 37% of the reduction of 3-PGA to GAP takes place in MC, meaning that 37% of the ATP and NADPH demand at this step in the CBB cycle is assigned to the MC. The malate/pyruvate shuttle mechanism, which transfers both CO₂ and reducing equivalents, is dominant. However, 29.4% of the CO₂ is transported via the aspartate/alanine mechanism, which transports CO₂ only (Figure 3). This means that only 70.6% of the possible capacity to move NADPH with the C₄ shuttle is used. Overall, the imposition of the energy budget parameters appears to create a realistic reference state with balanced energy budget for sucrose synthesis.

In Figure 3 we consider sucrose biosynthetic costs only. To simulate carbon assimilation more realistically, we need to calibrate the model to the overall photosynthetic capacity of a sorghum leaf and incorporate estimates on additional energetic costs due to photorespiration, CO₂ leaking back from BSC to MC and costs related to respiratory activity. As described in section “Materials and Methods,” to set the photosynthetic capacity we derived a typical daily CO₂ assimilation rate based on a survey of literature-reported maximal photosynthetic CO₂ fixation rates at mid-day, daily time course integrals as well as measured leaf dark respiration rates Accordingly, we estimate the total daily CO₂ assimilation rate for sucrose synthesis in a mature sorghum leaf to be 1.344 mol CO₂/m²/day, with a respiratory loss of 0.131 mol CO₂/m²/day. We next used these two rates in an iterative cycle that is outlined in Supplementary Figure 4 to set the model’s light uptake and maintenance fluxes as well as to adjust the energy budget parameters, beginning with the values for a_BS,M, f_LET,BS, and f_LET,M obtained above from the Yin and Struik model (Table 4, scenario 2).

In short, we first determined the minimal light levels for which the combined net daily CO₂ assimilation and respiratory loss (1.475 mol CO₂/m²/day) can be converted into sucrose. Second, while keeping these light fluxes constant, we raised the day-time maintenance cost so that the combined CO₂ rate was reduced back to the net daily CO₂ assimilation rate of 1.344 mol CO₂/m²/day. Assuming that maintenance respiratory loss during the 14 h day-time also applies to the night-time (10 h), the night-time ATP and NADPH drain reactions were set to 10/14 of the daytime rates. Third, to obtain a mature leaf reference state, the photo-assimilate export rate was maximized. The iterative cycle was repeated by adjusting f_LET,M until, across the two cell types, no ATP surplus or supplemental ATP production by mitochondrial oxidative phosphorylation could be detected. This means that at the final state the photosynthetic production of ATP and NADPH was balanced with the overall demands for synthesis of the exported photo-assimilate. In result, we obtained light uptake fluxes of 6.30 and 15.82 mol photons/m²/day for the BSC and MC, respectively. The daytime ATP and NADPH drain reactions were determined to be 0.369 and 0.123 mol/m²/day, respectively. The corresponding nighttime values for ATP and NADPH drain reactions were 0.264 and 0.088 mol/m²/day, respectively. These settings for nighttime maintenance cause a nighttime respiratory CO₂ loss that is 13.4% of the CO₂ uptake flux. Nighttime respiratory loss of non-stressed S. bicolor plants has been estimated before to be between 10 and 16% of the daily photosynthesis rate (Hodges et al., 1979). We therefore consider the respiratory loss of 13.4% to be realistic. The balanced C₄ energy budget was obtained after an adjustment of the value for f_LET,M from 0.933 (Table 4, scenario 2, sucrose) to 0.895, while a_BS,M and f_LET,BS were kept at their formerly determined values (Table 4, scenario 3, assimilate export). The predicted absorbed quantum yield for CO₂ uptake is 0.060 mol CO₂ / mol photons (Table 4, scenario 3, assimilate export), which is similar to the measured value for sorghum leaves of 0.061 mol CO₂ / mol photons (Ehleringer and Pearcy, 1983). Finally, this scenario was taken to define a mature leaf reference state with balanced energy budget based on which carbon partitioning between assimilate export and oil droplet deposition can be further characterized. All constraints that define this state are listed in Supplementary Table 2 and Figure 1 highlights the constraints related to light flux, photorespiration, CO₂ leakage and maintenance cost.

By default, TAG biosynthesis and accumulation can take place in all four sub-models of iTJC1414x4. Given the above reference state, we tested the efficiency of conversion of CO₂ into TAG in all four sub-models (Supplementary Figure 5). Overall, TAG production can take place in the BSC- and MC-day models at the same rate. In the dark models maximal TAG production is reduced by less than 1%, relative to the maximal day rate (Supplementary Figure 5). If energy budget parameters are applied (Table 4), the photosynthetic energy budgets are balanced for sucrose synthesis, while for TAG synthesis, a substantial ATP surplus is detected. For scenario 2 (Table 4), the ATP surplus values of 2.044 and 1.097 for TAG production in the BSC and MC, respectively, amount to 27 and 14.5% of the photosynthetic ATP production. For scenario 3 (Table 4) these percentages are 17.3 and 8.5%, respectively. ATP surplus under TAG synthesis does not occur when light flux distribution is freely adjustable (Table 4, scenario 1), which demonstrates that the ATP surplus is a consequence of application of the energy budget parameters. The model prediction of ATP surplus suggests that, unless there are mechanisms to flexibly adjust light energy input or to discharge the surplus in ATP in another process, high rates of TAG synthesis in a leaf might be stalled due to energy cofactor imbalances. To further explore model predictions of TAG biosynthesis, a scenario of daytime TAG accumulation in BSCs is further detailed in Figure 4. As seen for sucrose synthesis (Figure 3), CO₂ from the environment is moved from MC to BSC via the malate/pyruvate shuttle mechanism with minor participation of the aspartate/alanine shuttle. In contrast to the sucrose synthesis scheme, a part of the pyruvate produced by NADP-ME in the BSC serves as a precursor to acetyl-CoA and fatty acid biosynthesis (Figure 4). The same amount of pyruvate that is removed from the pyruvate/malate shuttle cycle is replenished in MCs by conversion of 3-PGA to phosphoenolpyruvate. Altogether, a scheme emerges by which pyruvate destined to fatty acid synthesis is generated from 3-PGA via phosphoenolpyruvate under participation of the 3-PGA/triose phosphate shuttle and the C₄ shuttles, while omitting the pyruvate kinase and PPDK steps (Figure 4, blue trace). The scheme is cost-saving relative to direct conversion of 3-PGA to pyruvate in BSC, since generating phosphoenolpyruvate in the C₄ cycle with PPDK comes at a net cost of 2 ATP. However, the overall savings per unit pyruvate entering fatty acid synthesis is only one ATP since using pyruvate kinase in the BSC would gain one ATP (Figure 4).

The amount of TAG that is accumulated in a leaf will depend on the fraction of photo-assimilate that can be re-directed into oil droplet biosynthesis by a metabolic engineering approach and by the number of days the accumulation can take place before the material is harvested. Imposition of a range of carbon partitioning ratios is demonstrated in Figures 5A,B. Carbon allocation is shifted between photo-assimilate export and oil droplet deposition while starch is synthesized at a constant rate of 4.8 g/m²/day (Figure 5A), which reflects a constant loss of daily net carbon assimilation by nighttime respiration. As carbon allocation toward oil droplets increased from 0 to 100% at the expense of photo-assimilate synthesis, the TAG mass flux as well as the carbon molar flux into TAG deposition did not increase linearly, but with a slightly downward bending trend (Figure 5). This indicates that net carbon fixation decreases as TAG deposition increases. Supplementary Figure 6 confirms that, as carbon allocation to TAG increased, the simulated net CO₂ uptake rate decreased from 1.32 to 0.97 mol carbon/m²/day, while the energy inputs (photon uptakes) remain constant. This means that, from the same light energy inputs, 72.4% of the amount of CO₂ that is fixed into sucrose can be fixed into TAG. Given a predicted daily TAG accumulation rate, the number of days until a desired amount of TAG per dry weight is reached can be calculated based on Equation 3 (see section “Materials and Methods”). Figure 5B shows the effect of change in carbon allocation on the predicted number of days until 1, 5, or 20% TAG (w/dw) are reached in mature leaves. For example, if 5% of net carbon assimilation was allocated to oil droplets, the time taken until a 20% TAG level is reached was predicted to be 18 days (Figure 5B). If 20% of net CO₂ assimilation was allocated to oil droplets, then this level could be reached within 5 days. Altogether, it appears that if only a small shift in carbon allocation from assimilate toward oil droplet production can be achieved, substantial TAG accumulation should be possible within a month.

Leaf lipase activities (e.g., TAG hydrolysis) and the capacity for oxidative degradation of fatty acids by peroxisomal β-oxidation can result in futile cycles of synthesis and degradation which diminish the rate of net TAG accumulation. In iTJC1414x4, cyclic synthesis and degradation of fatty acids and of TAG can be imposed onto the model by the numerical constraints shown in Equations 4, 5. We define cycling as the multiplication factor by which the rate of synthesis of fatty acids or TAG exceeds the rate at which fatty acids or TAG itself are deposited in TAG. Fivefold cycling of fatty acids (Equation 5), for example, means that the de novo synthesis of fatty acids is five times the amount of fatty acid that is processed and stored in TAG. All fatty acids that are produced but not stored can only enter the peroxisomal process of β-oxidation. Besides fatty acid cycling, TAG cycling means that newly synthesized TAG is hydrolyzed into glycerol and free fatty acids, which then can be re-utilized to synthesize TAG. This requires glycerol and free fatty acids to be activated to glycerol phosphate and CoA esters, both of which are ATP-driven processes (reactions “GlycerolKinase,” “PalmitoylCoASynthesis,” “OleoylCoASynthesis,” “LinoleoylCoASynthesis,” and “LinolenoylCoASynthesis”). Figure 4 compares a scenario of TAG synthesis from CO₂ without cycling to a scenario where fivefold FA cycling takes place. At fivefold cycling, fatty acid biosynthesis is increased by 2.7-fold while the net CO₂ uptake and TAG deposition rates decrease by 62 and 53%, respectively. It is also notable that the CBB cycle flux is predicted to markedly increase under the fivefold cycling condition (Figure 4). Degradation of fatty acids by β-oxidation transfers electrons onto oxygen and NAD⁺, and the produced NADH can contribute to the overall balance in reducing equivalents. Carbon from β-oxidation is predicted to be recovered as malate, and in the pyruvate/malate shuttle mechanism its production makes up for the reduced synthesis and transport of malate from MCs to BSCs (Figure 4).

Assuming all net CO₂ was allocated toward oil droplets, we simulated the impact on the rate of TAG storage by levels of increasing futile cycling (Figures 5C,D). As seen in Figure 5C, as TAG cycling was increased up to about 50-fold, the rate of TAG deposition decreased only slightly. Notably, fatty acid cycling more strongly affects TAG deposition than TAG cycling. As can be deduced from the graphs in Figure 5C, a 50% reduction of the maximal rate of TAG deposition is obtained for about fivefold cycling of fatty acids, while about 100-fold TAG cycling is needed to reduce TAG deposition in the same way. This difference can be in part explained by the difference in re-synthesis cost. Re-synthesis of oleoyl-CoA from acetyl-CoA requires 9 ATP and 17 reducing equivalents, while re-cycling of one free oleic acid into TAG requires activation of the free fatty acid to the CoA ester by acetyl-CoA synthetase (EC 6.2.1.1) in coordination with adenylate kinase (EC 2.7.4.3), which comes at a net cost of 2 ATP. The effect of cycling on the daily net TAG deposition rate can be used to predict TAG accumulation over time (Equation 3). The time needed to reach a 20% TAG target increases as FA or TAG cycling increase (Figure 5D). These scenarios were considered again for a 100% allocation of CO₂ net assimilation to oil droplets. As before, the time to reach target levels of TAG was much less substantially affected by TAG cycling than by FA cycling (Figure 5D).

Both FA and TAG futile cycling can be imposed in combination with a carbon partitioning scenario. Figure 6 shows a factorial plot for simulating the three variables for up to 250-fold cycling for fatty acids and TAG. The time duration to reach 20% leaf TAG is shown in color coding. In agreement with the analysis in Figures 5, 6 shows how FA cycling has generally more substantial effects than TAG cycling across a range of carbon partitioning ratios. Although we show up to 250-fold futile cycling, model simulations with high cycling rates might be unlikely to occur in reality due to physiological limitations in flux capacity. To estimate limitations on fatty acid and TAG cycling fluxes, we predicted fatty acid biosynthesis rates when total leaf biomass is produced (Supplementary Table 3) under conditions of light flux, photorespiration, CO₂ leakage and maintenance cost as defined above for the mature leaf reference state (Supplementary Table 2). This simulation predicted that fatty acids, while contributing to biomass formation, are newly formed at a rate of 0.11 mol/m²/day, which is equivalent to 12% of the net CO₂ fixation rate. We now supposed that 10 times that rate (1.1 mol/m²/day) might define the upper limit in flux capacity for growing tissue or for mature tissue engineered to accumulate TAG. Figure 7 shows imposition of this flux capacity limit for discrete slices of Figure 6 at 10, 5, and 2% net carbon allocation to oil droplets. Simulation results where TAG and fatty acid biosynthesis operate below the upper limits are demarcated. For example, for 5% carbon allocation to oil droplets, maximal about 20-fold TAG and 25-fold FA cycling are feasible and within the demarcated area, the time to reach 20% TAG per leaf dry weight varies between 18 and 28 days (Figure 7B).

The simulations on carbon allocation to TAG (Figures 5–7) were based on rates that can be expected when S. bicolor leaves receive full sunlight. Realistically, there are weather dependent reductions in sunlight, effects of temperature and draft conditions, effects of shading within the canopy, as well as reductions in metabolic capacity due to onset of senescence. Of these conditions we decided here to briefly explore the effect of shading, without detailed modeling of canopy effects. At 1 m depths into the canopy, Z. mays leaves have been measured to receive roughly 70% less light than fully exposed leaves (Dohleman and Long, 2009). Consequently, in our shade simulations, we constrained the MC and BSC to receive 70% less light, thus to reduce the total influx for light from 22.12 to 6.63 mol photons/m²/day. With these shade conditions, we simulated how the rate of TAG storage was affected if 0–100% of net CO₂ assimilation was allocated to oil droplets, and there was up to 1,000-fold futile cycle activity. At 100% carbon allocation to oil droplets and without futile cycling, TAG could accumulate at 1.35 g/m²/day and reach the 20% yield target in less than 2 weeks (Supplementary Figure 7). If only 30% of net carbon assimilation is allocated to oil droplets or if there is 200-fold TAG cycling, then 20% TAG (w/dw) would be reached within a month.