key: cord-013503-tjotro5h
authors: Herrmann, Helena A.; Dyson, Beth C.; Vass, Lucy; Johnson, Giles N.; Schwartz, Jean-Marc
title: Flux sampling is a powerful tool to study metabolism under changing environmental conditions
date: 2019-09-02
journal: NPJ Syst Biol Appl
DOI: 10.1038/s41540-019-0109-0
sha: 
doc_id: 13503
cord_uid: tjotro5h

The development of high-throughput ‘omic techniques has sparked a rising interest in genome-scale metabolic models, with applications ranging from disease diagnostics to crop adaptation. Efficient and accurate methods are required to analyze large metabolic networks. Flux sampling can be used to explore the feasible flux solutions in metabolic networks by generating probability distributions of steady-state reaction fluxes. Unlike other methods, flux sampling can be used without assuming a particular cellular objective. We have undertaken a rigorous comparison of several sampling algorithms and concluded that the coordinate hit-and-run with rounding (CHRR) algorithm is the most efficient based on both run-time and multiple convergence diagnostics. We demonstrate the power of CHRR by using it to study the metabolic changes that underlie photosynthetic acclimation to cold of Arabidopsis thaliana plant leaves. In combination with experimental measurements, we show how the regulated interplay between diurnal starch and organic acid accumulation defines the plant acclimation process. We confirm fumarate accumulation as a requirement for cold acclimation and further predict γ–aminobutyric acid to have a key role in metabolic signaling under cold conditions. These results demonstrate how flux sampling can be used to analyze the feasible flux solutions across changing environmental conditions, whereas eliminating the need to make assumptions which introduce observer bias.

High-throughput technologies have resulted in a rapid increase in available 'omic data sets. 1 Large-scale metabolic networks constructed using these data integrate known and predicted metabolic pathways. 2 These large-scale networks can be constrained using experimental data and the system behavior can be analyzed using metabolic modeling. Metabolism describes a cellular phenotype under given conditions, and changes in metabolite concentrations and reaction fluxes can be used to assess a cellular response to changing environmental conditions.

With the existence of large-scale metabolic networks comes the need to have appropriate modeling techniques available for their analysis. The majority of techniques for analyzing large-scale metabolic networks fall within the paradigm of constraint-based modeling (CBM). 3 CBM imposes stoichiometric constraints on the metabolic reactions and analyzes the possible flux solutions at steady state. Because the system is assumed to be at steady state, even genome-scale models can be solved at little computational expense.

Two of the most widely used forms of CBM are flux balance analysis (FBA) and flux variability analysis (FVA). 4 The key feature of FBA and FVA is that they compute the steady state of a model using an objective function. The objective function defines a reaction that is to be maximized or minimized when solving the system under the set constraints. A typical objective is "maximum biomass production", whereby essential macromolecules such as proteins and lipids are defined at known ratios or quantities in an outgoing reaction of the metabolic system. 5, 6 FBA computes single steady-state solutions, which satisfy the objective. However, often multiple solutions exist and their range can be computed using FVA. 7 FVA provides no indication as to whether all singlepoint solutions within the range are feasible and which solutions are the most likely. In order to reduce the number of feasible solutions of FBA and to further constrain the feasible flux ranges returned by FVA, it is common practice to introduce multiple objective functions. 8 However, defining one or multiple objective function(s) intrinsically introduces an observer bias as to what the main "goal" of the cell is, in the context of the analysis. 9 Although biomass production as estimated by FBA was seen to match experimental data in Escherichia coli, 10 this may not be an appropriate objective when studying short-term environmental changes. In the green algae Chlamydomonas reinhardtii, a 13 Cmetabolic flux analysis was shown to contradict the assumptions of a prior FBA analysis by demonstrating that maximum biomass and maximum ATP production cannot stand alone as cellular objectives. 11 Optimal growth conditions, to which most objective functions are tailored, are an exception in natural environments. 5 Evolutionarily, metabolism is most likely optimized for overall robustness across many conditions, rather than a single conditionspecific objective. Bacillus subtilis mutants, which outperform the wild-type in terms of biomass production in control conditions have been found experimentally; the wild-type, however, is more robust to both environmental and genetic perturbations and therefore holds an evolutionary advantage. 12 Understanding the optimization process that allows an organism to tolerate changing environmental conditions is of particular interest in crop sciences, where increasing yields will need to be achieved, despite rapidly changing environmental conditions. Owing to their sessile nature, plants are exposed to frequent and sometimes extreme environmental fluctuations. Plant metabolism can therefore be presumed to hold an inherent robustness to changing environmental conditions. Furthermore, plant metabolism is highly intricate as it includes multi-cellular autotrophic and heterotrophic tissues with complex cellular compartmentalization. Analyses of plant metabolism and of plant metabolic strategies can therefore be considered a comprehensive case study of metabolism in general. Genetic modifications do not always produce a desired effect owing to network robustness: 13, 14 in crop sciences, optimization for increased yield in control conditions rather than in changing environments is a likely explanation for discrepancies between laboratory and in field results. 15, 16 If we wish to use CBM techniques to assess network robustness and phenotypic plasticity, we must be able to capture all alternative solutions and the probability with which they occur (the solutions space) of a metabolic network across different conditions. To calculate the exact properties of the solution space, we can use mathematical techniques such as convex analysis and vertex enumeration; 17 however, owing to their computational intensiveness, these methods are only efficient when applied to small, simple networks 18 or genome-scale models with one or more objective function. 19 When wanting to analyze genomescale metabolic networks without the use of an objective function, such approaches cannot practically be applied. Sampling of feasible flux solutions provides a realistic alternative.

