Modeling the effect of daytime duration on the biosynthesis of terpenoid precursors

Terpenoids are valued chemicals in the pharmaceutical, biotechnological, cosmetic, and biomedical industries. Biosynthesis of these chemicals relies on polymerization of Isopentenyl di-phosphate (IPP) and/or dimethylallyl diphosphate (DMAPP) monomers, which plants synthesize using a cytosolic mevalonic acid (MVA) pathway and a plastidic methyleritritol-4-phosphate (MEP) pathway. Circadian regulation affects MVA and MEP pathway activity at three levels: substrate availability, gene expression of pathway enzymes, and utilization of IPP and DMAPP for synthesizing complex terpenoids. There is a gap in understanding the interplay between the circadian rhythm and the dynamics and regulation of the two pathways. In this paper we create a mathematical model of the MVA and MEP pathways in plants that incorporates the effects of circadian rhythms. We then used the model to investigate how annual and latitudinal variations in circadian rhythm affect IPP and DMAPP biosynthesis. We found that, despite significant fluctuations in daylight hours, the amplitude of oscillations in IPP and DMAPP concentrations remains stable, highlighting the robustness of the system. We also examined the impact of removing circadian regulation from different parts of the model on its dynamic behavior. We found that regulation of pathway substrate availability alone results in higher sensitivity to daylight changes, while gene expression regulation alone leads to less robust IPP/DMAPP concentration oscillations. Our results suggest that the combined circadian regulation of substrate availability, gene expression, and product utilization, along with MVA- and MEP-specific regulatory loops, create an optimal operating regime. This regime maintains pathway flux closely coupled to demand and stable across a wide range of daylight hours, balancing the dynamic behavior of the pathways and ensuring robustness in response to cellular demand for IPP/DMAPP.

1 Introduction

Terpenoids are a family of molecules with more than 22,000 different natural products (Harborne et al., 1991). Many family members have crucial biological functions. For example, in plants, they work as hormones (gibberellin, abscisic acid, etc.), photosynthetic pigments (chlorophyll, phytol, and carotenoids), electron carriers (ubiquinone, plastoquinone), mediators of the assembly of polysaccharides (polyprenyl phosphates) and structural components of membranes (phytosterols). They are also used for other purposes, such as antibiotics, herbivore repellents, toxins and pollinator attractants (Mcgarvey and Croteau, 1995).

Plants synthesize terpenoids from two metabolic precursors: Isopentenyl di-phosphate (IPP) and dimethylallyl diphosphate (DMAPP). Two compartmentally separated pathways synthesize these precursors (Figure 1). The mevalonic acid (MVA) pathway converts acetyl-CoA (Ac-CoA) to IPP and DMAPP. This pathway is mostly cytosolic, with a couple of reactions taking place in the peroxisome. IPP and DMAPP are then used in the synthesis of phytosterols and ubiquinone (Mcgarvey and Croteau, 1995). The enzyme 3-hydroxy-3-methylglutaril-CoA reductase (HMGR) is a key enzyme in the regulation of the MVA pathway (Schaller et al., 1995). The methyleritritol-4-phosphate (MEP) pathway is compartmentalized in plastids and is responsible for the production of carotenoids, lateral chains of chlorophylls, plastoquinone, abscisic acid (ABA) and tocopherols (vitamin E, precursors and derivatives) (Eisenreich et al., 2001).

Figure 1

Figure 1. Representation of the two terpenoid biosynthesis pathways plus the ectopic pathway, the MVA pathway (left, cytosol and peroxisome) and the MEP pathway (right, plastid). DXR, DXP reductoisomerase; MCT, 2-C-methyl-D-erythrtle 4-phosphate cytidylyl transferase; CDP-ME, 4-(Citidine 5’-difosfo)-2-C-methyl-D-eritritol; CMK, 4-difosfocitidil-2-C-methyl-D-erythrtol kinase; CDP-MEP, 2-Fosfo-4-(cytidine 5’- diphospho)-2-C-methyl-D-eritritol; MDS, 2-C-methyl-D-eritritol 2,4-cyclodifosphate synthase; MEcPP, 2-C-methyl-D-eritritol 2,4-cycdiphosphate; HDS, 4-hydroxy-3-methylbut-2-en-1-il diphosphate synthase; HMBPP, 4-hydroxy-3-methylbut-2-in-1-il diphosphate; HDR, 4-hydroxy-3-methylbut-2-en-1-il diphosphate reductase; IDI, isopentenyl diphosphate Delta-isomerase; PhyPP, phytyl diphosphate.

While metabolite tracing indicates that each of the two pathways is responsible for the production of a subset of terpenoid compounds downstream, there is evidence of crosstalk between them (Hemmerlin et al., 2003, 2012; Hemmerlin, 2013), with some of intermediates in both pathways diffusing between the cytosol and the plastid (Bick and Lange, 2003; Hemmerlin et al., 2003, 2006; Laule et al., 2003). The first intermediate of the MEP pathway, DXP, can diffuse between the plastid and the cytoplasm (Hemmerlin et al., 2003; Page et al., 2004; Lange et al., 2015). At the level of IPP and DMAPP, this exchange was measured to occur mainly in the plastid-to-cytoplasm direction, promoted by a one-way symport system (Bick and Lange, 2003; Dudareva et al., 2005). The direction of this metabolic exchange between cellular compartments may depend on physiological state and species. There is lack of convincing evidence that other intermediates of both pathways can diffuse between the two compartments (Hemmerlin, 2013), and it would be interesting to understand what the effect of losing this exchange might have on the production of IPP/DMAPP in each compartment.

Several studies used mathematical modeling to elucidate how using synthetic biology to modify terpenoid metabolism might lead to changes in the regulatory and dynamic behavior of that metabolism. When it comes to microorganisms, for example, a kinetic model of the MVA pathway in E. coli using parameters from the literature correctly predicts expression and inhibition changes (Weaver et al., 2015). Petri nets were also used to model the integration of both pathways in yeast (Baadhe et al., 2012). Another example is a mathematical model of the MEP pathway in the malaria parasite P. falciparum (Singh and Ghosh, 2013). This ODE 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.

Understanding the regulation of terpenoid biosynthesis in plants is an important biological issue that has relevance also for synthetic biology (Harborne et al., 1991; Mcgarvey and Croteau, 1995; Hemmerlin et al., 2012; Liao et al., 2016; Tetali, 2019; Zhou and Pichersky, 2020). One of the master regulators of metabolism in plants is the circadian rhythm imposed by Earth’s rotation. Mathematical models were also used to study the effect of that rhythm on the dynamics of terpenoid precursor biosynthesis in the MEP pathway of peppermint leaves (Rios-Estepa et al., 2010) and Arabidopsis (Pokhilko et al., 2015; Neiburga et al., 2023). These three models focus on and include metabolites downstream of the precursors IPP and DMAPP, and do not include the intermediates upstream.

Measurements of volatile terpenoid emissions in plants show circadian oscillations with peaks during the day and a general decrease at night (Loivamäki et al., 2007; Zeng et al., 2017; Zheng et al., 2017; Picazo-Aragonés et al., 2020; Mu et al., 2022). These oscillations are driven by light regulation and internal circadian clocks in the MEP pathway. There is evidence of light regulated DXR expression and internal circadian clock regulation downstream via regulation of isoprene synthase. However, the internal clock oscillations are best maintained when coordinated by light cycles (Loivamäki et al., 2007). Gene expression analysis has shown that MEP pathway genes DXS, DXR, CMK, MCT, MDS, HDS, HDR and IDI present circadian oscillation, as well as many other downstream genes involved in carotene, tocopherol and other phytohormone biosynthesis, their expression levels rising during the day and decreasing through the night (Covington et al., 2008; Zheng et al., 2017).

The key genes for circadian regulation that are conserved across many plant species are CCA1, LHY, TOC1 and ZTL (Nagel and Kay, 2012; McClung, 2013). CCA1 and LHY inhibit TOC1 and vice versa, while ZTL tags TOC1 for degradation under blue light, which favors the morning components CCA1 and LHY. These morning components then upregulate MEP pathway genes (Alabadí et al., 2001; Yon et al., 2017; Picazo-Aragonés et al., 2020).

Currently, the use of mathematical models to explore the dynamics and regulatory mechanisms within the MVA pathway, and more critically, the intricate interplay between the MEP and MVA pathways, remains in its embryonic stages (Basallo et al., 2023). To our knowledge (Basallo et al., 2023), provide the only example where this is done, creating and analyzing a comprehensive model that delineates these pathways in plants, accounting for all intermediates leading to the crucial precursors IPP and DMAPP. This lack of attention may stem from the predominant focus on the MEP pathway for synthesizing chemical species of interest in plants, while MVA-derived products, equally significant, have been predominantly studied in microorganisms. As such, integrating both pathways into mathematical models that can be used as tools to study the integrated dynamics and regulation of both pathways stands as a paramount necessity for comprehensive understanding and exploration in plant biochemistry and metabolism.

In this paper, we adapted the mathematical model of the MVA and MEP pathways in Oryza sativa (rice) (Basallo et al., 2023) to account for the regulation of the kinetic activity in those pathways by the circadian rhythm. We then used the model to investigate how the changes in that rhythm over the year and at different latitudes affect IPP and DMAPP biosynthesis. Finally, we performed a set of in silico experiments where we removed circadian regulation from different parts of the model to investigate how the various regulatory loops contribute to the dynamic behavior of IPP and DMAPP biosynthesis. This enabled us to examine how different regulatory designs for the network (the network genotype) influence the adaptation of the dynamic behavior of IPP/DMAPP biosynthesis (the network phenotype) to the changes in daylight hours that occur due to the Earth’s circadian rhythm.

2 Materials and methods

2.1 Mathematical modeling formalism

To model the biosynthesis of IPP/DMAPP, we employed systems of ODEs. The saturating formalism was employed as the mathematical framework to depict the flux dynamics (Sorribas et al., 2007; Alves et al., 2008). This formalism approximates the kinetics of any given reaction to parameters that have biochemical interpretations in enzyme kinetics. In this formalism, we approximate the rate of a reaction in an inverse space at an operating point by:

$\begin{array}{ll}v\approx \frac{V{\prod }_{i=1}^m{x}_i^{}}{{\prod }_{i=1}^m(K_i+x_i^{})+{\prod }_{b=1}^p(x_b^{}+K_b)}& (1)\end{array}$

V parameters represent apparent saturation rate constants for the reactions. K_i 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.

2.2 Mathematical models for the MVA and MEP pathways

We adapted the model published in (Basallo et al., 2023), and summarized in Tables 1–3. There, the pathways are modeled using the canonical reaction set for each pathway (shown in Figure 1), extracted from KEGG and the literature consensus. Exchange of IPP and DMAPP between the cytosol (produced by the MVA pathway) and the plastid (produced by the MEP pathway) is modeled by considering that, under physiological conditions, IPP and DMAPP mostly flow from the plastid into the cytosol (Bick and Lange, 2003). This is implemented by setting the flow from the plastid to the cytosol to be ten times the rate of the import to the cytosol. Bick and Lange (Bick and Lange, 2003) also reported that other pathway intermediates were not actively transported between the two compartments. Table 3 summarizes all reactions of material interchanged between plastid and cytosol.

