Extended Utilization of Constraint-Based Metabolic Model in a Long-Growing Crop

The constraint-based rMeCBM-KU50 model of cassava storage root growth was analyzed to evaluate its sensitivity, with respect to reaction flux distribution and storage root growth rate, to changes in model inputted data and constraints, including sucrose uptake rate-related data—photosynthetic rate, total leaf area, total photosynthetic rate, storage root dry weight, and biomass function-related data. These mainly varied within ±90% of the model default values, although exceptions were made for the carbohydrate (−90% to 8%) and starch (−90% to 9%) contents. The results indicated that the predicted storage root growth rate was highly affected by specific sucrose uptake rates through the total photosynthetic rate and storage root dry weight variations; whereas the carbon flux distribution, direction and partitioning inclusive, was more sensitive to the variation in biomass content, particularly the carbohydrate content. This study showed that the specific sucrose uptake rate based on the total photosynthetic rate, storage root dry weight, and carbohydrate content were critical to the constraint-based metabolic modeling and deepened our understanding of the input–output relationship—specifically regarding the rMeCBM-KU50 model—providing a valuable platform for the modeling of plant metabolic systems, especially long-growing crops.


Introduction
The complexity of plant metabolism hinders experimental studies to unravel the fate of metabolic substrates derivation for biomass production, and it is even harder to unwind the underlying regulations. As multicellular organisms, plants are composed of cells with several subcellular compartments, defined as organelles; as a result, metabolic processes are fragmented into parts that are exposed to heterogeneous surrounding environments [1]. A constraint-based modeling (CBM) approach with flux balance analysis (FBA) was employed to deepen the understanding of plant metabolism through simulation. The CBM approach enables the study of biological systems and characterization of metabolic network behavior by flux analysis under steady-state assumptions, integrating biochemical, genetic, and genomic information [2]. Constraints and the biomass objective function, which are the heart of the CBM approach for identifying biologically relevant flux solutions, were established based on organism-specific and/or condition-specific experimental measurements [3][4][5][6][7]. In case of

The rMeCBM-KU50 Model
The rMeCBM-KU50 model is a constraint-based metabolic model of carbon assimilation in cassava storage roots. It was originally developed to simulate metabolism during starch accumulation in the storage roots of Kasetsart 50 (KU50), a high starch-yielding cultivar that is widely grown in Southeast Asia [4]. The rMeCBM-KU50 model is composed of 468 reactions related to 393 metabolites, which are partitioned into 3 subcellular compartments, i.e., the cytosol, mitochondria, and plastid. The details of the rMeCBM-KU50 model (i.e., reactions, metabolites, and biomass function) can be found in the supplementary information of the study by Chiewchankaset et al. [4]. The rMeCBM-KU50 model successfully mimicked the measured storage root growth rate of KU50 cultivar, using the pre-determined specific sucrose uptake rate of 0.0548 mmol sucrose gDW −1 SRs day −1 , calculated based on the photosynthetic capacity (Equation (1)), and biomass components (92.07% carbohydrates, 5.78% proteins, 1.69% fibers, and 0.45% lipid, on dry weight basis). The optimal growth-associated ATP maintenance (S GAM ) of the rMeCBM-KU50 model was 9.8 mmol ATP gDW −1 SRs [4].
Speci f ic sucrose uptake rate mmol sucrose gDW −1 SRs day −1 = P N is photosynthetic rate (µmolCO 2 m −2 s −1 ), LA is total leaf area per plant (m 2 ), SR DW is the dry weight of cassava storage roots (g), F (= 0.37) is a faction of photosynthetic sucrose translocated to the root of the cassava-estimated based upon the proportion of sucrose to the total soluble carbohydrates in shoot and root parts on the dry weight basis [4]-and 7.2 is the constant value for unit conversion.

