Metabolic Flow Analysis

Introduction

Based on preliminary experiments, we hope to identify an appropriate gene knockout approach to create a Nicotiana benthamiana plant chassis with low chlorogenic acid. Thus, we plan 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 will 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 intend to use the GSMM for Nicotiana tabacum as a starting point and adapt it to develop a model for Nicotiana benthamiana.

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

1.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.

2.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.

Substitute 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} $$

$$\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.

3.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 can 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.

table

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 Caffeoylmalic 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 is in a steady state and predict 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 techniques

We found that after knocking out the HQT pathway, chlorogenic acid cannot be synthesized at all, allowing the accumulation of more metabolic precursors in the model, which aids 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 is a highly conserved pathway. Moreover, the plant cannot survive normally after the complete knockout of HQT. Therefore, we consider partial knockout HQT by reducing the flux. After revising the model, we obtained the following results.

Fig.3 Chlorogenic acid production under revised knockout techniques

We observed that in cases where HQT is partially knocked out or both HQT and C3'H metabolic pathways are knocked out, the accumulation of chlorogenic acid is relatively low. Among these scenarios, knocking out both HQT and C3'H pathways results 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 does not significantly reduce chlorogenic acid content. Therefore, partial knockout of HQT might be a better experimental target.

Shikimate pathway

(a) Caffeoylmalic acid production simulation

Due to the large accumulation of chlorogenic acid in organisms and its crucial role in the shikimic acid 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 Caffeoylmalic acid by adding Caffeoyl-CoA and Malate as substrates in the model.

(b) Resveratrol production simulation

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.

Aiming for efficient synthesis of resveratrol, we selected the VvSTS gene (AAT84992, 151bp) derived from grapes and introduced overexpression mutants YlARO3K225L or YlARO3D154N to inhibit the binding of l-phe with DAHP synthase.

Non-shikimate pathway -- Crocin production simulation

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 crocetin, 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 crocetin yield under different knockout conditions to understand the effects on its production.

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 reduces 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 caffeoylquinic 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 predict the accumulation curves of caffeoylmalic acid, resveratrol, and crocetin over time in Versatobacco, we have constructed dynamic reaction bounds and ordinary differential equations based on the FBA results. The simulations of the accumulation curves 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.

Shikimate Pathway

After constructing the gene knockout models, we introduced two synthetic pathways for caffeoylmalic 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:

$$ \text{proc.ub} = \begin{cases} \begin{array}{c} \frac{BF}{1 + e^{-r (t - t_0)}}, & \quad t < THR \\ \frac{BF}{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} $$

$$ inhibit = \frac{C}{k \cdot PBF} \tag{8} $$

Assumptions

• Assuming that t=0 marks the beginning of the plant growth phase, there will be an initial accumulation of caffeoylmalic acid and resveratrol, set at 0.2 and 0.2, respectively.
• Since the transfection method is transient, a significant accumulation of metabolites is expected between 48-96 hours.

Variables

Here is a table!!!

Parameters

Here is a selectable table!!!

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 group model, we found that the amplification parameter \( k \) plays a major role in controlling the final steady-state concentration. When \( k \) is within a certain range, the results are closer to the wet-lab experimental data (fig.8A). Parameter \( r \), which controls the curve growth rate, can induce a significant S-shaped growth trend when \( r \) 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

We used fixed parameters (see Parameters) to simulate all gene knockout models. It was observed that the growth rate and final accumulation of caffeoylmalic acid and resveratrol accumulation curves are positively correlated with the corresponding balanced fluxes in FBA, and both exhibit S-shaped metabolite growth following the start of the plant growth phase (fig.9). We mainly focus on the product accumulation curve of HQT knockout plants compared to the control group, and this visualization result will be presented in the final results analysis.

Non-shikimate Pathway

To explore the application potential of Versatobacco in alternative pathways, we have inserted the crocetin 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 described earlier. 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

• Assume that t=0 marks the beginning of the plant growth phase, so there will be an initial accumulation of Crocetin, including Crocin I and Crocin II set as 0.1, 0.1.
• Given that the transfection method is transient, there will likely be a significant accumulation of metabolites between 48-96 hours.

Variables

Here is a table!!!

Parameters

Here is a selectable table!!!

Result

In the synthesis of crocin, ########

Results Analysis

To validate the potential of Versatobacco in producing beneficial biochemical compounds, we concentrate on the synthesis of caffeoylmalic acid and resveratrol via the shikimate pathway, as well as the synthesis of crocetin through the non-shikimate pathway in HQT-knockout plants.

The results show that HQT-knockout plants exhibit a significantly enhanced ability to synthesize target metabolites through the shikimate pathway, with a substantial increase in product accumulation compared to the control group (fig.12). This is mainly because the original HQT step diverted a large amount of Caffeoyl-CoA and p-Coumaroyl-CoA. By simulating the knockout and limiting its flux to a minimal range, we can achieve greater accumulation of metabolic coenzymes, leading to the synthesis of more desired products.

What's more, in the synthesis of crocin, ########