We created three types of mutant rice lines using the procedures described in Pérez et al. (2022). Mutant Type I had exogenous HMGR; Mutant Type II had exogenous HMGS, HMGR and MVK; and Mutant Type III had exogenous HMGS, HMGR, MVK, PMK and MVD. We placed these six transgenes (BjHMGS, tHMGR, CrMK, CrPMK, and CrMVD) in three independent expression cassettes driven by endosperm-specific promoters. A transit peptide inserted at the beginning of the five enzymes in the MVA pathway directs them to the plastid. See Pérez et al. (2022) and Supplementary Appendix S1 for the full details.

To create the mutant lines, we bombarded seven-day-old mature zygotic rice embryos (Oryza sativa cv. EYI105) with gold particles coated with the transformation vectors. We recovered transgenic plantlets and regenerated and hardened them off in soil. Genomic DNA was isolated from the callus and leaves of regenerated plants to confirm presence of the BjHMGS, tHMGR, CrMK, CrPMK, CrMVD and OsWR1 through PCR (Pérez et al., 2022). We recovered 12 independent type I mutant lines, 10 independent type II mutant lines, and 12 independent mutant lines.

We analyzed all rice lines after 12 weeks of growing in soil.

For each line we measured chlorophyll levels and analyzed the cytokinins trans-zeatin, zeatin riboside and isopentenyl adenine (iP); the gibberellins GA1, GA3 and GA4; the auxin indole-3-acetic acid (IAA); ABA; salicylic acid; jasmonic acid; and the ethylene precursor 1-aminocyclopropane-1-carboxylic acid (ACC) as described in Pérez et al. (2022) and Supplementary Appendix S1 .

We also measured gene expression for BjHMGS, tHMGR, CrMK, CrPMK, CrMVD, and the endogenous MVA and MEP pathway genes OsHMGS, OsHMGR, OsMK, OsPMK, OsMVD, OsDXS, OsDXR, OsMCT, OsCMK, OsMDS, OsHDS, OsHDR and OsIPPI using qRT-PCR as described in Appendix S1.

We counted the number of leaves, measured the height of the plants from the base of the stem to the maximum extension of the flag leaf, and measured the length and maximum width of the last expanded leaf as described in Pérez et al. (2022). We used a SPAD meter at six points on the last expanded leaf to quantify leaf chlorophyll. We multiplied length and maximum width of the last expanded leaf by a correction factor of 0.75 to estimate leaf area.

All experimental measurements are provided in Supplementary Data S1 and described in the results section of Appendix S1 .

We used ordinary differential equation systems to model the biosynthesis of IPP/DMAPP. The mathematical formalism used to describe the flux dynamics is the saturating formalism (Sorribas et al., 2007; Alves et al., 2008). This formalism allows us to approximate the kinetics of any given reaction using a rational expression, where parameters have physical interpretations that are analogous to those found in classical enzyme kinetics rate expressions. In this formalism, we approximate the rate of a reaction in an inverse space at an operating point by:

V parameters represent apparent saturation rate constants for the reactions. Ki parameters represent apparent binding constants for the substrate(s) or inhibitor(s) of the reaction. While no activators were considered in our model, these can also be included using this formalism.

We modelled the wild type IPP and DMAPP production (i.e., the endogenous MVA and MEP pathways), using the canonical reaction set for each pathway, shown in Figure 1 . We modelled the kinetics of each process, as well as those for the exchange fluxes of IPP and DMAPP between cytoplasm and plastid, using the rate expressions in Supplementary Tables S1 – S3 . We assume that the organism is able to maintain homeostasis of Acetyl-CoA and Acetoacetyl-CoA.

To model type I mutants, we added the reaction that transforms HMG-CoA_pl into MVA_pl to the plastid, as well as the cytoplasm-plastid exchange reactions for these two metabolites ( Supplementary Tables S3 – S4 ).

We extended the model for type I mutants to create the model for type II mutants. We included the reactions catalyzed by HMGS and MVK to the plastid compartment ( Supplementary Table S4 ). Corresponding compartment exchange reactions for the substrates and products of these enzymes are also added ( Supplementary Table S3 ).

