Integration of relative metabolomics and transcriptomics time-course data in a metabolic model pinpoints effects of ribosome biogenesis defects on Arabidopsis thaliana metabolism

Ribosome biogenesis is tightly associated to plant metabolism due to the usage of ribosomes in the synthesis of proteins necessary to drive metabolic pathways. Given the central role of ribosome biogenesis in cell physiology, it is important to characterize the impact of different components involved in this process on plant metabolism. Double mutants of the Arabidopsis thaliana cytosolic 60S maturation factors REIL1 and REIL2 do not resume growth after shift to moderate 10 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\circ }\hbox {C}$$\end{document}∘C chilling conditions. To gain mechanistic insights into the metabolic effects of this ribosome biogenesis defect on metabolism, we developed TC-iReMet2, a constraint-based modelling approach that integrates relative metabolomics and transcriptomics time-course data to predict differential fluxes on a genome-scale level. We employed TC-iReMet2 with metabolomics and transcriptomics data from the Arabidopsis Columbia 0 wild type and the reil1-1 reil2-1 double mutant before and after cold shift. We identified reactions and pathways that are highly altered in a mutant relative to the wild type. These pathways include the Calvin–Benson cycle, photorespiration, gluconeogenesis, and glycolysis. Our findings also indicated differential NAD(P)/NAD(P)H ratios after cold shift. TC-iReMet2 allows for mechanistic hypothesis generation and interpretation of system biology experiments related to metabolic fluxes on a genome-scale level.


Scientific Reports
| (2021) 11:4787 | https://doi.org/10.1038/s41598-021-84114-y www.nature.com/scientificreports/ One possibility to investigate the feedback of ribosome biogenesis defects on metabolism is the characterization of reaction fluxes. Metabolic fluxes depend, in part, on the metabolite pools 3 . They also depend on the enzymatic setup of a cell, which is in turn governed by gene regulatory and signalling networks that affect protein activity 4 . However, determination of metabolic fluxes is a tedious and labour-intensive task [5][6][7] . A targeted analysis that predicts relevant fluxes for hypothesis generation based on integration of available high-throughput data sets from systems biology studies may streamline the planning of such time-consuming experimental flux studies.
In this regard, constraint-based approaches have proved as a valuable tool for hypotheses generation regarding flux distributions and their differential behaviour. For instance, the simplest of these approaches, flux balance analysis (FBA), can predict steady-state fluxes in bacteria at exponential growth 8 . In general, metabolic fluxes of a system are predicted under the assumption that this system operates in steady-state and optimizes an objective (e.g. biomass yield). If feasible, the resulting mathematical approach often results in a non-unique flux distribution. To this end, constraints defined through integration of high-throughput data can reduce the solution space of feasible flux distributions [9][10][11] . Such approaches have been shown to result in a more accurate prediction which is closer to the actual physiological state 12 . Despite the availability of methods that integrate high-throughput data, their full potential has yet to be realized 13 .
Of particular interest are approaches which allow integration of relative metabolite levels, since these datasets are easier to obtain in contrast to absolute metabolite concentrations used in thermodynamic flux balance analysis 14 as well as approaches that use time-series data (e.g. TREM-Flux 15 , uFBA 16 , and dFBA 17 ). iReMet-flux 18 is the only constraint-based approach to date that can integrate relative metabolite levels to investigate differential flux behaviour between two scenarios. It relies on a mass-action-like description of reaction rates (i.e. flux). In contrast to uFBA, iReMet-Flux does not require data on absolute quantification of metabolite levels and therefore allows for a broader application due to the availability of relative metabolomics data. In contrast to TREM-Flux, it does not assume a linear scaling with the change of metabolite levels between two time points. In addition, iReMet-Flux differs from a recent approach in which the relative metabolomics data are integrated on a qualitative level (i.e. increases or decreases) 14 . Similar to the objective on which MOMA is based 19 , iReMet-flux minimizes the flux differences between two scenarios, but does not rely on pre-calculated flux distributions for a reference scenario. Additionally, iReMet-flux allows for the integration of relative enzyme levels either by direct usage of quantitative or qualitative proteomics data, or via gene expression ratio that can serve as a proxy 10,20,21 . However, if employed to time-series data, it does not account for the magnitude of possible flux changes between time steps. To address this problem, we extended iReMet-flux to account for temporal changes, while keeping the possibility of multi-level high-throughput data integration.
Here, we aimed to develop a novel constraint-based approach, termed TC-iReMet2, that facilitates the integration of relative metabolite and transcript levels while accounting for temporal change of physiological parameters. We used TC-iReMet2 to investigate differential flux behaviour of A. thaliana Col-0 wild type and reil1-1 reil2-1 double mutant plants before and after cold shift. Finally, we provided directly testable hypotheses about the impact of REIL-mediated deficiency in ribosome biogenesis on metabolism.