Flux sampling generates a sequence of feasible solutions (called a chain) that satisfy the network constraints, until the entire solution space is analyzed. Enough samples need to be generated for the samples to provide an accurate representation of the feasible solution space. 18 A chain of samples is said to have converged once it contains enough samples to give an accurate representation of the solution space. 20 Flux sampling provides information both on the range of feasible flux solutions (similar to FVA) but also on their probability. Importantly, unlike FBA, flux sampling does not require (but also does not exclude the option for) an objective function to be specified. Therefore, flux sampling methods hold great potential for analyzing optimization strategies that are not defined by clear objectives such as a simple biomass reactions.

When plants are exposed to a change in environmental conditions, such as temperature, which last only for a few days, optimizing biomass production during this period is arguably of secondary importance; sustaining metabolic function with minimal cost may be a higher priority to plants. For example, the allocation of carbon into different transient storage compounds accumulated to maintain cellular processes, has been shown to change when plants are exposed to environmental stresses. [21] [22] [23] [24] [25] Starch, malate, and fumarate are the three major carbon storage compounds which accumulate during the day in leaves of the model plant Arabidopsis thaliana. 21, 26, 27 Increased cytosolic fumarate accumulation is a known cold response of A. thaliana leaves and has been linked to an increased photosynthetic capacity that sustains metabolism in cooler temperatures. 24 Evidently, tight regulation of carbon partitioning is required for successful cold acclimation. Given the large number of reactions and pathways involved in linking primary carbon assimilation to its downstream storage products starch, malate, and fumarate, CBM seems appropriate. Of the CBM methods available, flux sampling allows us to gain a detailed understanding of the solution space and the interdependence of the different carbon stores under different temperature conditions, without imposing the constraint of an objective function.

Multiple large-scale metabolic networks of the model plant A. thaliana have been constructed. [28] [29] [30] [31] [32] Here, we used three of them to formally assess for the efficiency of existing and easily accessible flux sampling algorithms: the coordinate hit-and-run with rounding (CHRR), 33 the artificially centered hit-and-run (ACHR) 34 and the optimized general parallel (OPTGP) 35 algorithms. We identified the most efficient sampling method based on runtime and convergence, and applied it to study plant acclimation to cold. We experimentally measured diurnal CO2 uptake and organic carbon accumulation of A. thaliana in control and cold conditions. By constraining a leaf metabolic model to the two conditions and using an appropriate flux sampling algorithm, we were able to explore inherent metabolic robustness to temperature and predict the metabolic changes required to support a photosynthetic acclimation response to cold.

Although flux sampling has previously been applied as a technique for studying the solution space of metabolic networks, 29, 36, 37 this will, to our knowledge, be the first time that the available algorithms are formally compared with one another in the context of metabolic modeling, and that flux sampling is applied to study network robustness across changing environmental conditions.