We extended the model for type II mutants to create the model for type III mutants. We added the metabolite MVPP to the plastid by including the reaction that synthesize it. We also added the reaction that transforms MVPP_pl into IPP_pl, ( Supplementary Table S4 ). The plastid-cytoplasm exchange flux of MVPP_pl is included as well ( Supplementary Table S3 ).

Under physiological conditions, IPP and DMAPP mostly flow from the plastid into the cytosol (Bick and Lange, 2003). We implemented this observation by assuming that metabolites flow from the plastid to the cytosol at ten times the rate of the import reaction from the cytosol. Bick and Lange (Bick and Lange, 2003) also reported that other pathway intermediates were not actively transported between the two compartments. Other studies confirm this observation (for example, Wright et al. (2014)). However, those same studies show that in mutants overexpressing DXS, there is a second pool of MEcPP outside the chloroplast. In addition early intermediates of the MVA pathway can be found in the plastid space (Schneider et al., 1977).

Introducing an MVA pathway into the plastid as we did, may cause changes in the flux of MVA intermediates between the plastid and the cytoplasm. As such, we allowed for the possibility that HMG-CoA, MVA, MVP and MVPP enter and leave the plastid, albeit at very slow rates. Supplementary Table S3 summarizes all reactions of material interchanged between plastid and cytosol.

Each metabolite has its own differential equation in the model. The kinetic rate function, f j , for each process that produces a metabolite M appears as a positive term in the differential equation that determines the dynamic behavior of that metabolite. Similarly, the kinetic rate function, f k , for each process that consumes a metabolite appears as a negative term in the differential equation that determines the dynamic behavior of that metabolite:

For the wild type and for each mutant type we assemble a type-specific system of ordinary differential equations (ODEs) that describes the dynamic behavior of all metabolites in the system. As such, we have four different ODE model types. These four models describe the dynamic behavior of all metabolites in the system in a type-specific manner.

Supplementary Tables S5 presents the basal kinetic constants for each reaction in the four model types. Supplementary Table S6 presents the concentrations for the independent variables of the four model types.

To personalize the models and make them line specific we need to weigh the rate constants of the relevant type-specific model by the variations in enzyme activity of the individual line of interest. To do so we searched the literature for information about the correlation between changes in gene expression and enzyme activities in the MEP and MVA pathways. As we found no such information, we modeled variations in the enzyme activities of the mutant lines as described in Comas et al. (2016): changes in gene expression with respect to the WT are assumed to be proportional to changes in protein activity. This is the simplest possible assumption about the relationship between changes in gene expression and changes in enzyme activity.

We implement this assumption in the models by explicitly considering the enzymes that catalyze each reaction in the rate expressions. As V m a x ≈ k c a t E n z y m e , the model for the WT sets the enzyme activity to be 1 (the basal state). As we model mutant lines, we assume that changes in gene expression are proportional to changes in enzyme activity. Thus,

In the end we obtained one line-specific model for the wild type rice, 12 line specific models for type I mutant lines, 10 line specific models for type II mutants, and 12 line specific models for type III mutants.

While many biological phenomena are rhythmic, overall, biological systems survive because they can achieve homeostasis (Wang et al., 2022). In other words, metabolism remains buffered and stable against normal environmental fluctuations. This occurs for the MEP pathway (Wright et al., 2014). Mathematically this situation is described by a stable steady state. It is well known that, when modeling homeostatic situations, lack of stability is a good diagnostic tool for model incompleteness (Savageau, 1975; Grimbs et al., 2007; Schmidt et al., 2011; Voit, 2013). For example, Ni and Savageau (1996a; 1996b) used this type of diagnostic tool to predict regulatory interactions that could stabilize a model of the red blood cell metabolism. Because of this, we perform stability analysis of the line-specific rice models in order to identify possible model improvements that can stabilize unstable steady states.

An efficient way to assess stability is by calculating the eigenvalues of the Jacobian matrix of the ODE system, which are complex numbers (Voit, 2013). If the real parts of all eigenvalues are negative, the system is stable. Otherwise, the system is unstable. The Jacobian matrix is constructed by taking the partially derivatives of the right-hand side of the ODEs ( f i ) with respect to each state variable ( x j ), as shown in Eq. 4.

