to encode GAPDH enzymes (per PPDB annotations and the presence of Gp_dh_N and
Gp_dh_C domains), three (GRMZM2G162529, GRMZM2G463280, GRMZM2G026024) appeared
encode to encode genuine phosphoribulokinases (per PPDB annotations and the
presence of PRK domains), and two appeared to encode CP12-type regulatory
proteins, with no obvious evidence for any individual protein sharing more
than one of these roles. The regulatory role of CP12 does involve forming a
complex with PRK and GAPDH, but this reduces, rather than enhancing or
enabling, their individual activities [126]. We removed the PRK associations
of the GAPDH and CP12 genes from our model. PPDB assigned these three GAPDH
genes to a plastidic location based on experimental evidence, so we associated
them with those reactions exclusively (removing associations with the
cytosolic instances of EC 1.2.1.13 and/or EC 1.2.1.12.)

##  10\. Biomass equation

We developed a biomass equation following that used in [21]. Our calculations
are based on supplementary file S4 of that paper222Specifically,
journal.pone.0021784.s004.xls, as downloaded from the PLoS One web site 20
November 2013, in particular sheet 2, ‘Biomass_rxn’.

That sheet derives a biomass equation corresponding to the production of one
gram of plant dry weight, based on literature data on biomass composition; the
description is divided into subreactions forming (e.g.) ‘nitrogenous
compounds’, ‘lignin’, etc., which then participate in an overall biomass
reaction.) The units of the stoichiometric coefficients are mmol.

We have adopted most of the biomass composition assumptions of Saha et al
wholesale, with gratitude for their efforts in compiling this data from the
literature. However, we have made some minor adjustments, resulting in a
different overall stoichiometry for biomass production.

###  10.1. Fatty acids

Saha et al represent the total lipid/fatty acid contribution to biomass as a
pool of triglycerides in proportions apparently based on a maize oil
measurement and thus probably reflective of seed triglyceride composition.

We substitute measurements of the fatty acid content of mature maize leaf
membrane lipids [127] and write a biomass sub-reaction which consumes the
relevant free fatty acids (rather than their derivatives in the form of
triacylglycerols, membrane lipids, etc.,) as shown in Table 1.

Fatty acid | CornCyc compound | mol. wt. (g/mol) | mole fraction
---|---|---|---
palmitic | PALMITATE | 255.42 | 0.104
palmitoleic | CPD-9245 | 253.4 | 0.056
stearic | STEARIC_ACID | 283.47 | 0.011
oleic | OLEATE_CPD | 281.46 | 0.044
linoleic | LINOLEIC_ACID | 279.44 | 0.132
linolenic | LINOLENIC_ACID | 277.43 | 0.646
Table 1. Fatty acid proportions in biomass.

Weighting the molecular weights by the mole fractions, we find one mole of
fatty acid in appropriate proportions weighs 272.4 g. Dividing the mole
fractions by the overall molar weight and multiplying coefficients by 1000 to
convert to millimoles, we arrive at the final equation:

> 0.382 PALMITATE \+ 0.206 CPD-9245 \+ 0.04 STEARIC_ACID \+ 0.162 OLEATE_CPD
> \+ 0.485 LINOLEIC_ACID \+ 2.372 LINOLENIC_ACID = fatty_acids_biomass

where the left-hand side represents 1 g.

Fractions add to less than 1.0 because we ignore trace (mole fraction $\leq$
0.01) amounts of C14:0 and C20:0 fatty acids. Note that the leaf fatty acid
composition is known to change along the developmental gradient, so specifying
any single composition is an approximation; see [128].

###  10.2. Hemicellulose

We adopted the hemicellulose production reaction as is, using the species
added to the model for this purpose, ‘polysaccharide_[sugar]_unit’. The
resulting equation is:

> 0.548 polysaccharide_arabinose_unit \+
> 1.248 polysaccharide_xylose_unit  \+
> 0.301 polysaccharide_mannose_unit  \+
> 0.144 polysaccharide_galactose_unit \+
> 3.254 polysaccharide_glucose_unit \+
> 0.166 polysaccharide_galacturonate_unit \+
> 0.166 polysaccharide_glucuronate_unit = hemicellulose_biomass.

###  10.3. Total carbohydrates

We recalculated the stoichiometries of the carbohydrate-producing reaction to
account for the differing molecular weight of our representation of cellulose
(‘CELLULOSE_monomer_equivalent’, effectively a glucose molecule), account for
the fact that one unit of hemicellulose represents one gram, not one
(milli)mole, and express pectin in terms of polysaccharide_galacturonate_unit,
reflecting a belief that UDP is released in the formation of pectin from UDP-
D-galacturonate, rather than retained in the polymer [129].

