Integration of crop growth model and constraint-based metabolic model predicts metabolic changes over rice plant development under water-limited stress

Rice is a major staple food worldwide and understanding its metabolism is essential for improving crop yield and quality, especially in a changing climate. Constraint-based modelling is an established method for studying metabolism at a systems level, but one of its limitations is the difficulty in directly integrating certain environmental factors, such as water potential, to the model for predicting metabolic changes in response to environmental changes. Here, we developed a framework to integrate a crop growth model and an upgraded diel multi-organ genome-scale metabolic model of rice to predict the metabolism of rice growth under normal and water-limited conditions. Our model was able to predict distinct metabolic adaptations under water-limited stress compared to normal condition across multiple developmental stages. Our modelling results of dynamic changes in metabolism over the whole-plant growth period highlighted key features of rice metabolism under water-limited stress including early leaf senescence, reduction in photosynthesis and significant nitrogen assimilation during grain filling.


IN TROD UCTION
Rice (Oryza sativa) is one of the world's most important food crop and highly popular in Asia, where it serves as the cereal for major nutritional intake (Fresco 2005). Around 50 % of the world's population depend on rice for their daily nourishment where 60 % of their daily caloric intake comes from it (Khush 2003;Gnanamanickam 2009). In this regard, the improvement in quality and yield of rice grain is one of the most important areas in the field of rice research (Khush 2003). Several breakthroughs have been made including the improvement of lodging resistance (Ookawa et al. 2010), improved water use efficiency (Karaba et al. 2007), vitamin-enriched 'golden rice' (Beyer et al. 2002) etc. In addition, there are multiple ongoing research projects that aim to improve photosynthesis (Long et al. 2006;Hibberd et al. 2008), including the engineering of C 4 photosynthesis in rice (Zhu et al. 2010). In the decades to come, improvement in rice cultivation is extremely important in feeding billions of people (Hibberd et al. 2008). In addition, rice is also a target plant to be used as biofuel feedstock given that a large quantity of rice residues remains as waste (Abbas and Ansumali 2010;Ranjan and Moholkar 2013). To improve the quality and quantity of rice production, it is important to have a holistic view on metabolism during the development of the rice plant including seed production. This will enrich our understanding of how carbon-and nitrogen-rich compounds are transported into seeds for the synthesis of major biomass components (such as protein and carbohydrate) and cofactors (including vitamin B2, B3, B5, B6 and folate). In light of the change in global climate (Reynolds 2010), it is especially important to understand the interactions between plant metabolic processes and climate-driven factors in order to strategically and efficiently target and engineer the naturally evolved metabolic processes to meet global food demand (Dubey et al. 2016). Naturally evolved biological processes are generally efficient but sometimes their behaviours challenge our intuition. For example, an increase in RuBisCO content in rice leaf does not necessarily contribute to enhanced photosynthesis, rather ribulose-1,5-bisphosphate carboxylase/oxygenase (RubiSco) can be used as a nitrogen storage that is mobilized as nitrogen source for the seeds during grain filling (Murchie et al. 2002).
In the context of improving metabolism in plants, previous work on the application of multi-tissue genome-scale metabolic model (GSM) coupled with dynamic flux balance analysis (FBA) demonstrated the importance of using modelling to predict metabolic reprogramming during plant development (Shaw and Cheung 2018) that can help to design better biotechnological strategies. However, one limitation of flux balance models is the difficulty in directly incorporating climate-related factors such as temperature, humidity, water availability etc. Climate change can directly or indirectly change local environmental growth conditions of crop plants (Mall et al. 2017), with one important factor being drought (Dai 2011). Hence, it is important to integrate climate-related factors in studying metabolic adaptations to climate change. A recent modelling study managed to capture the effects of temperature and humidity by coupling a 24-time step diel flux balance model with a biophysical gas-exchange model (Töpfer et al. 2020). In this study, we used the approach of coupling a flux balance metabolic model with a mechanistic crop growth that takes environmental factors as inputs.
Equation-based mechanistic models of plant growth that can calculate the carbon assimilation rates of canopy structure and the production of organ-specific biomass over the timescale of a few months are currently available (Van Diepen et al. 1989;Humphries and Long 1995). However, these models ignore the detailed metabolic dynamics as the primary objective of crop growth models is to predict the product yield driven by certain environmental factors. These factors include CO 2 and O 2 concentrations, soil condition, time of planting, temperature, humidity, water availability and altitude etc., which are input parameters that determine the outputs of crop growth model, many of which are important for agronomy, e.g. carbon assimilation and biomass yield. On the other hand, FBA-based modelling framework focuses on internal cellular behaviour by predicting the activities of metabolic processes such as nutrient and energy conversion mechanisms. Given the different strengths of the two modelling frameworks, the integration of crop growth models and metabolic models allows us to study dynamic metabolic changes during plant growth under different environmental and physiological conditions. For example, groundwater is an important factor that can impact potential yield (Boling et al. 2007), which is in decline in several regions around the world for use in irrigation (Rodell et al. 2009;Famiglietti 2014;Wada et al. 2014). Here we integrated the two different kinds of plant models to simulate rice growth under normal and limited groundwater conditions based on a day-wise FBA framework. This enables us to predict organ-specific metabolic fluxes dynamically during different stages of plant growth and to analyse carbon-nitrogen allocation and partitioning across different organs. In a previous study, a multiscale integration approach was used by integrating static organ models with dynamic whole-plant growth model for barley (Grafahrend-Belau et al. 2013). Here, we proposed a similar approach by integrating an FBA framework with a crop growth model with customizable primary input options including sowing time, soil and incorporation of stress. As a case study we focused on analysing the impact of water-limited stress on rice metabolism over multiple growth stages including grain filling. The carbon-nitrogen remobilization during the grain-filling stage determines the weight and quality of rice grain (Zhang et al. 2016). Analysis of metabolism over different growth stages of whole-plant development and its comparison between control and stress conditions allows us to further understand plant metabolic plasticity and adaptations to a changing climate, and potentially identify targets for improvement.