In addition to steady state stability, another tool for model diagnostic is steady state robustness. Reasonable models generate steady states that are robust, and have low sensitivity to parameter changes (Savageau, 1975; Voit, 2013). Sensitivity measures how much a dependent variable or output changes when a parameter is altered (Comas et al., 2016). Parameters with high sensitivities tend to identify where information may be incomplete or inaccurate.

As such, we performed a sensitivity analysis to identify which steps of the pathway could have additional regulation that we were ignoring. We calculated logarithmic, or relative, steady-state parameter sensitivities, which measure the “relative change in a system variable (X) that is caused by a relative change in a parameter (p)” (Voit, 1991):

This sensitivity analysis generates a matrix of sensitivities for each line. Each element of the matrix (of dimensions m × n ) represents the sensitivity S i , j of metabolite M i to parameter p j . To facilitate visualization of sensitivity analysis results and comparison between lines, we compressed the sensitivity analysis matrices for each line in two ways.

First, to see how much a line is sensitive to a certain parameter over all metabolites (variables of the system) we calculate the following index (Comas et al., 2016):

In other words, we calculate the size (or Euclidean norm) of the vector whose components are the sensitivity of each metabolite to parameter p j , normalized by the number of metabolites in each mutant line (for example, models for Type I lines have 16 metabolites, while those for Type III have 18). We use the Euclidean norm of the sensitivity vectors as a way to represent aggregate sensitivities with a single metric to facilitate visual representation and analysis. Further, and because the number of metabolites increases from type I to type II and from type II to type III models, we make S ' j comparable between models by normalizing it by the number of metabolites considered in the model. Second, to see how sensitive a metabolite is to all parameters in a line we calculate the following index:

As with Eq 6, we calculate the size (or Euclidean norm) of the vector whose components are the sensitivity of the same metabolite   M i to each parameter p j , normalized by the number of parameters in each mutant line. We use the Euclidean norm of the sensitivity vectors as a way to represent aggregate sensitivities with a single metric to facilitate visual representation and analysis. Further, and because the number of parameters increases from type I to type II and from type II to type III models, we make S ' j comparable between models by normalizing it by the number of parameters considered in the model.

To investigate if we could use plant hormone levels as predictors of dynamic behavior in IPP/DMAPP biosynthesis in seeds we performed correlation analysis between hormones and metabolites, as well as genes, as described in Section 1.9 of Appendix S1. Significant effects were then included in the ODE models using one of two possible formalisms:

where M i is the concentration of metabolite i in the model, H_j is the level of hormone j at twelve weeks, and α and g are constants. We chose between the two alternatives in the following way. First we adjust a linear model of Log[ M i ] as a function of Log[ H j ] . A combination of low adjusted R 2 and high | g i j | suggests a potentially strong influence of the hormone levels on metabolite concentrations (high d | g i j | ) over a small range of hormone levels (low adjusted R 2 ). In this situation, we assumed a saturation effect and used Eq 9 to model hormone influence on metabolite production and consumption. Otherwise, we used Eq 8, as it uses a smaller number of parameters and minimizes the possibility of overfitting the model to the data. The threshold for selecting the one or the other formalism was set at 0.5 for the ratio | R a d j 2 / g | . If | R a d j 2 / g | > 0.5 we use the power law formalism. The right-hand side of the equations modify the ODEs by multiplying the terms that involve production/consumption of the metabolite and involving those enzymes whose genes levels also correlate to hormone levels, in a way that makes the observed correlations affect the production and consumption rates of the metabolites.

For a more detailed procedure, see the Supplementary Appendix S1 .

We note that, when hormone levels were below the experimental detection threshold, we reverted the kinetic expression presented in Supplementary Table S7 to the original model, using a piece-wise approximation to solve the differential equations.