It is not clear what form the ‘mannose’ referred to by Penning de Vries et al
should be assumed to take, as free mannose is not found in plants under most
circumstances (see, e.g., [130, 131, 132].) Here we somewhat arbitrarily
choose mannose-6-phopshate.

Component | Species in model | unit wt (mg) | wt fraction | units/g product
---|---|---|---|---
Ribose | RIBOSE | 150.053 | 0.010 | 0.067
Glucose | GLC | 180.063 | 0.050 | 0.278
Fructose | FRU | 180.063 | 0.020 | 0.111
Mannose | MANNOSE-6P | 258.120 | 0.010 | 0.039
Galactose | GALACTOSE | 180.063 | 0.010 | 0.056
Sucrose | SUCROS | 342.116 | 0.050 | 0.146
Cellulose | CELLULOSE_monomer_equivalent | 180.160 | 0.400 | 2.220
Hemicellulose | hemicellulose_biomass | 1000.000 | 0.400 | 0.400
Pectin | polysaccharide_galacturonate_unit | 193.130 | 0.050 | 0.259
Table 2. Carbohydrate species in biomass.

Table 2 shows the calculation, resulting in the equation:

> 0.067 RIBOSE \+ 0.278 GLC \+ 0.111 FRU \+ 0.039 MANNOSE-6P \+ 0.056
> GALACTOSE \+ 0.146 SUCROSE \+ 2.220 CELLULOSE_monomer_equivalent \+ 0.400
> hemicellulose_biomass \+ 0.259 polysaccharide_galacturonate_unit =
> carbohydrates_biomass.

###  10.4. Organic acids

We adopt this reaction as is. In the terminology of our model, the resulting
equation is:

> 0.556 OXALATE \+ 0.676 GLYOX \+ 1.515 OXALACETIC_ACID \+ 0.746 MAL \+ 1.562
> CIT \+ 1.724 CIS-ACONITATE = organic_acids_biomass.

###  10.5. Protein and free amino acids

We adopt these reactions as is. In the terminology of our model, the resulting
equations are:

> 1.15 L-ALPHA-ALANINE \+ 0.0959 ARG \+ 0.414 L-ASPARTATE \+ 0.0313 CYS \+
> 1.53 GLT \+ 0.0445 GLY \+ 0.0915 HIS \+ 0.465 ILE \+ 1.51 LEU \+ 5.71e-05
> LYS \+ 0.123 MET \+ 0.314 PHE \+ 0.762 PRO \+ 0.612 SER \+ 0.175 THR \+
> 0.00409 TRP \+ 0.244 TYR \+ 0.25 VAL = protein_biomass

and

> 0.624 L-ALPHA-ALANINE \+ 0.319 ARG \+ 0.418 L-ASPARTATE \+ 0.231 CYS \+
> 0.378 GLT \+ 0.740 GLY \+ 0.358 HIS \+ 0.424 ILE \+ 0.424 LEU \+ 0.380 LYS
> \+ 0.373 MET \+ 0.337 PHE \+ 0.483 PRO \+ 0.529 SER \+ 0.467 THR \+ 0.272
> TRP \+ 0.307 TYR \+ 0.475 VAL = free_aa_biomass.

###  10.6. Lignin

We adopt this reaction as is. In the terminology of our model, the resulting
equation is:

> 2.221 COUMARYL-ALCOHOL \+ 1.851 CONIFERYL-ALCOHOL \+ 1.587 SINAPYL-ALCOHOL =
> lignin_biomass.

###  10.7. Nucleic acids

We adopt this reaction as is (though note that, as discussed above, nucleotide
triphosphates are not necessarily the appropriate best representation for
polymerized nucleic acids). In the terminology of our model, the resulting
equation is:

> 0.247 ATP \+ 0.239 GTP \+ 0.259 CTP \+ 0.258 UTP \+ 0.255 DATP \+ 0.247 DGTP
> \+ 0.268 DCTP \+ 0.259 TTP = nucleic_acids_biomass.

###  10.8. Nitrogenous compounds

We use the same nitrogenous compound weight fraction breakdown, but
recalculate the stoichiometric coefficients accounting for the fact that the
protein biomass, free amino acid biomass, and nucleotide biomass species each
represent one gram, so that the appropriate stoichiometric coefficients of
those species for the production of one total gram of nitrogenous compounds
are simply the weight fractions; see Table 3.

Component | Species in model | unit wt (mg) | wt fraction | units/g product
---|---|---|---|---
Amino acids | free_aa_biomass | 1000.000 | 0.100 | 0.100
Proteins | protein_biomass | 1000.000 | 0.870 | 0.870
Nucleic acids | nucleic_acids_biomass | 1000.000 | 0.030 | 0.030
Table 3. Nitrogenous biomass breakdown.

The resulting equation is