Rice GSM
A published rice (O. sativa) GSM  was manually curated to a complete mass-and charge-balanced model (Os2384; see Supporting Information-File S1). Apart from chloroplast (plastid in the present version), mitochondrion and cytosol in the previous version of the model, peroxisome and vacuole were also included in the updated version of the rice GSM. As we adapted the model to a diel multi-organ GSM for day-wise plant growth study, several organspecific biomass-producing reactions were introduced in Os2384, especially for seed metabolism (see next sections). Moreover, reactions for degradation of amino acids, chlorophyll (Chl), and sugars were also included in the upgraded model for modelling senescence. Os2384 includes 2049 reactions, including 18 environmental transport exchange reactions and 54 biomass export reactions, 1799 metabolites (in comparison to 1736 reactions and 1484 metabolites in the previous version) and 2384 genes. The new metabolic reactions and transporters added were based on the information from previous publications (Cheung et al. 2013;Chatterjee and Kundu 2015;Moreira et al. 2019;Shaw and Cheung 2019), literature evidence described in the following sections and metabolic databases including RiceCyc (archive.gramene.org/pathway/ricecyc.html) and PlantCyc (pmn. plantcyc.org) to enable the required metabolic functionalities needed for this study, including biomass degradation and seed biomass synthesis. Inter-compartmental (inter-organelle) transport reactions were renamed with appended suffixes '_pc', '_mc', '_xc' and '_vc' to represent exchange to/from cytosol from/to plastid, mitochondrion, peroxisome and vacuole, respectively. Also, along with maintenance cost consideration through an ATP-consuming reaction (ATPase in the cytosol), NADPH oxidation (NADPHOxid) was considered in four different compartments (cytosol, mitochondrion, plastid and peroxisome) (Cheung et al. 2013;Moreira et al. 2019;Shaw and Cheung 2019). The final GSM have 851 reactions updated, mostly to correct for charge balance. Our model's results from MEMOTE (Lieven et al. 2020), a standardized testing tool for GSMs, showed high scores (>96 %) for stoichiometry consistency, mass balance, charge balance and metabolite connectivity [see Supporting Information-File S2], with identified errors related to external metabolites. To validate the model functionalities after the model upgrade, phototrophic and heterotrophic biomass production and energetic consistency were re-validated using the same approach as our previous GSMs (Moreira et al. 2019;Shaw and Cheung 2019). In brief, a unit flux to all biomass components was fixed for production and photon influx, mineral nutrient uptake and gas exchanges were allowed with all other organic energy sources set to zero to validate phototrophic condition. To validate heterotrophic condition, we used the same constraints except that the photon flux was blocked and glucose was allowed as the energy (and carbon) source. The model was tested for energy conservation, and it was shown that the model cannot produce energy (ATP), reducing power (NADPH) and biomass without any external input source. The model was also tested to make sure the absence of transhydrogenase cycles that transfer reducing power from NADH to NADPH without any external energy source.

