Accepted Article

Received Date : 19-Sep-2014 Accepted Date : 19-Dec-2014 Article type

: Technical Advance

Integration of transcriptomics and metabolomics data Chlamydomonas’ metabolic response to rapamycin treatment

specifies

Sabrina Kleessen3#, Susann Irgang2,#, Sebastian Klie3, Patrick Giavalisco4,*, and Zoran Nikoloski1,* 1

Systems Biology and Mathematical Modeling Group, Max-Planck Institute of Molecular Plant Physiology, Am Mühlenberg 1, 14476 Potsdam, Germany 2 Genes and Small Molecules Group, Max-Planck Institute of Molecular Plant Physiology, Am Mühlenberg 1, 14476 Potsdam, Germany 3 Targenomix GmbH, Am Mühlenberg 11, 14476 Potsdam, Germany 4 Experimental Systems Biology Group, Max-Planck Institute of Molecular Plant Physiology, Am Mühlenberg 1, 14476 Potsdam, Germany

#

These authors contributed equally to this work

*Corresponding authors: Zoran Nikoloski Max-Planck Institute of Molecular Plant Physiology Am Mühlenberg 1 14476 Potsdam phone: +49 331 567 8630 fax: +49 331 567 8136 [email protected]

Patrick Giavalisco Max-Planck Institute of Molecular Plant Physiology Am Mühlenberg 1 14476 Potsdam phone: +49 331 567 8632 fax: +49 331 567 8236 [email protected]

Running title: Model-driven integration of metabolomics data This article has been accepted for publication and undergone full peer review but has not been through the copyediting, typesetting, pagination and proofreading process which may lead to differences between this version and the Version of Record. Please cite this article as an 'Accepted Article', doi: 10.1111/tpj.12763 This article is protected by copyright. All rights reserved.

modeling/

data

integration/

Accepted Article

Keywords: constraint-based metabolomics/ rapamycin

Chlamydomonas

reinhardtii/

Summary Flux phenotypes predicted by constraint-based methods can be refined by inclusion of heterogeneous data. While recent advances facilitate the integration of transcriptomics and proteomics data, purely stoichiometry-based approaches for prediction of flux phenotypes by considering metabolomics data are lacking. Here we propose a constraint-based method, termed TREM-Flux, for integrating time-resolved metabolomics and transcriptomics data. We demonstrate the applicability of TREM-Flux to dissect the metabolic response of Chlamydomonas reinhardtii to rapamycin treatment by integrating the expression levels of 982 genes and the content of 45 metabolites obtained from two growth conditions. The findings pinpoint cysteine and methionine metabolism to be most affected by the rapamycin treatment. Our study shows that the integration of time-resolved unlabeled metabolomics data in addition to transcriptomics data can specify metabolic pathways involved in the system’s response to a studied treatment. Introduction With the advent of high-throughput technologies there has been a surge of data from systems biology studies. These data are increasingly used to predicting the condition-specific behavior of (bio)chemical systems and components (Kitano, 2002; Westerhoff et al., 2009). However, they also pose the challenge of reconciling the observed behavior of components with complex physiological traits to support and further refine the underlying mechanisms.

Aside from the classical statistical approaches, modern methods for in silico analysis usually combine data with mechanistic models to arrive at (steady-state) reaction fluxes as integrative phenotypes (Masakapalli et al., 2013; Varma and Palsson, 1994a). For instance, the prominent approaches based on flux balance analysis (FBA) (Varma and Palsson, 1994a) have allowed the prediction of steady-state fluxes with the assumption that metabolic networks optimize an objective (e.g., biomass production) (Gianchandani et al., 2010; Orth et al., 2010; Raman and Chandra, 2009; Varma and Palsson, 1994b).

However, the steady-state assumption along with the employed optimization principles often does not suffice to obtain a unique flux phenotype. Integration of high-throughput data aims at reducing the solution space of feasible flux phenotypes (Reed, 2012; Blazier and Papin, 2012).

This article is protected by copyright. All rights reserved.

Accepted Article

For instance, under the assumption that gene expression is proportional to fluxes of corresponding reactions, transcriptomics data can provide additional constraints to improve the accuracy of flux predictions (reviewed in (Blazier and Papin, 2012)). However, these approaches rely on the steady-state assumption, inherent to FBA-driven methods, which precludes the analysis of the dynamics of flux (re)distributions and metabolite levels.

Additional constraints obtained from metabolite levels result in predictions of time-resolved flux phenotypes closer to the actual physiological state. For instance, dynamic FBA (DFBA) offers an alternative approach to predict time-resolved metabolite and flux profiles (Mahadevan et al., 2002; Luo et al., 2006; Luo et al., 2009; Kleessen and Nikoloski, 2012; Varma and Palsson, 1994b). However, the existing applications of DFBA still require some knowledge of enzyme kinetics together with their parameters, and have thus been mostly applied to simple metabolic networks (e.g., cycles and linear pathways) of small size. In addition, applications of DFBAbased methods to genome-scale models are currently hampered by the size of the resulting instances as well as by the lack of optimization platforms which scale well with the size of the models (Kleessen and Nikoloski, 2012). Development of methods for integration of additional constraints obtained from metabolite levels in the constraint-based framework holds the promise to further reduce the solution space for flux distributions, leading to predictions of time-resolved flux phenotypes closer to the actual physiological state. Concomitantly, such data integration will allow for more precise predictions, reflected in the smaller flux ranges compatible with the constraints. In comparison to genes and their expression levels, the quantified levels of metabolites, acting as substrates and products of biochemical reactions, are more closely related to reaction fluxes. Therefore, understanding the extent to which metabolomics data may contribute to determining the flux state of an organism is paramount for future developments of this field. Despite these opportunities of using metabolomics data, there exist only a few methods allowing integration of metabolite levels into large-scale constraint-based modeling approaches for the purpose of network reconstruct and prediction of phenotypes. However, the few existing approaches allowing integration of metabolite levels into large-scale constraint-based modeling approaches are either parameterdependent (e.g., Integrative Omics-Metabolic Analysis (IOMA) (Yizhak et al., 2010) and into mass-action stoichiometric simulation (MASS) models (Jamshidi and Palsson, 2010)), consider only the presence/absence of the metabolite to refine models (e.g., Gene Inactivation Moderated by Metabolism, Metabolomics, and Expression GIM3E (Schmidt et al., 2013)), or involve nonlinear optimization due to thermodynamics constraints (e.g., thermodynamics-based metabolic flux analysis (TMFA) (Henry et al., 2007)), all in a given steady state.

This article is protected by copyright. All rights reserved.

Accepted Article