> 0.100 free_aa_biomass \+ 0.870 protein_biomass \+ 0.030
> nucleic_acids_biomass = nitrogenous_biomass.

###  10.9. Inorganic materials

We ignore these entirely, as they play no other role in the model. (Note that
even in iRS1563 the two species involved, potassium and chloride, participate
only in source and sink reactions.)

###  10.10. Total biomass reaction

We drop the inorganic materials term (note that weight fractions now add to
0.95) and recalculate the stoichiometric coefficients, accounting for the fact
that the component biomass subspecies each represent one gram; see Table 4.

Component | Species in model | unit wt (mg) | wt fraction | units/g product
---|---|---|---|---
Nitrogenous compounds | nitrogenous_biomass | 1000.000 | 0.230 | 0.230
Carbohydrates | carbohydrates_biomass | 1000.000 | 0.565 | 0.565
Lipids | fatty_acids_biomass | 1000.000 | 0.025 | 0.025
Lignin | lignin_biomass | 1000.000 | 0.080 | 0.080
Organic acids | organic_acids_biomass | 1000.000 | 0.050 | 0.050
Table 4. Breakdown of total biomass.

The final equation is

> 0.230 nitrogenous_biomass \+ 0.565 carbohydrates_biomass \+ 0.025
> fatty_acids_biomass \+ 0.080 lignin_biomass \+ 0.050 organic_acids_biomass =
> total_biomass.

Saha et al additionally incorporate an ATP cost in their overall biomass
reaction, based on that used in an earlier Arabidopsis model (AraGEM [41])
Combining this ATP hydrolysis with a sink of total biomass, we arrive at the
overall equation for biomass production and growth (CombinedBiomassReaction):

> 1.0 total_biomass \+ 30.0 ATP \+ 30.0 WATER = 30.0 ADP \+ 30.0 Pi \+ 30.0
> PROTON

###  10.11. Protonation

Throughout, note that the molecular weights of species in our model may differ
somewhat from those used in the iRS1563 table because of differing assumptions
about protonation. The practical consequences of this difference should be
limited.

###  10.12. Oxalate

Early drafts of the model could not produce oxalate. CornCyc indicates its
production as resulting only from ascorbic acid catabolism with concomitant
production of L-threonate. Recent reviews suggest this is the primary pathway
of oxalate production in plant species which form calcium oxalate crystals,
with the threonate ultimately being oxidized to tartrate [133, 134, 135],
though the pathways of production of soluble oxalate are less clear [134]. We
found little immediate evidence that tartrate (or threonate) is formed in
maize leaves at levels comparable to that of oxalate, or of pathways which
could further metabolize the tartrate.

Of the three reactions in iRS1563 which could produce oxalate, only one has an
associated gene: oxalate carboxylase (oxalate = formate + CO2); KEGG R00522
(EC4.1.1.2). The gene, ‘ACG37538’, may correspond to GRMZM2G103512, whose best
Arabidopsis hit is AT1G09560.1 (germin-like protein 5); it may thus be more
likely to be an oxalate-consuming oxidase [136] than an oxalate carboxylase,
though no function was computationally predicted for GRMZM2G103512 in CornCyc.

We decided the available information did not allow us to accurately model
oxalate production in maize. However, to retain the iRS1563 biomass equation
and ensure that mass and elemental balance was preserved, we allowed
production of oxalate from oxaloacetate by oxaloacetase (EC 3.7.1.1; PlantCyc
OXALOACETASE-RXN, [137]). This simple reaction has been observed in fungi
[138] but is considered unlikely to be widespread in plants [134].

##  11\. Plasmodesmatal transport reactions

Species allowed to be exchanged between cell types through the plasmodesmata
included:

  * •

carbon dioxide and oxygen;

  * •

known C4 cycle metabolites alanine, aspartate, malate, PEP, and pyruvate;

  * •

the Calvin cycle intermediates glyceraldehyde 3-phosphate and
3-phosphoglycerate;

  * •

photorespiratory metabolites glycerate, glycolate, serine, and glycine;

  * •

nutrients sucrose, phosphate, nitrate, ammonia, sulfate and magnesium;

  * •

glutamate and 2-ketoglutarate;

  * •

and cysteine and glutathione [139].

The inclusion of compounds involved in NAD-ME C4 or C3-C4 intermediate
photorespiratory carbon concentrating mechanism is not meant to suggest such a
system is necessarily active in maize but merely reflects our knowledge that
significant transport of those species between mesophyll and bundle sheath can
occur under at least some circumstances.