Table 1

Table 1. MVA pathway reactions that were considered in the model.

Table 2

Table 2. MEP pathway reactions that were considered in the model.

Table 3

Table 3. Exchange of MVA and MEP intermediates between the cytosol and the plastid.

We modeled the kinetics of each step, as well as those for the exchange fluxes of IPP and DMAPP between cytoplasm and plastid, using the rate expressions in Tables 1-3. We assume that the organism maintains homeostasis of Acetyl-CoA and Acetoacetyl-CoA. Table 4 presents the kinetic constants for each reaction, retrieved from (Basallo et al., 2023) and from the primary literature. Table 5 collects the concentrations for the independent variables. Hereafter, we refer to this model as Model A. The regulatory structure of each model is summarized in Table 6.

Table 4

Table 4. Kinetic Parameters.

Table 5

Table 5. Concentration of independent Variables (Albe et al., 1990).

Table 6

Table 6. List of models and the regulation modules they contain.

2.3 Stability analysis

We assess stability of the steady states 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 (${\text{x}}_j$), as shown in Equation 2.

$\begin{array}{ll}\mathit{J}\mathbf{=}{\mathit{D}}_{\mathit{x}}\mathit{f}\mathbf{=}{\mathit{f}}_{\mathit{x}}\mathbf{=}\frac{\mathbf{\partial }{\mathit{f}}_{\mathit{i}}}{\mathbf{\partial }{\mathit{x}}_{\mathit{j}}}\mathbf{=}\left(\begin{array}{cccc}\frac{\mathbf{\partial }{\mathit{f}}_{\mathbf{1}}}{\mathbf{\partial }{\mathit{x}}_{\mathbf{1}}}& \frac{\mathbf{\partial }{\mathit{f}}_{\mathbf{1}}}{\mathbf{\partial }{\mathit{x}}_{\mathbf{2}}}& \mathbf{\cdots }& \frac{\mathbf{\partial }{\mathit{f}}_{\mathbf{1}}}{\mathbf{\partial }{\mathit{x}}_{\mathit{n}}}\\ \frac{\mathbf{\partial }{\mathit{f}}_{\mathbf{2}}}{\mathbf{\partial }{\mathit{x}}_{\mathbf{1}}}& \frac{\mathbf{\partial }{\mathit{f}}_{\mathbf{2}}}{\mathbf{\partial }{\mathit{x}}_{\mathbf{2}}}& \mathbf{\cdots }& \frac{\mathbf{\partial }{\mathit{f}}_{\mathbf{2}}}{\mathbf{\partial }{\mathit{x}}_{\mathit{n}}}\\ \mathbf{⋮}& \mathbf{⋮}& \mathbf{\ddots }& \mathbf{⋮}\\ \frac{\mathbf{\partial }{\mathit{f}}_{\mathit{n}}}{\mathbf{\partial }{\mathit{x}}_{\mathbf{1}}}& \frac{\mathbf{\partial }{\mathit{f}}_{\mathit{n}}}{\mathbf{\partial }{\mathit{x}}_{\mathbf{2}}}& \mathbf{\cdots }& \frac{\mathbf{\partial }{\mathit{f}}_{\mathit{n}}}{\mathbf{\partial }{\mathit{x}}_{\mathit{n}}}\end{array}\right)& (2)\end{array}$

2.4 Sensitivity analysis

Logarithmic sensitivity analysis of the system was performed by 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):

$\begin{array}{ll}\overline{\mathit{S}}(\mathit{X}\mathbf{,}\mathit{p})\mathbf{=}\frac{\mathbf{\partial }\mathit{X}/\mathit{X}}{\mathbf{\partial }\mathit{p}/\mathit{p}}\mathbf{=}\frac{\mathbf{\partial }\mathit{log}\mathit{X}}{\mathbf{\partial }\mathit{log}\mathit{p}}& (3)\end{array}$

We also calculated the steady state sensitivities of the aggregate input flux ($S_{ij}^{+}$) and aggregate output flux ($S_{ij}^{-}$) of each metabolite to the model parameters using (Equation 4).

$${\dot{\mathit{M}}}_{\mathit{i}}\mathbf{=}{\mathbf{\sum }}^{\mathbf{ }}\mathit{f}{(\mathit{M}(\mathit{p})\mathbf{,}\mathit{p})}_{\mathit{m}}\mathbf{-}{\mathbf{\sum }}^{\mathbf{ }}\mathit{ɡ}{(\mathit{M}(\mathit{p})\mathbf{,}\mathit{p})}_{\mathit{n}}$$

$\begin{array}{ll}{\mathit{S}}_{\mathit{i}\mathit{j}}^{\mathbf{+}}\mathbf{=}\frac{\mathbf{\partial }{\mathbf{\sum }}^{\mathbf{ }}\mathit{f}{({\mathit{M}}_{\mathit{S}\mathit{S}}\mathbf{,}\mathit{p})}_{\mathit{m}}}{\mathbf{\partial }{\mathit{p}}_{\mathit{j}}}& (4)\end{array}$

$${\mathit{S}}_{\mathit{i}\mathit{j}}^{\mathbf{-}}\mathbf{=}\frac{\mathbf{\partial }{\mathbf{\sum }}^{\mathbf{ }}\mathit{g}{({\mathit{M}}_{\mathit{S}\mathit{S}}\mathbf{,}\mathit{p})}_{\mathit{n}}}{\mathbf{\partial }{\mathit{p}}_{\mathit{j}}}$$

In Equation 4 $f{(M(p),p)}_m$ represents the mathematical function that describes production flux $f{(.,.)}_m.$ This function depends on a set of metabolites $M(p)$ and parameter vector $p$. Similarly, $g{(M(p),p)}_n$ represents the mathematical function that describes consumption flux $g{(.,.)}_n.$ M_SS is the steady state value of metabolite M_i and p is the set of parameters involved in the reactions, while p_j is the specific parameter with respect to which the aggregated sensitivity S_ij is calculated.

2.5 Modeling the circadian rhythm

To model the effect of the circadian oscillation on the dynamics of Model A we adapted the approach used by (Pokhilko et al., 2014, 2015), described by Equation 5:

$\begin{array}{ll}\mathit{L}(\mathit{t})=\mathbf{0.5}\left(\left(\mathbf{1}+\mathit{t}\mathit{a}\mathit{n}\mathit{h}\left(\frac{\mathit{t}-\mathit{p}\mathit{e}\mathit{r}\mathit{i}\mathit{o}\mathit{d}·\mathit{F}\mathit{l}\mathit{o}\mathit{o}\mathit{r}(\frac{\mathit{t}}{\mathit{p}\mathit{e}\mathit{r}\mathit{i}\mathit{o}\mathit{d}})-\mathit{d}\mathit{a}\mathit{w}\mathit{n}}{\mathit{T}}\right)\right)-\left(\mathbf{1}+\mathit{t}\mathit{a}\mathit{n}\mathit{h}\left(\frac{\mathit{t}-\mathit{p}\mathit{e}\mathit{r}\mathit{i}\mathit{o}\mathit{d}·\mathit{F}\mathit{l}\mathit{o}\mathit{o}\mathit{r}(\frac{\mathit{t}}{\mathit{p}\mathit{e}\mathit{r}\mathit{i}\mathit{o}\mathit{d}})-\mathit{d}\mathit{u}\mathit{s}\mathit{k}}{\mathit{T}}\right)\right)+\left(\mathbf{1}+\mathit{t}\mathit{a}\mathit{n}\mathit{h}\left(\frac{\mathit{t}-\mathit{p}\mathit{e}\mathit{r}\mathit{i}\mathit{o}\mathit{d}·\mathit{F}\mathit{l}\mathit{o}\mathit{o}\mathit{r}(\frac{\mathit{t}}{\mathit{p}\mathit{e}\mathit{r}\mathit{i}\mathit{o}\mathit{d}})-\mathit{p}\mathit{e}\mathit{r}\mathit{i}\mathit{o}\mathit{d}}{\mathit{T}}\right)\right)\right)& (5)\end{array}$

$L(t)$ can have values between 0 (no light) and 1 (maximal diurnal light intensity). The Floor function returns the greatest integer less than or equal to the input value. We set the period of oscillation to 24 h, the duration of the day. While we do not do so, we note that this function could also be used to model complex circadian situations, for example with dissimilar or very long dusk and dawn periods. To model the positive effect of the presence of light on a variable X_i or a parameter p_j we use Equation 6:

$${\mathit{X}}_{\mathit{i}}(\mathit{t})\mathbf{=}\mathit{L}(\mathit{t})\mathbf{·}{\mathit{X}}_{\mathit{i}\mathbf{0}}$$

$\begin{array}{ll}{\mathit{p}}_{\mathit{j}}(\mathit{t})\mathbf{=}\mathit{L}(\mathit{t})\mathbf{·}{\mathit{p}}_{\mathit{j}\mathbf{0}}& (6)\end{array}$

In contrast, to model the negative effect of the presence of light on a variable X_i or a parameter p_j we use Equation 7:

$\begin{array}{ll}{\mathit{X}}_{\mathit{k}}(\mathit{t})\mathbf{=}(\mathbf{1}\mathbf{-}\mathit{L}(\mathit{t}))\mathbf{·}{\mathit{X}}_{\mathit{k}\mathbf{0}}& (7)\end{array}$

$${\mathit{p}}_{\mathit{l}}(\mathit{t})\mathbf{=}(\mathbf{1}\mathbf{-}\mathit{L}(\mathit{t}))\mathbf{·}{\mathit{p}}_{\mathit{l}\mathbf{0}}$$

2.6 Modeling the effect of latitude and seasonality on the number of hours of light per day

Parameters dawn and dusk control the time of sunrise and sunset respectively. By setting dawn to 0 h, we can control daytime length through changing only the dusk parameter. The default value we use for dusk is 12 h. This leads to $L(t)$ having 12h of light and 12h of dark, which is the situation close to the equator. Increasing dusk above 12h makes $L(t)$ have more than 12h of daylight within the 24h period of the oscillation. Decreasing dusk below 12h makes $L(t)$ have less than 12h of daylight within the 24h period of the oscillation. In this way, dusk is a proxy for the effect of geographic latitude on the number of daylight hours within a circadian oscillation.

Parameter T emulates twilight duration. At T close to 0 h, the function approximates a step function (square wave) while progressively higher values of T lead to a less steep transition between darkness and full light levels (sinusoidal wave).

The number of daylight hours in a day changes with the seasons and with latitude. In the peak of the northern hemisphere winter, there are 0h of daylight in the north pole and 24h of daylight in the south pole. Supplementary Figure S1 illustrates how the L function can model the number of full daylight, twilight and night hours at different latitudes and during the duration of a single year.

2.7 Modeling the effects of circadian oscillations on the dynamics of the MVA and MEP pathways

There are three regulatory circadian modules we consider in our model: regulation of substrate production for the MVA and MEP pathways, regulation of gene expression in the pathways, and regulation of IPP/DMAPP consumption after they are synthesized.