We used a form of forward stepwise regression (Efroymson, 1960) to investigate how the different phenotypic variables might be predicted from hormones levels, gene expression and metabolite concentrations. We analyzed the following plant phenotypic characteristics: Height, number of Leaves, Leaf Length, Leaf Width, and Chlorophyl levels. We split experimental data according to mutant type, so that the analysis and model building was performed three times, one for each mutant type. We investigate phenotype as a function of the predictor variables gene expression, hormone levels, and metabolite concentrations.

The first step of the regression analysis was building independent linear models with one predictor variable.

The second step of the regression analysis was to select the predictor variables that had a significant ( α = 0.05 ) effect on the phenotype and whose model had an adjusted R² greater than 0.2. If only one model has a significant effect, we would choose that one. If more than one predictor variable has a significant effect in explaining the predicted variable, we chose the model for the variable with the highest adjusted R². If the adjusted R² is similar between models, we chose the model with the lowest AICc (AIC corrected for small sample sizes) score. The lower this score, the lower the chance that a model over fits the observations. If the AICc is similar for more than one model, we selected the predictor variable with the highest adjusted R². At this stage, we have a one variable model.

The third step of the analysis was to create models with all possible combinations of predictor variables where one of the elements of the pair is the predictor variable selected in step one. We then selected the best two variable models as described in the previous paragraph, while making sure they are not collinear.

We repeated steps two and three and stopped when adding a new variable did not improve the explanatory power of the model. Thus, for a given set of significant predictors { x 1 , x 2 ,   … , x n } , the model would be:

We used forward stepwise regression instead of the more traditional multilinear modelling approach that starts from Eq. 10 and eliminates all variables that have no effect because the number of data points is smaller than the number of parameters to fit to the full multilinear model.

The full mathematical description of IPP/DMAPP biosynthesis consists of 14, 16, 17 and 18 differential equations for the wild type (WT), Mutant Type I (MT-I), Mutant Type II (MT-II) and Mutant Type III (MT-III), respectively. Eq 11 shows the overall ODE systems for the four model types:

Table 1 provides the concentrations of each metabolite estimated from the model in the WT line. The system can achieve homeostasis (stable steady state in mathematical nomenclature). Stable steady states have negative real parts for the eigenvalues of the system’s Jacobean matrix ( Table 2 ). The model also estimates that pathway substrates (HMG-CoA and DXP) and end-products (DMAPP and IPP) concentrations are, in general, larger than those of intermediate metabolites, which is another hallmark of a well-behaved biosynthetic pathway (Alves and Savageau, 2000).

To understand how perturbations in parameters may affect the ability of the system to maintain homeostasis, we calculated the logarithmic sensitivity of the steady state Jacobian eigenvalues to each parameter of the model ( Supplementary Table S8 ). The model has over eighty parameters and eigenvalues have sensitivities that are above one (in absolute values) to thirty of them. The parameters to which more eigenvalues are sensitive concentrate in reactions r2 (HMG-CoA_cyt → MVA_cyt), r3 (MVA_cyt → MVP_cyt), r4 (MVP_cyt → MVPP_cyt), and r6 (IPP_cyt → DMAPP_cyt) of the MVA pathway and reactions r10 (Glyceraldehyde-3-P + Pyruvate → DXP) and r18 (IPP_pl → DMAPP_pl) of the MEP pathways. This suggests that modifying isoprenoid biosynthesis could destabilize the physiological steady states of the plant.

Sensitivity analysis identifies the parameters to which the various variables of the model are most sensitive, as described in (Sorribas et al., 2007; Alves et al., 2008). A high sensitivity of a variable to a parameter indicates that small changes in that parameter might lead to big changes in the variable.

Plausible models of biological systems have low sensitivities to most parameters (Savageau, 1976; Kitano, 2007). The logarithmic sensitivity analysis of the dependent concentrations of the WT model with respect to each parameter of the model we performed shows that our model fits this quality criterion. Only 51 out of 728 sensitivities are larger than one ( Supplementary Table S9 ). DMAPP and IPP are the metabolites with the highest sensitivities. High sensitivities are well known to identify the parts of a system that need to be modeled in more detail when additional information becomes available (Savageau, 1976; Kitano, 2007). This is consistent with the fact that we modeled IPP and DMAPP usage only through simple sink reactions, without considering any metabolic and regulatory details.

