Metabolic Flow Analysis

Introduction

Based on preliminary experiments, we hoped to identify an appropriate gene knockout approach to create a Nicotiana benthamiana plant chassis with low chlorogenic acid. Thus, we planed to use Constraint-Based Reconstruction and Analysis (COBRA) to analyze and compare the effects of different gene knockdown methods on chlorogenic acid reduction and evaluate their impact on other metabolic pathways. This would provide valuable guidance for subsequent wet-lab experiments.

GSMM Adjustments

Currently, there is no genome-scale metabolic model (GSMM) available for Nicotiana benthamiana. However, a GSMM exists for Nicotiana tabacum. Given the differences between these species and their distinct cultivation conditions, we intended to use the GSMM for Nicotiana tabacum as a starting point and adapted it to develop a model for Nicotiana benthamiana.

To enhance the accuracy of the model, we adjusted growth conditions, such as controlling photosynthesis and CO₂ absorption, and refined the optimization targets accordingly.

Photosynthesis

We assumed that only the leaf performs photosynthetic carbon reactions. Using literature data (Payyavula, Raja S et al, 2015), we obtained an average photosynthetic CO₂ amount of \(20 \mu mol/m^2/s\). The average leaf area is \(20,000 cm^2\), the fresh leaf weight is \(1700g\), and the dry weight is \(10\%\) of the fresh weight. We converted these units to \(mmol/gDW/h\) to standardize with the model units.

$$\small \begin{align} v_{photosynthesis} &=\frac{20\ \mu mol/m^2/s}{1000}\cdot 3600\cdot\frac{2\ m^2}{1700\ gFW\cdot 10\%}\\ &\approx 0.847mmol/gDW/h\end{align} \tag{1} $$

Finally, we used 0.80-0.88 as the model parameter.

CO₂ absorption

(a) Adjust CO₂ absorption to match photosynthesis

Using the bisection method, we adjusted the CO₂ absorption to bring the photosynthesis rate close to the target value of \(0.8 mmol/gDW/h\). Ultimately, we obtained a CO₂ absorption rate of \(39.7 mmol/gDW/h\), resulting in a photosynthesis flux of \(0.8155 mmol/gDW/h\). We considered this result to be preliminarily reasonable.

(b) Formula calculation verification

We assumed standard atmospheric pressure and uniform gas diffusion, using Fick's law to roughly estimate the upper limit of the reaction flux.

$$ v=A \cdot D \cdot \frac{\Delta C}{d} \tag{2} $$

$$ \Delta C=\frac{C_e-C_i}{V_m} \tag{3} $$

where \(A\) is the leaf area, \(D\) is the CO₂ gas diffusion coefficient, \(\Delta C\) is the concentration difference, \(V_m\) is the molar volume of the gas, \(d\) is the leaf thickness, \(C_e\) is the extracellular CO₂ concentration, and \(C_i\) is the intercellular CO₂ concentration. The dry weight is \(10\%\) of the fresh weight.

Substituted the data from the literature (Payyavula, Raja S et al, 2015) into the calculation:

$$\small \begin{align} v_{CO_2-absorption^{'}}&=2m^2\cdot 2\cdot{10}^{-5}m^2/s\cdot\frac{\left(300-210\right)\cdot{10}^{-6}mol}{22.4\cdot{10}^{-3}m^3}\cdot\frac{1}{0.2\cdot{10}^{-3}m}\\ &\approx 0.804\ mol/s \end{align} \tag{4} $$

Then converted this value into the corresponding units in the model, namely \(mmol/gDW/h\):

$$\small \begin{align} v_{CO_2-absorption}&=v_{CO_2-absorption^{'}}\cdot{10}^3\cdot 3600\cdot \frac{1}{1700gFW\cdot 10\%}\\ &\approx 34.03mmol/gDW/h \end{align} \tag{5} $$

Considering deviations in actual conditions, this value is generally consistent with the fitted value. Thus, we set the CO₂ absorption limit to 30-39.

Biomass

Considering the metabolic and accumulation processes of biomass during plant growth, we believe that targeting a single pathway for growth simulation is not reasonable. Instead, using a series of linear combinations as the biomass growth target could better simulate the actual growth process. Based on the data from the article (Yu, Jing et al, 2023), we constructed a linear combination of metabolic reaction accumulates and used this as the optimization objective for our model. The specific data are as follows.

Metabolite Production Simulation

In the previous section, we developed a tobacco GSMM model that fits the actual situation well. Next, we need to simulate and analyze several major products to explore the effects of different gene knockout methods on their accumulation. Based on an examination of metabolic pathways and preliminary experiments, we considered chlorogenic acid as a byproduct of the product metabolism and aimed to reduce its accumulation to achieve a higher flux of target product synthesis.

Here, we will analyze the synthesis of chlorogenic acid as well as phaselic acid, resveratrol, and crocin.

Methods

We employed Flux Balance Analysis (FBA), which relies on the steady-state assumption of the metabolic network where the production and consumption rates of all metabolites are balanced over time. FBA involves constructing a stoichiometric matrix and applying constraints to reaction fluxes. Through linear programming, FBA optimizes objective functions such as the organism's growth rate or the yield of metabolic products, thereby predicting the optimal flux distribution within the metabolic network (Heirendt Laurent et al, 2019).

Fig.1 The schematic of FBA

We assumed that the metabolic network was in a steady state and predicted the distribution of metabolic fluxes by optimizing a linear combination designed based on the specific accumulation ratios of biological components. This allowed us to predict and compare the accumulation rates of target products.

Chlorogenic acid production simulation

We observed the shikimate pathway in the model and compared the differences in chlorogenic acid metabolism under various gene knockout conditions. In the control group, chlorogenic acid, being a byproduct, accumulated in larger amounts, which is detrimental to the synthesis and accumulation of other secondary metabolites. By employing gene knockout techniques, we compared the accumulation of chlorogenic acid under the knockout of several key metabolic pathways. The results are shown below.

Fig.2 Chlorogenic acid production under different knockout strategies

We found that after knocking out the HQT pathway, chlorogenic acid could not be synthesized at all, allowing the accumulation of more metabolic precursors in the model, which aided in the synthesis of target metabolites.

However, during the experiments, we found that genes related to the HQT pathway appeared multiple times, indicating that this was a highly conserved pathway. Moreover, the plant could not survive normally after the complete knockout of HQT. Therefore, we considered partial knockout HQT by reducing the flux. After revising the model, we obtained the following results.

Fig.3 Chlorogenic acid production under revised knockout strategies

We observed that in cases where HQT was partially knocked out or both HQT and C3'H metabolic pathways were knocked out, the accumulation of chlorogenic acid was relatively low. Among these scenarios, knocking out both HQT and C3'H pathways resulted in extremely low chlorogenic acid content, which may significantly affect the normal growth of the plant. Additionally, in cases of partial HQT knockout, further knockout of other genes did not significantly reduce chlorogenic acid content. Therefore, partial knockout of HQT might be a better experimental target.

Target compound production simulation

(a) Addition of Missing Reactions

Phaselic acid

Due to the large accumulation of chlorogenic acid in organisms and its crucial role in the shikimate pathway, it directly affects the synthesis of related secondary metabolites. A decrease in its content can lead to the accumulation of other metabolites related to the upstream substance, Caffeoyl-CoA. Therefore, we chose to introduce the HCT-M gene to increase the yield of phaselic acid by adding Caffeoyl-CoA and Malate as substrates in the model.

Resveratrol

In the shikimate pathway, another upstream substance of chlorogenic acid, p-Coumaroyl quinic acid, can be generated from p-Coumaroyl-CoA, which is also an upstream substance for many valuable metabolites such as lignin and anthocyanins. Therefore, we chose p-Coumaroyl-CoA as a reactant to synthesize resveratrol, which has a promising market outlook, as another target product.

To achieve efficient synthesis of resveratrol, we intended to introduce the VvSTS gene and the YIARO3 gene. Consequently, we incorporated relevant reactions into the model to accurately simulate the actual conditions.

Crocin

In addition to investigating the shikimic acid pathway, it is crucial to assess the impact of gene knockouts within the non-shikimic acid pathway. Accordingly, we selected crocin, a high-value compound in the MEP pathway, for detailed analysis.

We incorporated four pertinent genes, GjCCD4a, GjALDH2C3, GjUGT94E13, and GjUGT74F8, into our model and added the associated missing reactions. We then evaluated the crocin yield under different knockout conditions to understand the effects on its production.

(b)Fitting Model with Wet Lab Data

We found that partial knockout of HQT may be a promising option in our preliminary analysis. To further optimize the model, we utilized the LC-MS/MS results of key compounds from the mutant (HQT knockout) and wild-type plants obtained from wet experiments, adjusting the bounds of relevant reactions in our model to fit the experimental data. The specific data are shown in the table below.

Table 2: Data from Wet Lab

Results

Through flux balance analysis, we simulated the yield of each substance under different knockout scenarios.

Based on our results, we found that, compared to other knockout methods, partial knockout of the HQT gene significantly reduced the accumulation of the by-product chlorogenic acid while preserving normal plant life functions. This approach also redirects metabolic flux, leading to increased production of phaselic acid and resveratrol from the shikimic acid pathway, as well as crocin from the non-shikimic acid pathway. Therefore, targeting the HQT gene appears to be the most effective strategy for constructing an optimized plant chassis.

Dynamic Flux Balance Analysis

To further simulate the accumulation of phaselic acid, resveratrol, and crocin over time in Versatobacco, we constructed dynamic reaction bounds and ordinary differential equations based on the results from FBA and wet lab. The simulations for these three substances not only serve as a validation for the successful construction of the chassis, but will also be used to identify the optimal timing for achieving the best economic benefits, supporting our ideas both in the wet lab and human practice!

Specifically, here we concentrate on the synthesis of phaselic acid and resveratrol via the shikimate pathway, as well as the synthesis of crocin through the non-shikimate pathway in HQT-knockout plants.

Shikimate Pathway

After constructing the gene knockout models, we introduced two synthetic pathways for phaselic acid and resveratrol at the Caffeoyl-CoA and p-Coumaroyl-CoA steps, respectively. We also inserted a reaction for the net output of metabolites secreted from intracellular to extracellular compartments, which we name “proc” (Fig.7).

The upper-bound flux of the proc is dynamically constrained, which first increases with plant growth, then gradually decreases after reaching the maximum production rate:

$$ proc.ub = \begin{cases} \begin{array}{c} \frac{PBF}{1 + e^{-r (t - t_0)}}, & \quad t < THR \\ \frac{PBF}{1 + e^{-r (THR - t_0)}} - d \cdot (t - THR), & \quad t \geq THR \end{array} \end{cases} \tag{6} $$

Then, we constructed ordinary differential equations (ODEs) for the accumulation concentrations of the metabolites:

$$ \frac{dC}{dt} = D2W \cdot PF \cdot (1 - inhibit) \tag{7} $$

Since plants maintain a stable internal environment, the inhibition rate of the generating reactions will increase as the production of metabolites accumulates. Eventually, we will observe that the products accumulate slowly within plants, similar to a linear change:

$$ inhibit = \min(\frac{C}{k \cdot PBF}, i) \tag{8} $$

Assumptions

• Assuming that t=0 marks the beginning of the plant growth phase, there will be an initial accumulation of phaselic acid and resveratrol, set at 0.18, 0.14 in control and 0.48, 0.22 in HQT-knockout models, respectively.
• Since the transfection method is transient, a significant accumulation of metabolites is expected between 48-96 hours.
• The analysis subject is consistent with FBA, targeting 1700gFW of leaves (approximately 2m²).

Variables

Parameters

Senesitivity Analysis

In the above model, sensitivity analysis is required for parameters 1-4 to validate the model's robustness and determine suitable ranges.

Here, using the control model, we found that the amplification parameter \( k \) plays a major role in controlling the final concentration, which should be used to fit the wet-lab experimental results (Fig.8A). Parameter \( r \), controlling the curve growth rate, can induce a significant S-shaped growth trend when it exceeds 0.4 (Fig.8B). Parameters \( t_0 \) and \( r \) together control the timing of the maximum growth rate (Fig.8C). The decline parameter for the proc reaction flux, \( d \), also inhibits the final accumulation of metabolite concentrations, and its inhibitory effect is related to the threshold setting of \( THR \) (Fig.8D).

Result

Then we tried to simulate the accumulation curves of control and HQT-knockout models. It was observed that the growth rate and final accumulation of phaselic acid and resveratrol accumulation curves are positively correlated with the corresponding balanced fluxes through FBA, and both exhibit S-shaped metabolite growth following the start of the plant growth phase (Fig.9). What's more, the HQT-knockout models exhibited a stronger ability to synthesize metabolites compared to the control models, with phaselic acid being 2.5 to 3.7 times higher and resveratrol 1.6 to 1.8 times higher. The dynamic changes further confirmed the significance of knocking out HQT in the shikimate pathway.

Non-shikimate Pathway

To explore the application potential of Versatobacco in alternative pathways, we have inserted the crocin biosynthesis pathway and are now constructing a local dynamic system for dFBA simulation. Similarly, we have named the reaction for the net output of metabolites secreted from intracellular to extracellular compartments as “proc”, and the reaction where Crocin II converts to Crocin I as “R_15” (Fig.10).

Here, we discuss two substances: Crocin I and Crocin II. The mathematical models for the dynamic upper bounds of proc and metabolite accumulation are consistent with those in shikimate pathway. Our focus is primarily on dynamically constraining the upper bound of the R_15 reaction:

$$ R15.ub = \max \left( \frac{\phi_1 \cdot C_2}{\phi_2 + C_2}, \phi_3 \right) \tag{9} $$

Assumptions

• Assuming that t=0 marks the beginning of the plant growth phase, there will be an initial accumulation of crocin. Crocin I and Crocin II are the primary final products, with a ratio of approximately six to one (Xie, 2023), set at 0.12, 0.02 in control and 0.16, 0.03 in HQT-knockout models.
• Given that the transfection method is transient, there will likely be a significant accumulation of metabolites between 48-96 hours.
• The analysis subject is consistent with FBA, targeting 1700gFW of leaves (approximately 2m²).

Variables

Parameters

Result

Compared to the control group, HQT-knockout plants exhibit better synthesis of crocin. However, the enhancement is not as significant as that seen with phaselic acid and resveratrol. Therefore, based on this simulation result, we believe that Versatobacco also has certain potential in the synthesis of metabolites through the exogenous pathway, but the benefits are not as substantial as those of the direct products from the shikimate pathway.

Results Analysis

Here, we present a comparison graph of metabolite synthesis between the control group and HQT-knockout plants.

In terms of why HQT-knockout plants exhibit a significantly enhanced ability to synthesize target metabolites through the shikimate pathway, we think that the original HQT step diverted a large amount of Caffeoyl-CoA and p-Coumaroyl-CoA. That means, by limiting the flux of HQT step to a negligible range close to zero, we can achieve greater accumulation of metabolic coenzymes, leading to the synthesis of more desired products. And our simulation results support all our ideas!

Then again, HQT-knockout models also show advantages in the synthesis of crocin, illustrating that MEP pathway is related to the pathway we discussed before. However, the mechanism is not clear through dry lab results, which guides wet lab to further explore it.