Results
Formulation of TC-iReMet2. We propose Time Course Integration of Relative Metabolite and Transcript levels (TC-iReMet2) that estimates fluxes based on the integration of time-course data on relative metabolite and transcript levels. The key feature of TC-iReMet2 is that it accounts for the possible magnitude of flux changes between time points and thus could provide a more accurate explanation of flux rerouting over time. We show that TC-iReMet2 can be applied to study flux redistributions in pathways in a large-scale metabolic network of A. thaliana. Unlike genome-scale metabolic networks 22 , we refer to large-scale models as those reconstructed following a bottom-up approach 23 .
Similar to other constraint-based approaches, TC-iReMet2 uses a stoichiometric matrix S of the considered metabolic model. The rows of the stoichiometric matrix correspond to metabolites, and columns stand for reactions. The integer entries denote the molarity of a product (positive entry) or a substrate (negative entry) in a reaction, ensuring mass and charge conservation. In the following, we assume that the investigated metabolic network contains P reactions and n metabolites, and that its functioning is compared between two experimental scenarios, denoted by A and B (e.g. mutant and wild type) over to time points t + 1 and t. Furthermore, we denote by p 1 the number of irreversible reactions and by p − p 1 the number of reversible reactions.
Under mass action kinetics, a flux through an irreversible reaction i, 1 ≤ p 1 ≤ p 1 , can be formally described by v i = k i E i n j=1 (x j ) |S ji | , where x j denotes the concentration of metabolite j, S ji denotes the stoichiometric coefficient with which a metabolite j enters a reaction i as a substrate, E i denotes the enzyme concentration and k i denotes the reaction specific rate constant. Note that this expression can be written equally for scenario A: where the rate constant k i is the only unchanged parameter ( k A i = k B i ) -as it summarizes the key property of the same enzyme. Therefore, the relationship of a single flux between two scenarios can be written as: www.nature.com/scientificreports/ To simplify the notation, we will refer to the ratio of metabolite levels of j as r j =  Determining the entirety of metabolite and enzyme concentrations is not possible with the existing  technologies 24,25 . For metabolite ratios, only a small portion of the metabolome, and hence metabolite ratios, can be quantified. To account for the case that a metabolite ratio cannot be measured, general upper and lower boundaries for metabolite ratios are introduced. If the ratio of metabolite j is experimentally quantified, it is indicated by χ(r j ) = 1 and otherwise by χ(r j ) = 0.
In the absence of enzyme ratios, we use the Gene Protein Reaction (GPR) rules of metabolic models to approximate enzyme ratios using transcriptomic data. The GPR roles are defined by a set of Boolean expressions that describe which genes encode an enzyme. For example, gene products encoding for isoenzymes or isoforms are linked by an OR operator. Conversely, protein subunits that must be present simultaneously to form an active enzyme are linked by an AND operator. In case of an enzyme encoded by one gene, the enzyme concentration is approximated by its expression value. For each reaction that is catalyzed by a complex requiring multiple genes, the enzyme concentration is set to the minimum expression value of gene products connected by the AND operator. For the OR operator, the sum of expression values for the respective genes is used. These rules were applied to each reaction in both scenarios, fractioned and assigned as the corresponding enzyme ratio. Therefore, an enzyme ratio is represented by a ratio of gene expression levels following the GPR rules. Equivalently to metabolite ratios, if a GPR rule for reaction i is defined, it is indicated by H(q i ) = 1 and for reactions without a defined GPR rule, by H(q i ) = 0.
In this setup, we only consider constraints for irreversible reactions, since more than 80% of reactions that are assumed to follow mass-action-like kinetics (this excludes artificial and transport reactions) are irreversible in the analyzed model of A. thaliana. This has been verified by performing flux variability analysis at a fixed flux through the biomass reaction, to specify that 80% of reactions operate in only one direction 18 . A ratio constraint for reaction i is included if not only the enzyme ratio, but also at least one of the substrate ratios corresponding to that reaction is available. For metabolites or enzymes whose ratios could not be determined we use the extremal values found at that specific time point. Let F(i) denote the set of substrates of reaction i. Additionally, let the set of irreversible reactions with at least one experimentally quantified metabolite ratio and approximated enzyme ratio be denoted by max m for unmeasured transcript ratios. Furthermore, to account for enzymes that are substrate saturated and in turn would lead to infeasibilities due to metabolite ratio constraints, slack variables ε i were introduced to relax the strict ratio constraints. To minimize these relaxations a weighting of the summed slack variables of ǫ = 0.01 was used. Hence, a ratio constraint was formulated as follows: Similarly, a ratio constraint for the biomass reaction can be formulated. To this end, a time-point specific biomass fraction, denoted by κ t+1 , can be calculated. First, the maximum biomass yield, denoted by opt, is calculated for both scenarios via FBA. A biomass fraction κ t+1 between both scenarios is then determined by using proxies for biomass (for a detailed description see Methods -Parameterizing the objection and of TC-iReMet2 and estimating fractions of biomass yield). We fix the biomass reaction of scenario B to its respective value derived from FBA. In contrast, biomass flux in scenario A is fixed to a fraction κ t+1 of its optimal yield. Lower and upper bounds are specified as deviations, denoted by δ , of the calculated fraction. Therefore, biomass fluxes for both scenarios can be constrained as follows: Furthermore, we assume that: (i) the metabolic network to operate in quasi-steady state at every time point. Hence, Sv A = Sv B = 0 , where v A and v B denote the flux distributions of scenarios A and B respectively. (ii) the biological system aims to maintain an optimal state given by the enzymatic setup. This assumption is captured by making sure that the flux distributions between the two scenarios at a given time point t + 1 are as close as pos- (iii) the physiological state at time t + 1 depends on the physiological state at time t. We model this assumption by accounting for the magnitude of possible physiological changes by assuring that the difference of flux distributions between time points is as small as possible, i.e.
respectively. This magnitude obviously depends on the difference between time points, where the magnitude of possible flux changes increases with time. To this end, we introduce weighting factors to minimize the difference of flux distributions between scenarios at the current time point, weighted by α , as well as for differences between prior time points for scenario A, weighted by β , and scenario B, weighted by γ.
In summary, the TC-iReMet2 approach is cast as a quadratic program (QP) as follows:

