The plant leaf is the main site of photosynthesis. This process converts light energy and inorganic nutrients into chemical energy and organic building blocks for the biosynthesis and maintenance of cellular components and to support the growth of the rest of the plant. The leaf is also the site of gas–water exchange and due to its large surface, it is particularly vulnerable to pathogen attacks. Therefore, the leaf's performance and metabolic modes are inherently determined by its interaction with the environment. Mathematical models of plant metabolism have been successfully applied to study various aspects of photosynthesis, carbon and nitrogen assimilation and metabolism, aided suggesting metabolic intervention strategies for optimized leaf performance, and gave us insights into evolutionary drivers of plant metabolism in various environments. With the increasing pressure to improve agricultural performance in current and future climates, these models have become important tools to improve our understanding of plant–environment interactions and to propel plant breeders efforts. This overview article reviews applications of large-scale metabolic models of leaf metabolism to study plant–environment interactions by means of flux-balance analysis. The presented studies are organized in two ways — by the way the environment interactions are modelled — via external constraints or data-integration and by the studied environmental interactions — abiotic or biotic.

The plant leaf is the main photosynthetic organ where light energy is converted into chemical energy to support the maintenance and growth of the leaf and the rest of the plant. The leaf is also the site of gas-exchange (O2, CO2 and water vapour) through the stomata and different forms of photosynthesis and leaf architecture — C3, C4 and crassulacean acid metabolism (CAM) — have evolved to cope with different environments. However, given its large surface area the plant leaf is also vulnerable to pathogen attacks. Due to these various factors, the leaf is a main determinant of plant growth and health. Therefore, gaining a better understanding of leaf–environment interactions is key to developing strategies to better equip our crop plants for productivity requirements in current and future climates. Computational models of leaf metabolism can aid this quest by explaining metabolic, architectural, and evolutionary aspects to guide the engineering of more productive and resistant crop species. In this review, I focus on flux-balance models of leaf metabolism and highlight their applications to study plant–environment interactions.

Flux-balance models rely on the stoichiometry of the metabolic network under consideration and an optimality assumption based on which metabolic steady-state fluxes can be predicted [1,2]. In the past decade, plant stoichiometric network models have been reconstructed for a large range of species and to varying extent, covering cell suspension models, multi-cell type, and whole-plant models [3–6]. Several pipelines for the automated reconstruction of plant metabolic network models are available [7–10]. The optimality assumption or ‘objective function’ depends on the metabolic system under consideration and reflects its function in the plant. For instance, for a growing leaf the objective function could be the maximization of biomass production and for a mature leaf the export of photosynthetic assimilates to the phloem [11]. Other optimality criteria have also been applied, such as the minimization of internal fluxes or photon usage efficiency [12,13]. The latter objectives are often chosen when metabolic outputs, such as growth rates, have been experimentally determined and can be applied as constraints. Additionally, external constraints, such as limited nutrient or light availability, can be modelled by constraining the respective uptake fluxes. Depending on these constraints the model then predicts optimal metabolic steady-state flux modes which can then be analyzed with respect to changes between conditions and complement experimental observations. See Figure 1 for a schematic representation.

### Schematic representation of flux-balance modelling.

Figure 1.
Schematic representation of flux-balance modelling.

The top part illustrates the metabolic system under consideration. (Left) External constraints, such as exchange and uptake rates and biochemical properties constrain the allowed metabolic fluxes. (Middle) The metabolic network coverts input compounds and generates a metabolic output via a series of biochemical reactions. (Right) This metabolic output can be measured as e.g. growth rate, biomass composition or flux patterns. The bottom part illustrates the mathematical representation of the metabolic system under consideration. (Left) External constraints as well as experimental data can serve to constrain the lower and upper flux bounds (vmin ≥ v ≥ vmax). (Middle) The metabolic network and fluxes are represented as a stoichiometric matrix, S and the metabolic flux vector, v, respectively. The (usually) linear optimization problem optimizes the biological objective, e.g. growth or phloem output, in a metabolic steady-state (S × v = 0). (Right) Result of this optimization is a set of flux distributions captured in v. These optimal solutions can then be subjected to further analysis.

Figure 1.
Schematic representation of flux-balance modelling.

The top part illustrates the metabolic system under consideration. (Left) External constraints, such as exchange and uptake rates and biochemical properties constrain the allowed metabolic fluxes. (Middle) The metabolic network coverts input compounds and generates a metabolic output via a series of biochemical reactions. (Right) This metabolic output can be measured as e.g. growth rate, biomass composition or flux patterns. The bottom part illustrates the mathematical representation of the metabolic system under consideration. (Left) External constraints as well as experimental data can serve to constrain the lower and upper flux bounds (vmin ≥ v ≥ vmax). (Middle) The metabolic network and fluxes are represented as a stoichiometric matrix, S and the metabolic flux vector, v, respectively. The (usually) linear optimization problem optimizes the biological objective, e.g. growth or phloem output, in a metabolic steady-state (S × v = 0). (Right) Result of this optimization is a set of flux distributions captured in v. These optimal solutions can then be subjected to further analysis.

Close modal

As an extension to this, various data-integrative approaches have been proposed to model condition-specific metabolic fluxes or to extract context-specific metabolic networks. Most of these approaches employ transcript, protein or metabolite abundances to constraint the flux boundaries of the respective reactions [14–17].

While any flux-balance model of leaf metabolism is inherently coupled to the environment via exchange reactions, such as light influx, CO2 and O2 exchange, and nutrient uptake; the focus here will be on studies that investigate the leaf's metabolic response to changes in those parameters or which integrate data collected from leaf material exposed to different environments. The highlighted studies are divided into two groups — environmental constraint-driven and data-integrative — based on how the environment-specificity is achieved. Moreover, the studies are organized according to the modelled interaction, i.e. abiotic or biotic. A schematic overview of the various modelling approaches is shown in Figure 2 and summarized in Table 1.

### Representation of environment-coupled models of leaf metabolism.

Figure 2.
Representation of environment-coupled models of leaf metabolism.