Construction of multi-organ flux balance model
Os2384 was used to build a three-organ diel model of rice including leaf, stem and seed [see Supporting Information-File S3]. The reconstruction process of the multi-organ model was similar to our previous studies (Shaw and Cheung 2018;Moreira et al. 2019;Shaw et al. 2019). Reaction and metabolite ids were duplicated with suffixes '_Leaf ', '_Stem' and '_Seed' to represent leaf-, stem-and seedspecific individual organ models, respectively, under '_Light' (day) and '_Dark' (night) phases (Fig. 1). The organ-specific models were connected by a common pool (CP) for exchangeable metabolites, representing the phloem and xylem. The light and dark phases were connected through the storage reactions (Shaw and Cheung 2018). In this study, the leaf model represents leaf blade, and the stem model represents straw fraction of culm and leaf sheath. Leaf sheath can be considered as a non-photosynthetic tissue as its photosynthetic capacity is very small (Yoshida 1981). Thus, photon uptake was only allowed in the leaf. Uptake of nitrate (NO 3 − ), ammonia (NH 3 ), phosphate (P i ), sulfate (SO 4 2− ) and magnesium (Mg 2+ ) was allowed from CP that can be distributed to leaf, stem and seed by their respective organ transport reactions ('_transport') ( Fig. 1). Figure 1. Schematic description of rice multi-organ GSM. Internal biosynthetic and degradation-specific processes are shown with solid and dashed arrows, respectively. The rice multi-organ GSM contains leaf-, stem-and seed-specific compartments connected by a CP. Amino acid (AA) and sucrose (Suc) transport was allowed between the three organs through CP. 'BiossSyn' and 'BiomassDeg' indicate single biomass reactions for the synthesis and degradation of organ-specific biomass, respectively, in the respective organs.
Twenty different amino acids and sucrose were allowed to transport between leaf, stem and seed through CP (Fig. 1). All organ-specific models can freely exchange CO 2 , O 2 , water and proton (H + ) to/ from the environment. Metabolite transport costs between the CP and organ-specific models were accounted for in our model by considering the energetic requirement for establishing the proton gradient required for metabolic transport. Each organ-specific model has its own biomass-synthesizing reaction ('_biomass_tx') in which the stoichiometry represents the proportion of biomass components specific to the corresponding organ [see Supporting Information- Table  S1]. Moreover, leaf and stem models also include reactions for biomass degradation ('_deg_tx'). We have also included reactions for the accumulation/export of individual biomass components ('_biomass') in all three organs to enable its ease of reuse in future studies. However, fluxes for all these individual biomass transporters were fixed to zero during the entire simulations (except for NCC, phytol, Mg 2+ and sulfite (SO 3 2− ) which were allowed to accumulate during leaf degradation). The overall set-up of the model allowed photosynthesis to occur in the leaf, storage of assimilates in the leaf and stem and remobilization of reserves in the seed during senescence for grain filling.

Biomass equations of leaf, stem and seed
The stoichiometry of components in the single biomass-synthesizing reactions for leaf, stem and seed was obtained from different sources in experimentally observed proportions for different rice organs. Leaf and stem biomass stoichiometries were taken from Poolman et al. (2013) and Lakshmanan et al. (2013), respectively. Other than the biomass components used by Poolman et al. (2013), palmitate, fructose and precursors for xylan synthesis were included for the stem fraction. Biomass stoichiometry for seed was obtained from USDA National Nutrient Database for Standard Reference (ndb.nal. usda.gov/ndb). Measurements suggest that husk is about 20 % of the weight of the whole grain and mainly formed by lignocellulosic biomass (Singh 2018). Thus, we included this estimate to represent the total whole grain or seed biomass stoichiometries [see Supporting Information- Table S1]. Seed biomass reactions include the components for protein, vitamin, carbohydrate, water-soluble (β-glucan) and water-insoluble fibres (lignin) and fatty acids. Stoichiometries of biomass components in all organs were normalized to represent 1 g of biomass production for one unit of flux ('Used stoichiometry' in Supporting Information- Table S1). Stoichiometries of biomass degrading reactions in leaf and stem were also normalized to represent 1 g of degrading products for one unit of flux. Detail of the biomass stoichiometries is given in Supporting Information-Table S1.  (Fig. 2). WOrld FOod STudies provides a stand-alone framework to obtain crop production potential (i.e. changes in biomass of different organs) under different environmental conditions (such as water availability), location, soil composition and sowing dates. Some of the major applications of WOFOST include analysing 'yield risk and inter-annual yield variability', 'yield variability over soil types, or over a range of agro-hydrological conditions' , 'relative importance of growth determining factors', 'effects of climate change' and 'regional assessments of crop yield potential' (www.wur. nl/en/show/A-gentle-introduction-to-WOFOST.htm). WOrld FOod STudies simulates plant growth based on eco-physiological processes where different model components, such as light interception, assimilation rates, transpiration rates, partitioning etc., are connected with modules of weather and soil, which works together to provide outputs on organ formation and growth potential with a fixed time step of 1 day (de Wit et al. 2019).