Application of TC-iReMet2 to data from the reil1-1 reil2-1 A. thaliana mutant. We employed
TC-iReMet2 to gain insights into the metabolic effects of the ribosome biogenesis defect that is caused by A. thaliana REIL deficiency. To this end, we compared predicted flux differences between Col-0 wild type and reil1-1 reil2-1 double mutant with deficiency in cytosolic 60S ribosome biogenesis. The REIL proteins are required for growth when plants are shifted to cold ( < 10 • C ) conditions, but not at optimal temperature ( ≃20 • C) 2 . The reil1-1 reil2-1 double mutant and wild type differ only slightly in size when grown at 20 • C . Young developing leaves of the mutants showed an acute tip and two basal serrations instead of the typical rounded leaves of the Col-0 wild type, and were similar to the pointed leaves phenotype of cytosolic ribosome mutants [26][27][28] . However, the pointed-leaf phenotype of the reil1-1 reil2-1 double mutant was no longer apparent after transfer to soil and at developmental stages < 1.10 29 that were analyzed in this study. When shifted to 10 • C (cold), both the mutant and the wild type stopped growing. Following seven days in the cold, the wild type resumed growth, while the mutant remained strongly growth-inhibited ( Fig. 1, Supplementary Table S1). The mutant survived at least four weeks after cold shift and maintained cellular integrity as was determined by electrolyte leakage assays of rosette leaves 30 . Growth parameters of wild type and reil1-1 reil2-1 were determined as proxies of relative biomass accumulation at day 0, day 1, days 7 and 21 after cold shift using morphometric data (see Methods -parameterizing the objective function). Along with the morphometric data, the relative changes of metabolite pools and transcripts were profiled 30 (see "Methods" section for details). The experimental setup and the availability of transcriptomics data and data on relative metabolite levels allowed the application of TC-iReMet2 29,31 to quantify the nominal and relative differences in metabolic fluxes of the wild type and the mutant (Supplementary Fig. S1). We refer to nominal changes as the sum of predicted flux differences, defined as the absolute value of difference between wild type and mutant flux, over all analyzed time points. The nominal changes may provide a skewed picture about the differences, particularly since the differences in fluxes between reactions in a given flux distribution differ in several orders of magnitude 20 . As a result, differences between fluxes that are anyhow small will be dominated by the differences between fluxes that take larger values. To remedy this issue, we also calculated the relative changes, defined as the sum of normalized flux differences over all analyzed time points, where the flux differences between wild type and mutant were normalized to their respective absolute maximum value over all time points. To apply TC-iReMet2 we used a bottom-up assembled model of A. thaliana, ArabidopsisCore 23 . This model consists of 549 reactions, of which 229 are transport reactions and artificial reactions representing growth (biomass) and non-growth-associated maintenance functions (NGAM 22 ).
Sum of predicted flux differences. The overall flux distance of wild type compared to mutant across all predicted reactions differed before cold shift, with the wild type having a higher overall flux ( Fig. 2A). This prediction was consistent with the slight growth advantage of the wild type at the optimized growth temperature (Fig. 1). The difference of fluxes between consecutive time points remained approximately constant during the common hibernation phase, up to day 7. When the wild type resumed growth in the cold, the overall predicted flux differences increased approximately 3-fold. When considering the sum of flux changes per time step for wild type (Supplementary Fig. S2) and mutant ( Supplementary Fig. S3), we find similar changes for the wild type and   www.nature.com/scientificreports/ the mutant at the steps from 0 days to 1 day and 1 day to 7 days, with an increase in the change from day 7 to day 21 (Fig. 2B). However, we observe that the changes between day 7 and 21 are considerably larger in the wild type in comparison to the mutant, in line with the resumed growth of the former in the cold. In the following, we identify the reactions and pathways which contribute most to these observed differences.
Analysis of differential reactions. We next considered the flux differences for each reaction in the metabolic model. Additionally, we investigated reactions displaying large changes in flux differences at early time points, as the most interesting to understand the changes in the metabolic network functionality in response to the cold shift.
K-means clustering of reaction behaviour. We focussed on differential behaviour of all reactions between mutant and wild type, excluding transport reactions and artificial reactions to avoid bias due to lack of gene association for these reactions. To this end, we applied K-means clustering to group reactions (Supplementary Table S2) with similar relative flux changes, where the number of clusters was determined by the silhouette index (Supplementary Fig. S4). As a result, we identified K = 7 clusters of reactions (Fig. 3), with a maximum silhouette index value of 0.78, based on the relative flux changes (Fig. 3A). For comparison, we also consider the K-mean clustering of the nominal flux changes (Fig. 3B). To provide an intuitive description of clusters as well as reaction behaviour over time, we introduce a three-character pattern consisting of Up (U), Down (D) and No changes (N) if the respective relative flux differences increased, decreased or stayed the same between two time points. Using this If we consider the top 10 reactions (Supplementary Table S3) with respect to relative and nominal changes directly after cold shift, we find H-serine dehydrogenase (HSerDHNADP_h (UDU), HSerDHNAD_h(DUD)), isocitrate dehydrogenase (iCitDHNADP_m (DDD), iCitDHNAD_m(UUD)) as well as 6-phosphogluconic dehydrogenase (6PGDHNAD_h(DUD)), glutamate dehydrogenase (GluDH1NADP_m(DUD)) and glutamate synthetase (GluSNAD_h(UDD)) conserved among both measures. All these reactions are redox reactions. Additionally, 6-phosphogluconic dehydrogenase (6PGDHNADP_h(UDU)), glutamate dehydrogenase (GluDH2NAD_m(DUD)) and glutamate synthetase (GluSNAD_h(UDD)) can only be found in the top 10 reactions of nominal changes. Conversely, malate dehydrogenase (MalDH2NADP_c(UNN)), fructose-biphosphate aldolase (SBPA_h(UDD)) and sedoheptulose-biphosphatase (SBPase_h(UDD)) can only be found in the top 10 reactions of relative changes.