(Left) Abiotic as well as biotic interactions between a plant leaf and the environment have been modelled using flux-balance analysis. (Middle) These interactions can occur at various sites of the leaf, e.g. light uptake through photosynthesis in the chloroplasts or gas-exchange through the stomata (marked A to H and listed on the right). (Right) Studies investigating abiotic constraints can be further classified based on environment-specific exchange constraints, such as light intensity or nutrient availability, data-integrative approaches, such as temporal response or accession specific outputs, or combinations of both. Data-integrative approaches have also been used to model biotic interactions.

Figure 2.
Representation of environment-coupled models of leaf metabolism.

(Left) Abiotic as well as biotic interactions between a plant leaf and the environment have been modelled using flux-balance analysis. (Middle) These interactions can occur at various sites of the leaf, e.g. light uptake through photosynthesis in the chloroplasts or gas-exchange through the stomata (marked A to H and listed on the right). (Right) Studies investigating abiotic constraints can be further classified based on environment-specific exchange constraints, such as light intensity or nutrient availability, data-integrative approaches, such as temporal response or accession specific outputs, or combinations of both. Data-integrative approaches have also been used to model biotic interactions.

Close modal
Table 1.
Environment-coupled flux-balance models of leaf metabolism
Author and yearSpeciesBiological questionModelling approach
Abiotic interactions — environmental constraint-driven models
Poolman et al. [18], Poolman et al. [11], Chatterjee et al. [20], Cheung et al. [19Rice, Arabidopsis Response to different light intensities Different light constraints
Simons et al. [21], Arnold and Nikoloski [22], Arnold et al. [25Maize, Arabidopsis Response to different nitrogen levels and sources Condition-specific biomass compositions
Shaw and Cheung [27Arabidopsis Resource partitioning in whole-plant model Dynamic FBA and different nutrient availability
Lakshmanan et al. [28], Chatterjee et al. [20], Yuan et al. [29], Shameer et al. [30Rice, Tomato, Generic CAM model Response to different CO2 levels Constraints on rubisco's vC/vO ratio
Mallmann et al. [33Flaveria genus Response to different CO2 levels Flux-balance model coupled to a kinetic model of photosynthesis
Blätke and Brautigam [34Generic C4 model Evolutionary drivers of C4 photosynthesis CCM-dependent rubisco population, cell type-specific light availability
Töpfer et al. [36Generic C3 — CAM model Water-saving flux modes in a C3 leaf 24 hour diel resolution, flux-balance model coupled to a biophysical model of gas–water exchange
Abiotic interactions — data-driven models
Töpfer et al. [44,64], Töpfer et al. [50Arabidopsis Response to changes in light and temperature Transcript and metabolomics data integration
Lakshmanan et al. [46Rice Response to changes in light Transcriptomics data integration
Liu et al. [65Arabidopsis Response to low and elevated CO2 Transcriptomics data integration
Bogart and Myers [48Maize Source to sink transition along the leaf Transcriptomics data integration, non-linear constraints
Nägele and Weckwert [49Arabidopsis Metabolite compartmentation in different accessions exposed to low temperature Metabolomics data integration
Sajitz-Hermstein et al. [51Arabidopsis High to low CO2 acclimation in wild type and photorespiratory mutants Metabolomics data integration
Biotic interactions — data-driven models
Botero et al. [56Potato Effect of pathogen attack on photosynthetic activity Transcriptomics data integration
Rodenburg et al. [58Tomato Nutrient exchange between host leaf and pathogen Transcriptomics data integration
Author and yearSpeciesBiological questionModelling approach
Abiotic interactions — environmental constraint-driven models
Poolman et al. [18], Poolman et al. [11], Chatterjee et al. [20], Cheung et al. [19Rice, Arabidopsis Response to different light intensities Different light constraints
Simons et al. [21], Arnold and Nikoloski [22], Arnold et al. [25Maize, Arabidopsis Response to different nitrogen levels and sources Condition-specific biomass compositions
Shaw and Cheung [27Arabidopsis Resource partitioning in whole-plant model Dynamic FBA and different nutrient availability
Lakshmanan et al. [28], Chatterjee et al. [20], Yuan et al. [29], Shameer et al. [30Rice, Tomato, Generic CAM model Response to different CO2 levels Constraints on rubisco's vC/vO ratio
Mallmann et al. [33Flaveria genus Response to different CO2 levels Flux-balance model coupled to a kinetic model of photosynthesis
Blätke and Brautigam [34Generic C4 model Evolutionary drivers of C4 photosynthesis CCM-dependent rubisco population, cell type-specific light availability
Töpfer et al. [36Generic C3 — CAM model Water-saving flux modes in a C3 leaf 24 hour diel resolution, flux-balance model coupled to a biophysical model of gas–water exchange
Abiotic interactions — data-driven models
Töpfer et al. [44,64], Töpfer et al. [50Arabidopsis Response to changes in light and temperature Transcript and metabolomics data integration
Lakshmanan et al. [46Rice Response to changes in light Transcriptomics data integration
Liu et al. [65Arabidopsis Response to low and elevated CO2 Transcriptomics data integration
Bogart and Myers [48Maize Source to sink transition along the leaf Transcriptomics data integration, non-linear constraints
Nägele and Weckwert [49Arabidopsis Metabolite compartmentation in different accessions exposed to low temperature Metabolomics data integration
Sajitz-Hermstein et al. [51Arabidopsis High to low CO2 acclimation in wild type and photorespiratory mutants Metabolomics data integration
Biotic interactions — data-driven models
Botero et al. [56Potato Effect of pathogen attack on photosynthetic activity Transcriptomics data integration
Rodenburg et al. [58Tomato Nutrient exchange between host leaf and pathogen Transcriptomics data integration

Overview of studies which analyze plant–environment interactions by coupling mathematical models of leaf metabolism to the environment. The presented approaches are organized by abiotic and biotic interactions and by the applied modelling approach — environmental constraint-driven and data-integrative.

### Environmental constraint-driven models

#### Response to different light intensities

In one of the earliest studies in 2013, Poolman et al. [18] analyzed a genome-scale model representing a developing rice leaf cell to identify changes in reaction fluxes in response to changes in light influx. One of their findings put forward a beneficial role for the energy-consuming photorespiratory pathway under supraoptimal light conditions. The photorespiratory pathway recycles toxic intermediates which are formed due to the dual activity of ribulose-1,5-bisphosphate carboxylase-oxygenase (rubisco), the main enzyme to convert atmospheric carbon dioxide to energy-rich molecules. With increasing temperature rubisco's ratio of carboxylation to oxygenation (vC/vO) shifts towards O2-fixation. This side reaction forms byproducts which are then recycled in the energy-wasting process of photorespiration. The model predicted that the photorespiratory pathway was spontaneously activated at supraoptimal light intensities and thereby dissipated excess light energy to prevent cell damage or overreduction.

A year later the same authors extended this study by exploring flux rearrangements during the transition from a growing to a mature leaf [11]. This was modelled as a transition from synthesizing leaf biomass to export of nutrients to the phloem as a function of light intensity. The authors found that the predicted flux patterns were relatively insensitive to the specific metabolic function of the leaf and attributed flux differences to changes in energy and redox requirements for the different metabolic outputs.

Cheung et al. [19] used the effect of light intensity on leaf metabolism together with a cost-weighted flux-balance formulation to study alternative metabolic pathways in Arabidopsis. The authors simulated C-limited conditions by fixing cellular output and maintenance costs and varied photon influx. They then simulated flux distributions for sets of randomly chosen flux weighting factors using flux minimization as the objective function. Their results emphasized the potential contribution of alternative fluxes in rebalancing and consuming ATP and NADPH, such as the chlorophyll and xanthophyll pigment cycles and various futile cycles, in different environmental and physiological conditions.

Later, Chatterjee et al. [20] used a different rice leaf model to study its responses to different light intensities and varying vC/vO ratios of rubisco, thereby mimicking increased photorespiration under drought stress, normal and suppressed photorespiration. Their comprehensive environment-scan combined with cost-weighted flux-balance analysis revealed different metabolic modes for maintaining redox and ATP balance including flux rearrangements across compartments and shifts in transport reactions. These findings highlighted alternative routes and metabolic flexibility for stress adaptation in leaf metabolism.

#### Response to different nitrogen levels and sources

Plant species, such as maize, sugarcane, and sorghum perform C4 photosynthesis. This specialized form of photosynthesis involves a carbon-concentration mechanism in which initial CO2 fixation and internal re-fixation by rubisco are separated in two adjacent cell types — the bundle sheath and mesophyll cells. This separation increases both nitrogen (N) and water use efficiency compared with C3 plants. To investigate the metabolic impact of N availability in maize Simons et al. [21] presented a C4 model capturing the interactions between bundle sheath and mesophyll cells. The authors used condition-specific biomass compositions and explored the effect of two knockdown mutants of glutamine synthetase (GS) — an enzyme essential for N assimilation. The authors used transcriptomic and proteomic data to introduce regulatory constraints that represented the metabolism of wild type and GS mutants under N-complete and -deficient conditions. This way, between nine and 100 reaction fluxes per condition were constrained. Using these condition-specific models they determined flux-sum ranges as a flux measure for the reactions associated with either the production or consumption of a given metabolite. They achieved up to 90% accuracy when comparing their model predictions with measured metabolite levels. Furthermore, the authors identified a set of genes for further testing which the model's predictions had related to changes in biomass formation in the considered conditions.

In the same year, Arnold and Nikoloski [22] presented an Arabidopsis core model equipped with three biomass compositions to simulate carbon-limited, N-limited, and optimal growth conditions. This core model and two previous Arabidopsis models [23,24] were later used to in silico analyze the effect of N supply (nitrate/ammonium) on individual amino acid synthesis costs under autotrophic/heterotrophic and day/night growth conditions [25]. The authors quantified the synthesis costs of amino acids in terms of ATP demand and found these costs to be highly dependent on the environment — most amino acid costs for night conditions were higher than those for heterotrophic day conditions and the costs for autotrophic conditions were higher than the costs for heterotrophic conditions. The study further confirmed $NH4+$ uptake to be cheaper than $NO3−$ uptake to meet the plant's N demand due to the extra cost for nitrate reduction.

Shaw and Cheung used their previously developed diel Arabidopsis [26] model to generate a sophisticated dynamic multi-tissue, day–night modelling framework in which they studied the effect of both N-rich and -limiting conditions on Arabidopsis leaf and root growth [27]. Their framework demonstrated the power of integrated whole-plant analysis to uncover optimal growth strategies. More specifically, their model showed that it is energetically most efficient to store nitrate taken up into the root during the night in the vacuole and to then transport and fix it in the leaf during the day when light energy is available.

#### Response to different CO2 levels

Studies that investigated the effect of different CO2 levels on leaf metabolism and particularly photorespiration include models of rice [20,28], tomato [29] and generic models of CAM photosynthesis [30]. Typically, constraints on rubisco's vC/vO ratio were implemented to model the leaf's metabolic response to changing CO2 levels or temperature (as vC/vO is temperature-dependent, see also section ‘Response to different light intensities’).

Of particular sophistication is a study by Mallmann et al. which combined a kinetic model of photosynthesis [31] and a stoichiometric model [32] to explore the evolution of C4 photosynthesis in the genus Flaveria [33]. Here, environment-specificity of the model was achieved by using the output of the kinetic model of C3–C4 intermediate photosynthesis to constrain key C4 parameters, such as net CO2 uptake, rubisco's vC/vO ratio in mesophyll and bundle sheath cells, CO2 leakage from the bundle sheath, PEP-carboxylase activity in the mesophyll, the activity of NADP-malic enzyme (ME) in the bundle sheath, plasmodesmatal flux of glycine and serine, and decarboxylation by the glycine decarboxylase complex in the genome-scale stoichiometric model of C4 photosynthesis. This way, the authors found that C2 photosynthesis, an intermediate state where CO2 is relocated from the mesophyll to the bundle sheath cells via the photorespiratory intermediates glycolate and glycerate, caused a N missbalance and that rebalancing N metabolism was a driving force for the evolution of C4 photosynthesis.

#### Spatial and temporal constraints as metabolic drivers

In a related study, Blätke and Bräutigam [34] studied the evolutionary trajectory of C4 metabolism by examining selective pressures that potentially have led to the occurrence of this C-fixation mechanism in a two-celled stoichiometric model representing the C4-typical mesophyll-bundle sheath cell Kranz anatomy. In contrast with earlier studies of C4 metabolism [21,32,33] in their model C4 photosynthesis was not enforced by fixing a priori defined flux patterns between the two cell types but the C4 syndrome was allowed to emerge under a set of different input constraints. The authors employed an artifice in which they approximated the CO2 concentration-dependent changes in rubisco's vC/vO ratio by modelling two rubisco populations in the bundle sheath cells — the native rubisco, which performs both carboxylation and oxygenation at a typical C3 plant ratio of 3 : 1 [35], and a carbon-concentration mechanism-dependent rubisco population which only catalyzed the carboxylation of ribulose 1,5-bisphosphate and used exclusively CO2 released in the bundle sheath cells. Using this setup, they tested several scenarios and found that high photorespiration and N limitation drove the emergence of C4 flux patterns in the model. The model also predicted that light availability and distribution across the leaf section, i.e. between mesophyll and bundle sheath cells could play a role in the evolutionary choices of possible decarboxylation enzymes NAD-ME, NADP-ME, and PEP-carboxykinase. Finally, the model predicted that C2 photosynthesis was optimal under particular conditions and supported the hypothesis of C2 photosynthesis as a stable intermediate state between C3 and C4 photosynthesis.

In our recent work, we focussed on studying the tradeoff between water-saving and leaf productivity and the analysis of alternative CAM-like flux modes in a C3 metabolic network [36]. In contrast with C4 photosynthesis, in CAM photosynthesis initial and re-fixation of CO2 are not spatially but temporally separated between day and night. This way CO2 can be taken up through the stomata at colder and more humid night-time hours, thereby minimizing water loss through transpiration. To capture this behaviour, we coupled a biophysical model of gas–water exchange to a time-resolved diel leaf metabolic model. This way, we were able to model the effect of three main abiotic parameters — temperature, relative humidity, and light intensity — on leaf metabolism. We used this environment-coupled model to study the emergence of water-saving flux modes on the Pareto frontier for water-saving and leaf productivity. Our analysis revealed that vacuolar storage capacity is a main determinant of the extent of CAM. Additionally, we identified the mitochondrial enzyme isocitrate dehydrogenase (ICDH) and an isocitrate-citrate-proline-2-oxo-glutarate cycle as a potential contributor to initial carbon fixation at night. Model analysis across a wide range of environmental parameters showed that the water-saving potential of CAM-like flux modes strongly depends on the environment and that under certain conditions CAM with night-time carbon fixation by ICDH could reach 11% total water saving for the conditions tested.

### Environment-specific leaf models through data integration

The approaches highlighted so far all have in common that the environment-specificity was achieved by constraining input parameters, such as light intensity, CO2 or N uptake, the light gradient in the leaf cross section, the gas–water exchange or few reaction fluxes based on experimental data. Most of the fluxes were left unconstrained and could vary between a generic lower and upper boundary. In contrast, several approaches that rely on the integration of Omics data and typically constrain several hundreds to thousands of reaction fluxes have been applied to plant metabolic models and resulted in environment-specific models and flux mode predictions. The majority of the methods employ transcriptomics data to constrain flux boundaries and can be grouped based on whether they employ binarized (e.g. PROM [37], iMAT [37,38], MADE [39], GIMME [40], AdaM [41]) or continuous gene-expression levels (e.g. E-Flux [42]) to constrain fluxes through the respective reactions.

#### Transcriptomics data integration

In an earlier study, we used an Arabidopsis metabolic network model accounting for both primary and secondary metabolism [43] and applied a transcriptomics data integrative approach to study changes in leaf metabolism in response to changes in temperature and light intensities [44,45]. The data integration was achieved by using a modified version of the E-Flux approach and enabled us to define three optimization-based indices to characterize different aspects of metabolic pathway behaviour in the context of the entire network. Our study highlighted pathways that showed differential behaviour with respect to a null model and their involvement in the different acclimation responses.

Lakshmanan et al. [46] investigated the effect of four different light treatments and darkness on rice leaf metabolism by combining a rice genome-scale model with transcriptome profiles. The authors performed flux sampling based on the E-Flux method and analyzed the up- and down-regulation of individual metabolic pathways. Additionally, they compared their modelling results to previously measured metabolomics data [47]. Using these combined approaches, they found that photosynthesis and secondary metabolism were up-regulated in blue light and reserve carbohydrates degradation was pronounced in the dark. The analysis further identified phytohormones, such as abscisate, ethylene, gibberellin, and jasmonate as keymarkers of light-mediated regulation and helped elucidate the transcriptional control of red and blue light signals.

Another data-integrative study explored Arabidopsis's response to low and elevated CO2 by generating condition-specific models based on the integration of transcriptomics data using the iMat approach. Liu et al. showed that the simulated CO2-fixation at different CO2 concentrations was consistent with measured CO2-assimilation curves. Additionally, they predicted post-transcriptionally regulated reactions across different CO2 concentrations and that low CO2 stress requires stronger metabolic adjustment than elevated CO2 conditions.

While the above-outlined approaches incorporated time-resolved and/or compartment-specific data, Bogart and Myers [48] used spatially resolved data to study the transition from source to sink tissue along the leaf gradient. Their analysis predicted metabolic fluxes by integrating spatially resolved transcriptomics and enzyme activity data with a genome-scale model of maize leaf metabolism. Additionally, they modelled the non-linear relationships between CO2 and O2 levels and reaction rates in the C4 system. This resulted in a non-linearly constrained model describing mesophyll and bundle sheath cells in 15 segments of the developing maize leaf. The analysis successfully recaptured results from radiolabeling experiments and the observed base-to-tip transition between C-importing and -exporting tissue along the leaf axis. Thus, this approach laid the foundation for studying the response of various heterogeneous metabolic systems to environmental and biochemical perturbations by means of flux-balance analysis.

#### Metabolomics data integration

Nägele and Weckwerth [49] proposed a metabolomics data integrative approach which they applied to study metabolite compartmentation in leaves of different Arabidopsis accessions exposed to low temperature. The authors used a compartmentalized Arabidopsis network reconstruction [43] and compartment-specific metabolomics data to generate a reduced metabolic interaction matrix, i.e. a matrix that contains only measurable metabolites and their respective compartmentation. The extracted subcellular metabolic network comprised over 500 metabolic intermediates and interactions and was integrated with metabolite covariances of different Arabidopsis thaliana accessions. Analysis of these data-integrated networks revealed differences in the regulation of carbohydrate compartmentation in the three investigated accessions and highlighted differential regulation of the interconversion and transport of vacuolar glucose and fructose.

In a followup study to the above-mentioned analysis of Arabidopsis's response to changes in temperature and light intensities [44] we tested the informative value of the transcriptomics data-integrative approach with respect to inferences that could be made on the level of metabolites [50]. We observed that substrates in pathways which our transcript-data integrative study had classified as important for the adjustment processes showed less temporal fluctuations. Moreover, these pathways had on average fewer substrates than the average pathway investigated. This led us to the conclusion that the observed substrate robustness is an inducible genetic mechanism, both depending on the metabolic network structure and the specific environmental condition.

In an attempt to model flux re-routing upon external perturbations we developed a metabolomics data-integrative approach and applied it to predict reactions and pathways with altered fluxes in wild type Arabidopsis leaves and four photorespiratory mutants undergoing high-to-low CO2 acclimation [51,52]. The approach relied on relative metabolite abundances for at least two metabolic states and employed a mass-action like representation of differential fluxes. We found that the observed flux alterations in the knock-out mutants were mainly attributed to ATP synthase and the photosystem I, and fluxes through reactions from the Calvin–Benson–Bassham cycle, proline biosynthesis, N and redox metabolism and glycolysis. The study highlighted the power of integrated differential flux analysis to complement labor-intensive flux measurements which are usually restricted to small-scale networks.

In recent years, several studies have focussed on modelling plant-microbe interactions by analyzing, both beneficial symbiont relationships, such as interactions between the roots and N fixing soil bacteria [53,54] and plant–pathogen interactions. Two studies have focused on modelling the plant leaf–pathogen interactions between the solanaceous species potato (Solanum tuberosum) and tomato (Solanum lycopersicum) and the late blight causing oomycete Phytophthora infestans.

Botero et al. aimed at gaining a better understanding of the molecular basis of the previously identified decrease in photosynthetic activity of infected potato plants [55,56]. To this end, the authors reconstructed a genome-scale model of potato and integrated gene-expression data from three time-points after pathogen infection [57] using the maximization of biomass as an objective function. They predicted decreased photophosphorylation, the activity of the Calvin–Benson–Bassham cycle and starch synthesis as well as transiently increased photorespiration during the host–pathogen interaction and compared their results to experimental observations from the literature.

In contrast with this host-centric study, Rodenburg et al. studied the nutrient flux from the host to the pathogen during the infection [59]. Thus, the authors coupled a genome-scale model of tomato leaf metabolism [29] and P. infestans [58] by allowing all common metabolites to be exchanged between the plant and the pathogen. Using this combined plant-pathogen model they analyzed four scenarios in which they maximized biomass production of the pathogen and determined (i) the minimal set of nutrients P. infestans needed to import from the tomato leaf, (ii) and (iii) the minimum and maximum number of reactions used by P. infestans, (iv) a combination of (i) and (ii) where they determined the minimal nutrient uptake combined with minimal usage of pathogen reactions. The authors reasoned that the four investigated scenarios represented extreme cases and that the pool of imported nutrients would likely be a combination of the metabolite pools predicted in these scenarios. Additionally, the authors integrated dual-transcriptome time series data of a full late blight infection cycle on tomato leaves. Analysis of these contextualized models indicated that, as the infection progresses, P. infestans performed less de novo synthesis of metabolites and more scavenging from the tomato plant.

• Environment-coupled flux-balance models of leaf metabolism are a versatile means for multi-omics data integration and interpretation [17,60], they have closed gaps in our understanding of plant metabolism and led to the suggestion of engineering strategies to enhance plant's performance [36].

• Yet, one should be aware that flux-balance models are best suited to model plant systems in optimal growth conditions for which an objective function as well as a cellular maintenance costs can be readily defined. For non-optimal growth conditions determining an adequate objective function is less obvious and cellular maintenance costs might vary substantially and could account for a significant portion of the total energy budget [13].

• In the future, we will see continued efforts to integrate the here discussed environment-coupled flux-balance models into multi-scale modelling frameworks which can cover computational models of gene regulation, protein synthesis, metabolic pathways, and plant architecture up to the level of ecosystem models [61–63]. With the advancement and automatization of plant computational modelling I anticipate these models to play a crucial role in future agriculture research by guiding the development of crop species with improved yield and resistance in current and future climates.

The author declares that there are no competing interests associated with this manuscript.

N.T. did not receive specific funding for this work.

• CAM

crassulacean acid metabolism

•
• GS

glutamine synthetase

•
• ICDH

isocitrate dehydrogenase

•
• ME

malic enzyme

1
Orth
,
J.D.
,
Thiele
,
I.
and
Palsson
,
B.Ø
. (
2010
)
What is flux balance analysis?
Nat. Biotechnol.
28
,
245
248
2
Collakova
,
E.
,
Yen
,
J.Y.
and
Senger
,
R.S.
(
2012
)
Are we ready for genome-scale modeling in plants?
Plant. Sci.
191–192
,
53
70
3
Shaw
,
R.
and
Cheung
,
C.Y.M.
(
2019
)
Multi-tissue to whole plant metabolic modelling
.
Cell. Mol. Life Sci.
77
,
489
495
4
Nikoloski
,
Z.
,
Perez-Storey
,
R.
and
Sweetlove
,
L.J.
(
2015
)
Inference and prediction of metabolic network fluxes
.
Plant Physiol.
169
,
1443
1455
5
Sweetlove
,
L.J.
and
Ratcliffe
,
R.G.
(
2011
)
Flux-Balance modeling of plant metabolism
.
Front. Plant Sci.
2
,
38
6
de Oliveira Dal'Molin
,
C.G.
and
Nielsen
,
L.K.
(
2018
)
Plant genome-scale reconstruction: from single cell to multi-tissue modelling and omics analyses
.
Curr. Opin. Biotechnol.
49
,
42
48
7
de Oliveira Dal'Molin
,
C.G.
and
Nielsen
,
L.K.
(
2013
)
Plant genome-scale metabolic reconstruction and modelling
.
Curr. Opin. Biotechnol.
24
,
271
277
8
Seaver
,
S.M.D.
,
Henry
,
C.S.
and
Hanson
,
A.D.
(
2012
)
Frontiers in metabolic reconstruction and modeling of plant genomes
.
J. Exp. Bot.
63
,
2247
2258
9
Seaver
,
S.M.D.
,
Gerdes
,
S.
,
Frelin
,
O.
,
Lerma-Ortiz
,
C.
,
,
L.M.T.
,
Zallot
,
R.
et al (
2014
)
High-throughput comparison, functional annotation, and metabolic modeling of plant genomes using the PlantSEED resource
.
111
,
9645
9650
10
Seaver
,
S.M.D.
,
Lerma-Ortiz
,
C.
,
,
N.
,
Mikaili
,
A.
,
Sreedasyam
,
A.
,
Hanson
,
A.D.
et al (
2018
)
PlantSEED enables automated annotation and reconstruction of plant primary metabolism with improved compartmentalization and comparative consistency
.
Plant J.
95
,
1102
1113
11
Poolman
,
M.G.
,
Kundu
,
S.
,
Shaw
,
R.
and
Fell
,
D.A.
(
2014
)
Metabolic trade-offs between biomass synthesis and photosynthate export at different light intensities in a genome-scale metabolic model of rice
.
Front. Plant Sci.
5
,
656
12
Holzhütter
,
H.-G.
(
2004
)
The principle of flux minimization and its application to estimate stationary fluxes in metabolic networks
.
Eur. J. Biochem.
271
,
2905
2922
13
Cheung
,
C.Y.M.
,
Williams
,
T.C.R.
,
Poolman
,
M.G.
,
Fell
,
D.A.
,
Ratcliffe
,
R.G.
and
Sweetlove
,
L.J.
(
2013
)
A method for accounting for maintenance costs in flux balance analysis improves the prediction of plant cell metabolic phenotypes under stress conditions
.
Plant J.
75
,
1050
1061
14
Blazier
,
A.S.
and
Papin
,
J.A.
(
2012
)
Integration of expression data in genome-scale metabolic network reconstructions
.
Front. Physiol.
3
,
299
15
Åkesson
,
M.
,
Förster
,
J.
and
Nielsen
,
J.
(
2004
)
Integration of gene expression data into genome-scale metabolic models
.
Metab. Eng.
6
,
285
293
16
Estévez S
,
R.
and
Nikoloski
,
Z.
(
2014
)
Generalized framework for context-specific metabolic model extraction methods
.
Front. Plant Sci.
5
,
491
17
Töpfer
,
N.
,
Kleessen
,
S.
and
Nikoloski
,
Z.
(
2015
)
Integration of metabolomics data into metabolic networks
.
Front. Plant Sci.
6
,
49
18
Poolman
,
M.G.
,
Kundu
,
S.
,
Shaw
,
R.
and
Fell
,
D.A.
(
2013
)
Responses to light intensity in a genome-Scale model of rice metabolism
.
Plant Physiol.
162
,
1060
1072
19
Cheung
,
C.Y.M.
,
Ratcliffe
,
R.G.
and
Sweetlove
,
L.J.
(
2015
)
A method of accounting for enzyme costs in flux balance analysis reveals alternative pathways and metabolite stores in an illuminated Arabidopsis leaf
.
Plant Physiol.
169
,
1671
1682
20
Chatterjee
,
A.
,
Huma
,
B.
,
Shaw
,
R.
and
Kundu
,
S.
(
2017
)
Reconstruction of oryza sativa indica genome scale metabolic model and Its responses to varying RuBisCO activity, light intensity, and enzymatic cost conditions
.
Front. Plant Sci.
8
,
2060
21
Simons
,
M.
,
Saha
,
R.
,
Amiour
,
N.
,
Kumar
,
A.
,
Guillard
,
L.
,
Clément
,
G.
et al (
2014
)
Assessing the metabolic impact of nitrogen availability using a compartmentalized maize leaf genome-scale model
.
Plant Physiol.
166
,
1659
1674
22
Arnold
,
A.
and
Nikoloski
,
Z.
(
2014
)
Bottom-up metabolic reconstruction of Arabidopsis and Its application to determining the metabolic costs of enzyme production
.
Plant Physiol.
165
,
1380
1391
23
Poolman
,
M.G.
,
Miguet
,
L.
,
Sweetlove
,
L.J.
and
Fell
,
D.A.
(
2009
)
A genome-scale metabolic model of Arabidopsis and some of its properties
.
Plant Physiol.
151
,
1570
1581
24
de Oliveira Dal'Molin
,
C.G.
,
Quek
,
L.-E.
,
Palfreyman
,
R.W.
,
Brumbley
,
S.M.
and
Nielsen
,
L.K.
(
2010
)
AraGEM, a genome-scale reconstruction of the primary metabolic network in Arabidopsis
.
Plant Physiol.
152
,
579
589
25
Arnold
,
A.
,
Sajitz-Hermstein
,
M.
and
Nikoloski
,
Z.
(
2015
)
Effects of varying nitrogen sources on amino acid synthesis costs in Arabidopsis thaliana under different light and carbon-source conditions
.
PLoS One
10
,
e0116536
26
Cheung
,
C.Y.M.
,
Poolman
,
M.G.
,
Fell
,
D.A.
,
Ratcliffe
,
R.G.
and
Sweetlove
,
L.J.
(
2014
)
A diel flux balance model captures interactions between light and dark metabolism during day-night cycles in C3 and crassulacean acid metabolism leaves
.
Plant Physiol.
165
,
917
929
27
Shaw
,
R.
and
Cheung
,
C.Y.M.
(
2018
)
A dynamic multi-tissue flux balance model captures carbon and nitrogen metabolism and optimal resource partitioning during Arabidopsis growth
.
Front. Plant Sci.
9
,
884
28
Lakshmanan
,
M.
,
Mohanty
,
B.
and
Lee
,
D.-Y.
(
2013
)
Identifying essential genes/reactions of the rice photorespiration by in silico model-based analysis
.
Rice
13
,
20
29
Yuan
,
H.
,
Cheung
,
C.Y.M.
,
Poolman
,
M.G.
,
Hilbers
,
P.A.J.
and
van Riel
,
N.A.W.
(
2016
)
A genome-scale metabolic network reconstruction of tomato (Solanum lycopersicum L.) and its application to photorespiratory metabolism
.
Plant J.
85
,
289
304
30
Shameer
,
S.
,
Baghalian
,
K.
,
Cheung
,
C.Y.M.
,
Ratcliffe
,
R.G.
and
Sweetlove
,
L.J.
(
2018
)
Computational analysis of the productivity potential of CAM
.
Nat. Plants
4
,
165
171
31
Von Caemmerer,
S.
(
2000
)
Biochemical Models of Leaf Photosynthesis
,
CSIRO Publishing, Clayton, Victoria, Australia
32
de Oliveira Dal'Molin
,
C.G.
,
Quek
,
L.-E.
,
Palfreyman
,
R.W.
,
Brumbley
,
S.M.
and
Nielsen
,
L.K.
(
2010
)
C4GEM, a genome-scale metabolic model to study C4 plant metabolism
.
Plant Physiol.
154
,
1871
1885
33
Mallmann
,
J.
,
Heckmann
,
D.
,
Bräutigam
,
A.
,
Lercher
,
M.J.
,
Weber
,
A.P.M.
,
Westhoff
,
P.
et al (
2014
)
The role of photorespiration during the evolution of C4 photosynthesis in the genus Flaveria
.
Elife
3
,
e02478
34
Blätke
,
M.-A.
and
Bräutigam
,
A.
(
2019
)
Evolution of C4 photosynthesis predicted by constraint-based modelling
.
Elife
8
,
e49305
35
Gutteridge
,
S.
and
Pierce
,
J.
(
2006
)
A unified theory for the basis of the limitations of the primary reaction of photosynthetic CO2 fixation: was Dr. Pangloss right?
103
,
7203
7204
36
Töpfer
,
N.
,
Braam
,
T.
,
Shameer
,
S.
,
Ratcliffe
,
R.G.
and
Sweetlove
,
L.J.
(
2020
)
Alternative CAM modes provide environment-specific water-saving benefits in a leaf metabolic model
.
Plant Cell
32
,
3689
3705
37
Chandrasekaran
,
S.
and
Price
,
N.D.
(
2010
)
Probabilistic integrative modeling of genome-scale metabolic and regulatory networks in Escherichia coli and Mycobacterium tuberculosis
.
107
,
17845
17850
38
Shlomi
,
T.
,
Cabili
,
M.N.
,
Herrgård
,
M.J.
,
Palsson
,
B.Ø
and
Ruppin
,
E.
(
2008
)
Network-based prediction of human tissue-specific metabolism
.
Nat. Biotechnol.
26
,
1003
1010
39
Jensen
,
P.A.
and
Papin
,
J.A.
(
2011
)
Functional integration of a metabolic network model and expression data without arbitrary thresholding
.
Bioinformatics
27
,
541
547
40
Becker
,
S.A.
and
Palsson
,
B.O.
(
2008
)
Context-specific metabolic networks are consistent with experiments
.
PLoS Comput. Biol.
4
,
e1000082
41
Töpfer
,
N.
,
Jozefczuk
,
S.
and
Nikoloski
,
Z.
(
2012
)
Integration of time-resolved transcriptomics data with flux-based methods reveals stress-induced metabolic adaptation in Escherichia coli
.
BMC Syst. Biol.
6
,
148
42
Colijn
,
C.
,
Brandes
,
A.
,
Zucker
,
J.
,
Lun
,
D.S.
,
Weiner
,
B.
,
Farhat
,
M.R.
et al (
2009
)
Interpreting expression data with metabolic flux models: predicting Mycobacterium tuberculosis mycolic acid production
.
PLoS Comput. Biol.
5
,
e1000489
43
Mintz-Oron
,
S.
,
Meir
,
S.
,
Malitsky
,
S.
,
Ruppin
,
E.
,
Aharoni
,
A.
and
Shlomi
,
T.
(
2012
)
Reconstruction of Arabidopsis metabolic network models accounting for subcellular compartmentalization and tissue-specificity
.
109
,
339
344
44
Töpfer
,
N.
and
Niokoloski
,
Z.
(
2013
)
Large-scale modeling provides insights into Arabidopsis's acclimation to changing light and temperature conditions
.
Plant Signal. Behav.
8
,
e25480
45
Caldana
,
C.
,
Degenkolbe
,
T.
,
,
A.
,
Klie
,
S.
,
Sulpice
,
R.
,
Leisse
,
A.
et al (
2011
)
High-density kinetic analysis of the metabolomic and transcriptomic response of Arabidopsis to eight environmental conditions
.
Plant J.
67
,
869
884
46
Lakshmanan
,
M.
,
Lim
,
S.-H.
,
Mohanty
,
B.
,
Kim
,
J.K.
,
Ha
,
S.-H.
and
Lee
,
D.-Y.
(
2015
)
Unraveling the light-specific metabolic and regulatory signatures of rice through combined in silico modeling and multi-omics analysis
.
Plant Physiol.
169
,
3002
3020
47
Jung
,
E.S.
,
Lee
,
S.
,
Lim
,
S.-H.
,
Ha
,
S.-H.
,
Liu
,
K.-H.
and
Lee
,
C.H.
(
2013
)
Metabolite profiling of the short-term responses of rice leaves (Oryza sativa cv. Ilmi) cultivated under different LED lights and its correlations with antioxidant activities
.
Plant Sci.
210
,
61
69
48
Bogart
,
E.
and
Myers
,
C.R.
(
2016
)
Multiscale metabolic modeling of C4 plants: connecting nonlinear genome-scale models to leaf-scale metabolism in developing maize leaves
.
PLoS One
11
,
e0151722
49
Nägele
,
T.
and
Weckwerth
,
W.
(
2013
)
A workflow for mathematical modeling of subcellular metabolic pathways in leaf metabolism of Arabidopsis thaliana
.
Front. Plant Sci.
4
,
541
50
Töpfer
,
N.
,
Scossa
,
F.
,
Fernie
,
A.
and
Nikoloski
,
Z.
(
2014
)
Variability of metabolite levels is linked to differential metabolic pathways in Arabidopsis's responses to abiotic stresses
.
PLoS Comput. Biol.
10
,
e1003656
51
Sajitz-Hermstein
,
M.
,
Töpfer
,
N.
,
Kleessen
,
S.
,
Fernie
,
A.R.
and
Nikoloski
,
Z.
(
2016
)
iReMet-flux: constraint-based approach for integrating relative metabolite levels into a stoichiometric metabolic models
.
Bioinformatics
32
,
i755
i762
52
Timm
,
S.
,
Mielewczik
,
M.
,
Florian
,
A.
,
Frankenbach
,
S.
,
Dreissen
,
A.
,
Hocken
,
N.
et al (
2012
)
High-to-low CO2 acclimation reveals plasticity of the photorespiratory pathway and indicates regulatory links to cellular metabolism of Arabidopsis
.
PLoS One
7
,
e42809
53
Pfau
,
T.
,
Christian
,
N.
,
Masakapalli
,
S.K.
,
Sweetlove
,
L.J.
,
Poolman
,
M.G.
and
Ebenhöh
,
O.
(
2018
)
The intertwined metabolism during symbiotic nitrogen fixation elucidated by metabolic modelling
.
Sci. Rep.
8
,
12504
54
diCenzo
,
G.C.
,
Tesi
,
M.
,
Pfau
,
T.
,
Mengoni
,
A.
and
Fondi
,
M.
(
2020
)
Genome-scale metabolic reconstruction of the symbiosis between a leguminous plant and a nitrogen-fixing bacterium
.
Nat. Commun.
11
,
2574
55
Restrepo
,
S.
,
Myers
,
K.L.
,
del Pozo
,
O.
,
Martin
,
G.B.
,
Hart
,
A.L.
,
Buell
,
C.R.
et al (
2005
)
Gene profiling of a compatible interaction between Phytophthora infestans and Solanum tuberosum suggests a role for carbonic anhydrase
.
Mol. Plant Microbe Interact.
18
,
913
922
56
Botero
,
K.
,
Restrepo
,
S.
and
Pinzón
,
A.
(
2018
)
A genome-scale metabolic model of potato late blight suggests a photosynthesis suppression mechanism
.
BMC Genomics
19
,
863
57
Gyetvai
,
G.
,
Sønderkær
,
M.
,
Göbel
,
U.
,
Basekow
,
R.
,
Ballvora
,
A.
,
Imhoff
,
M.
et al (
2012
)
The transcriptome of compatible and incompatible interactions of potato (Solanum tuberosum) with Phytophthora infestans revealed by DeepSAGE analysis
.
PLoS One
7
,
e31526
58
Rodenburg
,
S.Y.A.
,
Seidl
,
M.F.
,
de Ridder
,
D.
and
Govers
,
F.
(
2018
)
Genome-wide characterization of Phytophthora infestans metabolism: a systems biology approach
.
Mol. Plant Pathol.
19
,
1403
1413
59
Rodenburg
,
S.Y.A.
,
Seidl
,
M.F.
,
Judelson
,
H.S.
,
Vu
,
A.L.
,
Govers
,
F.
and
de Ridder
,
D.
(
2019
)
Metabolic model of the Phytophthora infestans-tomato interaction reveals metabolic switches during host colonization
.
mBio
10
,
e00454-19
60
Töpfer
,
N.
,
Seaver
,
S.M.D.
and
Aharoni
,
A.
(
2018
)
Integration of plant metabolomics data with metabolic networks: progresses and challenges
.
Methods Mol. Biol.
1778
,
297
310
61
Benes
,
B.
,
Guan
,
K.
,
Lang
,
M.
,
Long
,
S.P.
,
Lynch
,
J.P.
,
Marshall-Colón
,
A.
et al (
2020
)
Multiscale computational models can guide experimentation and targeted measurements for crop improvement
.
Plant J.
103
,
21
31
62
Marshall-Colon
,
A.
,
Long
,
S.P.
,
Allen
,
D.K.
,
Allen
,
G.
,
Beard
,
D.A.
,
Benes
,
B.
et al (
2017
)
Crops in silico: generating virtual crops using an integrative and multi-scale modeling platform
.
Front. Plant Sci.
8
,
786
63
Chew
,
Y.H.
,
Seaton
,
D.D.
and
Millar
,
A.J.
(
2017
)
Multi-scale modelling to synergise plant systems biology and crop science
.
Field Crops Res.
202
,
77
83
64
Töpfer
,
N.
,
Caldana
,
C.
,
Grimbs
,
S.
,
Willmitzer
,
L.
,
Fernie
,
A.R.
and
Nikoloski
,
Z.
(
2013
)
Integration of genome-scale modeling and transcript profiling reveals metabolic pathways underlying light and temperature acclimation in Arabidopsis
.
Plant Cell
25
,
1197
1211
65
Liu
,
L.
,
Shen
,
F.
,
Xin
,
C.
and
Wang
,
Z.
(
2016
)
Multi-scale modeling of Arabidopsis thaliana response to different CO2 conditions: From gene expression to metabolic flux
.
J. Integr. Plant Biol.
58
,
2
1l