Integration of crop growth model and constraint-based metabolic model
The rice multi-organ GSM provided a platform for integrating these day-specific biomass constraints for each organ from WOFOST (Fig. 2). In the constraint-based FBA framework, the metabolism of each day was modelled by fixing day-wise organ-specific growth rates from WOFOST and solving an FBA problem for each day. The parameters for generating the WOFOST predictions, which were subsequently used as constraints for the FBA simulation, can be found in Supporting Information- Table S2. In this study, two different simulations were conducted-normal plant growth (normal) and growth under water-limited condition (stress) (referred to as 'water-limited crop production without groundwater' in WOFOST). Groundwater is an important source of global freshwater that can vary due to the effect of climate change. The rising world population and the excessive use of groundwater in agriculture make it an important factor to be considered in modelling approaches to predict metabolic adaptations.
Leaves are the main source of remobilized N in panicles (Mae and Ohira 1981;Mae 1997) and about 70-90 % of panicle N comes from the remobilized source (Mae and Ohira 1981). In Japonica cultivar carbon remobilization provides about 14 % of grain carbohydrate (Yang et al. 2001). Here, both protein and sugar degradation were allowed from the leaf, whereas the stem was allowed to degrade only sugars [see Supporting Information- Table S1]. The end of the vegetative growth period was considered as the time when the seed starts to grow (start of the reproductive phase), the reproductive phase ends when remobilization of reserves begins (grain filling starts, leaf (and stem) started to degrade for 2 days in a row) and grain filling continues until the end of the simulation (Fig. 3). The ratio of biomass production (growth) between day and night phases was set to 70 %:30 % (estimated from net daytime growth and net night-time growth for Arabidopsis thaliana under 12-h photoperiod; Sulpice et al. 2014). The days in our simulations were assumed to be 12-h photoperiods. Given that the biomass outputs were fixed in our model simulations, photon uptake and CO 2 exchanges were kept as unconstrained during the entire simulations.
We used two types of biomass reactions in the leaf, one for the synthesis of biomass components ('biomass_tx' , LeafBiomassSyn) and another for C and N remobilization ('biomass_deg_tx' , LeafBiomassDeg) from stored proteins (and Chl; Fig. 1) and starch. As for stem, a reaction named StemBiomassSyn was used to synthesize stem biomass components, whereas sugar remobilization was modelled with the reaction named StemBiomassDeg. Leaf-and stem-specific absolute growth rates (AGRs; see Supporting Information- Table S3) (positive growth rates, obtained from WOFOST) integration to the respective BiomassSyn reactions was implemented until the leaf and stem started to degrade (negative growth rate). Beyond these points, the biomass degradation rates were used to set constraints for the respective BiomassDeg reactions (and BiomassSyn reactions were set to zero). Thus, we allowed protein and carbohydrate (glucose, sucrose and starch) remobilization from leaf and stem during grain filling. Seed AGRs from WOFOST were fixed to the seed biomass synthesis reaction (SeedBiomassSyn). All rates from WOFOST were scaled to a single plant per gram dry weight (DW) biomass considering an optimal density of ~250 000 plants per hectare (~25 plants per square metre; Salahuddin et al. 2009). The objective function for the FBA model was set to minimizing total absolute sum of fluxes reflecting efficient use of reactions (enzymatic machinery). This objective function also eliminates substrate cycles. Chlorophyll degradation pathway was included in the metabolic model based on known enzymatic activities involved in the degradation (Hörtensteiner 2006). Given the final destination of Chl degradation is the vacuole where Chl can be degraded to non-fluorescent Chl catabolites (NCCs) during senescence and not export to the rest of the plant as the N source (Hörtensteiner and Feller 2002), we only considered Chl to degrade and ultimately accumulate in vacuole as NCCs using a NCC transport reaction. On the other hand, RuBisCO serves as an important N remobilization source (Murchie et al. 2002) from senescing leaves along with other chloroplast proteins (Hörtensteiner and Feller 2002) during grain filling. Full elucidation of the process of RuBisCO and other chloroplast protein proteolysis is still underway (Hörtensteiner and Feller 2002;Feller et al. 2007;Ono et al. 2013). Thus, we assumed no interaction between proteolysis and other chloroplastic/cytosolic metabolic processes; therefore, net protein degradation was available directly as amino acids in the cytosol. Leaf biomass degradation included the degradation of Chl a and b. While the carbon from the degradation of Chl a and Chl b can feed into the rest of the metabolic system, N from Chl a and Chl b can only accumulate in the vacuole and cannot be used as N source for other metabolic processes. Amino acids and sucrose from leaf and stem degradation can be metabolized and imported into the seed. Thus, from the time leaf and stem started degrading, they can become a 'source' while the seed remained as the major 'sink' for the released metabolites. Therefore, the C, N and energy sources for seed development can include photosynthates directly from photosynthesis and the remobilized carbon and nitrogen products from leaf and stem. We used similar maintenance requirements for leaf and stem (7.27 mmol g DW −1 day −1 and 2.56 mmol g DW −1 day −1 through ATP and NADPH maintenance reactions at a ratio of 3:1; Cheung et al. 2013) according to their biomass (Shaw and Cheung 2018), and similarly seed maintenance was set to 0.1 mmol g DW −1 day −1 and 0.033 mmol g DW −1 day −1 for ATP and NADPH maintenance, respectively. Seed maintenance was estimated by temporally calculating terminal sucrose transport to the seed without any maintenance and a ~12 % raise in terminal sucrose demand in seed after addition of our estimated maintenance that match the experimental estimation that rice inflorescence with seed requires ~12 % additional sugars for maintenance (see Table 4 of De Vries et al. (1983)).