At the level of substrate production (Cockburn and McAulay, 1977), show that triose intermediates of glycolysis remain roughly constant over the circadian light cycle. In contrast, they also show that pyruvate in plants leaves can oscillate over the daylight cycle, changing over two-fold with respect to its average value (Pokhilko et al., 2014, 2015) assume that the circadian light cycle controls the availability of the triose G3P, which is a precursor of the MEP pathway. Similarly, the availability of Acetyl-CoA in the cytoplasm of plant leaves also changes over the circadian light cycle (Buchanan et al., 2002; Taiz et al., 2014; Lee et al., 2021). In the morning, acetyl-CoA levels rise, as photosynthesis becomes more active, reaching peak levels at midday. Acetyl-CoA availability decreases during the afternoon, remaining low as photosynthesis is absent and metabolic activity is reduced during the night. We follow the experimental evidence by Cockburn and McAulay (1977), combining it with the approach used by Pokhilko et al. (2014, 2015) and model the effect of circadian rhythms on the availability of MVA and MEP precursors using Equation 8:

$$\mathit{A}\mathit{c}\mathbf{–}\mathit{C}\mathit{o}\mathit{A}(\mathit{t})\mathbf{=}{\mathit{L}}_{\mathit{s}\mathit{u}\mathit{b}\mathit{s}\mathit{t}\mathit{r}\mathit{a}\mathit{t}\mathit{e}}(\mathit{t})\mathbf{·}\mathit{A}\mathit{c}\mathbf{-}\mathit{C}\mathit{o}{\mathit{A}}_{\mathit{H}}$$

$\begin{array}{ll}\mathit{P}\mathit{y}\mathit{r}\mathit{u}\mathit{v}\mathit{a}\mathit{t}\mathit{e}(\mathit{t})\mathbf{=}{\mathit{L}}_{\mathit{s}\mathit{u}\mathit{b}\mathit{s}\mathit{t}\mathit{r}\mathit{a}\mathit{t}\mathit{e}}(\mathit{t})\mathbf{·}\mathit{P}\mathit{y}\mathit{r}\mathit{u}\mathit{v}\mathit{a}\mathit{t}{\mathit{e}}_{\mathit{H}}& (8)\end{array}$

Here, Pyruvate _H and Ac-CoA _H represent the steady state values for pyruvate and acetyl-coA in model A (Table 5).

Model A modified through the addition of Equation 8 will be referred to hereafter as Model B. Model B reverts to Model A when L(t) = 1.

At the level of gene expression regulation, experimental evidence indicates that there is anti-phasic regulation of MVA and MEP pathway genes by daylight (Covington et al., 2008; Vranová et al., 2013; Atamian and Harmer, 2016; Jin et al., 2021). The data suggests that daylight activates MEP pathway genes (Jin et al., 2021) and deactivates MVA genes. To model this effect, we made protein activity time dependent and directly proportional to light level, using Equation 9:

$$\mathit{V}\mathit{m}\mathit{a}{\mathit{x}}_{\mathit{M}\mathit{V}\mathit{A}\mathbf{ }\mathit{i}}(\mathit{t})\mathbf{=}(\mathbf{1}\mathbf{-}\mathit{L}(\mathit{t}))\mathbf{·}\mathit{V}\mathit{m}\mathit{a}{\mathit{x}}_{\mathit{M}\mathit{V}\mathit{A}\mathbf{ }\mathit{i}\mathit{H}}$$

$\begin{array}{ll}\mathit{V}\mathit{m}\mathit{a}{\mathit{x}}_{\mathit{M}\mathit{E}\mathit{P}\mathbf{ }\mathit{j}}(\mathit{t})\mathbf{=}\mathit{L}(\mathit{t})\mathbf{·}\mathit{V}\mathit{m}\mathit{a}{\mathit{x}}_{\mathit{M}\mathit{E}\mathit{P}\mathbf{ }\mathit{j}\mathit{H}}& (9)\end{array}$

Index i represents each enzyme of the MVA pathway, while index j represents each enzyme of the MEP pathway. $\mathit{V}\mathit{m}\mathit{a}{\mathit{x}}_{\mathit{M}\mathit{V}\mathit{A}} {\mathit{i}}_{\mathit{H}}$ represents the basal value of $Vma{x}_{MVA i}$ in Model A (Table 5). Hereafter we will refer to the Model A modified with Equation 9 as Model C.

At the level of IPP/DMAPP consumption, and following (Pokhilko et al., 2014, 2015) and the experimental evidence (Loivamäki et al., 2007; Covington et al., 2008; Zeng et al., 2017; Zheng et al., 2017; Picazo-Aragonés et al., 2020; Mu et al., 2022), we model usage of IPP and DMAPP downstream of the MVA as being dependent of the circadian rhythm. In fact, the expression of genes from pathways that synthesize more complex terpenoids was observed to be coordinated to that of the genes from the MEP pathway (Jin et al., 2021). It is well established experimentally that plant emission of terpenoid oils, which are derived from IPP and DMAPP precursors is high during the light hours of the day and strongly decreases during nighttime. We model this effect by modifying the rate constants of reactions r20 and r21 according to Equation 10:

$${\mathit{k}}_{\mathbf{1}}(\mathit{t})\mathbf{=}\mathit{L}(\mathit{t})\mathbf{·}{\mathit{k}}_{{\mathbf{1}}_{\mathit{H}}}$$

$\begin{array}{ll}\mathit{k}{\mathbf{\text{'}}}_{\mathbf{1}}(\mathit{t})\mathbf{=}\mathit{L}(\mathit{t})\mathbf{·}{{\mathit{k}}^{\mathbf{\text{'}}}}_{\mathbf{1}\mathit{H}}& (10)\end{array}$

Here, k_1H and k’_1H are the values of k1 and k’1 used in Model A. Model A modified with Equation 10 will be referred to hereafter as Model D.

By combining Models B, C, and D we investigated how changes in the regulatory action of circadian rhythms might affect the dynamic behavior of the MEP and MVA pathways. This is the mathematical equivalent of mutating the genome to eliminate, create, or modify regulatory loops with the purpose of studying their effect. Thus, in Model BC we eliminate regulation of IPP/DMAPP consumption by the circadian rhythm, in Model BD we eliminate regulation of MVA and MEP gene expression by circadian rhythms, and in Model CD we eliminate regulation of pathway substrate production by circadian rhythms. Finally, we assemble the modifications of Models B, C, and D in a single model, which we refer to as Model E.

2.8 Model Implementation

The Mathematica code for the implementation of all models and Figures is provided in Supplementary Data Sheet S1.

3 Results

3.1 Quality analysis of the basal model for the MEP and MVA pathways

Equation 11 represents the joint basal mathematical model (Model A hereafter) of the MEP and MVA pathways adapted from (Basallo et al., 2023):

$$\dot{\mathbf{\text{HMG-CoAcyt}}}\mathbf{=}{\mathit{r}}_{\mathbf{1}}\mathbf{–}{\mathit{r}}_{\mathbf{2}}$$

$$\dot{{\mathbf{\text{MVA}}}_{\mathbf{\text{cyt}}}}\mathbf{=}{\mathit{r}}_{\mathbf{2}}\mathbf{–}{\mathit{r}}_{\mathbf{3}}$$

$$\dot{{\mathbf{\text{MVP}}}_{\mathbf{\text{cyt}}}}\mathbf{=}{\mathit{r}}_{\mathbf{3}}\mathbf{–}{\mathit{r}}_{\mathbf{4}}$$

$$\dot{{\mathbf{\text{MVPP}}}_{\mathbf{\text{cyt}}}}\mathbf{=}{\mathit{r}}_{\mathbf{4}}\mathbf{–}{\mathit{r}}_{\mathbf{5}}$$

$$\dot{{\mathbf{\text{IPP}}}_{\mathbf{\text{cyt}}}}\mathbf{=}{\mathit{r}}_{\mathbf{5}}\mathbf{+}{\mathit{r}}_{\mathbf{7}}\mathbf{+}{\mathit{r}}_{\mathbf{9}}\mathbf{–}{\mathit{r}}_{\mathbf{6}}\mathbf{–}{\mathit{r}}_{\mathbf{8}}\mathbf{–}\mathbf{4}{\mathit{r}}_{\mathbf{20}}$$

$\begin{array}{ll}\dot{{\mathbf{\text{IPP}}}_{\mathbf{\text{pl}}}}\mathbf{=}{\mathit{r}}_{\mathbf{8}}\mathbf{+}{\mathit{r}}_{\mathbf{16}}\mathbf{+}{\mathit{r}}_{\mathbf{19}}\mathbf{–}{\mathit{r}}_{\mathbf{9}}\mathbf{–}{\mathit{r}}_{\mathbf{18}}\mathbf{–}\mathbf{3}{\mathit{r}}_{\mathbf{21}}& (11)\end{array}$

$$\dot{{\mathbf{\text{DMAPP}}}_{\mathbf{\text{cyt}}}}\mathbf{=}{\mathit{r}}_{\mathbf{6}}\mathbf{+}{\mathit{r}}_{\mathbf{23}}\mathbf{–}{\mathit{r}}_{\mathbf{7}}\mathbf{–}\mathbf{2}{\mathit{r}}_{\mathbf{20}}\mathbf{–}{\mathit{r}}_{\mathbf{22}}$$

$$\dot{{\mathbf{\text{DMAPP}}}_{\mathbf{\text{pl}}}}\mathbf{=}{\mathit{r}}_{\mathbf{17}}\mathbf{+}{\mathit{r}}_{\mathbf{18}}\mathbf{+}{\mathit{r}}_{\mathbf{22}}\mathbf{–}{\mathit{r}}_{\mathbf{19}}\mathbf{–}{\mathit{r}}_{\mathbf{21}}\mathbf{–}{\mathit{r}}_{\mathbf{23}}$$

$$\dot{\mathbf{\text{DXP}}}\mathbf{=}{\mathit{r}}_{\mathbf{10}}\mathbf{–}{\mathit{r}}_{\mathbf{11}}$$

$$\dot{\mathbf{\text{MEP}}}\mathbf{=}{\mathit{r}}_{\mathbf{11}}\mathbf{–}{\mathit{r}}_{\mathbf{12}}$$

$$\dot{\mathbf{\text{CDP}}–\mathbf{\text{ME}}}\mathbf{=}{\mathit{r}}_{\mathbf{12}}\mathbf{–}{\mathit{r}}_{\mathbf{13}}$$

$$\dot{\mathbf{\text{CDP}}–\mathbf{\text{MEP}}}\mathbf{=}{\mathit{r}}_{\mathbf{13}}\mathbf{–}{\mathit{r}}_{\mathbf{14}}$$

$$\dot{\mathbf{\text{MEcPP}}}\mathbf{=}{\mathit{r}}_{\mathbf{14}}\mathbf{–}{\mathit{r}}_{\mathbf{15}}$$

$$\dot{\mathbf{\text{HMBPP}}}\mathbf{=}{\mathit{r}}_{\mathbf{15}}\mathbf{–}{\mathit{r}}_{\mathbf{16}}\mathbf{–}{\mathit{r}}_{\mathbf{17}}$$