Pathways enriched in reactions with highly altered fluxes across time points. Metabolic reac-
tions do not function in isolation, so analysis and interpretation of the prediction is best carried out in terms of pathways. To identify the pathways that are changed over time, we used the metabolic pathways as defined by the underlying A. thaliana model 23 Table S4). We inspected and considered as relevant those pathways that were enriched with reactions displaying large predicted flux differences between wild type and mutant (Fig. 4). A reaction was defined to exhibit large changes, if its absolute sum of flux changes across all time points was above the median of considered reactions present in the model (excluding transport and artificial reactions, as specified above). To identify pathways enriched with such reactions we used the Fishers exact test with significance threshold P considering multiple hypotheses correction following the Benjamini-Hochberg procedure (Supplementary Table S5). Considering nominal changes, we found five pathways to be enriched for reactions with large changes. These pathways, ordered by decreasing P-value, with p < 0.01, include: the Calvin-Benson-Cycle (CBC), photorespiration, gluconeogenesis, leucine synthesis, and in addition with < 0.05 , glycolysis. Considering relative instead of nominal changes, we found pathways with p < 0.01 to include the Calvin-Benson cycle, glycolysis, gluconeogenesis, and in addition with p < 0.05 , photorespiration.
Flux sampling analysis with quadratic constraints. We examined how specific these findings are by sampling the solution space given in optimal solution for each considered time point. As a sufficiently large www.nature.com/scientificreports/ enough sample size gives information about range of fluxes as well as their probability, it gives the means to explore for alternative solutions and so for the uniqueness of the solution. In this case for each considered time point the proposed approach (see Methods -Flux sampling for TC-iReMet2) did not find a solution after 1000 trials. Therefore, this analysis indicates that the findings are specific, in the sense that alternative optima are unlikely, and significant as there are no other possible flux distributions in optimal solution.