Software used for FBA simulations
All simulations were done in scobra (github.com/mauriceccy/scobra), an extension package for cobrapy (0.4.1) (Ebrahim et al. 2013), using the GNU Linear Programming Kit (GLPK) (www.gnu.org/software/glpk/) solver on Python 2.7 (www.python.org). The simulated flux distributions for both normal and water-limited conditions are provided in Supporting Information- Table S4.

Updated rice GSM and modelling framework
In this study, we present an updated rice GSM (Os2384) which was manually curated for mass-and charge-balancing with 2049 reactions, 1799 metabolites and 2384 gene associations. We paid special attention in the manual curation of metabolic reactions involved in multiple degradation pathways, in particular the Chl degradation pathway, to allow our genome-scale model to simulate metabolic processes involved in senescence.
A multi-organ version of the rice GSM was constructed, comprising leaf, stem and seed. To the best of our knowledge, this is the first rice multi-organ GSM with a seed model. The seed model has its own biomass equation and we have included an estimation of the cellular maintenance cost of the seed in this study (see Materials and Methods). While it is possible to include a root model in multiorgan models, as demonstrated in previous studies (Grafahrend-Belau et al. 2013;Shaw and Cheung 2018), a root model was not included in this study because of the lack of root growth rate prediction from the WOFOST simulations for rice growth. Instead, nutrients including nitrate (NO 3 − ), ammonia (NH 3 ), phosphate (P i ), sulfate (SO 4 2− ) and magnesium (Mg 2+ ) were available directly from the CP.
The day-wise FBA framework used in this study allowed us to investigate the metabolic changes involved in rice growth from seedling to seed development. Similar to the dynamic FBA approach previously used in studying the growth of barley (Grafahrend-Belau et al. 2013), the integration of the two models at different scales was mostly done through the setting of biomass constraints. The main distinctive feature of our study is the introduction of a general integration approach with a public stand-alone tool, WOFOST (Van Diepen et al. 1989;de Wit et al. 2019), which is a crop growth model that has the capability of incorporating weather data-a proxy for climate-environmental variables. In principle, the model integration approach can be applied with any models that can provide biomass constraints for the flux balance model. An additional distinctive feature of our model is the ability to model leaf and stem degradation (decrease in organ biomass over time) to analyse plant senescence process, allowing us to study the resource allocation during seed filling, which could depend on both newly assimilated and reallocated C and N.