Details about the mathematical form and parameter values of each flux function r_i are given in Tables 1-3. Here we investigate whether the model has a positive, biological reasonable, steady state that is stable and robust. Table 7 provides the steady state concentration of the variables in the Model. The concentrations are within normal metabolite concentration ranges. In addition, the steady state is stable, having negative real parts for all the eigenvalues of the system’s Jacobian matrix (Table 8). The model also predicts that the concentrations of pathway intermediates are lower than those of its substrates (HMG-CoA and DXP) and end-products (DMAPP and IPP), which is another hallmark of a well-behaved biosynthetic pathway (Alves and Savageau, 2000).

Table 7

Table 7. Predicted concentrations of metabolites (mM).

Table 8

Table 8. Eigenvalues for the steady state.

Sensitivity analysis identifies the parameters to which the 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 the value of that parameter might lead to large changes in the value of the variable.

To understand how sensitive the stability of the steady state is to perturbations in the parameters of the models, we calculated the logarithmic sensitivity of the steady state Jacobian eigenvalues to each parameter of the model (Table 9). The model has over eighty parameters. Eigenvalues have sensitivities above one (in absolute value) to thirty of those parameters. 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. The eigenvalues whose real part is closer to zero are most sensible to Vmax3, Vmax4, Vmax5, Km3, Ki31, and Km4, parameters from r3 (MVA_cyt → MVP_cyt), r4 (MVP_cyt → MVPP_cyt) and r5 (MVPP_cyt → IPP_cyt).

Table 9

Table 9. Logarithmic sensitivities of the eigenvalues to reaction parameters.

Plausible models of biological systems have low sensitivities to most parameters (Savageau, 1976; Kitano, 2007). The logarithmic sensitivity analysis of the dependent concentrations 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 0.5 and none is larger than 3 in absolute value (Table 10). DMAPP and IPP are the metabolites with the highest sensitivities. The parameters responsible for these high sensitivities are the maximum velocities of isomerization between DMAPP and IPP. In general, the parameters causing the highest sensitivities for each metabolite correspond to a reaction directly involved in producing or consuming that metabolite. In addition, the metabolites on the MEP pathway seem to share a high sensitivity to parameters from rate r10. DXS catalyzes this reaction, where G3P and pyruvate (the substrates of the pathway) produce DXP.

Table 10

Table 10. Logarithmic sensitivities of the concentrations to reaction parameters.

Parameters from r10 (G3P + Pyruvate → DXP), catalyzed by DXS, are overrepresented in the top 51 highest sensitivities (14 out of 51), particularly Vmax8 and Km8 and for MEP pathway intermediates (Table 10). The highest absolute values of DXS-related sensitivities correspond to MEcPP sensitivities. This is consistent with experimental results showing that DXS activity is a determinant of MEcPP levels (Wright et al., 2014).

When we analyze the sensitivity of the positive fluxes to parameters, we find that 25 out of 728 sensitivities are above 0.5 in absolute value. Similarly, 21 out of 728 sensitivities of the negative fluxes are above 0.5 (Table 11). There is no specific parameter being overrepresented, and the pattern is that fluxes are most sensitive to Vmax and Km parameters from one of the cognate parameters of each flux. Interestingly, the highest sensitivities of the IPP/DMAPP input fluxes are to parameters from the IDI enzyme activity, not MVD or HDR.

Table 11

Table 11. Logarithmic sensitivities of the aggregate flux to reaction parameters.

3.2 IPP and DMAPP levels are robust to circadian-dependent flux decrease

Model E implements regulatory effects of light on the biosynthesis of terpenoid precursors IPP and DMAPP. Thus, we analyze this model to understand how changes in the number of daylight hours affect IPP/DMAPP production.

First, we simulate the behavior of the system during two circadian oscillations. Comparing Supplementary Figure S2, where we represent pathway intermediates to Figure 2, where we show the dynamic curves for pathway end products illustrates how their concentration increases through the day and decreases during the night, over 48h at different locations on the globe and in different seasons.

Figure 2

Figure 2. Time course simulation of IPP and DMAPP concentrations throughout 48h at different latitudes and times of the year: equator/spring and fall equinoxes (dusk = 12h), middle latitudes (winter, dusk = 9h; summer, dusk = 15h) and near polar circle latitudes (winter, dusk = 3h; summer, dusk = 21h). T = 1h. Green lines – MEP pathway. Magenta lines – MVA pathway.

When there are approximately 12h of daylight, the concentrations of intermediate metabolites have small oscillations in both pathways throughout the 24h of the day (Supplementary Figure S2). This occurs consistently at or near the equator and in temperate latitudes during spring and fall. As the number of daylight hours per day increases, the relative amplitude of oscillation decreases. Similarly, as the number of daylight hours per day decreases, that amplitude increases. These changes are more pronounced for oscillations in the concentration of intermediate metabolites of the MVA pathway than for intermediates of the MEP pathway. The average daily concentrations of pathway intermediates are only slightly affected by the number of daylight hours. The picture is subtly reversed for the end products of the MVA and MEP pathways (Figure 2). The number of daylight hours per day has almost no influence in the relative amplitude of IPP and DMAPP concentration oscillations. In contrast, the daily average concentrations of IPP and DMAPP are directly correlated with the number of daylight hours.

Furthermore, the circadian rhythm significantly influences the system’s transient behavior. Cytosolic IPP and DMAPP steadily accumulate during daylight hours and decrease at night. In contrast, plastid IPP and DMAPP quickly reach quasi-steady state levels during daylight hours and gradually decrease at night. This behavior slightly changes during longer days, where quasi-steady states are reached for all four metabolite pools.

We further investigate how changing number of daylight hours between 0 and 24 affect IPP and DMAPP biosynthesis by simulating the dynamic behavior of the system when the number of daylight hours changes between 0 and 24. Then we measure the impact of those changes on concentration oscillations (Figures 3A, B) and production fluxes (Figures 3C–E) within the pathways. This experiment confirms that IPP concentration oscillations exhibit higher amplitudes than DMAPP (Figure 2). Additionally, we observe that the relative amplitude of IPP and DMAPP concentration oscillations remains relatively stable, except when the number of daylight hours approaches 0 or 24 (indicated by blue lines in Figure 3 panels), where the oscillation becomes a steady state with zero amplitude.

Figure 3

Figure 3. Relative amplitude of the circadian IPP concentration and flux oscillations as a function of the number of daylight hours (h). X-axis – number of daylight hours. Y-axis – minimum value/maximum value of a variable during a day. (A, B) Relative amplitude of IPP_cyt concentration, (C) Relative amplitude of the flux going through the MVA pathway to produce IPP, (D) Relative amplitude of the flux going through the MEP pathway to produce IPP and (E) Relative amplitude of the flux going through both pathways to produce IPP.

Overall, when plants experience between 3 and 21 daylight hours, the relative amplitude of IPP and DMAPP daily concentration oscillations remains roughly consistent (Figure 3A). The fluctuation in pathway-specific IPP production remains close to 100% regardless of daytime length (Figures 3C, D). However, overall production fluctuations can be significantly less than 100% if the daytime length exceeds 4h (Figure 3E), attributed to the anti-phasic regulation of pathway expression. Essentially, while flux through one pathway nearly ceases at night, flux through the other pathway almost stops during daylight hours. When considering metabolite replenishment through isomerization and compartment exchange, production fluctuation is at most around 25%, even under short days (Supplementary Figure S3).

3.3 The role of IPP-DMAPP exchange between the cytosol and the plastid

As mentioned in the introduction, while IPP and DMAPP are exchanged between the plastid and cytosol, there is little evidence that any other intermediate of either pathway also diffuses between compartments. Given that it is widely accepted that cytosolic and plastid IPP/DMAPP are used by the plant to synthesize quasi-orthogonal sets of more complex terpenoids, we were interested in understanding the effect of preventing the inter-compartmental diffusion of those metabolites. As such we removed that exchange from Model E, by setting diffusion rate parameters equal to zero. This significantly decreases the amplitude of the IPP/DMAPP concentration oscillations induced by the circadian rhythm, while maintaining the amplitude of the oscillations at the level of pathway intermediates (Supplementary Figure S4).

3.4 Contribution of the various regulatory modules towards the dynamic behavior of IPP and DMAPP biosynthesis during circadian light cycles

As described in section 2.7, three distinct regulatory modules connect the circadian rhythm to the regulation of MEP and MVA pathway activity. To assess the impact of these alternative regulatory modules, we investigate how the number of daylight hours per day influences the dynamics of IPP and DMAPP production in models B, C, D, BC, BD, and CD. We then replicate the analysis conducted for Model E.

3.4.1 Circadian regulation of pathway substrate availability

Model B focuses solely on the circadian regulation of pathway substrate concentrations. Under this framework, IPP and DMAPP levels remain relatively stable, showing resilience to changes in daylight hours (refer to Figures 3A, 4 and Supplementary Figure S5). Specifically, when exposed to more than 15h of daylight, IPP and DMAPP concentrations reach a quasi-steady state (as indicated by the green line in Figure 3). Oscillations become noticeable when daylight hours drop below 15, with the maximum relative amplitude of their concentration oscillations occurring at around 2h of daylight per day. Even under minimal daylight hours, slight oscillations still occur (as shown in Figure 4). However, the relative amplitude of concentration oscillations for pathway intermediates increases significantly compared to the fully regulated Model E (compare Supplementary Figures S2 and S5). A noteworthy difference between Model E and Model B lies in the relative amplitude of concentrations for intermediates of the MVA and MEP pathways. Contrary to Model E, in Model B, the relative amplitude of concentrations for MVA pathway intermediates is smaller than that for MEP pathway intermediates.

Figure 4

Figure 4. Model B Time course simulation of IPP and DMAPP concentrations throughout 48h at different latitudes and times of the year: equator/spring and fall equinoxes (dusk = 12h), middle latitudes (winter, dusk = 9h; summer, dusk = 15h) and near polar circle latitudes (winter, dusk = 3h; summer, dusk = 21h). T = 1h Green lines – MEP pathway. Magenta lines – MVA pathway.

Figure 3 illustrates the impact of varying daylight hours on the relative amplitude of IPP and DMAPP concentration oscillations. Days with more than 14h of daylight exhibit constant IPP and DMAPP concentrations. As daylight hours decrease, the effect of the circadian light rhythm on concentration oscillation amplitudes becomes similar to that observed in Model E (compare blue and green curves in Figure 3A). This shift results from significant decreases of up to 70% in metabolite concentrations during long nights compared to maximum daylight concentrations. It is also worth noting that Model B achieves nighttime IPP and DMAPP quasi-steady state concentrations in less than 1h while steadily accumulating these intermediates over daylight hours.

Circadian fluctuations in production via MVA and MEP pathways are observed on days shorter than 12h and 16h, respectively (refer to Figures 3C, D). Overall production of both pathways exhibits similar behavior (Figure 3E), and isomerization and compartment exchange dampen these oscillations.

3.4.2 Antithetic circadian regulation of MEP and MVA gene expression