Both run-time and convergence are fastest when using CHRR in MATLAB We compared the run-time and convergence of the CHRR, ACHR, and OPTGP algorithms using three metabolic models of A. thaliana. We tested the run-times of 500,000, 5,000,000, and 50,000,000 samples (S), of which 5000 were stored and the rest were discarded with constant measures of thinning. For S = 50,000,000, the CHRR algorithm was 2.5 times faster than the OPTGP and 5.3 times faster than the ACHR for the Arnold model (Fig. 2) . This difference in speed increases with model complexity, such that, for the Poolman model, the CHRR was 3.3 times faster than the OPTGP and 8.0 times faster than the ACHR (Fig. 2) . The OPTGP was run in two parallel processes; however, even when running it as a single process it is faster than the ACHR. Although we cannot exclude the fact that MATLAB may have a faster connection to Gurobi than Python, flux sampling is fastest when using the CHRR setup as available in the COBRA toolbox for MATLAB.

The number of reactions that did not satisfied the convergence criteria were assessed for each chain ( Table 1 ). The longer runtimes of ACHR and OPTGP implementations in Python are not outweighed by faster convergence. In fact, as the three pilot chains with varying thinnings show, the CHRR algorithm, across all model reactions, converges the fastest, with the lowest number of samples required for convergence, the least amount of autocorrelation, and the lowest discrepancy between chains (Table 1) . We confirmed this difference in convergence by inspecting trace and auto-correlations plots of individual reactions, such as those for the biomass reaction of the Arnold model shown in Fig. 2 (C). CHRR shows little dependence between consecutive samples even with a thinning of T = 100, whereas OPTGP and ACHR show low levels of autocorrelation only when T = 10,000.

Our results show differences in the outcomes when convergence is reached according to the different convergence diagnostics. All convergence diagnostics agree that the CHRR performs best (Table 1 ). According to the Raftery & Lewis and the IPSRF diagnostics, all of the flux samples of reactions in the Arnold and Poolman models converge in < 5000 samples with a thinning constant of 10,000 when using CHRR. The fact that such large numbers of samples are required for model convergence shows that, owing to the irregular solution shape of genome-scale metabolic networks, autocorrelation in chains is a common problem and must be overcome and tested for using appropriate convergence diagnostics. Analyzing samples that have not achieved convergence can lead to incorrect conclusions about the metabolic fluxes under study. Currently, many applications of the ACHR, OPTGP, and CHRR algorithms to biological networks do not report whether convergence has actually been achieved. 29, [37] [38] [39] Previous comparisons between the ACHR and CHRR algorithms 33 have been made using 15 different metabolic models but were based only on a single convergence diagnostic, the potential scale reduction factor (PSRF). 40 Notably, different convergence diagnostics test different features and may not always be in good agreement. Therefore, more than one diagnostic should be used to confirm that the sampling chain is likely to have reached convergence. 20, 41 The ACHR and OPTGP algorithms have been compared using five genome-scale metabolic models and three different convergence diagnostics, including the PSRF. 35 However, the PSRF assumes a normal distribution of solutions, 42 which is questionable given that sampling is most needed when the distribution of flux estimates is non-normal, 41 as is often the case in metabolic modeling.

CHRR-based flux sampling generates verifiable hypotheses concerning plant acclimation to cold When plants of A. thaliana are transferred from 20°C to 4°C, photosynthesis is inhibited. 24 Metabolism slows down in cooler temperatures and, in order to sustain normal metabolic functions, plants need to acclimate their CO 2 uptake, altering the concentrations of metabolic enzymes to achieve a new steady state. 24, 43, 44 After 7 days of cold, we observe that the A. thaliana wild-type Col-0 is able to achieve the same level of photosynthesis as measured in control conditions (Fig. 3) . The allocation of carbon to the three main carbon storage compounds, malate, fumarate, and starch, shifts in the cold as part of the metabolic acclimation response. After 7 days of acclimation both photosynthesis and transient carbon accumulation attain a new cold acclimated state. Most notably, after 7 days of cold treatment a larger proportion of carbon is partitioned into fumarate (Fig. 3 ).

Using the experimental data shown in Fig. 3 to constrain the CO 2 input and the malate, fumarate, and starch accumulation reactions using the Arnold model (please see methods for further details), we were able to compute converged flux sampling distributions for all reactions. We did so for both control and cold conditions, which allowed us to overlay the sampling distributions of reaction fluxes and to assess changes required in plant metabolic behavior for acclimation (Fig. 4) .