This succinct overview of the existing approaches indicates that currently there exists no method which allows the integration of quantitative unlabeled time-resolved metabolomics data with a large-scale metabolic network model to infer fluxes without inclusion of kinetic parameters. Such an approach should maintain the simplicity of constraint-based methods and their largely parameter-independent nature, while allowing for inclusion of additional data closer to the actual physiological state. Here we propose a method termed Time-Resolved Expression and Metabolite-based prediction of Flux values (abbreviated by TREM-Flux), facilitating the integration of time-resolved transcriptomics and metabolomics data into constraint-based modeling. The approach is intended to also demonstrate that the consideration of such metabolomics data may more precisely characterize the metabolic state of a biological system in comparison to the usage of transcriptomics (or proteomics) data. We used TREM-Flux to get insights in the systemic response of mixo- and photo-autotrophically grown Chlamydomonas reinhardtii (Chlamydomonas) and to identify specific metabolic functions (pathways) affected by to the anti-proliferative drug rapamycin. Rapamycin is a heterocyclic macrolide (Vézina et al., 1975) which, in eukaryotic cells, binds to the 12-kD FK506-binding protein (FKBP12) and the composed complex then inhibits the activity of an essential regulator of growth, namely, the target of rapamycin (TOR) kinase (Choi et al., 1996; Chen et al., 1995). TOR can regulate growth and development by activation or repression of global anabolic or catabolic processes (Laplante and Sabatini, 2012; Wullschleger et al., 2006). Contrary to plants (Dobrenel et al., 2011) but similar to human and yeast systems (Wullschleger et al., 2006), cell growth of Chlamydomonas is highly sensitive to rapamycin treatment (Crespo et al., 2005; Díaz-Troya et al., 2008). Therefore, it offers an ideal experimental system to explore the function of TOR in a photoautotrophic organism by testing how this drug affects Chlamydomonas metabolism. Constraint-based integration of highthroughput data may pinpoint pathways which are most affected by modulation of TOR, and can serve as candidates for further investigations. The following is a brief summary of our contributions which delineates the TREM-Flux approach from the contending alternatives: (1) we demonstrated that the inclusion of only transcript levels to derive transcript bounds under the (strong) steady-state assumption provides comparable precision of flux phenotypes as the inclusion of metabolite levels and transcripts without the steady-state principle; (2) we established that the steady-state principle, when combined with the transcript data, resulted in prediction of many inactive pathways that are not in correspondence with the changes observed on the level of metabolites; and (3) we showed that the predicted flux phenotypes are in agreement with the differential behavior of metabolites, thus providing the connection between the different systemic levels under two independent growth conditions for the investigated system. These findings hold equally well on data collected from two different growth conditions, photoautotrophic and mixotrophic. Altogether, alongside the generation and analysis of a rich multi-level data set, our study illustrates that the proposed

This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.

Accepted Article

over all measured metabolites in each compartment (frequently used in other studies, e.g., (Flamholz et al., 2013; Nöh and Wiechert, 2006)).

Formulation of TREM-Flux Since the predicted flux distribution for a time point metabolite levels measured at

is dependent on the difference between

and the consecutive time point

,

,

, the flux

distribution for the last measured time point cannot be predicted. Accordingly, for time points the flux distributions for time points ,

measured

, can be obtained as solutions

to the following linear program (Figure 1E):

where

I

II

(1)

III IV V

denote the vector of coefficients in the biomass reaction, used as a physiological

objective in the model,

denotes the total (non-compartmentalized) level of metabolite

min and max operators are taken over all measured metabolites,

and

derived bounds following E-Flux, as explained in the preceding section, and

, the

are transcriptthe number of

analyzed time points. The proposed approach (Figure 1E and Eq. (1)) assumes that fluxes scale linearly with the change in metabolite levels in the time period between two time points, denoted by . Moreover, is allowed to vary within the time-series experiment (i.e., the sampling does not have to be uniform). As a result, if the change between three or more time points is linear, then the flux values predicted on a finer time scale will not differ from those on a coarser time scale (over the same set of time points and with the same data used).

This article is protected by copyright. All rights reserved.

Accepted Article

Response of Chlamydomonas to rapamycin treatment Transcriptomics and metabolomics data were obtained from the cell-wall deficient Chlamydomonas strain CC503 grown under continuous light and ambient CO2 at 28°C in photoautotrophic and mixotrophic conditions with and without rapamycin treatment, referred to as treatment and control conditions, respectively (see Experimental Procedures). Absolute levels (contents) of 45 metabolites were available for eight time points after rapamycin addition, i.e., 0, 30, 60, 90, 120, 240, 360, 480 min (Data S1). Moreover, expression levels of 11455 genes were obtained for six time points (i.e., 0, 30, 60, 120, 240, 480 min) (Data S2). The time points corresponding to 90 and 360 min for the transcriptomics data were interpolated by using cubic splines to facilitate the integration with all available metabolomics data (see Experimental Procedures).

Application of TREM-Flux The proposed method was applied on a genome-scale metabolic model of Chlamydomonas for the same strain used in this study (Chang et al., 2011). The biomass reactions of the reconstruction for mixotrophic and photoautotrophic conditions were related to the cell counts (Data S3), measured over all considered time points and conditions (cf. Experimental Procedures). We applied the following two constraints: (i) the ratio of accumulated biomass between the two conditions equals the ratio of cell counts at each time point, and (ii) the biomass predicted at a given time point is not smaller than the increase in cell counts with respect to the previous time points. These constraints were included to the generic formulation of TREM-Flux, in Eq. (1) above, to provide more physiologically meaningful flux predictions. From the 2143 reactions in the reconstruction, 960 were associated with particular genes encoding the enzymes considered in the model. By using the transcriptomics data for the corresponding genes, the bounds of 398 out of these 960 reactions were affected in comparison to the default bounds included in the original model. The result of solving the optimization problem in Eq. (1) with time points is time-resolved flux distribution specifying the flux of every reaction in the first time points associated with optimal biomass accumulation. However, it is well known that the unique optimal value of the biomass objective in problem in Eq. (1) is not necessarily associated with a unique flux distribution. This issue is carefully investigated in the following section.

Flux variability and TREM-Flux Analyzing the alternative optima is crucial for interpreting the results from constraint-based modeling approaches (Lewis et al., 2012). To this end, flux variability was determined for each reaction, time point and condition following the application of TREM-Flux based on flux variability analysis (FVA) (Mahadevan and Schilling, 2003; Burgard et al., 2001) (Figure 2A). More specifically, we determined the lower and upper boundary for a reaction flux at a given This article is protected by copyright. All rights reserved.

Accepted Article

time point which guarantee the same optimal value for the biomass objective. Altogether, 2127 reactions (excluding the biomass reaction and the reactions whose lower and upper bounds coincided) were examined in this time-resolved FVA (tFVA, Figure 2B).

The magnitude of variation of reaction fluxes in alternative optima was used to quantify the capacity of the integration of the used data to narrow down the space of optimal flux distributions. To this end, the range of flux variation (i.e., maximum-minimum) was determined along with its contribution to the range from the bounds considered in TREM-Flux (i.e., ). This resulted in a percentage of variability for each reaction in the considered conditions and time points (Figure 2C). The percentages of variability for all reactions were divided into ten bins, each covering a range of ten percent, as shown in Figure 2D.