Model C only incorporates circadian regulation of gene expression in the MEP and MVA pathways, with an antithetic regulation pattern between the two pathways, as documented in previous studies (Covington et al., 2008; Vranová et al., 2013; Atamian and Harmer, 2016). Broadly, the behavior of IPP and DMAPP oscillations closely parallels that of model E: when plants experience between 3 and 21 daylight hours, the relative amplitude of daily concentration oscillations for IPP and DMAPP remains relatively stable. However, this amplitude is approximately half of that observed in Model E (compare blue and orange curves in Figure 3A). Additionally, a comparison between Figures 2, 5 and Supplementary Figures S2, S6 reveals that the dynamic behavior of pathway intermediate concentrations is similar between Models E and C. Moreover, variations in the fraction of daylight hours per day minimally affect these concentrations. Furthermore, variations in the fraction of daylight hours per day have minimal impact on these concentrations. In contrast, the relative amplitude of oscillations in the concentrations of pathway products IPP and DMAPP remains high when there are between 4h and 20h of daylight. Figure 3 demonstrates that the influence of daylight hours on the relative amplitude of the oscillations qualitatively mirrors that observed in model E. However, the depletion of IPP and DMAPP is approximately 20% smaller in this model compared to model E (compare orange and blue curves in Figure 3A). Interestingly, the system also demonstrates rapid adaptation to transitions between light and darkness, requiring approximately two hours to reach either a daylight quasi-steady state or a nighttime quasi-steady state. Model C also shows the same qualitative behavior as Model E in terms of flux circadian regulation. Global influx of IPP and DMAPP displays oscillations that are similar to those for concentration. The flux of material going through either the MVA or MEP pathways is reduced almost 100%, regardless of daytime length (Figures 3C, D). The aggregated fluxes show that the pathways take turns in producing IPP, as seen in Model E for long days (Figure 3E). In Model C, and because pathway substrate availability does not depend on the circadian rhythm, this type of dynamic behavior extends to short days as well.

Figure 5

Figure 5. Model C. Time course simulation of IPP and DMAPP concentrations throughout 48h at different latitudes and times of the year: equator/spring and fall equinoxes (dusk = 12h), middle latitudes (winter, dusk = 9h; summer, dusk = 15h) and near polar circle latitudes (winter, dusk = 3h; summer, dusk = 21h). T = 1h. Green lines – MEP pathway. Magenta lines – MVA pathway.

3.4.3 Circadian regulation of IPP and DMAPP consumption

Model D considers a situation where the circadian rhythm only regulates the activity of pathways that use IPP and DMAPP to synthesize more complex terpenoids. The red curve in Figure 3A shows that, under these conditions, the amplitude of IPP and DMAPP concentration oscillations undergoes sharp changes when the number of daylight hours is lower than 4. If the number of daylight hours is between 4 and 23, the oscillation amplitude experiences a slight increase. As the number of daylight hours goes up to 24, the oscillation amplitude diminishes, and the system reaches a steady state. Similar to Model C, Supplementary Figure S7 illustrates that pathway intermediate concentrations remain relatively stable across a broad range of daylight hours. Figure 6 shows that IPP and DMAPP concentrations oscillate in a pattern akin to that of Model E, while pathway intermediate levels maintain stability with low-amplitude oscillations, especially for MEP pathway intermediates. Throughout the day, both compartments witness accumulation of IPP and DMAPP, followed by rapid depletion at night. In Model D, maximum relative depletions occur for days with 23h of light.

Figure 6

Figure 6. Model D. Time course simulation of IPP and DMAPP concentrations throughout 48h at different latitudes and times of the year: equator/spring and fall equinoxes (dusk = 12h), middle latitudes (winter, dusk = 9h; summer, dusk = 15h) and near polar circle latitudes (winter, dusk = 3h; summer, dusk = 21h). T = 1h. Green lines – MEP pathway. Magenta lines – MVA pathway.

The flux of material going through the MVA pathway is approximately constant regardless of daytime length (Figure 3C). The flux of material going through the MEP pathway is reduced, at most, by 10% (see the behavior in Figure 3D). The same behavior is observed for the sum of both fluxes, but with slightly weaker oscillations (Figure 3E). Maximum oscillation amplitude happens at daytime length around 18h.

3.4.4 Pairwise combination of circadian regulatory modules

To test whether there is a synergistic effect of circadian regulation between different modules, we also created Models BC, CD, and BD. These model emulate situations where circadian regulation is lost in only one of the regulatory modules. In Model BC circadian regulation of substrate availability is lost. In Model BD circadian regulation of gene expression is lost. In Model CD circadian regulation of IPP/DMAPP utilization is lost. We apply the same analysis as with the previous models (Figure 3E). In general, the relationship between the amplitude of concentration oscillations and the number of daylight hours in a model with two active circadian regulation modules resembles the combined dependencies of models where each individual module is the sole active circadian regulator (compare Figures 3A, B). In addition, the amplitude of concentration and flux oscillations is always smaller than that observed for the fully regulated Model E.

4 Discussion

4.1 Cellular demand for IPP/DMAPP and the MVA and MEP pathways

IPP and DMAPP are the precursor monomers to terpenoids, a family of molecules that contains many chemicals with importance in biology, pharmacy, biotechnology, biomedicine and cosmetics, such as squalene, cholesterol, some vitamins and most plant hormones. Plants produce those monomers using two biosynthetic pathways: the MVA pathway in the cytosol, and the MEP pathway in the plastid. IPP and DMAPP are used as the building blocks for more complex terpenoids, ranging from protective molecules such as carotenoids to hormones such as strigolactones. A version of the MVA pathway was present in the last common ancestor of archaea and eukaryotes, while ancestral bacteria contained a version of the MEP pathway (Lombard and Moreira, 2011; Hoshino and Gaucher, 2018; Zeng and Dehesh, 2021). In plants, the current MEP pathway seems to have evolved from the ancestral MEP pathway present in the early symbiotic cyanobacteria that became the chloroplast (Lichtenthaler, 1999). While both pathways have a level of crosstalk and IPP and DMAPP can be exchanged between the cytosol and the plastid, the contribution of each pathway to the biosynthesis of complex terpenoids is not the same. The MVA pathway mainly provides flux for the biosynthesis of sesquiterpenes, sterols, polyprenols, and triterpenes, while the MEP pathway is the main provider for the biosynthesis of chlorophylls, tocopherols, quinones, carotenoids, monoterpenes and strigolactones, among others (Pérez et al., 2022). Both pathways provide flux and materials for cellular metabolism. As such one would expect that their regulation would be mainly demand-driven, that is, its flux should be mainly regulated by the cellular demand for the material produced by the pathways.

4.2 Regulatory design for the MVA and MEP pathways is consistent with design principles for demand driven pathways

Several hallmarks of demand-driven pathways exist (Savageau, 1972; Alves and Savageau, 2000; Bromig et al., 2020). In terms of regulation, demand-driven pathways are more efficiently regulated by negative feedback (Savageau, 1972; Alves and Savageau, 2000). The most efficient negative feedback configuration is created by overall feedback, where the product of the pathway inhibits the flux of the first reaction (Savageau, 1972; Alves and Savageau, 2000, 2001; Ye and Medzhitov, 2019). This is clearly the case with the MEP pathway (Figure 1). Intriguingly, this type of regulation is absent from the MVA pathway. One can speculate why this is so. A probable reason for the absence of overall feedback in the MVA pathway is that the synthesis of the final products of the cytosolic MVA pathway, IPP and DMAPP, occurs in the peroxisome, not in the cytosol (Vranová et al., 2013). This prevents overall negative feedback from these products to the initial step of the pathway. As such, the cascading inhibitory feedback observed for the MVA pathway, where each enzyme is inhibited by its own product, still conveys information about cellular demand for the end product backwards through the biosynthetic chain one step at a time. This feedback structure creates a pathway whose response time to changes in cellular demand is slower than that permitted by the overall feedback configuration (Alves and Savageau, 2000), but still responsive to cellular demand. Another signature behavior for demand-driven pathways is that the concentration of pathway intermediates is lower than that of the pathway final products. This is observed when the steady state of Model E is calculated (Table 7). Stability (Table 8) and robustness of the steady state to changes in the parameter values (Tables 10 and 11) are other signature of demand-driven pathways (Savageau, 1972; Alves and Savageau, 2000). Both pathways have stable steady states when subjected to constant light conditions, and that steady state is robust to change in model parameters. Additionally, we find that IPP and DMAPP levels are less sensitive to restrictions on carbohydrate availability than the levels of pathway intermediates.

4.3 Circadian regulation of supply, gene expression and demand for the MVA and MEP pathways

The circadian light cycle regulates MVA and MEP pathway activity at three levels: pathway substrate availability, expression of genes coding for pathway enzymes, and activity of pathways that use IPP and DMAPP for the synthesis of more complex terpenoids. Figure 2 and Supplementary Figure S2 show that some of the steady state design principles described in section 4.2 are also observed when pathway dynamic behavior is away from the steady state and driven by the circadian light cycle. On the one hand, IPP and DMAPP concentrations are orders of magnitude higher than those for most pathway intermediates for both the MVA and the MEP pathways (Supplementary Figure S2). On the other hand, the metabolic oscillation itself is insensitive over a wide range to changes between shorter and longer days. For days with between 4h and 20h of light, the amplitude of the oscillation remains surprisingly stable (Figure 3).

Understanding how each of the three regulatory levels contributes to that stability drove us to create models of the pathways were circadian regulation was only active for one of the levels. When the light cycle only regulates pathway substrate availability, we find that the concentration waves for pathway intermediates are much more sensitive to changes in the number of daylight hours in the circadian rhythm than for the fully regulated model (Supplementary Figure S5). However, the concentration of the products IPP and DMAPP remains very insensitive to the supply of pathway substrate (Figure 4). In fact, if the number of daylight hours is above 15, the oscillation is lost, and the pathway operates at or near a steady state. This further strengthens the argument that the structure and parameters of these pathways have been selected by evolution to be consistent with a demand driven pathway. Our model shows that MVA-produced IPP and DMAPP oscillate with smaller relative amplitudes than those for the concentrations of the intermediate metabolites of the pathway. The relative amplitude of the IPP and DMAPP concentration oscillations is approximately half of that observed when the circadian rhythm is considered to regulate the three levels of pathway activity, while remaining equally insensitive to changes in the number of daylight hours.

When the light cycle only regulates gene expression, concentrations of pathway intermediates are very insensitive to changes in the number of daylight hours (Supplementary Figure S6). 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). When this is the only type of circadian regulation acting on the MVA and MEP pathways, our models indicate that IPP/DMAPP concentrations are significantly less robust to fluctuations in enzyme activity than the concentrations of pathway intermediates (Figure 5). As such, circadian regulation of gene expression alone, would create a regulatory structure for the pathway that would be suboptimal for regulating a demand driven pathway. We also found that having antithetic regulation of the gene expression between the MVA and the MEP pathways leads to oscillations that have smaller amplitudes, making the dynamic flux going through the pathways less sensitive to changes between day and night.