We implement the models for the biosynthetic pathways in each mutant line using the same procedure as that for the WT line (sections 2.3 to 2.6). These systems do not reach homeostasis, having unstable steady states with a few intermediate metabolites accumulating over time. This suggests that the models are not plausible representations of the biological situation (Savageau, 1971; Savageau, 1975; Savageau, 1976; Voit, 2013).

Biological systems can stabilize steady states and reach homeostasis by adjusting the activity of enzymes in a pathway, for example through post-transcriptional regulation of protein levels and activity. We investigated whether emulating this type of adjustment would stabilize the steady states in the models.

First, we identified the metabolites that accumulated in each line, which were DXP, CDP-MEP, MEcPP, or combinations thereof. Reactions r11, r14, and r15 of Table S2 either produce or consume these metabolites. Using a minimal intervention policy, we scanned the values for the Vmax parameters of reactions Vmax9, Vmax12 and Vmax13 in order to identify the minimum change in those parameters that would stabilize the steady state of each mutant line.

To stabilize homeostasis in the models we scanned one-dimensional, two-dimensional and three-dimensional parameter spaces and found the values of Vmax that stabilized the steady state in each model. The sets of Vmax that made the model for each mutant line stable were stored in a candidate sets list and we chose the final parameter set as the one with minimum normalized Euclidian distance to the original set of parameter values for that line. We provide the list of stabilized parameter values for each line in Supplementary Data S1 – Model_stabilization. Supplementary Figure S1 shows that Type II and Type III mutants require larger changes in parameter values than Type I mutants.

We performed a sensitivity analysis of the stable homeostatic concentrations with respect to each enzyme parameter, in each mutant line ( Figure 2A ). We find that those concentrations are very robust, with 3% of all individual sensitivities being higher than 1 in absolute value in Type I mutants. This number goes down to 1.8% in Type II mutants and 1.4% in Type III mutants. The total number of individual sensitivities calculated for each line is higher than 900.

Globally, we found that, for each parameter p j , the aggregated sensitivity S ' j   of all metabolites to that parameter decreases in mutant lines with respect to the WT. The global sensitivities to each parameter decrease in the following order: WT>Type I lines> Type II lines> Type III lines. We also found that the aggregated sensitivity of each metabolite M i to all parameters also decreases in the same order ( Figure 2B ). Thus, our results suggest that post-transcriptional regulation of a small number of enzymes is sufficient to maintain homeostasis of IPP/DMAPP biosynthesis in each mutant line.

While having a line-specific, data-driven, model is a more accurate way of describing and predicting the behavior of each mutant line, these are less than helpful in predicting how a new mutant line of any of the three types will behave dynamically.

To create general, type-specific models that are more useful for predicting the dynamic behavior and characteristics of generic new mutant lines, we created a median experimental line for each mutant type. To do so we use the median gene expression activities for each gene in all lines of a given mutant type. Then, we follow the procedure described in methods to generate three new models, one per mutant type. Their steady state concentrations, stability, and sensitivity analysis in Supplementary Table S10 ; Figure 2 ; Supplementary Data S1 . Figure 2 also shows that these lines have sensitivity profiles that are similar to those of the individual mutant lines of the same type. Moreover, the models for the median lines of each mutant type have homeostatic behavior that is robust to mutations in enzyme parameters, which is a hallmark of a plausible model.

We also wanted to investigate whether early plant hormone levels in the plant might be a proxy for subsequent changes in the biosynthesis of IPP/DMAPP in seeds. To do so, we calculated how changing hormone levels could explain changes in gene expression and metabolite levels, as described in the methods section 2.11 and in Supplementary Appendix S1 . We are not assuming that the endogenous MEP, MVA, and exogenous MVA Pathways have an influence in hormone production, only measuring if there is a phenomenological correlation between variations in hormone levels and metabolites or levels of gene expression.