In order to demonstrate how the application of an objective function for an FBA analysis can lead to vastly different conclusions, we have overlaid FBA results for maximum biomass production (under the same model constraints as applied for the sampling) over the flux sampling distributions (Fig. 4) . This further emphasizes how an objective function, if inappropriate for the analyses under consideration, can be misleading.

Both sucrose export to other tissues and cytosolic pyruvate production as a precursor for the tricarboxylic acid (TCA) cycle are predicted to be unchanged. The model suggests that, given equal carbon assimilation in control and cold conditions, cellular maintenance and export functions are supported equally in both conditions on the time scale considered. Although sucrose export is difficult to measure experimentally, we observed the rate of respiration in the cold to be the same as in control conditions (Fig. 3) . Given that, for both conditions, the model shows equal fluxes from cytosolic pyruvate into the TCA cycle (Fig. 4) , which feeds directly into respiration, this model prediction is in agreement with our experimental data.

For fumarate (and other metabolites not shown) the model predicts a shift in cellular compartmentalization with temperature. Model results show an increase in fumarate export from the mitochondrion into the cytosol in the cold. The reverse is shown for fumarate export from the chloroplast, where it is produced via the breakdown of arginosuccinate (Fig. 4) . Leaves developed in the cold have increased cytoplasmic and decreased vacuolar volumes; 45 a reshuffling of metabolites across cellular compartments has therefore previously been proposed as an important temperature acclimation response. 46 Model results suggest an increase in flux from cytosolic malate to fumarate via cytosolic fumarase (FUM2) (Fig. 4) . This is consistent with previous experimental results from mutant studies 24 This reaction is thus essential for photosynthetic acclimation to cold.

Flux sampling distributions suggest a link between carbon and nitrogen metabolism to support cold acclimation The sampling distributions suggest a trade off between increased carbon compound accumulation and decreased amino-acid production (Fig. 4) , linking nitrogen and carbon metabolism. Synthesis of γ-aminobutyric acid (GABA), however, is predicted to increase in the cold. GABA has previously been reported to accumulate in response to environmental stresses, including cold treatment. 47 GABA has been suggested as a signaling molecule of the carbon to nitrogen status in plant leaves and evidence for its role in regulating nitrate uptake exists for both rapeseed and A. thaliana. [48] [49] [50] A. thaliana plants in the cold show increased nitrogen assimilation compared with those in control conditions. 51 If GABA is indeed involved in carbon:nitrogen signaling, 48 it may be counteracted by the increased accumulation of malate. Increased levels of malate have previously been shown to suppress nitrate reductase expression and activity in tobacco leaves. 52 Malate levels in Col-0 may be kept below a certain threshold, by redirecting carbon to fumarate, thereby supporting adequate nitrogen assimilation and an increased photosynthetic capacity. This hypothesis is supported by the observation that A. thaliana mutants, which show increased levels of malate and decreased levels of fumarate, grow significantly less well in highnitrogen conditions than Col-0. 22 Fumarate has few known metabolic functions in A. thaliana leaves, 26 and may thus serve the purpose of a carbon storage buffer in changing environmental conditions.

Based on run-time and platform, CHRR is faster than OPTGP, which is itself faster than ACHR. CHRR also converges faster than both OPTGP and ACHR. Users with unrestricted access to MATLAB are therefore recommended to use CHRR. For those who wish to work using an open-source platform, OPTGP is recommended over ACHR; in general, OPTGP converges faster than ACHR, has a shorter run-time and allows for parallel processes.

When running sampling algorithms, sets of flux samples are produced for each reaction in the model. Here, we have tested convergence for all reactions of the models using three different sample chains. We have highlighted the importance of checking for convergence using different diagnostics when analyzing an irregular solution space of large networks. If only a subset of reaction fluxes are of interest, only those distributions will have to be checked for convergence. Flux sampling provides a powerful tool for exploratory analyses assessing metabolic differences across different environmental conditions. Our results further confirm the notion that it is not possible to fully automate convergence analyses using a single diagnostic 20 and that results should confirmed via manual inspections of trace and autocorrelation plots.