Investigation of the first five time points showed only minor differences with respect to the variability between the growth conditions. However, for time points 6 and 7 (240-360 min after rapamycin treatment) there was an increased variability in flux predictions for rapamycin-treated cells. An explanation for this observation is that the rapamycin-treated cells do severely restrict their growth without fully switching off the carbon fixation. Therefore, the biological processes in the later time points were expected to be more severely affected by the treatment. This behavior resembled a starvation-like phenotype under full nutrition. Previously, it has been shown that in Arabidopsis thaliana this behavior results in accumulation of storage carbon (starch and lipids) as well as amino acids (Caldana et al., 2013).

Nevertheless, between 1045 and 1553 of the 2127 reactions were found in the bin denoting 0 - 10 percent variability within the modeled bounds over all time points. Altogether, remarkable restriction in variability was already obtained by the 45 metabolites from the 1063 included in the metabolic reconstruction (see below for comparative analyses).

tFVA-based comparison of TREM-Flux and E-Flux Here we aimed at demonstrating the effect of integrating metabolomics data in addition to gene expression levels in further constraining the predicted flux distributions. To this end, we compared TREM-Flux to a modification of the E-Flux approach to facilitate prediction of timeresolved flux distributions. In the time-resolved E-Flux approach, a steady state was assumed for each time point, and the sum of the objective over time was optimized. For comparison of the results from TREM-Flux, flux values were predicted for all but the last time point. Accordingly,

This article is protected by copyright. All rights reserved.

Accepted Article

the flux distributions for time points

,

were obtained by changing lines III and

IV or Eq. 1 with the steady state assumption

.

(2)

The percentage of variability was next analyzed in the same way as for TREM-Flux, as shown in Figure S1. The number of reactions in the first bin was slightly larger in E-Flux (1143-1636 reactions) compared to TREM-Flux (1045-1553 reactions). Moreover, the number of reactions in the last bin for E-Flux (190-311 reactions) was similar to the observed behavior for TREM-Flux (193-346 reactions). The last two time points demonstrated a behavior different from that for the other time points, as also observed in TREM-Flux. While the remaining bins were typically smaller in E-Flux, we nevertheless have to take into account that the first bin of the timeresolved E-Flux results includes the 526 blocked reactions arising due to the steady-state assumption (Table S1). In contrast, TREM-Flux does not result in any blocked reaction.

The blocked reactions in the modified version of E-Flux were enriched for galactose, histidine and beta-alanine metabolism as well as transport in flagella, Golgi and nucleus (see Table S2 for the full list). Moreover, cysteine and methionine metabolism, which represents an essential pathway in the assimilation of sulfur and a precursor of the synthesis of the redox-regulator glutathione, was also observed to be enriched for blocked reactions. In total six from the nine reactions participating in the cysteine and methionine metabolism were blocked, and, thus, the time-resolved E-Flux, assuming a steady state, cannot be employed to predict the effect of the rapamycin treatment on this metabolic function. This observation was especially interesting since it has been shown previously that the redox-balance of the rapamycin-treated cell might be affected, and the synthesis of the downstream product of cysteine, namely, glutathione, was upregulated in Arabidopsis (Ren et al., 2012; Caldana et al., 2013). Altogether, these findings suggest that the conclusions drawn only by integration of transcriptomics data together with the steady-state assumption may be misleading.

Effects of integrating metabolite and transcript levels quantified by tFVA To further quantify the effects of integrating the time-resolved metabolite and transcript profiles, we first determined the distribution of flux ranges in the alternative optima, obtained from tFVA, over all reactions and time points. The characteristics of the distributions from two different conditions/approaches can be used to compare the effects of the data integration on the specificity of the flux phenotype predictions. To this end, we used the mean, median, and the

This article is protected by copyright. All rights reserved.

Accepted Article

skewness for each individual distribution of flux ranges (Figure S2). In Figure 3, we only present the findings of the control under the mixotrophic growth condition in boxplots, as similar observation can be made for photoautotrophic growth (Table S3).

As summarized in Figure 3, we observed that the inclusion of more metabolites results in a smaller median and mean as well as larger skewness (see Table S3), supporting the general decrease of flux ranges obtained from TREM-Flux. We next sampled 40, 30, 20, and 10 metabolites from the 45 measured, and the values for the characteristics of the distributions shown in Figure 3A and Table S3 online are the averages over 10 samples. While the mean and median of the time-resolved E-Flux were marginally smaller in comparison to TREM-Flux, we note that this is largely due to the steady-state assumption and/or bias due to consideration of the blocked reactions. Therefore, we also presented the characteristics of the distributions from which blocked reactions were removed. Our claim about the effects of the steady-state assumption was further supported by the characteristics of the distributions obtained upon a simulated reduction of the range in which metabolite levels may change between two consecutive time points (i.e., from 10 to 0.1 mmol gDW-1). Similar behavior could be observed by using only the flux bounds included in the metabolic reconstruction (Figure 3B and Table S3).

Comparison of flux distributions In the previous sections we demonstrated that majority of reactions show flux values which are well constrained due to the integration of both metabolomics and transcriptomics data in the large-scale metabolic model. Therefore, we can use one of the flux distributions as a representative to conduct additional differential analysis as one would do in classical data analysis. There are two principles ways in which the time-resolved flux distributions can be used: (i) for each reaction one can determine the difference in flux between consecutive time points or between two conditions in a given time point, and can then rank the reactions based on these differences (Figure S3A); (ii) for each reaction, one can consider the flux profile of each reaction over time, and the comparison for the reaction’s behavior between two conditions can be obtained by applying similarity measures to the respective flux profiles (Figure S3B). In the following, we focus on these two types of comparisons for the representative flux distributions to dissect the response of Chlamydomonas to the rapamycin treatment. Comparison of flux distributions between investigated growth conditions In the following, for each time point, we compared representative flux distributions, predicted by using TREM-Flux, under control and treatment conditions for both mixotrophic and photoautotrophic growth (illustrated in Figure 4). Following the approach (i) above, given two consecutive time points, we determined the absolute difference of the corresponding flux This article is protected by copyright. All rights reserved.

Accepted Article

