Abstract
During the last decade genome sequencing has experienced a rapid technological development resulting in numerous sequencing projects and applications in life science. In plant molecular biology, the availability of sequence data on whole genomes has enabled the reconstruction of metabolic networks. Enzymatic reactions are predicted by the sequence information. Pathways arise due to the participation of chemical compounds as substrates and products in these reactions. Although several of these comprehensive networks have been reconstructed for the genetic model plant Arabidopsis thaliana, the integration of experimental data is still challenging. Particularly the analysis of subcellular organization of plant cells limits the understanding of regulatory instances in these metabolic networks in vivo. In this study, we develop an approach for the functional integration of experimental high-throughput data into such large-scale networks. We present a subcellular metabolic network model comprising 524 metabolic intermediates and 548 metabolic interactions derived from a total of 2769 reactions. We demonstrate how to link the metabolite covariance matrix of different Arabidopsis thaliana accessions with the subcellular metabolic network model for the inverse calculation of the biochemical Jacobian, finally resulting in the calculation of a matrix which satisfies a Lyaponov equation. In this way, different strategies of metabolite compartmentation and involved reactions were identified in the accessions when exposed to low temperature.
Introduction
The rapidly increasing knowledge about whole plant genome sequences represents a corner stone in the understanding of plant metabolism. Next-generation sequencing (NGS) technologies have been developed allowing for the fast and cheap production of huge sets of genome data sequences (Metzker, ). By functional genome annotation, the information about coded proteins is derived. Based on these postulated gene functions, enzymatic reactions can be predicted, e.g., representing metabolite interconversions or transport processes. Genome-scale metabolic reconstructions (GEMs) can be described as advanced functional annotations (De Oliveira Dal'molin and Nielsen, ). Here, the functional gene annotations are combined with a network topology which is derived by connecting the metabolic educts and products of the predicted reactions. The stoichiometric matrix N characterizes each metabolic interconversion and transition in a metabolic network and is typically organized in such way that rows of the matrix represent the metabolites while reactions, i.e., connections, represent the columns of the matrix. The stoichiometric coefficients of each reaction are then given in the matrix cells. GEMs have been published for multiple organisms (an overview is given in Collakova et al., ) and comprehensive protocols for all stages of the reconstruction process are available (Thiele and Palsson, ). Due to the huge metabolic coverage of GEMs, which, in principle, comprises all metabolic interactions known so far from genome annotation, strategies for genome-scale experiments are needed for efficient validation of model outputs. In this context, flux measurements and high-throughput measurements of the transcriptome, proteome and metabolome play a crucial role. Hence, technologies like transcriptomics, proteomics and metabolomics are central to systems biology approaches aiming at the system wide understanding of biological networks (Weckwerth, ,; Blazier and Papin, ; Nägele and Weckwerth, ). Recently, we connected data from metabolomics experiments to a simplified metabolic network structure of leaf primary metabolism in Arabidopsis thaliana to characterize metabolic shifts during cold exposure (Doerfler et al., ). This allowed us to differentiate short and long term metabolic response to low temperature and to identify key points of regulation such as the interface of primary and secondary metabolism mediated by the shikimic acid pathway. This model did not include any subcellular compartmentation of metabolism. Plant cells, however, show a high degree of compartmentation resulting in a complex biochemical network comprising metabolite interconversions as well as intracellular transport processes (Lunn, ). Therefore it is not surprising that physiological responses to a changing environment could be related to subcellular metabolic reprogramming (Knaupp et al., ; Schulze et al., ; Nägele and Heyer, ). To enhance the comprehensive biochemical and physiological output of studies analyzing plant-environment interactions a platform would be desirable focusing the linkage of experimental data with the underlying subcellular metabolic network structure.
Here, we present a workflow aiming at the development of such a platform. Based on a recently published metabolic network reconstruction accounting for subcellular organization of leaf metabolism in Arabidosis thaliana (Mintz-Oron et al., ), we derived a metabolic model structure including all metabolic intermediates which are accessible by experimental high-throughput measurements. The model comprises all subcellular compartments of leaf cells which can be robustly analyzed by a non-aqueous fractionation technique. For efficient data integration we present a strategy of mathematical modeling which we prove to be successful in providing a comprehensive overview of metabolic changes induced by a perturbed plant-environment interaction.
Material and methods
Adaptation of a genome-scale metabolic reconstruction model to a set of experimentally accessible data
The model adaptation was performed starting with the original metabolic reconstruction model for (juvenile) leaf metabolism, which was derived by Mintz-Oron and co-workers (Mintz-Oron et al., ). This reconstruction model comprises 7 compartments which are the cytoplasm, endoplasmic reticulum, golgi apparatus, mitochondrion, peroxisome, plastid and vacuole with a total of 2463 metabolic intermediates and 2769 reactions (Mintz-Oron et al., ). In a first step, we reduced the compartments in the model to the cytoplasm, plastid and vacuole which can experimentally be analyzed from the same sample using the non-aqueous fractionation (NAF) method (Nägele and Heyer, ). In a second step, we reduced the metabolic intermediates of these compartments to a set of metabolites which are experimentally accessible by a gas chromatography coupled to mass spectrometry (GC-MS) analysis. This reduction step was performed by comparison of metabolites in the model to metabolites in the Golm Metabolome Database for GC-MS based metabolite profiling (Hummel et al., ). Additionally, the model contains the plastidial starch pool as well as CO2. In the following step, all intermediates in the reduced model were (re)connected manually according to the reactions described in the original reconstruction model and no further reactions were added. Hence, the reduced model describes a subset of the metabolic connections in the original model with a lower degree of detail. A step of this reduction procedure is exemplarily shown in Figure S1. Finally, our reduced model contained 3 compartments, 524 metabolic intermediates and 548 metabolic interactions. The reduced model of leaf metabolism is provided in the supplements in Systems Biology Markup Language (SBML; File ‘Model_S3.xml’). Model reduction, connection and graphical evaluation was performed using the open source software COPASI (Version 4.8; http://www.copasi.org) (Hoops et al., ) and CellDesigner™ (Version 4.3; http://www.celldesigner.org/) (Funahashi et al., ).
Data integration and jacobian matrix calculation
A metabolic interaction matrix was derived from the reduced model describing all metabolic interactions in the model. Hence, the metabolic interaction matrix represents a simplified version of the stoichiometric matrix of the original metabolic reconstruction model. This metabolic interaction matrix was applied for inverse calculation of the Jacobian matrix of the metabolic model. The calculation procedure was based on an algorithm, which is implemented in the metabolomics toolbox COVAIN (Sun and Weckwerth, ), solving (Eq. 1) by applying a total least square optimization procedure. To test the reliability of our calculations we applied the inverse calculation on a data set on subcellular carbohydrate distribution from non-cold exposed and 7 days cold exposed Arabidopsis leaf samples which were published recently (Nägele and Heyer, ). This experimental data set contained 5 biological replicates for each metabolite concentration (Nägele and Heyer, ). The reduced model structure, i.e., the metabolic interaction matrix, was adapted to the metabolite pools which had been experimentally analyzed by Nägele and Heyer. The covariance matrix C was built from the experimental data on subcellular carbohydrate concentration and linked to the underlying biochemical system by the following equation (Steuer et al., ; Sun and Weckwerth, ): Here, J represents the Jacobian matrix and D is the fluctuation matrix. All diagonal entries of D were randomly drawn from a standard normal distribution. In general, the Jacobian matrix characterizes the local dynamics at a steady state condition. In context of a metabolic network, entries of the Jacobian J represent the elasticities of reaction rates to any change of the metabolite concentrations which are characterized by equation 2: N is the stoichiometric matrix, r represents the rates for each reaction and M represents metabolite concentrations. In our approach, we replaced the stoichiometric matrix N by the reduced metabolic interaction matrix. This results in a Jacobian matrix referring to the underlying stoichiometric simplification. Although N was changed, the term Jacobian matrix can still be used as it directly refers to the first partial derivative of (now simplified) metabolic functions to changes in metabolite levels. This corresponds to the general definition of the Jacobian matrix in context of the mean value theorem of vector functions, which is provided elsewhere (Strehmel et al., 2012). The Jacobian matrices Ja and Jb, which describe two different metabolic states, were calculated 105 times each. Medians of calculated Jacobians were normalized to the square of the interquartile distance in order to increase the median-to-noise ratio of the inverse calculations (Eq. 1). To compare two different metabolic states, we determined the absolute values of the differential Jacobian, dJij, abs, defining the relative change of the two normalized Jacobians Ja,norm and Jb,norm which are associated with different treatments or genotypes: All calculations of Jacobian matrices as well as statistical analyses (t-tests) were performed using the numerical software environment Matlab® (V7.12.0 R2011a).
Results
A subcellular metabolic network model comprising experimentally accessible intermediates from high-throughput analysis of one sample
The reduction process of the metabolic reconstruction network of leaf metabolism (Mintz-Oron et al., ) resulted in a metabolic network model comprising the compartments cytoplasm, plastid and vacuole with a total of 524 metabolic intermediates and 548 metabolic interactions. In addition to metabolite pools which are accessible by a GC-MS experiment, the model also comprises the plastidial starch pool and CO2. These intermediates are experimentally accessible by, for example, GC-MS analysis of starch hydrolysate, photometric assays (starch) and infrared gas analysis (CO2) as previously described (Wienkoop et al., ; Nägele and Heyer, ; Valledor et al., ).
While the metabolic network reconstruction of subcellular leaf metabolism resulted in a stoichiometric matrix derived from genome sequence information (Mintz-Oron et al., ), this stoichiometric matrix was changed during the reduction process by deleting and reconnecting components. Hence, the resulting matrix of the reduced model differs a lot from the original stoichiometric information and therefore it is termed as the metabolic interaction matrix. It describes all experimentally accessible metabolic interactions by superpathways. These superpathways implicitly describe all metabolic steps which are involved in a metabolic interaction. If all metabolic intermediates of a reaction or pathway are included in both the original and reduced model then the entries of the metabolic interaction matrix equal the entries of the stoichiometric matrix.
Application of the reduced metabolic interaction matrix to analyse regulation of carbohydrate compartmentation
To test the applicability of the reduced model structure to analyse subcellular metabolic interaction, we used the underlying metabolic interaction matrix for inverse calculation of Jacobian matrices to a recently published data set on subcellular carbohydrate compartmentation (Nägele and Heyer, ). We derived a specific metabolic interaction matrix which contained only those pools which were measured by Nägele and Heyer by further reduction of the metabolic network structure. Finally, the model described the pools of sucrose, raffinose, glucose, fructose and starch in the cyotsol, plastid and vacuole as well as their metabolic interactions. In addition to these measured pools we also included all direct metabolic interactions which were not experimentally analyzed, i.e., all metabolites which were connected to the measured pools by one further metabolic interaction. Including those interactions may be helpful for estimating the impact of modeling results with respect to the metabolite interconversions which are directly connected to the set of metabolites which were experimentally analyzed. Yet, this model extension does not affect the modeling results but is rather suggested to support the interpretation. The resulting metabolic interaction matrix describing all reactions and transports of the measured metabolites was used for the inverse calculation of a Jacobian matrix. Calculations were performed exemplarily on experimental data sets for three natural accessions of Arabidopsis thaliana, C24, Rschew (Rsch) and Tenela (Te) before (non-acc) and after a 7 day exposure to 4°C (7d acc) (Nägele and Heyer, ). Results of inverse calculations indicated a main difference between C24 and Te to exist in vacuolar hexose interactions as well as in the cytosolic and plastidial sucrose interaction (Figures 1 A,B; p < 0.05). Comparison of calculations for 7d acc samples of C24 and Rsch revealed a similar and still significant interaction pattern, but to a lower extent (Figure 1C; p < 0.05). For non-acc samples, the comparison between C24 and Rsch revealed no further signficant differences than the comparison of C24 and Te (data not shown). To exclude the possiblity that these differences in metabolic interactions represent artifacts from the model reduction process itself, the reactions in the reduced model structure were compared to the original metabolic reconstruction work (Table S2, yellow marked lines). In the reduced model, these interactions still correspond to the original model and were not a consequence of inserted reactions. Hence, they are directly linked to the enzymatic reactions described in the original metabolic network reconstruction.
Figure 1
Discussion
The development of genome-scale metabolic network models has become a central approach to approximate the topology of metabolic networks in vivo. While such networks have successfully been applied in several studies to analyse network properties (Poolman et al.,
Our reduced metabolic model of subcellular primary leaf metabolism in Arabidopsis thaliana intends to provide a platform aiming at the integration of as many as possible experimental data to metabolic network information. Although the experimental coverage of metabolic intermediates may significantly deviate from the metabolic coverage of the model, the model can be adjusted to the available experimental data set. We have exemplarily performed such an adjustment to an experimental data set on the central carbohydrate metabolism which was published previously (Nägele and Heyer,
Previous studies have shown that NAF coupled to high-throughput analysis enables the comprehensive characterization of a metabolic homeostasis (Klie et al.,
Figure 2

A workflow for deriving subcellular metabolic network structures according to experimental data sets. In the present study, the first step of deriving a genome-scale metabolic network reconstruction model was not performed, but a published reconstruction work was applied instead (Mintz-Oron et al.,
Statements
Author contributions
Thomas Nägele performed model programming, calculation, modeling and wrote the manuscript. Wolfram Weckwerth wrote the manuscript. Thomas Nägele and Wolfram Weckwerth performed the design of the study.
Acknowledgments
We would like to thank the members of the Department Ecogenomics and Systems Biology at the University of Vienna for critical and helpful discussions. We thank the EU-Marie-Curie ITN MERIT (GA 2010-264474) for financial support of Thomas Nägele. Finally, we would like to thank the reviewers for their comments and suggestions to improve the manuscript.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: http://www.frontiersin.org/Journal/10.3389/fpls.2013.00541/abstract
Figure S1Schematic overview of the model reduction process. In contrast to metabolite (B), metabolites (A) and (C) are accessible by a GC-MS measurement and are kept in the reduced model structure.
Table S2Reactions included in the reduced model (column B) are compared to reactions of the original reconstruction network (column E). Reactions which have been identified to be affected differentially during cold exposure of sensitive and tolerant Arabidopsis accessions are marked in yellow (referring to Figure 1).
Model S3Reduced model structure in format of Systems Biology Markup Language (xml-file).
References
1
BlazierA. S.PapinJ. A. (2012). Integration of expression data in genome-scale metabolic network reconstructions. Front. Physiol. 3:299. 10.3389/fphys.2012.00299
2
CollakovaE.YenJ. Y.SengerR. S. (2012). Are we ready for genome-scale modeling in plants?Plant Sci. 191–192, 53–70. 10.1016/j.plantsci.2012.04.010
3
De Oliveira Dal'molinC. G.NielsenL. K. (2013). Curr. Opin. Biotechnol. 24, 271–277. 10.1016/j.copbio.2012.08.007
4
DoerflerH.LyonD.NägeleT.SunX.FragnerL.HadacekF.et al. (2013). Granger causality in integrated GC-MS and LC-MS metabolomics data reveals the interface of primary and secondary metabolism. Metabolomics9, 564–574. 10.1007/s11306-012-0470-0
5
FeistA. M.PalssonB. O. (2008). The growing scope of applications of genome-scale metabolic reconstructions using Escherichia coli. Nat. Biotechnol. 26, 659–667. 10.1038/nbt1401
6
FunahashiA.MorohashiM.KitanoH.TanimuraN. (2003). CellDesigner: a process diagram editor for gene-regulatory and biochemical networks. BIOSILICO1, 159–162. 10.1016/S1478-5382(03)02370-9
7
GibonY.BlaesingO. E.HannemannJ.CarilloP.HöhneM.HendriksJ. H. M.et al. (2004). A robot-based platform to measure multiple enzyme activities in Arabidopsis using a set of cycling assays: comparison of changes of enzyme activities and transcript levels during diurnal cycles and in prolonged darkness. Plant Cell16, 3304–3325. 10.1105/tpc.104.025973
8
HoopsS.SahleS.GaugesR.LeeC.PahleJ.SimusN.et al. (2006). COPASI - a COmplex PAthway SImulator. Bioinformatics22, 3067–3074. 10.1093/bioinformatics/btl485
9
HummelJ.SelbigJ.WaltherD.KopkaJ. (2007). The golm metabolome database: a database for GC-MS based metabolite profiling, in Metabolomics, eds NielsenJ.JewettM. (Berlin; Heidelberg: Springer), 75–95. 10.1093/bioinformatics/bti236
10
KlieS.KruegerS.KrallL.GiavaliscoP.FluggeU. I.WillmitzerL.et al. (2011). Analysis of the compartmentalized metabolome - a validation of the non-aqueous fractionation technique. Front. Plant Sci. 2:55. 10.3389/fpls.2011.00055
11
KnauppM.MishraK. B.NedbalL.HeyerA. G. (2011). Evidence for a role of raffinose in stabilizing photosystem II during freeze-thaw cycles. Planta234, 477–486. 10.1007/s00425-011-1413-0
12
KruegerS.GiavaliscoP.KrallL.SteinhauserM.-C.BüssisD.UsadelB.et al. (2011). A topological map of the compartmentalized Arabidopsis thaliana leaf metabolome. PLoS ONE6:e17806. 10.1371/journal.pone.0017806
13
LunnJ. E. (2007). Compartmentation in plant metabolism. J. Exp. Bot. 58, 35–47. 10.1093/jxb/erl134
14
MayR. M. (1974). Stability and complexity in model ecosystems. Princeton, NJ: Princeton University Press.
15
MetzkerM. L. (2010). Sequencing technologies - the next generation. Nat. Rev. Genet. 11, 31–46. 10.1038/nrg2626
16
Mintz-OronS.MeirS.MalitskyS.RuppinE.AharoniA.ShlomiT. (2012). Reconstruction of Arabidopsis metabolic network models accounting for subcellular compartmentalization and tissue-specificity. Proc. Natl. Acad. Sci. U.S.A. 109, 339–344. 10.1073/pnas.1100358109
17
NägeleT.WeckwerthW. (2013). Eigenvalues of Jacobian matrices report on steps of metabolic reprogramming in a complex plant-environment interaction. Appl. Math. 4, 44–49. 10.4236/am.2013.48A007
18
NägeleT.HeyerA. G. (2013). Approximating subcellular organisation of carbohydrate metabolism during cold acclimation in different natural accessions of Arabidopsis thaliana. New Phytol. 198, 777–787. 10.1111/nph.12201
19
NägeleT.StutzS.HörmillerI. I.HeyerA. G. (2012). Identification of a metabolic bottleneck for cold acclimation in Arabidopsis thaliana. Plant J. 72, 102–114. 10.1111/j.1365-313X.2012.05064.x
20
NägeleT.WeckwerthW. (2012). Mathematical modeling of plant metabolismŕfrom reconstruction to prediction. Metabolites2, 553–566. 10.3390/metabo2030553
21
PoolmanM. G.MiguetL.SweetloveL. J.FellD. A. (2009). A genome-scale metabolic model of Arabidopsis and some of its properties. Plant Physiol. 151, 1570–1581. 10.1104/pp.109.141267
22
RohwerJ. M. (2012). Kinetic modelling of plant metabolic pathways. J. Exp. Bot. 63, 2275–2292. 10.1093/jxb/ers080
23
SchaberJ.LiebermeisterW.KlippE. (2009). Nested uncertainties in biochemical models. IET Syst. Biol. 3, 1–9. 10.1049/iet-syb:20070042
24
SchulzeW. X.SchneiderT.StarckS.MartinoiaE.TrentmannO. (2012). Cold acclimation induces changes in Arabidopsis tonoplast protein abundance and activity and alters phosphorylation of tonoplast monosaccharide transporters. Plant J. 69, 529–541. 10.1111/j.1365-313X.2011.04812.x
25
SteuerR. (2007). Computational approaches to the topology, stability and dynamics of metabolic networks. Phytochemistry68, 2139–2151. 10.1016/j.phytochem.2007.04.041
26
SteuerR.KurthsJ.FiehnO.WeckwerthW. (2003). Observing and interpreting correlations in metabolomic networks. Bioinformatics19, 1019–1026. 10.1093/bioinformatics/btg120
27
StrehmelK.WeinerR.PodhaiskyH. (2012). Numerik gewöhnlicher Differentialgleichungen. Wiesbaden: Vieweg+Teubner Verlag.
28
SunX. L.WeckwerthW. (2012). COVAIN: a toolbox for uni- and multivariate statistics, time-series and correlation network analysis and inverse estimation of the differential Jacobian from metabolomics covariance data. Metabolomics8, S81–S93. 10.1007/s11306-012-0399-3
29
ThieleI.PalssonB. O. (2010). A protocol for generating a high-quality genome-scale metabolic reconstruction. Nat. Protoc. 5, 93–121. 10.1038/nprot.2009.203
30
ValledorL.FuruhashiT.HanakA. M.WeckwerthW. (2013). Systemic cold stress adaptation of Chlamydomonas reinhardtii. Mol. Cell. Proteomics12, 2032–2047. 10.1074/mcp.M112.026765
31
WeckwerthW. (2011a). Green systems biology - from single genomes, proteomes and metabolomes to ecosystems research and biotechnology. J. Proteomics75, 284–305. 10.1016/j.jprot.2011.07.010
32
WeckwerthW. (2011b). Unpredictability of metabolism-the key role of metabolomics science in combination with next-generation genome sequencing. Anal. Bioanal. Chem. 400, 1967–1978. 10.1007/s00216-011-4948-9
33
WienkoopS.WeissJ.MayP.KempaS.IrgangS.Recuenco-MunozL.et al. (2010). Targeted proteomics for Chlamydomonas reinhardtii combined with rapid subcellular protein fractionation, metabolomics and metabolic flux analyses. Mol. Biosyst. 6, 1018–1031. 10.1039/b920913a
34
WormitA.TrentmannO.FeiferI.LohrC.TjadenJ.MeyerS.et al. (2006). Molecular identification and physiological characterization of a novel monosaccharide transporter from Arabidopsis involved in vacuolar sugar transport. Plant Cell18, 3476–3490. 10.1105/tpc.106.047290
Summary
Keywords
plant systems biology, genome-scale models, metabolomics, inverse calculation, mathematical modeling, subcellular compartmentation
Citation
Nägele T and Weckwerth W (2013) A workflow for mathematical modeling of subcellular metabolic pathways in leaf metabolism of Arabidopsis thaliana. Front. Plant Sci. 4:541. doi: 10.3389/fpls.2013.00541
Received
23 August 2013
Accepted
12 December 2013
Published
24 December 2013
Volume
4 - 2013
Edited by
Sonia Osorio, Málaga University, Spain
Reviewed by
Alexander R. van der Krol, Wageningen University, Netherlands; Kansuporn Sriyudthsak, RIKEN Center for Sustainable Resource Science, Japan
Copyright
© 2013 Nägele and Weckwerth.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Thomas Nägele, Department Ecogenomics and Systems Biology, University of Vienna, Althanstr. 14, 1090 Vienna, Austria e-mail: thomas.naegele@univie.ac.at
This article was submitted to Plant Systems Biology, a section of the journal Frontiers in Plant Science.
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.