Flux sampling is currently an under-utilized technique in the metabolic modeling of large-scale networks. Using cold acclimation of the model plant A. thaliana, we demonstrate how flux sampling can be used effectively to analyze alternative feasible solutions across multiple conditions whilst eliminating the need to make assumptions that introduce observer bias. Given short-term environmental changes, adequately sustaining basic metabolic functions with minimal resource investment may be a more-likely cellular objective on that time scale than, for example, maximizing growth. We therefore did not set an objective function for flux sampling but used four experimentally measured flux values (CO 2 input and fumarate, malate and starch accumulation) to constrain a leaf metabolic model. 28 We further demonstrated how these flux sampling results can lead to different conclusion than traditional FBA analyses. By constraining the model to both cold and control conditions we were able to select reactions that show different flux distributions across the two conditions. Our model highlights reactions that are essential to change with temperature (i.e., the flux distribution of the two temperature conditions do not overlap) such as the production of cytosolic fumarate via malate. 24 The model further demonstrates the properties of the flux distributions of GABA to differ in cold and control conditions, highlighting GABA as a plausible signaling molecule for supporting a shift in the nitrogen and carbon balance, required to sustain photosynthesis in the cold. Thus, through flux sampling, we were able to generate novel hypotheses about the roles of GABA, fumarate and malate in cold acclimation, which would have been unfeasible to detect using FBA and FVA methods.

By overlaying different FBA solutions onto flux sampling distributions obtained under condition-specific model constraints, FBA in combination with flux sampling, could, in future work, be used to determine plausible objective functions and help generate predictions about how cellular objectives might be changing in response to environmental changes.

Constraint-based reconstruction and analysis (COBRA) methods for genome-scale metabolic networks are integrated in the COBRA toolbox 53, 54 for the MATLAB programming language and the COBRApy package 55 for the open-source Python programming language. Three algorithms for flux sampling exist across the two platforms: CHRR (MATLAB), ACHR (MATLAB and Python), and OPTGP (Python). Further flux-sampling algorithms exist; 37,56-58 however, as they are not currently available in the COBRA packages, they are here not considered for comparison. 59, 60 1. The artificially centered hit-and-run (ACHR) sampler estimates the center of the solution space in a "warm-up" phase. This estimate is then continuously revised with further sampling. The center estimate is used to inform the direction of further sampling such that the full solution space is covered in fewer steps than in traditional hit-and-run sampling (where the direction of the next sample is chosen at random). 34 Although the Markovian nature of hit-and-run (i.e., the fact that each future sampling state is dependent only on the current sampling state) is lost in the ACHR, it overcomes the edge-trapping limitation of the standard hit-andrun algorithm (i.e., it no longer gets stuck at the bounds of a solution space if these are of an elongated shape, a frequent feature in metabolic models). 2. The optimized general parallel sampler (OPTGP) is argued to be an improvement on the ACHR, because from the warm-up point it generates multiple short chains from the estimated center and only considers the last point in the chain as a sample. 35 It thereby increases the randomness and efficiency with which the total solution space is explored. Furthermore, it allows for parallel sampling. Larger samples can thus be generated in shorter runtimes. 3. The coordinate hit-and-run with rounding (CHRR) algorithm starts with a pre-processing step that rounds the solution space to a more regular, convex shape, and therefore a Markov chain can be used to explore the rounded solution space without the limitation of edgetrapping. After sampling, the solutions are back-transformed to match the original solution space in order to obtain the true value of the sampled points. 33, 61 Assessing convergence

The aim of flux sampling is to generate enough consecutive samples (a long enough chain) in order to get an accurate depiction of the solution space. A sample chain is considered to have converged once it can be assumed that the sampled subset of solutions represent the properties of H.A. Herrmann et al.

the true solutions obtained from an infinite amount of samples (i.e., when the shape of the flux distribution no longer changes with more samples). 62 In flux sampling, an algorithm's efficiency is defined by the run-time and the number of samples required for apparent convergence to the true flux distribution. Here, we apply three different diagnostics in order to test for convergence:

1. Raftery & Lewis diagnostic: this diagnostic estimates the total number of samples, N max , required for a set of samples to achieve convergence, based on a given subset of samples (pilot chain). 63 We estimated N max required for convergence using the CODA package in R. 64 The diagnostic further returns a dependence factor, I, which is indicative of autocorrelation (the degree of dependence between consecutive samples). Chains with I > 5 are here considered to be problematic as this value suggests high dependence between samples or influential starting values (i.e., the chain was not run long enough). 2. Interval-based potential scale reduction factor (IPSRF): adapted from the original Gelman-Rubin diagnostic, PSRF, 40 the IPSRF is based on comparing the differences between consecutive samples (within sequence interval length) with the total differences observed between all samples (between sequence interval length). Because interval length, rather than variance of the samples (on which the original version is based), is considered, a normal distribution of the samples is no longer a requirement. 42 As with PSRF, the IPSRF Fig. 1 Summary of the general influx and outflux of the Arnold model, as set up in our analysis. Differences in pathways linking carbon assimilation and diurnal carbon storage were compared across the two temperatures (20°C and 4°C) in order to predict metabolic changes required for cold acclimation of A. thaliana 55 The linear programming solver used was Gurobi, version 8.0.0. The samplers were used in accordance with COBRA documentation using default parameter settings unless otherwise specified. Although an ACHR implementation is available both in the COBRA toolbox and in COBRApy, the comparisons made here are based on the Python implementation of the ACHR sampler because it is open-source and because the OPTGP, which is available only in Python, directly builds on the general parallel sampler (GP) of the ACHR. We did collect preliminary convergence and run-time data of the ACHR algorithm in MATLAB; however, because its convergence was evidently slower than CHRR, we have chosen to omit its MATLAB implementation from further analyses. OPTGP and CHRR algorithms were run using two parallel processes; there is currently no option to run ACHR as multiple processes.

To compare the performance of the samplers, three published genomescale A. thaliana metabolic models were obtained in SBML format and are from here on referred to by their first author: the Poolman, 32 Arnold, 28 and Dal'Molin 2 models. The Poolman model is based on the Aracyc database 66 and describes a non-compartmentalized heterotrophic culture using 1406 reactions. The Arnold model describes a compartmentalized, photoautotrophic system and is based on 549 manually curated A. thaliana specific reactions. The Dal'Molin model is a compartmentalized network reconstruction of 1601 reactions, applicable to both photosynthetic and non-photosynthetic tissues of plant metabolism. 2 The original model constraints 2,28,32 were used when comparing sampler performance across these models, such that the Arnold, Poolman, and Dal'Molin model had 270, 645, and 330 degrees of freedom respectively.

Samples that are close together within a sample chain can be autocorrelated (i.e., be similar to one another due to the way in which the algorithm works). In order to avoid the effects of autocorrelation and to ensure convergence, a large number of samples should be run and a technique called thinning can be applied. 67 Sample chains that are thinned store only every k th sample in the chain, where T = k is called the thinning constant. In order to compare autocorrelation of the three algorithms, we applied three different thinning constants (T = 100, 1000, 10,000) storing 5000 samples for each chain. Three replicate chains were run for each value of T. Convergence diagnostics were calculated for sample chains produced for each of the model reactions. Run-times were obtained on a personal laptop (7.2 GB RAM, i5-6200U CPU processor, 4 cores, 2.3 GHz capacity, Ubuntu 16.04.4 OS).

Carbon assimilation and respiration (CO 2 influx and outflux) by the A. thaliana wild-type Columbia-0 (Col-0) was measured using infrared gas analysis at 100 μmol m −2 s −1 light under control conditions (20°C) and after 1 week of cold treatment (4°C) as described previously. 24 Fumarate, malate and starch concentrations at the onset and at the end of the photoperiod were measured for both control and cold conditions using enzyme assays. 24 Averages of 3-4 replicated were calculated. Measurements taken across the different conditions were tested for significant differences using an unpaired t test (p < 0.05) assuming unequal variances, as implemented in R.

Given that previously published data confirm an approximately constant rate of accumulation of transient carbon storage products during the day, 21, 24 we subtracted the beginning of day metabolite concentrations from the end of day concentrations for each of the products in order to obtain a flux value for carbon storage over one photoperiod. These flux values were used to constrain the Arnold model (Fig. 1) . This reduced the degrees of freedom of the model to 266. Measurements were taken at control conditions (red) and after 7 days of cold treatment (blue), at the beginning of day (BOD) and at the end of day (EOD), using infrared gas analysis and enzyme assays. 24 The s.e.m. of the 3-4 replicates for each measurement is shown as error bars. The above data (excluding carbon outflux) were converted to mmol gFWÁDay and used to constrain the Arnold model as outlined in the method. Significant differences between measurements across the two temperature conditions are indicated by an asterisk