distributions, and ranked the reactions according to the resulting values (Figure S3A). We then examined the reactions present in the top 10% (i.e., 214 reactions) to determine enriched metabolic functions in these selected reactions. Altogether, the analysis of the flux distributions indicated that 30 out of 82 metabolic functions were present in the top 10% of reactions showing the highest absolute difference in fluxes between control and treatment at every time point. At significance level of 0.05, we observed that six of these 30 metabolic functions were associated to transport (Figure 4 and Table S4). This observation was in line with the fact that 17.5 % of reactions in the used model interconnect the 10 compartments. Of these, transport to chloroplast, mitochondria, and extracellular transport were determined to exhibit the most remarkable difference in more than three time points. This was not unexpected since the reactions and metabolic functions affected by the treatment were reactions which rely on transport of metabolic intermediates between these compartments. Our findings, indeed, indicated that reactions from the TCA cycle as well as the oxidative phosphorylation (mitochondrion), pyruvate metabolism, the carbon and CO2 fixation (chloroplast), glycolysis/gluconeogenesis (cytosol) were ranked in the top 10% with respect to changes in flux over most time points in both growth conditions. Moreover, three of these metabolic functions were also found as active over all conditions and time points; hence, the aforementioned metabolic functions were among the most dynamically behaving between the control and rapamycin treatment, and would be expected to result in the strongest effects on the changes of the metabolic state. According to previous studies in Chlamydomonas, but also Arabidopsis, yeast and humans, these observations are not unprecedented: TOR is known to function as a central sensor of nutrient availability and energy distribution (Wullschleger et al., 2006). As can be seen from the abovementioned metabolic functions, inhibition of TOR regulates the conversion of light energy into metabolic intermediates necessary for growth and development. Especially the TCA cycle and the oxidative phosphorylation have been shown to be negatively influenced by the repression of TOR in Arabidopsis thaliana (Ren et al., 2012; Moreau et al., 2012; Caldana et al., 2013). Next to the above mentioned metabolic functions the alanine and aspartate metabolism, the pentose phosphate pathway (3 time points), retinol metabolism (2 time points) and valine, leucine and isoleucine degradation were largely affected under mixotrophic growth, followed by a few (single) large changes in glutamate metabolism, valine, leucine and isoleucine biosynthesis and phenylalanine, tyrosine and tryptophan metabolism. In contrast, methionine metabolism, as well as glycine, serine and threonine metabolism, N-metabolism, and arginine and proline metabolism were specific to photoautotrophic growth. These findings were again to some extent in full agreement to previous results from research on the effects of TOR repression in eukaryotes. Especially the influence on glutamate and general amino acid and N-metabolism have been welldocumented (Laplante and Sabatini, 2012; Wullschleger et al., 2006). Repression of TOR has

This article is protected by copyright. All rights reserved.

Accepted Article

been initially defined to resemble N-starvation in yeast and to affect the amino acid homeostasis in the cell (Komeili et al., 2000). These early findings allowed connecting N-metabolism, especially the synthesis of glutamate and glutamine, to the TCA cycle and to the signaling between the different compartments (Loewith and Hall, 2011). Taken together these results suggested that similar metabolic functions of large expected absolute flux values were affected between the treatment and control in the two growth conditions. Especially the connection between amino acid metabolism, TCA cycle, and carbon fixation represent common pathways affected by the inhibition of TOR.

Time-resolved behavior of metabolic functions and concordance with differential behavior of metabolites Here we followed the approach (ii), above, to investigate the correlation between the timeresolved flux distributions under the control and treatment for every reaction. This approach would allow us to identify the reactions which demonstrate large negative correlations. We hypothesize that such reactions, showing opposite trends of the flux profiles between the considered conditions, are responsible for altering the pool sizes of the participating metabolites. However, as in the previous section, reactions can be grouped into a metabolic function; here, we used the median of the correlation values for the participating reactions to characterize the behavior of the entire metabolic function (Figure S3B) and determined those which are in line with the abovementioned hypothesis (see Figure 5 for illustration).

We found that the most pronounced negative values for the medians under mixotrophic growth occurred in: protein synthesis, ascorbate and aldarate metabolism, tyrosine metabolism, cofactor recycling, phenylalanine, tyrosine and tryptophan biosynthesis, sphingolipid metabolism and sulfur metabolism. Under photoautotrophic conditions, the metabolic functions with such behavior included: protein synthesis, sulfur metabolism, glutathione metabolism, sphingolipid metabolism, ascorbate and aldarate metabolism, nicotinate and nicotinamide metabolism and Oglycan biosynthesis (see Figure 5 and Table S5 for the full list).

The metabolic functions characterized with opposite trends between control and treatment for both growth conditions included (Figure 5): protein synthesis, sulfur metabolism, shingolipid metabolism, ascorbate and aldarate metabolism, nitrogen metabolism, oxidative phosphorylation as well as vitamin B6 metabolism.

This article is protected by copyright. All rights reserved.

Accepted Article

If the median of the correlation values is negative, we hypothesize that the respective metabolic function harbors the differentially behaving metabolites (i.e., the metabolites whose pool sizes change between the control and treatment). Indeed, analysis of differential behavior (see Experimental Procedures) indicated that 32 of the 45 metabolites show differential behavior in at least one of the growth conditions between treatment and control (see Table S6). The differentially behaving metabolites were enriched for cysteine and methionine metabolism under both mixotrophic and photoautotrophic growth, while under mixotrophic growth this was also the case of mitochondrial transport. Cysteine and methionine metabolism an essential pathway in the assimilation of sulfur, and the synthesis of the downstream product of cysteine, namely, glutathione, was shown to be affected in Arabidopsis upon rapamycin treatment (Ren et al., 2012; Caldana et al., 2013). Reassuringly, glutathione metabolism exhibited negative median correlation only under photoautotrophic growth (see Table S5).

Time-resolved behavior of metabolic functions and concordance with differential behavior of genes Moreover, at significance level of 0.01, we identified 5833 and 3286 out of 11455 genes as differentially expressed (DE) between rapamycin treatment and control condition under mixotrophic and photoautotrophic growth, respectively (see Experimental Procedures). This indicated the broad response to rapamycin, which renders it difficult to delineate the main effects of the treatment which could be investigated through further experiments.

For both growth conditions, we determined the metabolic functions for which the DE genes were enriched and the number of DE genes and the type of change (i.e., up/down-regulation) in the enriched metabolic function (Data S4 and S5). The enrichment for metabolic functions was carried out only for those DE genes which could be mapped onto the metabolic network model of Chlamydomonas (i.e., altogether 982 out of 11455 genes).

Under mixotrophic growth condition, we found that the DE genes were enriched for three metabolic functions: cysteine and methionine metabolism (including eight DE genes from the total of nine genes, denoted by 8/9), porphyrin and chlorophyll metabolism (32/53), and urea cycle and metabolism of amino groups (6/10). Five DE genes in cysteine and methionine metabolism were down-regulated, while 24 and 4 DE genes were up-regulated in the two metabolic functions: porphyrin and chlorophyll metabolism as well as urea cycle and metabolism of amino groups, respectively. Under photoautotrophic growth conditions, we found the DE genes to be enriched for six metabolic functions, including: biosynthesis of steroids (including 15 DE genes from total of 51, i.e., 15/51), glutamate metabolism (13/23), glycerophospholipid

This article is protected by copyright. All rights reserved.

Accepted Article

metabolism (15/73 genes), methane metabolism (2/2), nitrogen (N) metabolism (12/22) as well as phenylalanine, tyrosine, and tryptophan biosynthesis (13/27). Trends of up-regulation were typical for biosynthesis of steroids (11/15) and phenylalanine, tyrosine, and tryptophan biosynthesis (9/13); while down-regulation was predominant in glutamate metabolism (10/13), glycerophospholipid metabolism (10/15), methane metabolism (2/2), and N-metabolism (7/12). Analogous findings were obtained when considering the MapMan ontology (Thimm et al., 2004) (Data S6). Note that the mixture of up- and down-regulation of DE genes in the same metabolic function provides further challenges in the interpretation.

These findings further indicate that the analysis of the predicted flux distributions was in agreement with the results from enrichment analysis of DE genes with respect to metabolic functions and MapMan bins (see Tables S4, S5).