Discussion
Here, we proposed a computational approach, termed TC-iReMet2, and showed that it provides the means for time-resolved predictions of fluxes while keeping the simplicity of the constraint-based modelling framework and allowing for the integration of relative metabolomic, transcriptomic, and morphometric data. The findings of this study indicate that TC-iReMet2, a differential flux profiling method, can be used to identify differential fluxes between wild type and mutants over time. It is important to note that TC-iReMet2 uses the ratio of transcripts as a proxy for the ratio of enzyme abundance (following GPR rules). This is a strong assumption, knowing that post-translational modifications and translational efficiency have a large effect on both the abundances and ratios of proteins. However, since transcript ratios are used as one component of the constraints, such an approach provides a better coverage of metabolic networks than proteomics data 14 . With the advances in the proteomics profiling, TC-iReMet2 has the potential to provide further applications closer to the assumptions of the approach. Moreover, the enzyme kinetic assumed in TC-iReMet2 does not consider saturation effects no presence of regulators (e.g. activators or inhibitors) of enzyme activity. Inclusion of a saturation effect, like in Michaelis-Menten kinetic, would not allow casting the problem with only linear constraints, rendering application to large-scale networks computationally challenging. Similar problem arises when considering the inclusion of regulation, which is additionally problematic due to the lack of information on how the effect of the regulator is captured in the enzyme kinetic form used, particularly for plants 32,33 . One possible approach to overcome this issues is to use a power-law formalism 34 , which would allow the constraints to remain linear, at the cost of making assumptions about which regulators affect a reaction rate and with what strength. For this reasons we have decided that TC-iReMet2 is formulated based on mass-action-like kinetic, while allowing for discrepancy to model possible effects due to the mentioned saturation and regulation.
Applying TC-iReMet2 to the comparison of growth deficient reil1-1 reil2-1 double mutant to A. thaliana Col-0 wild type before and after cold shift strongly supported for the previously suggested theory stating that REIL mediated ribosome biogenesis deficiency feeds back into metabolism. Overall, we find that flux differences are more similar during hibernation phase, with strong flux redistributions occurring at day 21. This is evident from the data from the morphometric analysis (Fig. 1, Supplementary Table S6). Mutant and wild type plants grow similar but start to differ strongly between days 7 and 21 after temperature shift.
Even more important, TC-iReMet2 enables a way to compare wild type and mutant differential fluxes prior to cold shift and in the early hibernation phase. Thus, it allows for the analysis of the mutant system relative to wild type without being obscured by the effects of differential growth occurring between days 7 and 21 of the current experiments. When considering differential fluxes during the hibernation phase, we find that REIL mediated ribosome biogenesis deficiency might feed back into metabolism by altering the RuBisCO carboxylase to oxygenase ratio (Supplementary Table S7). Additionally, mutant associated deficiencies of the CBC and glycolysis fluxes combined with mutant-specific increase of all fluxes in the photorespiratory pathway support this hypothesis (Fig. 4). Overall the strongest mutant flux deficiencies appear to be in the RuBisCO (carboxylation), FNR, malate dehydrogenase and alanine transaminase reactions (Supplementary Table S5).
Predicted relative flux changes directly after cold shift appear to be small. However, one day after cold shift, the fluxes of reactions distributed across various pathways of central metabolism, including carbohydrate, organic acid and amino acid metabolism differ between mutant and wild type. What is common to these reactions is that they all require NAD(P) as a cofactor. These changes may indicate either an altered redox state of these cofactors or more likely differential use of NAD and NADP after cold shift in the mutant. For example, when considering the predicted inverse flux changes of the mitochondrial iCit dehydrogenase isozyme reactions, iCitDHNAD_m and iCitDHNADP_m, we can deduce in agreement with our metabolic model that the mutant switches to preferential use of NAD rather than NADP for this reaction. Inversely, a preferential use of NADP is predicted for 6-phosphogluconic dehydrogenase reactions, 6PGDHNADP_h and 6PGDHNAD_h, and for the H-serine dehydrogenase isozyme-reactions, HserDHNADP_h and HserDHNAD_h. Taken together with the additional indicated flux changes of glutamate synthetase (GluSNAD_h), and of the mitochondrial malate dehydrogenase (MalDHNAD_m), or glutamate dehydrogenases, GluDH1NADP_m and GluDH2NAD_m, we hypothesize that the reil1-1 reil2-1 double mutant defect is associated with a NAD/ NADP cofactor deregulation.
Generation of this hypothesis would not have been possible by stand-alone analysis of the transcriptome data alone. When we compare the results of TC-iReMet2 with a differential analysis of the transcriptomics data reaction per reaction following GPR rules, only the increase of isocitrate dehydrogenase (iCitDHNAD_m) flux in the mutant can be found overlapping with TC-iReMet2's predictions (Supplementary Table S8). This indicates that TC-iReMet2's integration of metabolomics and transcriptomics data provides added value compared to the sole analysis of either, thus, allowing new and additional support for hypothesis generation. Verification of these findings and hypothesis testing can be performed by subsequent studies and detailed quantification of NAD and NADP levels and their redox states under same and extended experimental set-ups. Altogether, the predictions from TC-iReMet2 suggest that altered use of NAD and NADP or of their redox state is an important mechanism by which REIL mediated ribosome biogenesis deficiency feedbacks into metabolism early after cold shift.
The current formulation of TC-iReMet2 has the potential to be further optimized, since the weighting parameters of the objective function are chosen based on the assumption that the consecutive increase between time points equals consecutive decrease of weighting. Validation of predictions gave robustness to this assumption. www.nature.com/scientificreports/ Yet, the usage of different weights could be considered based on other insights from independent physiological measurements. In addition, rather than using relative transcriptomics data, relative proteomic data 9 or enzyme activity measurements 35 could be integrated to provide more reliable predictions that are less influenced or obscured by post-transcriptional levels of regulation than transcriptome data. Therefore, TC-iReMet2 improves existing constraint-based approaches for differential flux prediction by accounting for possible temporal physiological change while also allowing for the integration of morphometric data.