iP is a proxy for changes in enzyme activity in all three mutant types. For all mutant types, iP correlates to changes in the activity of the early MEP pathway steps. In addition, for type I mutants, iP also correlates to changes in the activity of MVA pathway early steps. ABA is also a proxy influence in the early steps of the MEP pathway for all mutant types. Other hormones have a mutant specific effect. Table S7 summarizes the results for all hormones and presents the hormone dependency equations for each mutant type.

To validate the resulting multilevel models, we investigated if they could reproduce the correlation between experimental hormone levels and model metabolites in the following way.

For each mutant type and metabolite whose concentration is significantly affected by a given hormone, we took the median model described in Section 3.6 and calculated the concentration of the various metabolites as a function of hormone levels. Then, we calculated the correlation in the simulation plot and compared that correlation to the one computed when we plot concentration of the same metabolite vs experimental hormone levels. We summarize the results of this analysis in Table 3 . We find that the models maintain 48 out of 53 expected correlations between hormone levels and metabolite concentrations. This is consistent with the extended models being plausible multilevel descriptions of IPP/DMAPP biosynthesis in the three mutant types.

We estimated how the production of IPP/DMAPP changes across mutant lines and types, with the help of the models. Figure 3 summarizes the results. Overall, our models suggest that median metabolic flux going into the plastid’s IPP/DMAPP producing pathways increases in the following way: Type III>Type II>Type I. We dissected the production rates of IPP by HDR, IDI and MVD (both endogenous and ectopic, where applicable), and the exchange from plastid to cytosol. We also find that the mean production flux for IPP increases from Type I to Type II to Type III.

We used linear regression to investigate if hormone, gene expression, and metabolite levels can predict macroscopic plant phenotypes, such as plant height, number of leaves, leaf width, leaf length and chlorophyll levels, as described in methods. Table 4 shows that the subset of phenotype parameters that can be predicted from molecular data is different for different mutant types and provides the best fit models for each phenotype parameter and mutant type. Leaf length and chlorophyll content can be predicted for types I, II, and III. Height can be predicted for types I and III. Leaf width can be predicted for types II and III. The number of leaves can be predicted for mutant types I and II. Variations in the levels of the metabolites, genes and hormones shown in Table 4 are the best predictors for the variations measured in the macroscopic phenotype. The metabolic levels of DXP, MEcPP, HMG-CoA are useful in predicting phenotypical characters of leaves in all mutant types. DMAPP levels are useful in predicting chlorophyll content in type III mutants. Table 4 also includes the adjusted R² for each model, which is a measure of the percentage of variation in the phenotype that can be explained by the model. Models explain between 40% and 77% of the phenotypical variability, depending on the specific phenotype being measured and the mutant type.

The results from the previous section allow us to create a multilevel model, connecting metabolite concentrations, and hormone and gene expression levels to the macroscopic plant phenotypes. We summarize the multilevel model building process in Figure 4 . In Figure 5 we use the models from Table 4 and the experimental data for gene expression and hormone levels to calculate what is the expected value for the predicted phenotype, according to the relevant model. Then, we include the experimental determination for the same phenotype. We show that each type-specific model can semi-quantitatively predict macroscopic plant phenotype.

We further investigated if the type-specific models were accurate in predicting the phenotype of the other mutant types. We used the median model for each mutant type, feeding it with the experimental determinations for all the lines and measured how accurate the phenotype predictions were for each type. What we observed is that the type-specific models do not properly predict other mutant types. For example, Figure 6 shows that our type I model can only accurately predict chlorophyll levels for type I lines, failing to do so for mutant types II and III. Supplementary Table S11 shows that phenotype predictions are only accurate when made with the model for the correct mutant type. This suggests that prediction for new mutant types would require developing a data driven model for that mutant type.