Conclusions Understanding the metabolic response of a biological system exposed to different environmental stimuli is of great importance to obtain insights in the contribution of individual system components and pathways to the physiology of the system. While gene expression data only provide putative insights in the physiological response of a system, metabolomics data and the corresponding flux profiles reflect the ultimate physiological response. Here we propose an optimization-based approach, termed TREM-Flux, which can be used to integrate time-resolved transcriptomics and quantitative metabolomics data and to predict the corresponding flux. Application of TREM-Flux in a genome-scale model reconstruction pinpoints the metabolic functions involved in the time-resolved metabolic response, as illustrated in the case of rapamycin treatment in Chlamydomonas reinhardtii.

Our findings point out the added value of considering metabolite levels in genome-scale models by obtaining more accurate predictions of the active physiological state. The predicted flux distributions can be used to determine differences in particular metabolic functions between control and treatment under different growth conditions. The integration of additional data, such as enzyme-activities and protein levels, could further reduce the solution space in future analysis. To this end, experiments will have to be carefully planned to account for differences in time scales on which cellular process take place with sampling at high enough resolution to allow better insights into the changes of the levels of molecular components. Although the present study relies on data and a genome-scale model from green algae, we believe that any other organisms could be analogously interrogated if data sets and network reconstructions of similar quality are available.

This article is protected by copyright. All rights reserved.

Accepted Article

Experimental Procedures Model and data We used the genome-scale metabolic reconstruction of Chlamydomonas reinhardtii (strain CC503) from Chang et al. (Chang et al., 2011). After exclusion of external metabolites 2143 reactions and 1636 metabolites remain, from which 1063 metabolites were unique. Functional annotation in the metabolic network was performed based on the gene models from the Augustus update 5 of Chlamydomonas genome assembly version JGI v4.0, while the microarray chip is based on Augustus 10.2. The mapping was performed by using the Algal Functional Annotation tool (http://pathways.mcdb.ucla.edu/algal/id_conversion.html).

Chlamydomonas reinhardtii cultivation and metabolite extraction The Chlamydomonas reinhardtii strain CC503 used in this study is a cell–wall-deficient strain, generated from the wild type strain 137C+ in 1965 (Harris, 2001) and the genomic sequence was released in 2007 (Merchant et al., 2007). This strain, which can be grown on solid plates and in liquid culture, was pre-cultured in 100 ml cultures under mixotrophic and photoautotrophic conditions, using either TRIS-Acetate-Phosphate (TAP) or TRIS-Phosphate (TP) medium, respectively (Harris, 2001). Pre-cultures were grown on a rotary shaker with a speed of 120 rpm at 24°C at 30μE m-2 continuous light until a cell density of 5x106 cells ml-1 was reached.

We used 800 ml fermenters containing either TAP- or TP-medium inoculated with a starting culture density of 1x106 cells ml-1. The fermenters, illuminated with 200 μE, were constantly stirred and bubbled with 3% CO2 enriched air. The temperature was kept constant at 28°C using a water sheath and connected water bath. The pH was checked regularly to be stable at ~7. The cell counts and the cell diameter within each growing culture were determined from three independent measurements per sample, using a Beckman Coulter particle counter (Z2, Beckman Coulter, Krefeld, Germany) according to the manufacturer’s instructions.

For the sampling, 5-10 million cells were taken from the fermenter using a sterile syringe connected via sterile glass tubing to the fermenter system. The collected cells were pelleted by centrifugation for 5 min using 2000 xg and immediately frozen in liquid nitrogen. Samples for each condition, mixotrophic and photoautotrophic growth with and without 10 µM rapamycin, were taken at different eight time points (0, 30, 60, 90, 120, 240, 360, 480 min after rapamycin addition) in six replicates. The metabolites from these samples were extracted using a methyltert-butyl ester:methanol (3:1, v:v) extraction and the polar phase of these metabolites was, after

This article is protected by copyright. All rights reserved.

Accepted Article

appropriate derivatization, used to generate GC-ToF-MS metabolic profiles (Caldana et al., 2013; Giavalisco et al., 2011).

From the total number of detected metabolites 45 compounds were re-measured from commercially available standards and quantified by using serially diluted response curves (Fiehn et al., 2000; Kleessen et al., 2012) and the estimated absolute concentrations over the above mentioned time points were calculated using the mean of the replicates (Data S1). To fit the unit considered in the reconstruction we assumed the dry weight composition (0.024 and 0.020 chlorophyll in one fraction of dry weight in photoautotrophic and mixotrophic conditions, respectively) from Boyle et al. (Boyle and Morgan, 2009) to convert nmol g-1 chlorophyll from the measurements to mmol gDW-1. Furthermore, the RNA of three biological replicates, which was extracted using an RNAEasy Kit (QiaGen), was pooled and hybridized to an Agilent 4x44K custom-made array (Toepel et al., 2011). Samples from all different conditions for the microarray experiment were taken at six time points (0, 30, 60, 120, 240, 480 min). The intensities were first quantile-normalized (Bolstad et al., 2003) and the expression of each gene was quantified by the log2-transformed median of the three probes (Data S2). Missing time points for the transcriptomics data were interpolated by cubic splines to facilitate integration with the metabolite levels measured (Boor, 2001). All of the quantified metabolites and 982 of the genes could be mapped to the metabolic reconstruction. In addition, cell counts were obtained for the three replicates using a Beckman Coulter particle counter according to the manufacturer’s instructions. The final cell number and the average cell diameter were estimated from three independent counts per sample (Data S3). Integration of transcriptomics data We used the log2 transformed gene expression values (Data S2) to constraint the fluxes. We followed Colijn et al. (Colijn et al., 2009) for constructing the constraint vectors and . In this approach the allowed flux boundaries are constrained by integrating gene expression data in a continuous way (i.e., without imposing any thresholds). As described in detail in (Colijn et al., 2009), a linear mapping function was used between the measured expression data and the flux boundaries. Thereby, for each reaction which is catalyzed by only one gene, was set to the measured expression value in the respective condition. The upper bound for a reaction catalyzed by a protein complex was set to the minimum expression values of the associated genes. In addition, a reaction which can be catalyzed by any gene, from a given set, was set to the sum of their expression values. For instance, given three reactions, , for which is catalyzed by gene ); and

by

; and