Growth patterns of rice leaf, stem and seeds
The changes in biomass of leaf, stem and seed over the growth period under normal and water-limited (stress) conditions are shown in Fig. 3. The overall growth period can be divided into three stages, vegetative (seedling and up to seed starts to emerge), reproductive (seed development by fresh photosynthates) and grain filling (seed development by fresh and remobilized sources). Under both normal and stress conditions, plant growth goes through all three stages. However, water-limited stress had lower growth rates and shorter growth period (Fig. 4) which results in lower overall biomass than normal condition (Fig. 3).
Under water-limited condition, continuous lower growth rates during the vegetative phase (Fig. 4) led to lower leaf biomass (Fig. 3) that in turn reduces photon uptake capacity (Fig. 5A) and N and C assimilation ( Fig. 5B and C), which resulted in a lower amount of reserve for grain filling (Fig. 3).
Photon utilization for plants under water-limited stress is comparable to that under normal growth during the vegetative stage. However, the effect of water limitation becomes apparent after the vegetative stage where leaf growth (Fig. 3), photon utilization (Fig. 5A) and carbon assimilation (Fig. 5C) flatten off under water-limited stress, in contrast to the continued leaf growth and increase in photon uptake under normal condition beyond the vegetative stage (Fig. 5A). Another key difference between normal and water-limited conditions is the timing of beginning of seed growth. Under normal condition, significant amount of seed filling occurs relative early (50-60 days after sowing (DAS)), when there is still positive leaf growth. However, under waterlimited condition, significant seed growth only occurs when the leaf Figure 4. Growth rates under normal and water-limited conditions plotted as 10-day moving averages. Positive values represent growth; negative values represent degradation. Leaf, Stem and Seed represent leaf, stem and seed organ-specific growth and degradation (_deg), respectively. Separate plots for growth and degradation rates are provided in Supporting  Information-Fig. S1. Vertical lines are the same end of stage labels as Fig. 3. and stem started to senesce (after 80 DAS). This means the carbon source for seed filling under water-limited condition is primarily from the senescence products of leaf and stem, whereas a significant amount of carbon from photosynthesis can be used for seed growth in normal condition. The differences in the pattern of resource allocation and metabolism between the two conditions were explored by looking at the predictions from the multi-organ GSM.

Investment in stem growth supports final seed biomass under water-limited stress.
Plant metabolism is known to change as different tissues and organs go through their developmental stages (Grafahrend-Belau et al. 2013). Growth under normal condition follows a sigmoidal growth pattern where seeds receive fresh photoassimilates in the reproductive stage and sufficient reallocated compounds from leaf (and stem) degradation later during the grain filling (Fig. 3). Among the changes in leaf, stem and seed biomass under water-limited condition, reduction in total seed yield was more severely affected than leaf and stem (Fig. 3). Less N and C assimilation ( Fig. 5B and C), transport and storage throughout the water-limited stressed growth period resulted in reduced grain-filling potential (Figs 5D and 6B) which led to reduced seed biomass. Data from the WOFOST model suggest that grain filling under water-limited condition starts earlier (Day 50) than the normal condition (Day 65) which is a result of early onset of senescence that lasts longer following a short reproductive stage (Fig. 4). A sharp spike in seed growth under water-limited condition arises when stem degradation and its contribution to grain filling starts (Day 85; Fig. 4).
Water-limited condition results in a drop in N and C assimilation after a certain growth period ( Fig. 5B and C) that can influence Chl content (Mafakheri et al. 2010), and consequently photon capture capacity (Fig.  5A). However, the spike in seed growth under water-limited condition was coupled with an increase in sucrose transport (Fig. 5D) along with leaf and stem degradation (Fig. 4) to support the majority of the seed biomass production (Fig. 3), suggesting a shift in resource allocation to support seed growth.

Amino acids as nutrient source for the seed.
Under water-limited conditions, grain filling at the later stage of seed development was predicted to require significant amount of fresh N assimilation, while only minimal amount was required under normal condition (Fig. 5B). Based on the flux balance model, arginine (ARG), glycine (GLY) and serine (SER) were predicted to be the preferred source of transported amino acids to the seed under both normal and stress conditions ( Fig. 6A and B). Interestingly, these three amino acids are not most abundant amino acids in the seed biomass (glutamate being the highest followed by aspartate and leucine). In developing soybean cotyledons, ARG was found to be one of the major constituents of the total free amino acid pool that is important for N nutrition (Micallef and Shelp 1989). Our model predicted that GLY, ARG and SER are the most efficient sources of nitrogen to the seed for the synthesis of other nitrogen-containing compounds with respect to the applied objective function of minimizing total flux.