Model Sensitivity Analysis
The sensitivity of the inputted data for the rMeCBM-KU50 model was investigated through a simple perturbation and response method, similarly to other various model analysis studies, including constraint-based metabolic models [26,[28][29][30][31][32]. The individual inputted data were systematically varied throughout their range of feasible values, while observing the alteration of the simulation from the default prediction. The impacts of the experimental measurements and constraints on the model simulation were assessed based upon the simulated storage root growth rates and the predicted carbon fluxes, as demonstrated in Figure 1. Two groups of model data were investigated, consisting of (i) substrate-uptake related data (the specific sucrose uptake rate, photosynthetic rate, total leaf area, storage root dry weight, and total photosynthetic rate), and (ii) biomass function-related data (the composition of biomass components), as shown in Figure 1a,b, respectively. The values of these model data were thoroughly varied within ±90% the range of the original values [4], to cover a full range excluding carbohydrate and starch contents. The carbohydrate and starch contents, respectively, were varied in the ranges of −90% to 8% and −90% to 9% to ensure positive values (see details in Supplementary Materials File S1).

Figure 1.
Overall framework of rMeCBM-KU50 model assessment to changes of (a) substrate-uptake related data (i.e., the specific sucrose uptake rate, photosynthetic rate, total leaf area, storage root dry weight, and total photosynthetic rate) and (b) the biomass function-related data (i.e., composition of biomass components). Each data variation, the sensitivity analysis of model prediction examined the changes in the predicted growth rate, as well as the changes of metabolic flux direction and partitioning. Inactive fluxes refer to zero-flux reactions (gray arrows). Active fluxes refer to non-zero flux reactions (black arrows-with unchanged direction; and red arrows-with changed direction). The thickness of the arrows indicates the flux magnitude. Overall framework of rMeCBM-KU50 model assessment to changes of (a) substrate-uptake related data (i.e., the specific sucrose uptake rate, photosynthetic rate, total leaf area, storage root dry weight, and total photosynthetic rate) and (b) the biomass function-related data (i.e., composition of biomass components). Each data variation, the sensitivity analysis of model prediction examined the changes in the predicted growth rate, as well as the changes of metabolic flux direction and partitioning. Inactive fluxes refer to zero-flux reactions (gray arrows). Active fluxes refer to non-zero flux reactions (black arrows-with unchanged direction; and red arrows-with changed direction). The thickness of the arrows indicates the flux magnitude. For each group of data, the sensitivity on model prediction was assessed based upon the simulated storage root growth rate ( Figure 1, left panel) and the predicted metabolic fluxes (Figure 1, right panel). The sensitivity of the simulated growth rate (Figure 1, left panel) was determined by the percentage of changes, as shown in Equation (2).
where ε is the percentage of changes of the simulated specific storage root growth rate (v v ) from the default prediction (v o ). The sensitivity of the predicted metabolic fluxes (Figure 1, right panel) was evaluated both upon the changes in direction of individual flux, and upon flux partitioning at each metabolic branch point where a metabolite was converted into a product by multiple reactions (production and consumption reaction fluxes) [35]. The percentage of production and consumption reaction flux partitioning was calculated and, respectively, denoted as positive and negative signs in this study. The details of the metabolic branch points of the rMeCBM-KU50 model are provided in the Supplementary Materials (File S2). The percentage of data variation, where the direction and/or partitioning of each predicted reaction flux began to deviate from the original prediction, was denoted here as a breakpoint. All simulations were performed using the COBRA Toolbox 2.0.5 [36] in MATLAB (The MathWorks, Inc., version R2015a, Natick, MA, USA).

Results
The sensitivity of a prediction to variances in model input data determines the capability of the model to simulate real systems. Hence, the quality of the experimental measurements and constraints used for the prediction is critical to the reliability of the predictions. Accordingly, the rMeCBM-KU50 model was evaluated according to the sensitivity of predicted results, based upon (i) substrate-uptake related data (i.e., the specific sucrose uptake rate, photosynthetic rate, total leaf area, storage root dry weight, and total photosynthetic rate), and (ii) biomass function-related data (i.e., the compositional ratios of biomass components). The values for the observed input data were systematically varied within a ±90% range of default experimentally determined optimal values for cassava storage root growth [4] to ascertain the model sensitivity, considering values that have been reported for cassava [37][38][39][40][41][42][43][44].

The Impacts of Substrate-Uptake Related Data on Model Prediction
The rMeCBM-KU50 model primarily utilizes sucrose as the main carbon substrate to simulate the growth of cassava storage roots. The specific sucrose uptake rate employed in the model was calculated from experimentally determined data on the photosynthetic rate, total leaf area, and storage root dry weight of cassava. Details of the range of values used for the sensitivity analysis, with respect to sucrose uptake related data, are shown in Table 1.
Results for specific sucrose uptake rates are shown in Figure 2. The photosynthetic rate, total leaf area, and total photosynthetic rate were positively related to the specific sucrose uptake rates, in contrast to storage root dry weight (Figure 2a). The specific sucrose uptake rate was mostly affected by the storage root dry weight-changing by approximately +850% and −50% of the original value at SR DW of −90% and +90%, respectively-followed by the total photosynthetic rate (Figure 2b). The photosynthetic rate and total leaf area both showed little effects.  [4]. b The range of data variation within ± 90% of the original value in the rMeCBM-KU50 model. c The total photosynthetic rate was the product of the total leaf area and photosynthetic rate. d The specific sucrose uptake rate was the product of the total leaf area and photosynthetic rate divided by storage root dry weight.     The sensitivity of the model-regarding the specific storage root growth rate prediction, to specific sucrose uptake rate-was analyzed based on changes to the photosynthetic rate, total leaf area, total photosynthetic rate, and storage root dry weight ( Figure 3). The experimentally derived values used as the baseline for evaluating the sensitivity of the model are marked as red crosses, and the gray areas on the left panel represent the range of possible values, based on the literature. The results showed a linear, positive response between the predicted specific storage root growth rate and the specific sucrose uptake rate-in relation to variations in photosynthetic rate, total leaf area, and total photosynthetic rate-but that for storage root dry weight was nonlinear (Figure 3a-d). Moreover, it was found that the sensitivity of the predicted growth rate to specific sucrose uptake rate in relation to variations in total photosynthetic rate was multiplicative to the effect of the individual dependent factors (i.e., the variations in photosynthetic rate and total leaf area). The percent error of predicting the specific storage root growth rate from the specific sucrose uptake rate by varying the photosynthetic rate and total leaf area to within ±90% of the original values ranged from +90% to −88.89%, while that for total photosynthetic rate ranged from +260% to −97.78%, and storage root dry weight ranged from −46.67% to +880% ( Figure A1a in Appendix A). Thus, the storage root dry weight and total photosynthetic rate were considered important data for determining the specific storage root growth rate.
The sensitivity study covered the entire range of expected values in cassava. The photosynthetic rate was varied from 13 to 35 µmolCO 2 m −2 s −1 to represent values that have been reported under different growth and environmental conditions. Cassava photosynthetic rates of 13-24 µmolCO 2 m −2 s −1 and 20-35 µmolCO 2 m −2 s −1 are typical under greenhouse or growth chamber conditions [39,40] and field conditions [41], respectively. The total leaf area of cassava is around 1.24-3.38 m 2 depending on variety, shade level, water availability, temperature, and plant age [38,42]. The average storage root dry weight of cassava, based on over 100 cultivars, under rain-fed field conditions, ranges from 0.01 to 1.40 kg plant −1 [43]. This range supported the representativeness of the rMeCBM-KU50 model for mimicking the metabolism of cassava storage roots, as derived by proper model data (marked as red crosses in Figure 3).
To analyze the sensitivity of carbon flux distribution to changes in the specific sucrose uptake rate, the direction and partitioning of the carbon flux were examined. Details of the flux direction and partitioning are provided in Figures A2-A5 and Figures A6-A8 in Appendix A, respectively (see the details in Supplementary Materials File S3-S4). The flux direction was identical for all the model variables analyzed, as indicated by the blue-colored active fluxes with unchanged direction. Also, analysis of the flux partitioning revealed the variable data considered were largely comparable, except for the total photosynthetic rate from −60% variation downward, when the flux partitioning became different from the original value (corresponding to a decrease in the specific sucrose uptake rate of 84%), herein referred to as the breakpoint (Figures 2b and A6 in Appendix A). These results indicated that the carbon flux prediction was relatively robust to changes in the specific sucrose uptake rate.

The Impacts of Biomass Function-Related Data on Model Prediction
The biomass function is typically integrated into the model as an objective function, which is used to draw resources from the metabolic network and define cellular components along with their functions [27]. The default biomass function in the rMeCBM-KU50 model was retrieved from literature and consisted of carbohydrates, proteins, fibers, and lipids. Cassava comprises mainly of starch, which could account for up to 90% of the storage root dry weight, depending on variety, growth condition, and age [45]. Therefore, we sought to determine if the starch content could solely be used to define the biomass characteristics of the cassava storage root model. The sensitivity of the model prediction was analyzed by varying the total carbohydrate, starch, protein, fiber, and lipid compositions ( Table 2). The prediction of the specific storage root growth rate was examined against the variation of biomass components represented in Figure 4. The default storage root growth rate in the rMeCBM-KU50 model, marked as red crosses, were used as the baseline for sensitivity analysis. The results showed that the predicted storage root growth rate was more robust to changes in biomass components (Figure 4a-e) with a percent error of 18.89% ( Figure A1b in Appendix A), in comparison to that of specific sucrose uptake rate (800%; Figure A1a in Appendix A). The predicted storage root growth rate increased as the carbohydrate and starch contents were increased, relative to the default values ( Figure 4a,d), but decreased as the protein content was increased (Figure 4c). Variations in the lipid and fiber contents resulted in a minimal response (Figure 4b,e).
In contrast to the growth rate, the simulated carbon flux distribution was sensitive to the composition of biomass components. The biomass components-related variation in flux distribution could be inferred from the positions of their breakpoints (red circles in Figure 4) relative to the default values. The flux distribution was greatly influenced by the carbohydrate content with a decrease of 1% and an increase of 4%, relative to the default value, resulting in changes in flux distribution. The protein, fiber, lipid, and starch contents had less effect on the predicted fluxes, retaining similar predictions as the original values until a breakpoint of around −50%.
In addition, the impact of the biomass was studied, based on the direction of the individual reaction fluxes and the partitioning of each metabolic branch point. The results suggested that the direction of predicted fluxes was highly associated with the composition of biomass ( Figures A9-A12 in Appendix A and Supplementary Materials File S5). The variations in the direction of metabolic fluxes and breakpoints, due to changes in the biomass content, are shown in Figure 5. Each breakpoint contained at least one reaction flux exhibiting a directional change relative to the default value, highlighted in red. The differences in flux direction were found in similar reactions (red-colored areas in the black-dashed rectangle; Figure 5) within the starch and sucrose biosynthesis pathway, pentose phosphate pathway, and respiration pathway in cytosol and plastid, and amino acid biosynthesis pathway in the cytosol.    Further analysis of the predicted flux partitioning showed that it was also highly influenced by the changes in biomass contents (Figures 6 and A13-A15 in Appendix A, and Supplementary Materials File S6). Reactions with a similar flux partitioning as the original model were color-coded similarly, whereas, those indicating a breakpoint were color-coded differently ( Figure 6). The results indicated that the predicted flux direction and flux partitioning had similar sensitivity, as indicated by the shared breakpoints (−1% and +4% for carbohydrates, −50% and +6% for proteins, −50% for fibers, −60% for lipids, and −70% for starch). The metabolic fluxes in the plastid were more sensitive to variations in biomass content than those in the cytosol. The sensitive reaction fluxes were primarily found in the pentose phosphate pathway, followed by the starch and sucrose biosynthesis pathway. D-xylulose-5-phosphate and D-erythrose-4-phosphate branch points in the pentose phosphate pathway showed relatively high sensitivity, both in the cytosol and plastid compartments ( Figure 6II,VI), reflecting their primary role as metabolic precursors for energy balance [46]. Metabolic fluxes in the starch and sucrose biosynthesis pathway in the cytosol were highly influenced by the protein content ( Figure 6I), whereas, carbohydrate-related biomass components were more important in the plastid ( Figure 6V). The starch and sucrose biosynthesis pathway and pentose phosphate pathway play important roles in carbon utilization for storage root growth, as demonstrated by the rMeCBM-KU50 model [4].

Discussion
The constraint-based modeling approach has been adopted to help unravel the complexity of plant metabolism. The reliability of this approach is very much dependent on the quality of data used and its representativeness of real systems, among others. This study investigated the sensitivity of the rMeCBM-KU50 model to inputted data and constraints used for the prediction, which included (i) specific sucrose uptake-related data (photosynthetic rate, total photosynthetic rate, total leaf area,

Discussion
The constraint-based modeling approach has been adopted to help unravel the complexity of plant metabolism. The reliability of this approach is very much dependent on the quality of data used and its representativeness of real systems, among others. This study investigated the sensitivity of the rMeCBM-KU50 model to inputted data and constraints used for the prediction, which included (i) specific sucrose uptake-related data (photosynthetic rate, total photosynthetic rate, total leaf area, and storage root dry weight) and (ii) biomass function-related data. The results showed the specific storage root growth rate was highly dependent on the specific sucrose uptake rate (Figure 3), whereas, the carbon flux distribution was more sensitive to the proportion of biomass components (Figures 4-6).
Our analysis showed that the storage root dry weight and total photosynthetic rate had a remarkable effect on the specific sucrose uptake rate and also, the specific storage root growth rate prediction (Figure 3). The effects of the photosynthetic rate and total leaf area were comparatively minimal, despite being important indicators of plant growth [47]. These results suggest storage root dry weight and total photosynthetic rate can reliably be used to model the storage root growth rate in cassava under different conditions, an approach, which if implemented could enable researchers to channel limited resources and focus on the most relevant plant physiological traits in their quest to understand plant systems.
The cellular biomass composition, which is the proportion of constituent macromolecular compounds such as carbohydrates and proteins, is indicative of the metabolic products and perhaps, more importantly, reflects the physiological conditions, among others, to which the cells are exposed. An accurate representation of the biomass is, therefore, important for analyzing the metabolic flux [32]. Sensitivity analysis showed that the rMeCBM-KU50 model was more sensitive to variations in the biomass function than it was to changes in the sucrose uptake-related data. The biomass content had a greater effect on the metabolic flux prediction than the sucrose uptake-related data (Figures 3 and 4). Also, the carbon flux prediction by the rMeCBM-KU50 model was most sensitive to the carbohydrate content, which accounts for up to 90% of cassava storage roots biomass on a dry weight basis. Similar observations were reported in the yeast model [31] and oilseed rape model [28]. The biomass of yeast predominantly contains protein, 35-65% by weight [31], whereas, the main biomass component of oilseed rape is oil mostly in the form of triacylglycerol, accounting for over 40% by weight [28]. However, the flux distribution in Arabidopsis remained unchanged when the biomass function was varied, probably because the biomass components are more comparable [30].
Taken together, our study demonstrated that the rMeCBM-KU50 model is a well representative model of carbon assimilation metabolism in cassava storage roots. The model was able to predict the storage root growth rate under typical growing conditions (Figures 3 and 4), and showed sensitivity to the specific sucrose uptake rate, total photosynthetic rate, storage root dry weight, and carbohydrate content. Based on the findings, the specific sucrose uptake rate, carbohydrate content, and starch content could reliably be used as the important data for model prediction and used as the guideline of experimental design for model development.

Conclusions
A sensitivity analysis was performed to determine how sensitive the rMeCBM-KU50 model is to changes in input data and constraints, enabling us to understand and quantify the impact of the range of variables on the model outcome. By mostly varying the substrate uptake-related data and biomass function within ±90% of the model default values, we aimed to better understand the input-output relationship and identify variables that are most critical to the modeling process under specific conditions, allowing for a focused, wider application of the model. The results showed the total photosynthetic rate, storage root dry weight, and carbohydrate content had predominant impacts on the rMeCBM-KU50 model. The predicted storage root growth rate was highly dependent on the specific sucrose uptake rate through total photosynthetic rate, and storage root dry weight. The carbon flux distribution, flux direction and flux partitioning, was most sensitive to changes in biomass content, particularly carbohydrate, the major biochemical compound in cassava storage root. The sensitivity of reaction flux direction to variations in sucrose uptake rate-related data, File S4: The sensitivity of reaction flux partitioning (%) to variations in sucrose uptake rate-related data, File S5: The sensitivity of reaction flux direction to variations in biomass function-related data, File S6: The sensitivity of reaction flux partitioning (%) to variations in biomass function-related data.
Author Contributions: T.S. and S.K. supervised the project; P.C. performed all computational modeling and prepared the figures; T.S. and P.C. conducted the part of analyzed the results and drafted the manuscript; T.S., S.K., and P.C. discussed, interpreted the results, wrote the article, and approved the final manuscript.

Acknowledgments:
The authors would like to thank The National Center for Genetic Engineering and Biotechnology (BIOTEC, NSTDA) for P.C. post-graduate scholarship. This project was supported by National Science and Technology Development Agency (NSTDA) and National Research Council of Thailand (NRCT) (P-13-50395 and P-15-50974). Furthermore, we acknowledge the financial support provided by KMUTT through the KMUTT 55th Anniversary Commemorative Fund and the International Strategic Output and Outcome Fund (the outbound student and researcher). We also thank the Biological Modeling Laboratory (BML) members and KMUTT for useful suggestions.

Conflicts of Interest:
The authors declare no conflict of interest.     The percentage changes of predicted growth rate with respect to the original data. Each radar graph, the radial axes is the error percentage of model simulation to measured storage root growth rate. LA; total leaf area per plant; P N , photosynthetic rate; P N, total , total photosynthetic rate per plant; SR DW , storage root dry weight per plant; C, carbohydrate content; F, fiber content; L, lipid content; P, protein content; S, starch content.

Appendix A. Supplemental Figures
Processes 2019, 7, 259 19 of 37 growth rate. LA; total leaf area per plant; PN, photosynthetic rate; PN, total, total photosynthetic rate per plant; SRDW, storage root dry weight per plant; C, carbohydrate content; F, fiber content; L, lipid content; P, protein content; S, starch content.     partitioning (breakpoint). The production and consumption fluxes at each metabolic branch point were respectively denoted as positive and negative signs. LA; total leaf area per plant; PN, photosynthetic rate; PN, total, total photosynthetic rate per plant; SRDW, storage root dry weight per plant.