by (i.e.,

or

(i.e., and

and

catalyze isozymes of the enzyme catalyzing

build a protein complex catalyzing

upper flux bounds included in the network reconstruction for each reaction

This article is protected by copyright. All rights reserved.

), let

be the

. Thus, the upper

Accepted Article

bounds modified by the transcripts , where

are set to

;

and

denotes the transcript level for gene . In contrast to the original E-

Flux method, if the gene catalyzing a transport reaction is known the bound was also constrained by their expression value(s). Finally, we normalized the maximum of the obtained values, over all considered conditions and time points, to 1000. This was done to relate the transcript bounds to the general bounds of the reconstruction. In detail, this is formalized for each reaction, , condition,

and time point, , by

the lower bounds,

, were set to

. Finally, for irreversible and

for reversible

reactions, respectively.

Integration of cell counts Cell counts were included to further constrain the system. To this end, the mean of the cell count from the three replicates was integrated into the modeling. The cell counts were measured in number of cells per ml at eight time points. (0, 30, 60, 90, 120, 240, 360 and 480 min after rapamycin treatment, Data S3). For the mixotrophic and photoautotrophic condition the ratio between cell counts in control and rapamycin was calculated for each time point. In addition the difference in cell counts between the time points for each condition was included. The resulting linear programming formulation of TREM-Flux applied in predicting the time-resolved flux phenotypes was:

This article is protected by copyright. All rights reserved.

(3 )

Accepted Article

where

corresponds either to photoautotrophic or mixotrophic condition, and

all other notation follows the generic formulation of TREM-Flux. The control and treatment of each growth condition is denoted by and , respectively. The coefficients for the biomass reaction in the objective function were obtained from the original publication (Chang et al., 2011) (Supplementary Table S10 therein).

To integrate cell counts in the time-resolved variant of E-Flux, the predictions of flux distributions over time were obtained in a single linear program for control and rapamycin treated cells in the analyzed growth conditions given by:

This article is protected by copyright. All rights reserved.

(4)

Accepted Article

where the notation follows the generic TREM-Flux formulation (see Eq. 3).

Time-resolved flux variability analysis We used a modified version of flux variability analysis (Mahadevan and Schilling, 2003; Burgard et al., 2001) to determine the ranges for each reactions in alternative optimal solutions of TREM-Flux. For each reaction , , where is the number of reactions and for each ,

where

, the variability of the fluxes is calculated by:

(5)

is the obtained value for the biomass reaction calculated previously from Eq.

3, and all other notation is equivalent to that of the generic TREM-Flux formulation. The first constraint is substituted by the steady-state condition in the flux variability analysis without inclusion of metabolite data.

Comparison of flux ranges in different modeling variants. The distribution of flux ranges over all time points and reactions in alternative optima, obtained from tFVA, of different modeling variants was evaluated. 17 different modeling variants were investigated. TREM-Flux and E-Flux were formulated as explained earlier. In addition, we randomly sampled 40, 30, 20, and 10 metabolites from the 45 measured. This sampling was repeated 10 times and the average flux range for each reactions and time point over all samples were used to obtain the further analyzed distribution. Moreover, varying ranges of metabolite changes were considered. Therefore, Eq. 3 was reformulated to

This article is protected by copyright. All rights reserved.

(6 )

Accepted Article

where

defines the considered changes, varying from 10 to 0.1 mmol gDW -1 for the different

modeling variants. The notation follows the general TREM-Flux application (see Eq. 3). This formulation was further modified, therefore, the flux bounds obtained from the transcript data ( , , , and ) were replaced by the flux bounds included in the metabolic reconstruction (Chang et al., 2011).

Distribution of flux variability within modeled bounds For the analysis of the predicted flux ranges, we focused on these reactions for which

,

i.e., they are not constraint to constants by the given boundaries. This was the case for 2127 out of the 2143 reactions. Furthermore, the range of flux variation was obtained by the absolute value of the difference between the flux bounds for each reaction ( in the investigated time points and conditions. In addition, we calculated the range of the bounds considered in the constraint-based modeling approach, by the absolute value of the difference of and ( ). These were then used to determine the percentage of the range obtained from tFVA from the flux bound range in the model (Figure 2C).

This article is protected by copyright. All rights reserved.

Accepted Article

Enrichment analysis Enriched metabolic functions were determined by using a hypergeometric test at significance level of 0.05. Adjusted p-values were obtained by the Benjamini and Hochberg procedure (Benjamini and Hochberg, 1995) via the multtest R package (Pollard et al., 2011).

Comparative analysis of flux distributions between different conditions For each time point, we conducted comparison of a single growth condition (i.e., mixotrophic and photoautotrophic) between the two scenarios (i.e., control and rapamycin treatment). For each comparison the absolute difference of the predicted flux distributions for the particular time pointed was calculated for a schematic presentation. Furthermore, the obtained differences were sorted in decreasing order and the top ranked 214 reactions, representing 10% of all reactions, respectively, were selected for further analysis. For the selected reactions, the enriched metabolic functions were obtained as detailed above.

Correlation of time-resolved behavior of metabolic functions For each of the 82 metabolic functions, the Kendall tau rank correlation coefficient between time-resolved fluxes for each reaction was determined in control and rapamycin treatment under both mixotrophic and photoautotrophic conditions. The behavior of a metabolic function was then quantified by the median of the Kendall tau correlations over all participating reactions of this function.

Analysis of differential gene expression and metabolite contents To conduct the analysis of differential expression, the expression levels of genes were normalized relative to the first time point in each scenario/growth condition. Since a single replicate for gene expression was available for each time point and condition, all sampled time points per condition were pooled for statistical analysis. To test the effect of the rapamycin treatment, differential expressed (DE) genes were identified by conducting multiple moderated ttests with the experimental condition as a factor (i.e., control and rapamycin treatment) for each growth (i.e., photoautotrophic and mixotrophic), separately (Smyth, 2005). In the case of metabolomics data, replicated time points were summarized by the median of the respective metabolite contents. Differentially behaving (DB) metabolites were then determined in the same manner as DE genes, without further data transformations. Finally, DE genes and DB

This article is protected by copyright. All rights reserved.

Accepted Article

metabolites are selected using FDR-corrected p-value at significance level of 0.01. For the DE genes, the enriched metabolic functions as well as MapMan bins were obtained based on a hypergeometric distribution test at significance level of 0.05. In addition, enriched metabolic functions were obtained for the DB metabolites based on the same test statistic. Adjusted pvalues were obtained by the Benjamini and Hochberg procedure (Benjamini and Hochberg, 1995) via the multtest R package (Pollard et al., 2011).

Implementation The pre-processing and integration of the data as well as the optimization was performed in MATLAB with the optimization platform TOMLAB v.7.9 using the CPLEX solver with default parameter and accuracy settings. Further analysis of the predictions was conducted in the R 2.15.2 environment for statistical computing.

Acknowledgments The authors would like to thank Prof. Dr. Lothar Willmitzer for discussions on the manuscript.

Supporting Information Supporting Figures Figure S1. Distribution of flux variability in time-resolved E-Flux within model bounds. Figure S2. Characteristics of the distributions of flux ranges from tFVA. Figure S3. Schematic representation of differential analysis of flux distributions. Supporting Tables Table S1. Results from tFVA for applying time-resolved E-Flux. Table S2. Enriched metabolic functions of blocked reactions. Table S3. tFVA comparison of three scenarios.

Table S4. Comparative analysis of flux distributions between different conditions. Table S5. Index for metabolic functions. This article is protected by copyright. All rights reserved.

Accepted Article

Table S6. Differentially behaving measured metabolites. Supporting Data Data S1. Absolute metabolites for the four scenarios. Data S2. Gene expression levels for the four scenarios. Data S3. Cell counts for the four scenarios. Data S4. Differentially expressed genes. Data S5. Enriched metabolic functions for differentially expressed genes. Data S6. Enriched MapMan bins for differentially expressed genes.

References Benjamini, Y. and Hochberg, Y. (1995) Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J. R. Stat. Soc. Ser. B, 57, 289 – 300. Blazier, A.S. and Papin, J.A. (2012) Integration of expression data in genome-scale metabolic network reconstructions. Front. Physiol., 3, 299. Bolstad, B.M., Irizarry, R.A., Astrand, M. and Speed, T.P. (2003) A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics, 19, 185–193. Boor, C. de (2001) A Practical Guide to Splines, New York: Springer. Boyle, N.R. and Morgan, J.A. (2009) Flux balance analysis of primary metabolism in Chlamydomonas reinhardtii. BMC Syst. Biol., 3, 4. Burgard, A.P., Vaidyaraman, S. and Maranas, C.D. (2001) Minimal reaction sets for Escherichia coli metabolism under different growth requirements and uptake environments. Biotechnol. Prog., 17, 791–7. Caldana, C., Li, Y., Leisse, A., Zhang, Y., Bartholomaeus, L., Fernie, A.R., Willmitzer, L. and Giavalisco, P. (2013) Systemic analysis of inducible target of rapamycin mutants reveal a general metabolic switch controlling growth in Arabidopsis thaliana. Plant J., 73, 897–909.

This article is protected by copyright. All rights reserved.

Accepted Article

Chang, R.L., Ghamsari, L., Manichaikul, A., et al. (2011) Metabolic network reconstruction of Chlamydomonas offers insight into light-driven algal metabolism. Mol. Syst. Biol., 7, 518. Chen, J., Zheng, X.F., Brown, E.J. and Schreiber, S.L. (1995) Identification of an 11-kDa FKBP12-rapamycin-binding domain within the 289-kDa FKBP12-rapamycin-associated protein and characterization of a critical serine residue. Proc. Natl. Acad. Sci., 92, 4947– 4951. Choi, J., Chen, J., Schreiber, S.L. and Clardy, J. (1996) Structure of the FKBP12-Rapamycin Complex Interacting with Binding Domain of Human FRAP. Science (80-. )., 273, 239– 242. Colijn, C., Brandes, A., Zucker, J., et al. (2009) Interpreting expression data with metabolic flux models: predicting Mycobacterium tuberculosis mycolic acid production. J. A. Papin, ed. PLoS Comput. Biol., 5, e1000489. Crespo, J.L., Díaz-Troya, S. and Florencio, F.J. (2005) Inhibition of target of rapamycin signaling by rapamycin in the unicellular green alga Chlamydomonas reinhardtii. Plant Physiol., 139, 1736–49. Díaz-Troya, S., Florencio, F.J. and Crespo, J.L. (2008) Target of rapamycin and LST8 proteins associate with membranes from the endoplasmic reticulum in the unicellular green alga Chlamydomonas reinhardtii. Eukaryot. Cell, 7, 212–22. Dobrenel, T., Marchive, C., Sormani, R., Moreau, M., Mozzo, M., Montané, M.-H., Menand, B., Robaglia, C. and Meyer, C. (2011) Regulation of plant growth and metabolism by the TOR kinase. Biochem. Soc. Trans., 39, 477–81. Fiehn, O., Kopka, J., Dörmann, P., Altmann, T., Trethewey, R.N. and Willmitzer, L. (2000) Metabolite profiling for plant functional genomics. Nat. Biotechnol., 18, 1157–61. Flamholz, A., Noor, E., Bar-Even, A., Liebermeister, W. and Milo, R. (2013) Glycolytic strategy as a tradeoff between energy yield and protein cost. Proc. Natl. Acad. Sci. U. S. A., 110, 10039–44. Gianchandani, E.P., Chavali, A.K. and Papin, J. a (2010) The application of flux balance analysis in systems biology. Wiley Interdiscip. Rev. Syst. Biol. Med., 2, 372–82. Giavalisco, P., Li, Y., Matthes, A., et al. (2011) Elemental formula annotation of polar and lipophilic metabolites using (13) C, (15) N and (34) S isotope labelling, in combination with high-resolution mass spectrometry. Plant J., 68, 364–76. Harris, E.H. (2001) Chlamydomonas as a model organism. Annu. Rev. Plant Physiol. Plant Mol. Biol., 52, 363–406.

This article is protected by copyright. All rights reserved.

Accepted Article

Henry, C.S., Broadbelt, L.J. and Hatzimanikatis, V. (2007) Thermodynamics-based metabolic flux analysis. Biophys. J., 92, 1792–805. Jamshidi, N. and Palsson, B.Ø. (2010) Mass action stoichiometric simulation models: incorporating kinetics and regulation into stoichiometric models. Biophys. J., 98, 175–85. Kitano, H. (2002) Systems biology: a brief overview. Science (80-. )., 295, 1662–1664. Kleessen, S., Araújo, W.L., Fernie, A.R. and Nikoloski, Z. (2012) Model-based Confirmation of Alternative Substrates of Mitochondrial Electron Transport Chain. J. Biol. Chem., 287, 11122–31. Kleessen, S. and Nikoloski, Z. (2012) Dynamic regulatory on/off minimization for biological systems under internal temporal perturbations. BMC Syst. Biol., 6, 16. Komeili, A., Wedaman, K.P., O’Shea, E.K. and Powers, T. (2000) Mechanism of metabolic control. Target of rapamycin signaling links nitrogen quality to the activity of the Rtg1 and Rtg3 transcription factors. J. Cell Biol., 151, 863–78. Laplante, M. and Sabatini, D.M. (2012) mTOR signaling in growth control and disease. Cell, 149, 274–93. Lewis, N.E., Nagarajan, H. and Palsson, B.O. (2012) Constraining the metabolic genotypephenotype relationship using a phylogeny of in silico methods. Nat. Rev. Microbiol., 10, 291–305. Loewith, R. and Hall, M.N. (2011) Target of rapamycin (TOR) in nutrient signaling and growth control. Genetics, 189, 1177–201. Luo, R., Wei, H., Ye, L., et al. (2009) Photosynthetic metabolism of C3 plants shows highly cooperative regulation under changing environments: a systems biological analysis. Proc. Natl. Acad. Sci. U. S. A., 106, 847–52. Luo, R.Y., Liao, S., Tao, G.Y., Li, Y.Y., Zeng, S., Li, X.Y. and Luo, Q. (2006) Dynamic analysis of optimality in myocardial energy metabolism under normal and ischemic conditions. Mol Syst Biol, 2, 2006 0031. Mahadevan, R., Edwards, J.S. and Doyle, J.F. (2002) Dynamic Flux Balance Analysis of Diauxic Growth in Escherichia coli. Biophys J, 83, 1331–1340. Mahadevan, R. and Schilling, C.H. (2003) The effects of alternate optimal solutions in constraint-based genome-scale metabolic models. Metab. Eng., 5, 264–276. Masakapalli, S.K., Kruger, N.J. and Ratcliffe, R.G. (2013) The metabolic flux phenotype of heterotrophic Arabidopsis cells reveals a complex response to changes in nitrogen supply. Plant J., 74, 569–82.

This article is protected by copyright. All rights reserved.

Accepted Article

Merchant, S.S., Prochnik, S.E., Vallon, O., et al. (2007) The Chlamydomonas genome reveals the evolution of key animal and plant functions. Science (80-. )., 318, 245–50. Moreau, M., Azzopardi, M., Clément, G., et al. (2012) Mutations in the Arabidopsis homolog of LST8/GβL, a partner of the target of Rapamycin kinase, impair plant growth, flowering, and metabolic adaptation to long days. Plant Cell, 24, 463–81. Nöh, K. and Wiechert, W. (2006) Experimental design principles for isotopically instationary 13C labeling experiments. Biotechnol. Bioeng., 94, 234–51. Orth, J.D., Thiele, I. and Palsson, B.Ø. (2010) What is flux balance analysis? Nat. Biotechnol., 28, 245–8. Pollard, K.S., Gilbert, H.N., Ge, Y., Taylor, S. and Dudoit, S. (2011) multtest: Resamplingbased multiple hypothesis testing. Raman, K. and Chandra, N. (2009) Flux balance analysis of biological systems: applications and challenges. Br. Bioinform, 10, 435–449. Reed, J.L. (2012) Shrinking the metabolic solution space using experimental datasets. J. A. Papin, ed. PLoS Comput. Biol., 8, e1002662. Ren, M., Venglat, P., Qiu, S., et al. (2012) Target of rapamycin signaling regulates metabolism, growth, and life span in Arabidopsis. Plant Cell, 24, 4850–74. Schmidt, B.J., Ebrahim, A., Metz, T.O., Adkins, J.N., Palsson, B.Ø. and Hyduke, D.R. (2013) GIM3E: condition-specific models of cellular metabolism developed from metabolomics and expression data. Bioinformatics, btt493–. Smyth, G.K. (2005) Limma: linear models for microarray data. In R. Gentleman, V. J. Carey, W. Huber, R. A. Irizarry, and S. Dudoit, eds. Bioinformatics and computational biology solutions using R and Bioconductor. Springer, pp. 397–420. Thimm, O., Bläsing, O., Gibon, Y., et al. (2004) MAPMAN: a user-driven tool to display genomics data sets onto diagrams of metabolic pathways and other biological processes. Plant J., 37, 914–39. Toepel, J., Albaum, S.P., Arvidsson, S., Goesmann, A., Russa, M. la, Rogge, K. and Kruse, O. (2011) Construction and evaluation of a whole genome microarray of Chlamydomonas reinhardtii. BMC Genomics, 12, 579. Varma, A. and Palsson, B.O. (1994a) Metabolic Flux Balancing: Basic Concepts, Scientific and Practical Use. Nat. Biotechnol., 12, 994–998.

This article is protected by copyright. All rights reserved.

Accepted Article

Varma, A. and Palsson, B.O. (1994b) Stoichiometric flux balance models quantitatively predict growth and metabolic by-product secretion in wild-type Escherichia coli W3110. Appl. Environ. Microbiol., 60, 3724–31. Vézina, C., Kudelski, A. and Sehgal, S.N. (1975) Rapamycin (AY-22,989), a new antifungal antibiotic. I. Taxonomy of the producing streptomycete and isolation of the active principle. J. Antibiot. (Tokyo)., 28, 721–6. Westerhoff, H. V, Winder, C., Messiha, H., Simeonidis, E., Adamczyk, M., Verma, M., Bruggeman, F.J. and Dunn, W. (2009) Systems biology: the elements and principles of life. FEBS Lett., 583, 3882–90. Wullschleger, S., Loewith, R. and Hall, M.N. (2006) TOR signaling in growth and metabolism. Cell, 124, 471–84. Yizhak, K., Benyamini, T., Liebermeister, W., Ruppin, E. and Shlomi, T. (2010) Integrating quantitative proteomics and metabolomics with a genome-scale metabolic network model. Bioinformatics, 26, i255–60.

Figure legends Figure 1. Schematic representation of TREM-Flux. Time-resolved transcriptomics and metabolomics data (i.e., contents) are mapped to genome-scale metabolic reconstruction. (A) The data integration of gene expression data is performed based on the E-Flux method (Colijn et al., 2009). (B) Following the E-Flux method the flux bounds included in the metabolic network (min/max flux model) are replaced by the bounds from the transcript levels, and . (C) The integration of the time-resolved absolute metabolite levels assumes a transient state which may entail a change in the content of a metabolite between consecutive time points. (D) For a metabolite present in more than one compartment the sum of rows in the stoichiometric matrix is determined. (E) The generic TREM-Flux method combines the integration of metabolomics and transcriptomics data, optimizing a given objective over time. Figure 2. Time-resolved flux variability analysis (tFVA). Schematic representation of (A) flux variability analysis performed for each time point and condition, (B) reactions of zero, constant (non-zero) and varying fluxes over all time points for each condition (given in orange, blue and purple, respectively) and of (C) calculations of the percentage of flux variability within the model bounds. (D) Results of the distribution of flux variability within the bounds of flux predictions. The range of flux variation (maximum-minimum) in the alternative optima is obtained by tFVA. This range is then compared with range from the transcriptomics-based bounds ( and ) to obtain the percentage of variability. The results are presented for

This article is protected by copyright. All rights reserved.

Accepted Article

applying tFVA with TREM-Flux integrating metabolomics and transcriptomics data. The percentages are divided into ten bins, each covering 10 percent. The results are shown separately for each time point ( - ) for all investigated growth conditions: control (green) and rapamycin treatment (red) under mixotrophic growth as well as control (blue) and rapamycin treatment (yellow) under autotrophic growth. Figure 3. tFVA comparison for the control condition under mixotrophic growth. (A) The distributions of flux ranges for the control condition of mixotrophic growth are shown for (A) 12 different modeling variants including TREM-Flux, E-Flux, the sampling of different numbers of metabolites as well as a varying range of metabolite changes; and (B) 5 modeling variants with a varying range of metabolite changes or assuming a steady state by using only the flux bounds included in the metabolic reconstruction. The range of metabolite changes in the modeling variants which do not include any measured metabolites are modeled by the constraint (see Experimental Procedures for details).

Figure 4. Visualization of differential metabolic functions. Enriched metabolic functions for the comparison of flux distributions between investigated growth conditions are presented. Given two consecutive time points, the absolute difference of the corresponding flux distributions is determined, and the reactions are ranked according to the resulting values (see Figure S3A). The reactions present in the top 10% are examined and an enrichment analysis is conducted to determine enriched metabolic functions in these selected reactions. The presented enriched metabolic functions are obtained based on the hypergeometric test (adjusted p-value 0.05, see Experimental Procedures). The time point and condition wise results are presented in Table S4). The image was created by modifying a MapMan Template (Thimm et al., 2004). Figure 5. Visualization of time-resolved behavior of metabolic functions. A schematic overview of the median correlation values characterizing the time-resolved behavior of metabolic functions is presented. A metabolic function having a negative median for both investigated growth conditions is presented in red. If a metabolic function has only a negative correlation median under mixotrophic conditions and positive (or zero) in photoautotrophic it is marked in green (orange for the reversal of conditions). A metabolic function with a positive median in both conditions is shown in blue. The image was created by modifying a MapMan Template (Thimm et al., 2004).

This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.