When the light cycle only regulates IPP and DMAPP utilization, concentrations of pathway intermediates remain insensitive to changes in the number of daylight hours (Supplementary Figure S7). The amplitude of the IPP and DMAPP concentration oscillations is similar to that observed when the circadian rhythm is considered to regulate the three levels of pathway activity (Figure 6). However, it is more sensitive to changes between longer and shorter days, making the concentrations and fluxes for the system very responsive to changes in cellular demand for IPP and DMAPP.

When the light cycle influences only two of the three regulatory modules in the pathway, an additive regulatory effect is observed. The system’s sensitivity to variations in daylight duration is intermediate between the fully regulated model and the models in which only a single module is affected by the circadian rhythm. Similarly, the amplitude of concentration oscillations falls between that of the fully regulated model and the models with regulation limited to a single module, whether it is substrate availability, gene expression, or product utilization.

Taken together, these results suggest that the three levels of circadian regulation, plus the MVA- and MEP-specific regulatory inhibition loops contribute differently to creating an operating regime that maintains pathway flux strongly coupled to demand and insensitive to changes over a wide range of daylight hours. Inhibitory feedback stabilizes the pathway product concentrations, when circadian rhythms change pathway supply availability, at the cost of amplifying concentration oscillations of pathway intermediates. While circadian regulation acts only on gene expression or on demand for IPP/DMAPP, the same inhibitory feedback creates product concentration oscillations with bigger amplitudes and decreases the amplitude of the oscillations for pathway intermediates. Thus, circadian regulation of gene expression or of demand for IPP/DMAPP alone would create a pathway whose dynamic response is suboptimal to demand for the pathways’ final products. It is the joining of the three regulatory modules that balances the dynamic behavior of the pathway, making it as robust as possible to the cellular demand for pathway products.

4.4 Limitations of this work

Here we present what we believe are the main limitations of this study and discuss how they could affect the robustness and applicability of its findings.

One important limitation of this study is the black box manner in which we model the production of substrate for the MVA and MEP pathways. This is a limitation that is shared with other modeling studies (Rios-Estepa et al., 2010; Pokhilko et al., 2015; Weaver et al., 2015; Basallo et al., 2023; Neiburga et al., 2023), and we believe that this valid simplification facilitates the analysis of the intrinsic dynamic and regulatory behavior of the MVA and MEP pathways.

It is also worthy of note that a direct comparison of our modeling results with the behavior or real circadian mutants is complex. The circadian clock regulates carbohydrate metabolism, secondary metabolism, and volatile production in plants through complex transcriptional networks and signaling pathways that are not yet fully understood. Mutations in circadian clock genes have pleiotropic effects on all three of these processes. To determine the independent impact of clock mutations on each module, it would be necessary to modify the promoters of the genes involved so that their activity becomes independent of circadian regulation. However, the circadian-dependent regulatory regions within these promoters have not yet been fully identified, which complicates this task, especially given the large number of genes involved. Additionally, accurately measuring changes in the intermediates of the MEP and MVA pathways remains technically challenging. Still, the characterization of circadian-dependent, gene specific, promoter regions for terpenoid metabolism genes is possible (Loivamäki et al., 2007) and it is predictable that such experiments could be undertaken in the near future.

Another limitation is the fact that we use a similar approach to model the consumption of IPP and DMAPP out of the two pathways for the biosynthesis of more complex terpenoids. The biosynthesis of these and other terpenoid final products is also treated in this paper as a black box that draws flux from the MEP and MVA pathways in a way that is dependent on the circadian rhythm. This is valid for many volatile terpenoids, as well as carotenoids and other phytohormones (Loivamäki et al., 2007; Covington et al., 2008; Zeng et al., 2017; Zheng et al., 2017; Picazo-Aragonés et al., 2020; Mu et al., 2022). Still, to more accurately understand how the biosynthesis of these more complex terpenoids affects the dynamics of the MEP and MVA pathways, additional research is needed. This research should develop, analyze, and integrate detailed models of the biosynthetic pathways for those terpenoids with the MVA and MEP pathway models. That development and integration should take into account that certain terpenoids are synthesized using material that is drawn mainly from only one of the pathways (Rodríguez-Concepción et al., 2004; Yang et al., 2012; Qiao et al., 2021; Wu et al., 2021; Chandrasekaran et al., 2022; Liu et al., 2024). MVA-derived isoprenoid end products in plants are sterols and cytokinins that modulate membrane architecture, plant growth and development, and brassinosteroids that work as steroid hormones. In contrast, MEP-derived end products include photosynthesis-related isoprenoids (carotenoids and the side chains of chlorophylls, plastoquinones, and phylloquinones), gibberellins and abscisic acid hormones, and root volatile monoterpenes.

An additional limitation of this work is that the sensitivity analysis we performed, while comprehensive in identifying parameters to which the model is highly sensitive, is differential. This approach may not fully capture the dynamic complexities of biological systems and may miss higher order interactions between simultaneous changes in more than one parameter. We also note that, while the study provides insights into how different regulatory modules influence the effect of the circadian rhythm and latitude on IPP and DMAPP production, it does not fully address the potential interactions between these modules and other environmental or physiological factors that might affect the circadian rhythm. The models used to simulate these regulatory effects (Models B, C, D, and their combinations) are simplified representations and may not capture all the nuances of circadian regulation in real biological systems. For example, the antithetic regulation of gene expression in Model C and the regulation of consumption pathways in Model D may not fully represent the intricate feedback mechanisms present in vivo.

4.5 Conclusions

Several conclusions can be drawn from this work. Because the dynamic behavior of our model is robust, this allows us to conclude that its use to simulate physiological situations is likely to be appropriate. Our analysis also concludes that the feedback inhibition of enzymes by pathway intermediates and end products is compatible with a situation where the regulation of the flux going through those pathways is significantly driven by the cellular demand for their end product. Finally, we can also conclude that the three regulatory modules at which the circadian rhythm affects IPP and DMAPP production interact to make that production less sensitive to the seasonal changes in the number of daylight hours observed at different latitudes on our planet. Finally, the antithetic regulation in gene expression contributes to buffer the global production of IPP and DMAPP against shifts between day and night.

Thus, even with very limited quantitative information available, mathematical models can point at relevant features of these pathways and propose scenarios for future experimental exploration that may facilitate the modification of IPP and DMAPP production through synthetic biology efforts.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Author contributions

OB: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Conceptualization. AL: Writing – review & editing, Methodology, Formal analysis. AS: Writing – review & editing, Formal analysis. AM-S: Writing – review & editing, Formal analysis. EV: Writing – review & editing, Formal analysis. EM: Writing – review & editing, Formal analysis. AE: Writing – review & editing, Formal analysis. RA: Writing – review & editing, Writing – original draft, Visualization, Supervision, Software, Methodology, Investigation, Formal analysis, Conceptualization.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. PROSTRIG, an ERANET project from FACEJPI (PCI2019-103382, MICIUN), partially funded this project. This project ended in 2022. AL received Ph. D. fellowship funding from the European Union’s H2020 research and innovation programme under Marie Skłodowska-Curie grant agreement No. 801586. OB received a Ph. D. fellowship from AGAUR (2022FI_B 00395). AL and OB both received 1-year Ph. D. extension fellowships from IRBlleida.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2024.1465030/full#supplementary-material

Supplementary Figure 1 | Modeling circadian light cycles with the L function. (A) Approximate daylighthours at the equator. (B) Approximate daylighthours at the polar circle during the peak of winter. (C) Approximate daylighthours at the polar circle during the peak of summer.

Supplementary Figure 2 | Time course simulation of the system throughout 48h at different latitudes and times of the year: equator/spring and fall equinoxes (dusk = 12h), middle latitudes (winter, dusk = 9h; summer, dusk = 15h) and near polar circle latitudes (winter, dusk = 3h; summer, dusk = 21h). T = 1h.

Supplementary Figure 3 | Model E. Oscillation amplitude of IPP and DMAPP production (normalized to the maxima) for different values of T (Twilight) and different values of dusk (Daytime).

Supplementary Figure 4 | Model E without diffusion of IPP/DMAPP between compartments. Time course simulation of the system throughout 48h at different latitudes and times of the year: equator/spring and fall equinoxes (dusk = 12h), middle latitudes (winter, dusk = 9h; summer, dusk = 15h) and near polar circle latitudes (winter, dusk = 3h; summer, dusk = 21h). T = 1h.

Supplementary Figure 5 | Model B. Time course simulation of the system throughout 48h at different latitudes and times of the year: equator/spring and fall equinoxes (dusk = 12h), middle latitudes (winter, dusk = 9h; summer, dusk = 15h) and near polar circle latitudes (winter, dusk = 3h; summer, dusk = 21h). T = 1h.

Supplementary Figure 6 | Model C. Time course simulation of the system throughout 48h at different latitudes and times of the year: equator/spring and fall equinoxes (dusk = 12h), middle latitudes (winter, dusk = 9h; summer, dusk = 15h) and near polar circle latitudes (winter, dusk = 3h; summer, dusk = 21h). T = 1h.

Supplementary Figure 7 | Model D. Time course simulation of the system throughout 48h at different latitudes and times of the year: equator/spring and fall equinoxes (dusk = 12h), middle latitudes (winter, dusk = 9h; summer, dusk = 15h) and near polar circle latitudes (winter, dusk = 3h; summer, dusk = 21h). T = 1h.

References

Alabadí, D., Oyama, T., Yanovsky, M. J., Harmon, F. G., Más, P., Kay, S. A. (2001). Reciprocal regulation between TOC1 and LHY/CCA1 within the Arabidopsis circadian clock. Science 293, 880–883. doi: 10.1126/SCIENCE.1061320/SUPPL_FILE/1061320S1_THUMB.GIF

PubMed Abstract | Crossref Full Text | Google Scholar

Albe, K. R., Butler, M. H., Wr1ghtt, B. E. (1990). Cellular concentrations of enzymes and their substrates. J Theor Biol. 143 (2), 163–95. doi: 10.1016/s0022-5193(05)80266-8

PubMed Abstract | Crossref Full Text | Google Scholar

Alves, R., Savageau, M. A. (2000). Effect of overall feedback inhibition in unbranched biosynthetic pathways. Biophys. J. 79, 2290–2304. doi: 10.1016/S0006-3495(00)76475-7

PubMed Abstract | Crossref Full Text | Google Scholar

Alves, R., Savageau, M. A. (2001). Irreversibility in unbranched pathways: preferred positions based on regulatory considerations. Biophys. J. 80, 1174–1185. doi: 10.1016/S0006-3495(01)76094-8

PubMed Abstract | Crossref Full Text | Google Scholar

Alves, R., Vilaprinyo, E., Hernádez-Bermejo, B., Sorribas, A. (2008). Mathematical formalisms based on approximated kinetic representations for modeling genetic and metabolic pathways. Biotechnol. Genet. Eng. Rev. 25, 1–40. doi: 10.5661/bger-25-1

PubMed Abstract | Crossref Full Text | Google Scholar

Atamian, H. S., Harmer, S. L. (2016). Circadian regulation of hormone signaling and plant physiology. Plant Mol. Biol. 91, 691–702. doi: 10.1007/S11103-016-0477-4

PubMed Abstract | Crossref Full Text | Google Scholar