Methods
Flux sampling for TC-iReMet2. Uniform flux sampling provides an unbiased characterization of the solution space. When enough flux distributions are sampled, they can be used to analyze their probability distributions or the range of specific fluxes. In this setup, flux sampling was performed using a random walk algorithm (Hit and Run). A linear program with a single quadratic constraint is defined to find possible alternative solutions. For this, the solution space has to be defined. In addition to the defining constraints given by TC-iReMet2, a single quadratic constraint, due to quadratic objective function TC-iReMet2 is based on, has to be introduced. It fixes the value of the objective function to be the same as in optimal solution, forcing the optimization to find alternative optima. Here, z denotes the value of TC-iReMet2's objective function and z * the value found in optimal solution at the specific time point. Fluxes at the current time point are denoted by v t+1 , whereas v t denotes the flux distribution of the prior time point. We also allow for a small deviation, denoted by ζ , of the objective function at the optimum to counteract numerical problems. This way the solution space, containing all possible solutions at the time point specific optimum, is defined. To sample this space, steps are done as follows: 1. Select an initial point v 0 in solution space (here, we used v t+1 in optimal solution, derived from the main optimization problem of TC-iReMet2 as v 0 , since it must lie in the solution space).
2. Select a random direction v direction pointing in solution space.
3. Find the extreme point in solution space described by v t+1 = v 0 + v direction by solving a linear program with a quadratic constraint (due to the quadratic problem the main objective is based on): If there is no solution to the optimization problem given above, v direction does not point into solution space. As a consequence, steps 2. and 3. are repeated until a solution is found. 4. If there is a solution for the optimization program at step 3, a new point at the edge of the solution space can be determined to form a line segment with the initial point.  Numerical stability of TC-iReMet2. The multiplication of relative metabolite levels, substitution of unmeasured metabolite ratios with their respective minimum and maximum ratio value together with approximated enzyme ratios can lead to immense ratio constraints, which in turn could lead to numerical instabilities. Determining a maximum considered ratio constraint is therefore crucial to ensure numerical stability. To this end, we calculated the flux distributions allowing for a maximum ratio constraint ranging from 10 1 to 10 21 (includes the maximum possible ratio constraint in this setup) 10 times. Since each of those repeated measurements resulted in the same flux distribution, we used Pearson correlation to measure similarity between each flux distribution to its prior and successive one. Overall, correlations between flux redistributions were very high being above 0.9. The highest correlated region, while having a feasible solution at each considered time point, was detected when allowing for a maximum ratio constraint of 10 8 . This leads to 252, 252, 253 and 189 ratio constraints made at each considered time point. Hence, no ratio constraint exceeding 10 8 was considered in this setup to ensure numerical stability.
K-means clustering. We used R statistical programming languages implementation of the K-means algorithm with seven assumed clusters ( K = 7 ) and Euclidean distance as distance measure. The reason for selecting K = 7 is the following: We assumed there to be one cluster of no differences fluxes, one cluster displaying stronger wild type flux with a rise at 21 days and inversely the same for mutant. Two clusters of inverse behaviour where wild type or mutant flux is stronger at the first time points with a shift in sign at 21 days. Lastly, we assumed two clusters of consistent flux difference favouring wild type or mutant conserved over all time points. This line of reasoning was supported by the silhouette index analysis, which specify the number of clusters K = 7 to maximize the value of the index.