Given that the Arnold model is leaf specific and manually curated, it suits the purpose of studying photosynthetic plant acclimation. In order to do so we constrained the model with our experimental diurnal flux data for malate, fumarate, and starch accumulation as well as CO 2 influx. We set the cytosolic fumarase reaction, which produces cytosolic fumarate from malate, to be reversible. 24 Outgoing fumarate, malate, and starch reactions were added to the model in order to simulate diurnal carbon storage (Fig.  1 ). Diurnal accumulation of the metabolites was calculated by subtracting average beginning of day concentrations from average end of day concentration values. Metabolite accumulation and CO 2 influx were converted to mmol (gFW) −1 Day −1 . The resulting values were applied as model constraints. Upper and lower bounds were applied according to the calculated standard errors of three to four replicates as shown in Fig. 3 . To ensure convergence of all flux sampling distributions, 100,000 flux samples with a thinning of 10,000 were generated using the CHRR algorithm in the COBRA Toolbox. A Kruskal-Wallis test, as implemented in the SciPy Python package, version 0.19.1, was used to assess whether flux samples generated using either the cold or the control constrained model stemmed from the same distribution. 68 Reporting summary Further information on research design is available in the Nature Research Reporting Summary linked to this article.

The authors declare that all data supporting the findings of this study are available within the paper and online or can be generated using the publicly available source code (https://github.com/HAHerrmann/FluxSamplingPlantModels; https://doi.org/ 10.5281/zenodo.3239075).

Flux sampling distributions of key reactions linking photosynthetic input (CO 2 ) to transient carbon storage for control (red) and cold (blue) conditions. Condition-specific carbon assimilation and fumarate (FUM), malate (MAL), and starch accumulation were constrained according to experimentally measured results. Resulting fluxes of reactions including triose phosphate (TP), fructose 1,6bisphosphate (F6P), glyceraldehyde 3-phosphate (G3P), pyruvate (PYR), arginosuccinate (AS), arginine (ARG), phosphoenolpyruvate carboxylase (PEP), oxaloacetate (OAA), γ-aminobutyric acid (GABA), acetyl coenzyme A (ACoA), and sucrose (Suc) are shown. We overlaid FBA results for maximum biomass production (under the same model constraints as applied for the sampling) as vertical blue and red bars over the flux sampling distributions. Reactions for which the two distributions are significantly different (p < 0.001; Kruskal-Wallis) are marked with an asterisk in the top right corner

Multi-omics and metabolic modelling pipelines: challenges and tools for systems microbiology

AraGEM, a genome-scale reconstruction of the primary metabolic network in Arabidopsis

Constraint-based models predict metabolic and associated cellular functions

What is flux balance analysis?

The biomass objective function

Flux balance analysis of plant metabolism: the effect of biomass composition and model structure on model predictions

Methods and advances in metabolic flux analysis: a minireview

A multi-objective constraintbased approach for modeling genome-scale microbial ecosystems

Comparison and analysis of objective functions in flux balance analysis

Stoichiometric flux balance models quantitatively predict growth and metabolic by-product secretion in wild-type Escherichia coli W3110

Metabolic flux analysis of heterotrophic growth in Chalmydomonas reinhardtii

Large-scale in vivo flux analysis shows rigidity and suboptimal performance of Bacillus subtilis metabolism

Biological robustness

Phenotypic plasticity and robustness: evolutionary stability theory, gene expression dynamics model, and laboratory experiments

Food for thought: lower-than-expected crop yield stimulation with rising CO 2 concentrations

Crop yield gaps: their importance, magnitudes, and causes

Use of randomized sampling for analysis of metabolic networks

Monte Carlo sampling can be used to determine the size and shape of the steady-state flux space

Interplay between constraints, objectives, and optimality for genome-scale stoichiometric models

Convergence assessment techniques for Markov chain Monte Carlo

Coordination of carbon supply and plant growth

Arabidopsis has a cytosolic fumarase required for the massive allocation of photosynthate into fumaric acid and for rapid plant growth on high nitrogen

Acclimation of metabolism to light in Arabidopsis thaliana: the glucose 6-phosphate/phosphate translocator GPT2 directs metabolic acclimation

FUM2, a cytosolic fumarase, is essential for acclimation to low temperature in Arabidopsis thaliana

Mathematical modeling of diurnal patterns of carbon allocation to shoot and root in Arabidopsis thaliana

Fumaric acid: an overlooked form of fixed carbon in Arabidopsis

Analysis of Arabidopsis with highly reduced levels of malate and fumarate sheds light on the role of these organic acids as storage molecules

Bottom-up reconstruction of Arabidopsis and its application to determining the metabolic costs of enzyme production

A multi-tissue genomescale metabolic modeling framework for the analysis of whole plant systems

A diel flux balance model captures interactions between light and dark metabolism during day-night cycles in C3 and crassulacean acid metabolism leaves

Reconstruction of Arabidopsis metabolic network models accounting for subcellular compartmentalization and tissue-specificity

A genome-scale metabolic model of Arabidopsis and some of its properties

CHRR: coordinate hit-and-run with rounding for uniform sampling of constraint-based models

Direction choice for accelerated convergence in hitand-run sampling

optGpSampler: an improved tool for uniformly sampling the solution-space of genome-scale metabolic networks

Uniform sampling of steady-state flux spaces: means to design experiments and to interpret enzymopathies

Sampling the solution space in genome-scale metabolic networks reveals transcriptional regulation in key enzymes

Connecting extracellular metabolomic measurements to intracellular flux states in yeast

Genome-scale metabolic modeling elucidates the role of proliferative adaptation in causing the Warburg effect

Bayesian data analysis

Markov chain monte carlo convergence diagnostics: sa comparative review

General methods for monitoring convergence of iterative simulations

Carbon partitioning and export in transgenic Arabidopsis thaliana with altered capacity for sucrose synthesis grown at low temperature: a role for metabolite transporters

Altering flux through the sucrose biosynthesis pathway in transgenic Arabidopsis thaliana modifies photosynthetic acclimation at low temperatures and the development of freezing tolerance

Acclimation of Arabidopsis leaves developing at low temperatures. Increasing cytoplamic volumes accompanies increased activities of enzymes in the Calvin cycle and in the sucrose-biosynthesis pathway

Approximating subcellular organisation of carbohydrate metabolism during cold acclimation in different natural accessions of Arabidopsis thaliana

Metabolism of γ-aminobutyric acid during cold acclimation and freezing and its relationship to frost tolerance in barley and wheat

Putative role of γ-aminobutyric acid (GABA) as a long-distance signal in up-regulation of nitrate uptake in Brassica napus L

Closing the loop on the GABA shunt in plants: are GABA metabolism and signaling entwined? Front

Nitrate uptake and utilization is modulated by exogenous γ-aminobutyric acid in Arabidopsis thaliana seedlings

Source of nitrogen associated with recovery of relative growth rate in Arabidopsis thaliana acclimated to sustained cold treatment

Influence of malate and 2-oxoglutarate on the NIA transcript level and nitrate reductase activity in tobacco leaves

Quantitative prediction of cellular metabolism with constraint-based models: the COBRA Toolbox v2.0

Creation and analysis of biochemical constraint-based models using the COBRA Toolbox v.3.0

COnstraintsbased reconstruction and analysis for python

ll-ACHRB: a scalable algorithm for sampling the feasible solution space of metabolic networks

Non-stationary forward flux sampling

An ensemble evolutionary constraint-based approach to understand the emergence of metabolic phenotypes

The RAVEN toolbox and its use for generating a genome-scale metabolic model for Penicillium chrysogenum

A metabolic core model elucidates how enhanced utilization of glucose and glutamine, with enhanced glutamine-dependent lactate production, promotes cancer cell growth: the WarburQ effect

Uniform sampling of steady states in metabolic networks: heterogeneous scales and rounding

Markov chain Monte Carlo: an introduction for epidemiologists

How many iterations in the Gibbs sampler

CODA: convergence diagnosis and output analysis for

Evaluating the accuracy of sampling-based approaches to calculating posterior moments

AraCyc: a biochemical pathway database for arabidopsis

Are We There Yet? When to Stop a Markov Chain while Generating Random Graphs. International Workshop on Algorithms and Models for the Web-Graph, WAW: Algorithms and Models for the Web Graph

Use of ranks in one-criterion variance analysis

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/.