Baadhe, R. R., Mekala, N. K., Palagiri, S. R., Parcha, S. R. (2012). Development of petri net-based dynamic model for improved production of farnesyl pyrophosphate by integrating mevalonate and methylerythritol phosphate pathways in yeast. Appl. Biochem. Biotechnol. 167, 1172–1182. doi: 10.1007/S12010-012-9583-1/FIGURES/9

PubMed Abstract | Crossref Full Text | Google Scholar

Basallo, O., Perez, L., Lucido, A., Sorribas, A., Marin-Saguino, A., Vilaprinyo, E., et al. (2023). Changing biosynthesis of terpenoid percursors in rice through synthetic biology. Front. Plant Sci. 14. doi: 10.3389/FPLS.2023.1133299/BIBTEX

PubMed Abstract | Crossref Full Text | Google Scholar

Bick, J. A., Lange, B. M. (2003). Metabolic cross talk between cytosolic and plastidial pathways of isoprenoid biosynthesis: unidirectional transport of intermediates across the chloroplast envelope membrane. Arch. Biochem. Biophys. 415, 146–154. doi: 10.1016/S0003-9861(03)00233-9

PubMed Abstract | Crossref Full Text | Google Scholar

BRENDA Enzyme Database (2024). (Accessed January 25, 2024).

Google Scholar

Bromig, L., Kremling, A., Marin-Sanguino, A. (2020). Understanding biochemical design principles with ensembles of canonical non-linear models. PloS One 15, e0230599. doi: 10.1371/journal.pone.0230599

PubMed Abstract | Crossref Full Text | Google Scholar

Buchanan, B. B., Gruissem, W., Jones, R. L. (2002). Biochemistry & molecular biology of plants. Eds. Buchanan, B. B., Gruissem, W., Jones, R. L. (Hoboken, New Jersey: American Society of Plant Physiologists).

Google Scholar

Chandrasekaran, U., Byeon, S., Kim, K., Kim, S. H., Park, C. O., Han, A. R., et al. (2022). Short-term severe drought influences root volatile biosynthesis in eastern white pine (Pinus strobus L). Front. Plant Sci. 13. doi: 10.3389/fpls.2022.1030140

PubMed Abstract | Crossref Full Text | Google Scholar

Cockburn, W., McAulay, A. (1977). Changes in metabolite levels in kalanchoë daigremontiana and the regulation of Malic acid accumulation in crassulacean acid metabolism. Plant Physiol. 59, 455–458. doi: 10.1104/PP.59.3.455

PubMed Abstract | Crossref Full Text | Google Scholar

Cordoba, E., Salmi, M., León, P. (2009). Unravelling the regulatory mechanisms that modulate the MEP pathway in higher plants. J. Exp. Bot. 60, 2933–2943. doi: 10.1093/JXB/ERP190

PubMed Abstract | Crossref Full Text | Google Scholar

Covington, M. F., Maloof, J. N., Straume, M., Kay, S. A., Harmer, S. L. (2008). Global transcriptome analysis reveals circadian regulation of key pathways in plant growth and development. Genome Biol. 9, 1–18. doi: 10.1186/GB-2008-9-8-R130/FIGURES/6

PubMed Abstract | Crossref Full Text | Google Scholar

Dudareva, N., Andersson, S., Orlova, I., Gatto, N., Reichelt, M., Rhodes, D., et al. (2005). The nonmevalonate pathway supports both monoterpene and sesquiterpene formation in snapdragon flowers. Proc. Natl. Acad. Sci. U. S. A. 102, 933–938. doi: 10.1073/PNAS.0407360102/ASSET/F5C94EBD-80FE-4A87-A155-E7CC4E980CE3/ASSETS/GRAPHIC/ZPQ0030569670005.JPEG

PubMed Abstract | Crossref Full Text | Google Scholar

Eisenreich, W., Rohdich, F., Bacher, A. (2001). Deoxyxylulose phosphate pathway to terpenoids. Trends Plant Sci. 6, 78–84. doi: 10.1016/S1360-1385(00)01812-4

PubMed Abstract | Crossref Full Text | Google Scholar

Flores-Pérez, Ú., Sauret-Güeto, S., Gas, E., Jarvis, P., Rodríguez-Concepcióna, M. (2008). A mutant impaired in the production of plastome-encoded proteins uncovers a mechanism for the homeostasis of isoprenoid biosynthetic enzymes in arabidopsis plastids. Plant Cell 20, 1303–1315. doi: 10.1105/TPC.108.058768

PubMed Abstract | Crossref Full Text | Google Scholar

Guevara-García, A., San Román, C., Arroyo, A., Cortés, M. E., de la Gutiérrez-Nava, M. L., León, P. (2005). Characterization of the arabidopsis clb6 mutant illustrates the importance of posttranscriptional regulation of the methyl-d-erythritol 4-phosphate pathway. Plant Cell 17, 628–643. doi: 10.1105/TPC.104.028860

PubMed Abstract | Crossref Full Text | Google Scholar

Han, M., Heppel, S. C., Su, T., Bogs, J., Zu, Y., An, Z., et al. (2013). Enzyme inhibitor studies reveal complex control of methyl-D-erythritol 4-phosphate (MEP) pathway enzyme expression in Catharanthus roseus. PloS One 8(5), e62467. doi: 10.1371/JOURNAL.PONE.0062467

PubMed Abstract | Crossref Full Text | Google Scholar

Harborne, J. B., Tomas-Barberan, F. A., Phytochemical Society of Europe (1991). Ecological chemistry and biochemistry of plant terpenoids (Oxford: Clarendon Press).

Google Scholar

Hemmerlin, A. (2013). Post-translational events and modifications regulating plant enzymes involved in isoprenoid precursor biosynthesis. Plant Sci., 203–204. doi: 10.1016/J.PLANTSCI.2012.12.008

PubMed Abstract | Crossref Full Text | Google Scholar

Hemmerlin, A., Harwood, J. L., Bach, T. J. (2012). A raison d’être for two distinct pathways in the early steps of plant isoprenoid biosynthesis? Prog. Lipid Res. 51, 95–148. doi: 10.1016/J.PLIPRES.2011.12.001

PubMed Abstract | Crossref Full Text | Google Scholar

Hemmerlin, A., Hoeffler, J. F., Meyer, O., Tritsch, D., Kagan, I. A., Grosdemange-Billiard, C., et al. (2003). Cross-talk between the cytosolic mevalonate and the plastidial methylerythritol phosphate pathways in tobacco bright yellow-2 cells. J. Biol. Chem. 278, 26666–26676. doi: 10.1074/jbc.M302526200

PubMed Abstract | Crossref Full Text | Google Scholar

Hemmerlin, A., Tritsch, D., Hartmann, M., Pacaud, K., Hoeffler, J. F., Van Dorsselaer, A., et al. (2006). A Cytosolic Arabidopsis d-Xylulose Kinase Catalyzes the Phosphorylation of 1-Deoxy-d-Xylulose into a Precursor of the Plastidial Isoprenoid Pathway. Plant Physiol. 142, 441–457. doi: 10.1104/PP.106.086652

PubMed Abstract | Crossref Full Text | Google Scholar

Hoshino, Y., Gaucher, E. A. (2018). On the origin of isoprenoid biosynthesis. Mol. Biol. Evol. 35, 2185–2197. doi: 10.1093/molbev/msy120

PubMed Abstract | Crossref Full Text | Google Scholar

Jin, X., Baysal, C., Drapal, M., Sheng, Y., Huang, X., He, W., et al. (2021). The coordinated upregulated expression of genes involved in MEP, chlorophyll, carotenoid and tocopherol pathways, mirrored the corresponding metabolite contents in rice leaves during de-etiolation. Plants 10, 1456. doi: 10.3390/plants10071456

PubMed Abstract | Crossref Full Text | Google Scholar

Kitano, H. (2007). Towards a theory of biological robustness. Mol. Syst. Biol. 3, 137. doi: 10.1038/MSB4100179

PubMed Abstract | Crossref Full Text | Google Scholar

Lange, I., Poirier, B. C., Herron, B. K., Lange, B. M. (2015). Comprehensive assessment of transcriptional regulation facilitates metabolic engineering of isoprenoid accumulation in arabidopsis. Plant Physiol. 169, 1595–1606. doi: 10.1104/PP.15.00573

PubMed Abstract | Crossref Full Text | Google Scholar

Laule, O., Fürholz, A., Chang, H.-S., Zhu, T., Wang, X., Heifetz, P. B., et al. (2003). Crosstalk between cytosolic and plastidial pathways of isoprenoid biosynthesis in Arabidopsis thaliana. Proc. Natl. Acad. Sci. U. S. A. 100, 6866–6871. doi: 10.1073/pnas.1031755100

PubMed Abstract | Crossref Full Text | Google Scholar

Lee, C. P., Elsässer, M., Fuchs, P., Fenske, R., Schwarzländer, M., Millar, A. H. (2021). The versatility of plant organic acid metabolism in leaves is underpinned by mitochondrial malate–citrate exchange. Plant Cell 33, 3700. doi: 10.1093/PLCELL/KOAB223

PubMed Abstract | Crossref Full Text | Google Scholar

Liao, P., Hemmerlin, A., Bach, T. J., Chye, M.-L. (2016). The potential of the mevalonate pathway for enhanced isoprenoid production. Biotechnol. Adv. 34, 697–713. doi: 10.1016/J.BIOTECHADV.2016.03.005

PubMed Abstract | Crossref Full Text | Google Scholar

Lichtenthaler, H. K. (1999). THE 1-DEOXY-D-XYLULOSE-5-PHOSPHATE PATHWAY OF ISOPRENOID BIOSYNTHESIS IN PLANTS. Annu. Rev. Plant Physiol. Plant Mol. Biol. 50, 47–65. doi: 10.1146/annurev.arplant.50.1.47

PubMed Abstract | Crossref Full Text | Google Scholar

Liu, J., Yin, X., Kou, C., Thimmappa, R., Hua, X., Xue, Z. (2024). Classification, biosynthesis, and biological functions of triterpene esters in plants. Plant Commun. 5, 100845. doi: 10.1016/j.xplc.2024.100845

PubMed Abstract | Crossref Full Text | Google Scholar

Loivamäki, M., Louis, S., Cinege, G., Zimmer, I., Fischbach, R. J., Schnitzler, J. P. (2007). Circadian rhythms of isoprene biosynthesis in grey poplar leaves. Plant Physiol. 143, 540–551. doi: 10.1104/PP.106.092759

PubMed Abstract | Crossref Full Text | Google Scholar

Lombard, J., Moreira, D. (2011). Origins and early evolution of the mevalonate pathway of isoprenoid biosynthesis in the three domains of life. Mol. Biol. Evol. 28, 87–99. doi: 10.1093/molbev/msq177

PubMed Abstract | Crossref Full Text | Google Scholar

McClung, C. R. (2013). Beyond Arabidopsis: The circadian clock in non-model plant species. Semin. Cell Dev. Biol. 24, 430–436. doi: 10.1016/J.SEMCDB.2013.02.007

PubMed Abstract | Crossref Full Text | Google Scholar