Water-limited stress-specific metabolic responses.
As sessile organisms, plants are evolutionary well-adapted to respond metabolically under stress conditions. While many reactions in central metabolism are constitutively active for normal growth, there are other reactions which are only active in response to different stimuli (Wink 2006). Moreover, many key reactions in central carbon metabolism can also be compromised by reduction in flux under stressed conditions that can result in stress phenotypes such as stunted growth. To identify reactions that are differentially regulated between normal and waterlimited stress conditions, we looked at the reactions that were predicted to carry a non-zero flux (higher than 10 −6 mol per phase (12 h) to eliminate numerical errors) at different developmental stages under normal and stress conditions ( Fig. 7; see Supporting Information- Table  S5). Figure 7A shows that 30 and 8 reactions are unique to the vegetative and reproductive stages, respectively, whereas a comparatively large number of reactions (182) are unique to the grain-filling stage under normal condition. These 182 reactions are required to activate metabolism specific for grain filling under normal condition, whereas under water-limited condition the number is even higher (i.e. 407; Fig. 7B). Moreover, it is not surprising to see that there are very small number of reactions that are unique to the reproductive stage, the stage connecting between vegetative and grain filling, suggesting that gradual metabolic changes from vegetative stage to reproductive stage and from reproductive stage to grain-filling stage. The majority of the unique reactions under water-limited condition were found to be during the grainfilling stage (Fig. 7E), suggesting that there could be a large number of reactions that are activated in response to water-limited stress during the grain-filling stage. By inspecting the reactions that were unique to water-limited condition in the grain-filling stage [see Supporting Information- Table S5], we found many reactions involved in amino acid metabolism. It is worth noting that most of the unique reactions during the grain-filling stage under water-limited stress are in the leaf, suggesting that leaf metabolism is most impacted by water limitation during the grain-filling stage compared to stem and seed.
Our results suggest that the grain-filling stage is particularly important under water-limited stress. Under water-limited stress, there was a short reproductive stage and a long grain-filling stage (Fig. 3), which suggests that there could be major metabolic contributions during the grain-filling stage for seed development. This is also reflected in the larger number of unique reactions during the grain-filling stage under water-limited condition (407 from Fig. 7B) compared to normal condition (182 from Fig. 7A).
Among the three stages, the vegetative stage has the least amount of metabolic changes between normal and water-limited conditions ( Fig. 7C) compared to the reproductive stage and grain-filling stage ( Fig. 7D and E). This is consistent with the growth patterns (Fig. 3) where the least difference between normal and water-limited conditions was obtained from the WOFOST model. The majority of the 16 reactions unique to water-limited stress in the vegetative stage were reactions in the stem, with some reactions related to glutathione (glutathione synthase, S-formylglutathione hydrolase, thioredoxindisulfide reductase, etc.) and one-carbon compounds. It has been known that one-carbon metabolism is involved in the regulation of plant defence response (González and Vera 2019) and has a role to cope with oxidative stress (Gorelova et al. 2017). Glutathione metabolism is also known to play a role in plant defence system under environmental stress (Rennenberg and Brunold 1994) and is involved in the sulfur nutrition pathway. Our modelling results suggest that glutathione metabolism in the stem could be involved in water-limited stress during the vegetative stage.

A modelling framework integrating crop growth and metabolism
Rice, being a major staple food in the world, has been used to study different aspects of crop plants in multiple research directions (Paine et al. 2005;Jeong et al. 2010;Poolman et al. 2013). In the field of computational systems biology, genome-scale metabolic modelling has been an important tool for studying metabolic robustness in a holistic approach (Oberhardt et al. 2009;Poolman et al. 2009Poolman et al. , 2013Shaw and Cheung 2018). A mass-and charge-balanced GSM with gene-reaction associations can play an indispensable role in simulating the dynamic metabolic plasticity under perturbed conditions as well as normal metabolic throughput (Moreira et al. 2019;Shaw and Cheung 2019).
In this study, we curated an existing rice GSM  for mass-and charge-balancing, corrected for compartmentalization of metabolic reactions and added gene-reaction associations to the model. We also included the seed biomass and curated its biosynthetic pathways which will be an important resource for future metabolic studies on rice grain yield (Usui et al. 2016;Zhang et al. 2019). Moreover, the inclusion of degradation pathways for C and N sources in the model can help the study of resource reallocation during grain filling and plant metabolic functions in general. This study has presented a multi-organ diel version of the GSM of rice [see Supporting Information-File S3] that will allow future metabolic studies at a whole-plant level (Shaw and Cheung 2020). Figure 8 illustrates a potential use case of the multi-organ model in studying the pattern of metabolic interactions between the organs and with the environment at different developmental stages. A better understanding of resource partitioning and senescence can enhance productivity (Kumar et al. 2019). We can use our multi-organ model to study source-sink relationship, e.g. our model predicted the pattern of resource allocation from leaf to stem and seed, which were qualitatively similar between normal and water-limited conditions (Fig. 8).
In a changing climate, one of our priorities in the field of plant metabolic modelling is to study the metabolic adaptions of major crop plants to various climate-related factors at a whole-plant scale (Rotter and Van de Geijn 1999). One of the challenges has been the difficulty in directly relating climate-related parameters, such as groundwater potential, altitude, etc. that are used in crop growth models, to be used as constraints in a GSM. The modelling framework presented in this study illustrates how we can integrate the two modelling approaches to study metabolic adaptions to changes in climate. WOrld FOod STudies supplies the constraints based on environmental conditions and the constraint-based GSM uses these constraints to simulate metabolic phenotypes. This approach can be generalized to any species with available crop growth models and constraint-based metabolic models.