List of suppfigures 1  Phloem transport. Transport of nitrogen (upper panel)
and sulfur (lower panel) through the phloem in the best-fitting solution.
Dotted lines indicate minimum and maximum predicted values consistent with an
objective function value no more than 0.1% worse than the optimum.  List of
suppfigures 2  Photosystem II in mesophyll and bundle sheath. Dashed and
dotted lines indicate minimum and maximum predicted values consistent with an
objective function value no more than 0.1% worse than the optimum.  List of
suppfigures 3 Bundle sheath PEPC flux in the best-fitting solution. Dotted
lines indicate minimum and maximum predicted values consistent with an
objective function value no more than 0.1% worse than the optimum. List of
suppfigures 4  Summary of predictions for the gradient model using the least-
squares method without per-reaction scale factors. In eq. (8), $s_{i}=0$ for
all reactions $i$. (a) Sucrose and CO2 uptake rates (compare to figure 3a).
(b) Rates of carboxylation by PEPC and Rubisco (compare to figure 4b). (c)
Predicted rate for the reactions of the chlorophyllide A synthesis pathway
(compare to figure 6b.) (d) Predicted rates at the arogenate branch point
(compare to figure 6d). (e) Predicted oxygen and carbon dioxide levels in the
bundle sheath, with straight lines showing mesophyll levels (compare to figure
4d). (f) Distribution of correlation coefficients between data and predicted
fluxes for each reaction. (blue, this method; red, standard method.)
Correlation coefficients for reactions with zero predicted flux are taken to
be zero, resulting in the visible peak in the histogram. List of suppfigures 5
Summary of predictions for the gradient model using the E-Flux method. For
explanation of each panel, see S4_Figure.  List of suppfigures 6  Summary of
predictions for the gradient model using the E-Flux method with fixed biomass
composition. The biomass composition is fixed to that used by iRS1563, as
adapted (see Outline). For explanation of each panel, see S4_Figure. Note that
the chlorophyllide A synthesis pathway is blocked when the fixed biomass
composition is used.  List of suppfigures 7 Summary of predictions for the
gradient model with fixed biomass composition. For explanation of each panel,
see S4_Figure. Note that the chlorophyllide A synthesis pathway is blocked
when the fixed biomass composition is used. List of suppfigures 8  Predicted
biomass production rates in mesophyll and bundle sheath cells with fixed
biomass composition.

reaction name in model associated genes malate dehydrogenase (NADP)
MALATE_DEHYDROGENASE_NADP__RXN_chloroplast 1 alanine aminotransferase
ALANINE_AMINOTRANSFERASE_RXN 10 aspartate aminotransferase ASPAMINOTRANS_RXN 7
NAD-malic enzyme EC_1_1_1_39 2 NADP-malic enzyme (cytosol) MALIC_NADP_RXN 4
NADP-malic enzyme (chloroplast) MALIC_NADP_RXN_chloroplast 2 PEPCK
PEPCARBOXYKIN_RXN 6 PPDK PYRUVATEORTHOPHOSPHATE_DIKINASE_RXN_chloroplast 2
adenylate kinase ADENYL_KIN_RXN_chloroplast 6 pyrophosphatase
INORGPYROPHOSPHAT_RXN_chloroplast 2

List of supptables 1 Detailed parameters contributing to the effective PEP
regeneration rate: reactions in the genome-scale model which contribute to the
effective maximum PEP regeneration capacity, and the number of genes
associated with each. In addition to the reactions listed, transport
capacities of pyruvate, PEP, alanine, aspartate and malate across the
plasmodesmata and pyruvate, PEP, malate and oxaloacetate across the
chloroplast inner membrane could limit this rate; the model currently
associates no genes with these transport reactions. List of suppfigures 9
Predicted variable values in an FBA calculation that does not incorporate
expression data, compared to the best-fit and E-Flux methods. The FBA
calculation minimizes total flux while achieving the same total rate of CO2
assimilation as predicted at the tip of the leaf in the fitting results. Left
panel, FBA reaction rates vs. reaction rates predicted at the tip of the leaf
in the best-fitting solution; right panel, FBA reaction rates vs. reaction
rates predicted at the tip of the leaf by the E-Flux method. Axis limits
exclude a small number of reactions of particularly large flux. Fluxes in µmol
m-2 s-1.

## Index to additional supporting information files

Except as noted, these are available as arXiv ancillary files.

### S11 Model

iEB5204 in SBML format.

### S12 Model

iEB2140 in SBML format.

### S13 Model

iEB2140x2 in SBML format.

### S14 Protocol

Source code for the nonlinear constraint-based modeling package fluxtools.
Available at http://github.com/ebogart/fluxtools.

### S15 Protocol

Source code and input files for the calculations discussed above. Available at
http://github.com/ebogart/multiscale_c4_source.

### S16 Table

Predicted variable values along the leaf gradient.

### S17 Table

Upper and lower bounds on predicted values of selected variables along the
leaf gradient, from FVA calculations.