Mcgarvey, D. J., Croteau, R. (1995). Terpenoid metabolism. Plant Cell 7, 1015–1026. doi: 10.1105/tpc.7.7.1015

PubMed Abstract | Crossref Full Text | Google Scholar

Mu, Z., Llusià, J., Zeng, J., Zhang, Y., Asensio, D., Yang, K., et al. (2022). An overview of the isoprenoid emissions from tropical plant species. Front. Plant Sci. 13. doi: 10.3389/FPLS.2022.833030/BIBTEX

PubMed Abstract | Crossref Full Text | Google Scholar

Nagel, D. H., Kay, S. A. (2012). Complexity in the wiring and regulation of plant circadian networks. Curr. Biol. 22, R648–R657. doi: 10.1016/J.CUB.2012.07.025

PubMed Abstract | Crossref Full Text | Google Scholar

Neiburga, K. D., Muiznieks, R., Zake, D. M., Pentjuss, A., Komasilovs, V., Rohwer, J., et al. (2023). Total optimization potential (TOP) approach based constrained design of isoprene and cis-abienol production in A. thaliana. Biochem. Eng. J. 190, 108723. doi: 10.1016/J.BEJ.2022.108723

Crossref Full Text | Google Scholar

Page, J. E., Hause, G., Raschke, M., Gao, W., Schmidt, J., Zenk, M. H., et al. (2004). Functional analysis of the final steps of the 1-deoxy-d-xylulose 5-phosphate (DXP) pathway to isoprenoids in plants using virus-induced gene silencing. Plant Physiol. 134, 1401–1413. doi: 10.1104/PP.103.038133

PubMed Abstract | Crossref Full Text | Google Scholar

Pérez, L., Alves, R., Perez-Fons, L., Albacete, A., Farré, G., Soto, E., et al. (2022). Multilevel interactions between native and ectopic isoprenoid pathways affect global metabolism in rice. Transgenic Res. 31, 249–268. doi: 10.1007/s11248-022-00299-6

PubMed Abstract | Crossref Full Text | Google Scholar

Picazo-Aragonés, J., Terrab, A., Balao, F. (2020). Plant volatile organic compounds evolution: transcriptional regulation, epigenetics and polyploidy. Int. J. Mol. Sci. 21, 8956. doi: 10.3390/IJMS21238956

PubMed Abstract | Crossref Full Text | Google Scholar

Pokhilko, A., Bou-Torrent, J., Pulido, P., Rodríguez-Concepción, M., Ebenhöh, O. (2015). Mathematical modelling of the diurnal regulation of the MEP pathway in Arabidopsis. New Phytol. 206, 1075–1085. doi: 10.1111/NPH.13258

PubMed Abstract | Crossref Full Text | Google Scholar

Pokhilko, A., Flis, A., Sulpice, R., Stitt, M., Ebenhöh, O. (2014). Adjustment of carbon fluxes to light conditions regulates the daily turnover of starch in plants: a computational model. Mol. Biosyst. 10, 613–627. doi: 10.1039/C3MB70459A

PubMed Abstract | Crossref Full Text | Google Scholar

Qiao, Z., Hu, H., Shi, S., Yuan, X., Yan, B., Chen, L. (2021). An update on the function, biosynthesis and regulation of floral volatile terpenoids. Horticulturae 7, 451. doi: 10.3390/horticulturae7110451

Crossref Full Text | Google Scholar

Rios-Estepa, R., Lange, I., Lee, J. M., Markus Lange, B. (2010). Mathematical modeling-guided evaluation of biochemical, developmental, environmental, and genotypic determinants of essential oil composition and yield in peppermint leaves. Plant Physiol. 152, 2105–2119. doi: 10.1104/PP.109.152256

PubMed Abstract | Crossref Full Text | Google Scholar

Rodríguez-Concepción, M., Forés, O., Martínez-García, J. F., González, V., Phillips, M. A., Ferrer, A., et al. (2004). Distinct light-mediated pathways regulate the biosynthesis and exchange of isoprenoid precursors during arabidopsis seedling development. Plant Cell 16, 144–156. doi: 10.1105/tpc.016204

PubMed Abstract | Crossref Full Text | Google Scholar

Sauret-Güeto, S., Botella-Pavía, P., Flores-Pérez, U., Martínez-García, J. F., San Román, C., León, P., et al. (2006). Plastid cues posttranscriptionally regulate the accumulation of key enzymes of the methylerythritol phosphate pathway in Arabidopsis. Plant Physiol. 141, 75–84. doi: 10.1104/PP.106.079855

PubMed Abstract | Crossref Full Text | Google Scholar

Savageau, M. A. (1972). “The behavior of intact biochemical control systems,” in Current topics in cellular regulation. Eds. Horecker, B. L., Stadtman, E. R. (Cambridge, Massachusetts: Academic Press), 63–130. doi: 10.1016/B978-0-12-152806-5.50010-2

Crossref Full Text | Google Scholar

Savageau, M. A. (1976). Biochemical systems analysis : a study of function and design in molecular biology. (HArlow, England: Longman Higher Education), 379.

Google Scholar

Schaller, H., Grausem, B., Benveniste, P., Chye, M. L., Tan, C. T., Song, Y. H., et al. (1995). Expression of the hevea brasiliensis (H.B.K.) mull. Arg. 3-hydroxy-3-methylglutaryl-coenzyme A reductase 1 in tobacco results in sterol overproduction. Plant Physiol. 109, 761–770. doi: 10.1104/pp.109.3.761

PubMed Abstract | Crossref Full Text | Google Scholar

Singh, V. K., Ghosh, I. (2013). Methylerythritol phosphate pathway to isoprenoids: Kinetic modeling and in silico enzyme inhibitions in Plasmodium falciparum. FEBS Lett. 587, 2806–2817. doi: 10.1016/J.FEBSLET.2013.06.024

PubMed Abstract | Crossref Full Text | Google Scholar

Sorribas, A., Hernández-Bermejo, B., Vilaprinyo, E., Alves, R. (2007). Cooperativity and saturation in biochemical networks: a saturable formalism using Taylor series approximations. Biotechnol. Bioeng. 97, 1259–1277. doi: 10.1002/bit.21316

PubMed Abstract | Crossref Full Text | Google Scholar

Taiz, L., Zeiger, E., Møller, I. M., Ian, M., Murphy, A. S. (2014). Plant physiology and development. 6th ed. Eds. Taiz, L., Zeiger, E., Møller, I. M., Murphy, A. S. (Sunderland, Massachusetts: Sinauer), 20.

Google Scholar

Tetali, S. D. (2019). Terpenes and isoprenoids: a wealth of compounds for global use. Planta 249, 1–8. doi: 10.1007/S00425-018-3056-X/FIGURES/2

PubMed Abstract | Crossref Full Text | Google Scholar

Voit, E. O. (1991). Canonical nonlinear modeling : S-system approach to understanding complexity. (NY, NY: Van Nostrand Reinhold), 365.

Google Scholar

Voit, E. O. (2013). Biochemical systems theory: A review. ISRN Biomath. 2013, 1–53. doi: 10.1155/2013/897658

Crossref Full Text | Google Scholar

Vranová, E., Coman, D., Gruissem, W. (2013). Network analysis of the MVA and MEP pathways for isoprenoid synthesis. Annu Rev Plant Biol. 64, 665–700. doi: 10.1146/ANNUREV-ARPLANT-050312-120116

PubMed Abstract | Crossref Full Text | Google Scholar

Weaver, L. J., Sousa, M. M. L., Wang, G., Baidoo, E., Petzold, C. J., Keasling, J. D. (2015). A kinetic-based approach to understanding heterologous mevalonate pathway function in E. coli. Biotechnol. Bioeng. 112, 111–119. doi: 10.1002/BIT.25323

PubMed Abstract | Crossref Full Text | Google Scholar

Wright, L. P., Rohwer, J. M., Ghirardo, A., Hammerbacher, A., Ortiz-Alcaide, M., Raguschke, B., et al. (2014). Deoxyxylulose 5-phosphate synthase controls flux through the methylerythritol 4-phosphate pathway in arabidopsis. Plant Physiol. 165, 1488–1504. doi: 10.1104/PP.114.245191

PubMed Abstract | Crossref Full Text | Google Scholar

Wu, W., Du, K., Kang, X., Wei, H. (2021). The diverse roles of cytokinins in regulating leaf development. Hortic. Res. 8, 1–13. doi: 10.1038/s41438-021-00558-3

PubMed Abstract | Crossref Full Text | Google Scholar

Xie, Z., Kapteyn, J., Gang, D. R. (2008). A systems biology investigation of the MEP/terpenoid and shikimate/phenylpropanoid pathways points to multiple levels of metabolic control in sweet basil glandular trichomes. Plant J. 54, 349–361. doi: 10.1111/J.1365-313X.2008.03429.X

PubMed Abstract | Crossref Full Text | Google Scholar

Yang, D., Du, X., Liang, X., Han, R., Liang, Z., Liu, Y., et al. (2012). Different roles of the mevalonate and methylerythritol phosphate pathways in cell growth and tanshinone production of salvia miltiorrhiza hairy roots. PloS One 7, e46797. doi: 10.1371/journal.pone.0046797

PubMed Abstract | Crossref Full Text | Google Scholar

Ye, J., Medzhitov, R. (2019). Control strategies in systemic metabolism. Nat. Metab. 1, 947–957. doi: 10.1038/s42255-019-0118-8

PubMed Abstract | Crossref Full Text | Google Scholar

Yon, F., Kessler, D., Joo, Y., Cortés Llorca, L., Kim, S. G., Baldwin, I. T. (2017). Fitness consequences of altering floral circadian oscillations for Nicotiana attenuata. J. Integr. Plant Biol. 59, 180–189. doi: 10.1111/JIPB.12511/SUPPINFO

PubMed Abstract | Crossref Full Text | Google Scholar

Zeng, L., Dehesh, K. (2021). The eukaryotic MEP-pathway genes are evolutionarily conserved and originated from Chlaymidia and cyanobacteria. BMC Genomics 22, 137. doi: 10.1186/s12864-021-07448-x

PubMed Abstract | Crossref Full Text | Google Scholar

Zeng, L., Wang, X., Kang, M., Dong, F., Yang, Z. (2017). Regulation of the rhythmic emission of plant volatiles by the circadian clock. Int. J. Mol. Sci. 18, 2408. doi: 10.3390/IJMS18112408

PubMed Abstract | Crossref Full Text | Google Scholar

Zheng, R., Liu, C., Wang, Y., Luo, J., Zeng, X., Ding, H., et al. (2017). Expression of mep pathway genes and non-volatile sequestration are associated with circadian rhythm of dominant terpenoids emission in osmanthus fragrans lour. Flowers. Front. Plant Sci. 8. doi: 10.3389/FPLS.2017.01869/BIBTEX

PubMed Abstract | Crossref Full Text | Google Scholar

Zhou, F., Pichersky, E. (2020). More is better: the diversity of terpene metabolism in plants. Curr. Opin. Plant Biol. 55, 1–10. doi: 10.1016/J.PBI.2020.01.005

PubMed Abstract | Crossref Full Text | Google Scholar