Several modeling efforts focused on analyzing the biosynthetic and signaling dynamics of complex terpenoids in plants (Latowski et al., 2000; Bruggeman et al., 2001; Rios-Estepa et al., 2008; Rios-Estepa et al., 2010; Band et al., 2012; Beguerisse-Diaz et al., 2012; Pokhilko et al., 2013; Pokhilko et al., 2015; Allen and Ptashnyk, 2017; Nazareno and Hernandez, 2017; Dalwadi et al., 2018; Rizza et al., 2021). For example, Band et al. and Rizza et al. use compartmental modeling to study the biosynthesis and diffusion of gibberellins in root tips (Band et al., 2012; Rizza et al., 2021). Terpenoid signaling is often also a target for modeling. For example, Allen & Ptashnyk use models to study signaling interactions between brassinosteroid and gibberellin signaling pathways (Allen and Ptashnyk, 2017), and Nazareno & Hernandez do the same to study signaling interactions between of abscisic acid, ethylene and methyl jasmonate on stomatal closure (Nazareno and Hernandez, 2017). Modeling studies of the biosynthesis and regulation of terpenoid precursors in plants are less common. In fact, we are only aware of two such modeling effort for the MEP pathway in plants (Pokhilko et al., 2015) (Neiburga et al., 2023) and another in Plasmodium falciparum (Singh and Ghosh, 2013). The plant MEP model was used to study how circadian rhythms regulate the dynamics of the pathway in plants, while the P. falciparum model was used to investigate the regulation of the pathway and to predict the effects of genetic manipulations on the production of isoprenoids with the addition of in silico inhibitors. Regarding the MVA pathway, we found no model in plants. Still, this pathway was modeled in bacteria (Weaver et al., 2015, Dalwadi et al., 2018). These MVA models study the dynamics of the pathway in the context of introducing it in bacteria using synthetic biology. Finally, we know of only one other example where both pathways were modeled together in yeast, using Petri net-based dynamic modeling (Baadhe et al., 2012). Taking all this into account, creating models that can be used to study the dynamics and interactions of the MVA and MEP pathways in plants is an important goal, towards which this paper contributes.

Trusting the models normally requires that they are validated by comparison with experimental data that was not used to build them. This was one of the technical limitations in building line-specific models, as we needed most of the available data to estimate the parameter values for each line. As such, the measurements can be used for either model building or model validation, but not both simultaneously. To sidestep this problem, we used three approaches.

First, we used sensitivity and stability analysis as general model quality assessment tools (Savageau, 1975; Kitano, 2007; Voit, 2013) to both evaluate the quality of the line-specific models and identify the parameters that could be used to improve that quality. The WT model is of high quality, being robust to changes in parameters (>92% of parameter sensitivities lower than 0.5) and producing a stable steady state, with metabolite concentrations that are well within the accepted biological ranges (Albe et al., 1990). In contrast, the steady states for the original line specific mutant lines are unstable, leading to unbound accumulation of MEP pathway intermediates. This indicated that the models needed improvement. We hypothesized that the simplest reason for model instability could be consequence of a nonlinear relationship between changes in the expression of pathway genes and changes in enzyme activity. To test if the hypothesis is consistent with our data, we scanned the parameters that represent the enzyme activities that directly produce or use the accumulating metabolites in the model: Vmax9, Vmax12 and Vmax 13. By executing this procedure, we identified the sets of minimal changes to the values for these parameters that generated stable and robust steady states, with metabolic concentrations within known physiological ranges. We remark that other, more complex explanations might also be consistent with the experimental data. Still, Occam’s razor argues that the simplest explanation is the most likely one, in the absence of additional data (Borgqvist and Palmer, 2022; Piasini et al., 2023).

Second, we created “median” models for each of the three types of pathways. To do so, we pooled together all mutant lines of a given type and calculated the median of the pool for each variable. Then, we used that median to create the median model for each mutant type, in the same way we create line specific models. These median models also allowed us to estimate the metabolites for each individual line. Comparing the results with the line specific models shows that these median models have a similar behavior to the models of the individual lines of the same mutant type ( Figure 2 ). In addition, when we use the median models to predict the phenotype of the individual lines, the predictions have an error that is similar to that of the individual line models. This suggests that we can use the median model of a given mutant type to study newly created lines for that mutant type.