Integrated modelling framework allowed predictions of metabolic differences between normal and water-limited conditions
Plant growth and development depends on many environmental factors. Water-limited condition is known to have significant impact on plant growth and development that result in altered metabolic functions ( Jaleel et al. 2009). Drought can affect different developmental stages and can lead to large-scale changes in the expression of metabolic genes to cope with water-limited condition (Oladosu et al. 2019).
In this study, we were able to predict the dynamic metabolic changes over multiple stages of rice growth under normal and water-limited conditions, which showed a clear difference in metabolism between the two conditions. Our modelling results showed a nearly flat photon utilization after the vegetative growth under water-limited condition (Fig. 5A). Early leaf senescence (Fig. 3) could be one of the main reasons of stunted plant growth under water-limited condition as it led to reduced photosynthetic rate compared to normal condition ( Fig. 5A and C). Given that the only differences between the simulations of normal and water-limited conditions were the differences in biomass  constraints of each of the three organs to the flux balance model, our results suggest that the biomass constraints alone were sufficient to capture the differences in metabolic phenotypes of plants, e.g. patterns of resource allocation, between normal and stress conditions. Growth under water-limited stress started in a similar way as growth under normal condition up to the vegetative stage (Figs 3 and 4). Significant metabolic differences were observed between the two conditions only after the vegetative stage (Figs 5 and 7). Our model predictions suggest that significant nitrogen assimilation is necessary for grain filling at the late stage of seed development under water-limited condition but not under normal condition (Fig. 5B). Moreover, there were more grain filing-specific reactions under water-limited condition ( Fig. 7B and E), suggesting that there could be a specific set of metabolic adaptions during seed development in response to limited water availability.

CON CLUSION
This study provided an updated rice GSM and a multi-organ diel version of the model. We demonstrated the use of a day-wise FBA framework to integrate our multi-organ metabolic model with a mechanistic plant growth model, which allowed us to predict the metabolism of multiple rice organs over multiple stages of development under different environmental conditions. Our results showed clear metabolic differences between normal and water-limited conditions and revealed the patterns of nutrient assimilation and resource allocation between different organs. The modelling framework presented in this study can be applied to model the metabolism and growth of other crop plants under different climatic conditions, which could aid future research on metabolic analysis of crop growth and organ development under different climatic conditions.

SUPPORTIN G INFOR M ATION
The following additional information is available in the online version of this article- Figure S1. Growth and degradation rates under normal and water-limited conditions. Separate plots for growth and degradation are shown for both normal and water-limited conditions. File S1. Genome-scale metabolic model of rice (Oryza sativa). The model includes curated degradation pathways and gene-reaction associations. File S2. Genome-scale metabolic model (GSM) overview report. Metabolic model comprehensive overview based on scores computed by MEMOTE. File S3. Multi-organ genome-scale metabolic model of rice (Oryza sativa). The model includes leaf, stem and seed along with their respective biomass synthesis reactions. Table S1. Stoichiometry components of leaf, stem and seed biomass. Scaling and normalization procedures are also included. Table S2. Rice-specific soil and weather parameters used to obtain WOrld FOod STudies (WOFOST) output in the form of biomass yields under normal and water-limited conditions. Table S3. Leaf, stem and seed growth rates from WOrld FOod STudies (WOFOST). Table S4. Metabolic flux solutions from the multi-organ model. Table S5. Raw files of the Venn diagrams shown in Fig. 7.