Fishers exact test for enrichment analysis.
We used a right-tailed Fisher's exact test to determine the enrichment in regulated reactions of a pathway. To this end, we defined a reaction as regulated if its sum of relative or nominal differences across time points was above its corresponding median of all considered reactions, else we considered the reaction to be unregulated. Therefore, we tested the association between regulated reactions and pathways for both relative and nominal differences. This test was conducted with a significance level of 0.05 through Matlab's 'fishertest' function. The resulting P-values were corrected for multiple hypotheses testing following the Benjamini-Hochberg procedure.
Parameterizing the objective function of TC-iReMet2 and estimating fractions of biomass yield. The analyzed time points differed in scale and therefore constituted a good case to test the assumption that a subsequent flux distribution is dependent on the previous one, therefore allowing for different magnitudes of physiological change. To model the dependency of flux distributions between time points, we assumed a steady decrease in the dependence as the interval between the points increases. More specifically, the following weights were used at each time point depicted by Table 1. Weighting of flux distribution dependency of scenario A and B to a prior one is denoted by β and γ respectively. Difference of flux distributions between scenario A and B at analyzed time point was weighted by α.
To model the fraction of biomass yield in scenario A to the biomass yield in scenario B at each specific time point κ t+1 , we used four biomass proxy parameters (Supplementary Table S1), two diameter measurements of the A. thaliana rosette, diameter 1 and 2, the apparent planar leaf area, and leaf perimeter (i.e. the circumference). The morphometric parameters apparent planar leaf area and perimeter underestimate biomass accumulation, since rosette leaves could slightly overlap. We contrasted these estimates by using the sum of diameter 1 and 2 as proxies of biomass accumulation. The diameter may slightly overestimate biomass because only the longest leaves are considered. Accordingly, we integrated the four biomass proxy parameters by giving equal weight to each one of them (Supplementary Table S1). In detail, the morphometric parameter measurements were averages separately for the wild type and the mutant. Ratios of mutant and wild type were calculated per time point based on the averages of biomass proxies. Finally, the resulting ratios across the four biomass proxies were averaged to obtain time point specific biomass fractions κ t+1 . We allowed to a deviation δ of +/-0.05 from the biomass fraction κ t+1 as these calculations are estimates. In absence of biomass estimates, e.g. at day 1, we assumed a steady www.nature.com/scientificreports/ decrease of biomass fraction κ t+1 between day 0 and day 7. Therefore, the ratio at day 1 is a seventh closer to the biomass fraction κ t+1 of day 7 compared to day 0, resulting in a fraction of 0.77. All used biomass fractions κ t+1 are depicted in Table 1.
Transcriptomics and metabolomics data used. The transcriptomics data are already published and are uploaded to the Gene Expression Omnibus (https ://www.ncbi.nlm.nih.gov/geo/) and are available through accession number GSE101111. The metabolomics data are obtained from the Supplemental Table S1 of Beine-Golovchuk et al., 2018 30 . The morphometric data are obtained from Schmidt et al., 2013 2 . All data are included in the provided GitHub repository (https ://githu b.com/tcire met2/TC-iReMe t2) as well as in the Supplementary Tables to ensure easy access and reuse of the provided implementation.

Implementation and tools.
For implementation of TC-iReMet2 we used "MATLAB 2017b, The Math-Works" 36 in conjunction with the Tomlab optimization environment 37 . Statistical analysis and creation of figures was done with R 38 programming language and R's ggplot2 library 39 and "MATLAB 2017b, The MathWorks" 36 . The implementation is available at https ://githu b.com/tcire met2/TC-iReMe t2.

Code availibility
The transcriptome data are available from the Gene Expression Omnibus (https ://www.ncbi.nlm.nih.gov/geo/) through accession number GSE101111. The metabolome data are previously published among the supplemental data of Beine-Golovchuk and co-authors 2,3 .