Third, we used the multilevel models to predict plant phenotype and compare the results with the experimentally determined phenotype, achieving a prediction accuracy of up to 80% ( Supplementary Table S11 ). We summarize the process in Figure 4 . We make different assumptions about the relationship between variables in our modeling. We model the effect of changes in gene expression on the concentrations of pathway intermediates assuming a direct cause-effect relationship between changes in the expression of a gene and changes in the corresponding enzyme activity. We model the effect of enzyme activity on the flux going through the reaction catalyzed by that enzyme using traditional enzyme kinetics. In both cases this assumes a causal relationship between variables. In contrast, we assume that there might be a phenomenological relationship between changes in hormone levels and changes in gene expression and phenotypes and test for that relationship. When we find statistical evidence for that relationship, we include it in our models in different ways. The influence of hormone levels on gene expression is added to the metabolic pathway models through the use of approximation theory, in a way that is mathematically well justified (Salvador, 1996; Sorribas et al., 2007; Alves et al., 2008). The phenomenological influence of hormones, genes, and metabolites on phenotype was accounted for by using statistical linear models. As additional experimental studies become available, the phenomenological parts of the models will need to be adjusted in order to account for the knowledge generated by those experiments. Overall, the quality assessment steps we performed suggest that our models can be used as reasonable semi-quantitative prediction tools to help in better understanding the biology of isoprenoid biosynthesis modification in rice.

Posttranscriptional regulation is important for the proper functioning of the MEP and MVA pathways in plants (Laule et al., 2003; Guevara-García et al., 2005; Sauret-Güeto et al., 2006; Flores-Pérez et al., 2008; Xie et al., 2008; Cordoba et al., 2009; Han et al., 2013). Our modeling and analysis indicate that IPP/DMAPP production is robust to fluctuations in enzyme activity in WT rice ( Table 2 ; Supplementary Tables S8 , S9 ). Further, it suggests that posttranscriptional modulation of enzyme activity is important in stabilizing IPP/DMAPP production in mutant lines ( Figure 2 ; Supplementary Table S2 ). In all cases, stabilizing the steady state of a mutant line requires that the activity of a protein is upregulated with respect to the changes in gene expression for that protein ( Supplementary Data S1 ). This is fully consistent with the observation that, when compared to the WT Arabidopsis plants, changes in the activities of DXS and DXR proteins are bigger than the changes in expression of the corresponding genes (compare panels A and C of Figure 3 in Flores-Pérez et al. (2008)). Interestingly, DXR is one of the proteins flagged in our models as a potential stabilizing influence for the steady state of the rice mutants ( Supplementary Data S1 ).

Our results also suggest that variations in plant hormone levels can predict, to some extent, both plant phenotype ( Table S11 ) and the biosynthetic fluxes of IPP/DMAPP in the seeds ( Table S7 ). Further, variations in plant gene expression levels in combination with variations in plant hormones can improve phenotype predictions ( Tables 4 , S11 ). However, the more complex the genetic manipulation was, the less accurate the phenotype predictions become. While the median adjusted R² for the predictions is approximately 60% in mutant types I and I, this number goes down to 35% in type III mutants ( Table S11 ).

Pérez et al. (2022) also reported that plant development is more similar to that of the WT in mutant type III, followed by mutant type II, and finally I. Our modeling suggests an explanation for this observation. The analysis predicts that the global production of IPP/DMAPP is on average higher in the mutant types that have a more complete version of the exogenous MVA pathway in the plastid ( Figure 3 ). The average amount of flux going through the endogenous MEP and MVA pathways in Type III mutants is the most similar to that of the WT, followed by the flux going through the endogenous pathways in Type II mutants. The least similar flux to WT is that of Type I mutants ( Figure 3 ). Figure 3 also shows that the total amount of flux going through the IPP/DMAPP pools in mutant types I and II is similar to that of the WT. The flux going through IPP/DMAPP is larger in Type III mutants than in the WT rice. Taken together, these observations suggest that plants can distinguish between the IPP/DMAPP produced by each of the pathways. Too little flux going through the endogenous MEP pathway compromises the production of developmental hormones leading to plants with developmental defects. In conclusion, we believe that an iterative modeling-experimental process as the one presented here would be an effective way to identify which parts of each pathway are more sensitive to further manipulation, and which parts are more likely to be good targets for modification in order to increase the flux without disrupting the development of the plant.