Loading...
Dry Lab

Genome Scale Metabolic Model

Overview

Genome Scale Metabolic Models (GEMs) are computational representations of the metabolic networks of organisms including bacteria, archaea, and eukaryotes. They are composed of known metabolic reactions and metabolites from databases, and comprehensive annotated genome information of the target organism, in a way that enables mathematical modelling of metabolic processes (Passi et al., 2022). Genome Scale Metabolic Models often serve to predict metabolic fluxes by Flux Balance Analysis, as well as in silico metabolic engineering. By performing analyses using GEMs, researchers may identify crucial reactions and associated genes which may be engineered to optimize the production of the target substance. Our team aims to achieve high yields of Docosahexaenoic Acid through bioproduction using Y. Lipolytica, an oleaginous yeast widely used in production of oily substances, such as oleate and biodiesel. However, since the DHA production pathway is introduced into our chassis and has quite a heavy resource demand, there are several uncertainties that we must solve to optimize the yield of the polyunsaturated fatty acid:

  1. How will the addition of the pathway affect the original metabolism of Y. Lipolytica?
  2. How could we direct the metabolic flux towards DHA, while ensuring survival?
  3. What is the optimal process of production?

Using existing Genome Scale Metabolic Models of Y. Lipolytica with our own modifications, we attempt to address the questions and provide a direction for further wet lab experimental optimization of DHA yield.

Selection of Model

There are currently seven existing GEMs for Y. Lipolytica. Amongst these models, we believe that iYli21(Guo et al., 2022), the newest model, would be best suited for our purpose:

    1. High Accuracy:

    • iYli21 outperforms other models with an 85.7% accuracy in predictions regarding nutrient utilization

    2. Strain Suitability:

    • Our chosen strain of Y. Lipolytica, po1f, is a derivative of strain W29, the reference genome of iYli21. All previous models were built upon the genome of strain CLIB122.

3. Log & Modifications

To better prepare the model for our use, we made some modifications and notes to the existing model iYLI21. We created a log to track our modifications and list reactions and metabolites that are particularly relevant to production of DHA, such as fatty acid metabolism reactions.

    A. Deleting Duplications
    • i. There were initially several duplicated metabolites associated with different reactions. We identified the metabolites and replaced the duplicate:
      Metabolite Duplicate Name Compartment
      m60 m2037 S-adenosyl-L-methionine_C15H22N6O5S mi
      m200 m1855 malonyl-CoA_C24H38N7O19P3S cy
      m222 m2017 hydrogen peroxide_H2O2 pe
      m255 m2038 S-adenosyl-L-homocysteine_C14H20N6O5S cy
      m772 m2043 L-tryptophan_C11H12N2O2 cy/td>
      m878 m1963 N-acetyl-alpha-D-glucosamine 1-phosphate_C8H16NO9P cy
      m1954 m1959 3,4-dihydroxyphenylethyleneglycol cy
    B. Deleting Unnecessary Reactions
    • There are more than one biomass reactions in iYli21, the additional reactions were unnecessary, so we removed them:
      Reaction ID Name
      R1372 yeast 5 biomass pseudoreaction
      R1387 yeast 6 biomass pseudoreaction
      R1710 yeast 8 biomass pseudoreaction
      newBiom newBiom
    • C. Note relevant information
      • i. What are the compartments?
        • 1. According to NNU-iGEM-2022, the main compartments considered by our team are abbreviated as follows:
          Abbreviation Compartment
          C_cy Cytosol
          C_ex Extracellular
          C_mi Mitochondria
          C_pe Peroxisome
        ii. What are the lower/upper bounds?
        Name (Variable) Value (mmol/gDW/h)
        FB1N1000 -1000
        FB2N2 -2.43
        FB3N0 0
        FB4N8 7.8625
        FB5N1000 1000
      iii. Relevant metabolites and reactions:
      No. Name Compartment
      68 acetyl-CoA Cytosol
      200 malonyl CoA Cytosol
      40 NADPH Cytosol
      41 NADP+ Cytosol
      32 H2O Cytosol
    iv. iYli21 has several reactions related to lipid metabolism. Of which, we identified xPOOL_FA_LP as the main lipid accumulation reaction.

4. Flux Balance Analysis (Growth & Lipid Accumulation)

To gain an initial understanding of our chassis, we evaluated its growth and lipid flux:

    A. Biomass
  1. i. The biomass flux (growth rate) in iYli21 was experimentally validated in the original paper at certain uptake rates of glucose:
    Uptake Rate (mmol/gDW/h) Predicted Growth Rates (h⁻¹) Rates Determined Experimentally (h⁻¹)
    -0.61 0.031 0.047
    -0.64 0.036 0.048
    -2.43 0.28 0.26
  2. ii. Exchange reactions simulates substrate exchange with the external environment in iYli21, which includes essential nutrients such as carbon sources and nitrogen sources. However, quantification of nutrient uptake in the model does not depend on the concentration, but the uptake rates.
  3. iii. Through literature mining, with seven glycerol transporters and two channels, glycerol was identified to be more preferred by Y. lipolytica than glucose. We found that the uptake rate of glycerol in wild type Y. lipolytica may be 0.32g/gDW/h, multiplying by the molar mass of glycerol (92.09g/mol) and diving by 1000, it converts into approximately 3.47 mmol/gDW/h (Erian et al., 2022).
  4. iv. We set the lower bound of glucose exchange to zero and to -3.47 for glycerol exchange and observed the predicted growth rates. However, results were less ideal than the original settings (-2.43 mmol/gDW/h; Glucose), differing from the expected increase of biomass flux.
  5. Carbon Source Reaction ID Uptake Rate (mmol/gDW/h) Predicted Growth Rates (h-1)
    Glucose R1070 -2.43 0.2883
    Glycerol R1141 -3.47 0.2250

Addition of PUFA Synthase Pathway Reactions :

Since Y. lipolytica does not produce DHA, our DHA production occurs by an enzyme complex known as the PUFA synthase from microalgae Schizochytrium sp. 31. We did some literature mining to gain an understanding of the mechanism, then we wrote simplified balanced equations to represent each iterative cycle, consisting of elongation and sometimes reduction to C=C bonds (see DHA Yield Model: Pathway Derivation for more details).

B. Reactions Added:

See ‘GEM Reactions.pdf’

5. Flux Balance Analysis (DHA)

    A. What changed?
  • i. In the previous section, we created a list of reactions to observe. After the addition of PUFA synthase pathway reactions, we are interested in the effect these reactions had on those reactions, to obtain a brief idea of how DHA production would affect the cells’ metabolism.
  • ii. We ran the model without setting an objective function to observe any changes in selected reactions.
    Reaction ID Reaction Flux of Reaction Before Flux of Reaction After
    Biomass_C Biomass growth 0 0
    R1070 Glucose Exchange -2.43 -2.43
    R1016 Ammonium Sulfate Exchange -1.91 -1.948
    R2004 Acetyl CoA Production from Malonyl CoA 3.68E-25 2.967
    R2119 Acetyl CoA Production 0.2461 0
    R2121 Main Producing Reaction of Malonyl CoA and NADPH 1.03E-18 2.967
    xPOOL_FA_LP Fatty Acids Pool -0.001171 -0.001194

    After adding reactions, more nitrogen source was consumed, more acetyl-CoA, Malonyl CoA and NADPH were produced. Accumulation of fatty acids did not occur, and consumption of fatty acid increased (due to its negative value). However, the increased flux of acetyl CoA was from malonyl CoA instead of its original production reaction.

    Note: Acetyl CoA, Malonyl CoA and NADPH in compartment cytosol were targeted

  • iii. We ran the model setting the objective as biomass and DHA exchange separately to observe any changes in selected reactions.

    Nitrogen source uptake was increased by 67%, and less accumulated lipids were consumed. However, acetyl CoA, malonyl CoA, and NADPH production decreased. We ran further simulations to find that the reaction R2121 contributed to 98.50% of NADPH production before and only 41.80% after.

    B. DHA, DPAw6, and Biomass
  • i. DHA vs DPAw6
    • We set the DHA exchange reaction (R2313, Cytoplasm to Extracellular) as the objective function and optimized the model. Though we do not expect an efflux of DHA, this reaction must be set and targeted in iYli21 to achieve a theoretical yield. The resulting flux is 0.02773 mmol/gDW/h. Under this condition, we observed that the flux of biomass and DPAw6 are zero.

      We then set DPAw6 as the exchange reaction (R2314, Cytoplasm to Extracellular) as the objective function and optimized the model. The resulting flux is 0.02606 mmol/gDW/h. The biomass and DHA flux under this condition are again both zero.

      Objective Reaction ID Objective Reaction Flux of DHA exchange Flux of DPA exchange Biomass Flux
      R2313 DHA exchange 0.02773 0 0
      R2314 DPA exchange 0 0.02606 0
  • ii. Biomass vs DHA
  • 1. From the previous simulation, maximizing our introduced reactions lead to a biomass flux of 0. Further analysis was done to address the relationship between biomass and DHA production.

    Dynamic flux balance analysis (dFBA) was performed to evaluate the difference in biomass when biomass is the objective and when DHA is the objective.

    Figure: Comparison between dFBA data with different objective functions, fitted to a logistic curve and with r-squared evaluation.

    Figure: Comparison between dFBA data with different objective functions, fitted to a logistic curve and with r-squared evaluation.

    Data points after nutrients used up are replicates of the final data point, an assumption is made that death phase does not occur in this graph.

    The growth of Y. lipolytica is significantly slower and uses up nutrients at a much lower cell dry weight when DHA production is prioritized.

  • 2. The difference in DHA production with different set biomass fluxes was also evaluated, beginning from 0 and incrementing by 0.1 until there is no solution. Figure: Maximized flux of DHA and DPAw6 when biomass is set to a certain flux value.

    Figure: Maximized flux of DHA and DPAw6 when biomass is set to a certain flux value.

    Results shows that higher biomass lower bounds lead to less DHA production.

  • 3. Results from flux balance analysis (FBA; Figure 2.) and dynamic FBA (Figure 1.) above suggests that to optimize the production of DHA, growth of Y. lipolytica may be negatively affected, which is not ideal for industrial processes. Keeping cell growth rate low may be beneficial to increasing DHA production. Also considering the fact that lipid accumulation mainly occurs during stationary phase, stage control separating growth and DHA production may be ideal.

iii. DHA & DPAw6 quadratic objective:

An attempt was made to optimize both DHA and DPA exchange reactions at the same time by assigning a python dictionary containing {Reactions: Coefficients} to the objective. However, this produced the same result as simply setting the objective as DHA exchange. Manipulating the coefficients produced results that were identical to setting the objective as whichever has the higher coefficient, but when they had the same coefficient, DHA exchange was dominant. This suggests that the DHA production pathway may be dominant over DPA production.

Future Modifications

While iYli21 contains thousands of reactions and is arguably the current best genome scale metabolic model for Y. lipolytica, lots of further modifications are required for it to be customized to model the production of DHA. Below are future steps

    A. Evaluation of the glycerol metabolism pathways and determination of suitable uptake rates
  • i. Widely mentioned in literature, Y. lipolytica prefers glycerol over glucose. However, in iYli21, using glycerol as the sole carbon source produces results that contradict with experimental observations. In the future, more work would be done to evaluate the reactions associated with glycerol reactions, with the goal of producing results that agree with experimental observations.
  • B. Addition of pathways to model the peroxisomal degradation of DHA and DPA.
  • i. In cycle 1, we knock out peroxisomal gene PEX10 to inhibit degradation of DHA. In iYli21, this gene does not exist, so we were unable to perform gene deletion simulations. In the future, we hope to include degradation reactions for DHA and DPA, which is crucial to predict DHA production and accumulation in Y. lipolytica accurately.
  • C. Ratio of DHA : DPAw6
  • i. According to literature, the DHA to DPAw6 ratio produced by our chosen enzyme complex is 2.3 : 1. However, as discussed in 5.B.i. (DHA vs DPAw6), the model is currently not capable of producing non-zero DHA and DPAw6 flux values at the same time. We aim to improve the model such that it can depict the ratio accurately.
  • D. Prediction of Further Optimization Strategies by Upregulation/Downregulation of Reactions and Gene Expression
  • i. With more research and modifications, we would achieve a model with higher accuracy in predicting DHA production via the PUFA synthase. We hope to analyze and predict new directions for DHA yield optimization via various algorithms within the cobra toolbox.

8. METHOD & TOOLS

Software Documentation/Installation
MATLAB 2017a https://www.mathworks.com/academia/student-competitions/igem.html
Cobrapy https://cobrapy.readthedocs.io/en/latest/index.html
COBRA Toolbox https://opencobra.github.io/cobratoolbox/stable/installation.html
SBML Validator https://sbml.bioquant.uni-heidelberg.de/validator_servlet/
Matplotlib https://matplotlib.org/
Scipy https://scipy.org/
Pandas https://pandas.pydata.org/docs/index.html

Note: Check compatibility of solvers if using cobra in MATLAB

REFERENCES

Passi, A., Tibocha-Bonilla, J. D., Kumar, M., Tec-Campos, D., Zengler, K., & Zuniga, C. (2022). Genome-Scale Metabolic Modeling Enables In-Depth Understanding of Big Data. Metabolites, 12(1), 14. https://doi.org/10.3390/metabo12010014
Guo, Y., Su, L., Liu, Q., Zhu, Y., Dai, Z., & Wang, Q. (2022). Dissecting carbon metabolism of Yarrowia lipolytica type strain W29 using genome-scale metabolic modelling. Computational and Structural Biotechnology Journal, 20, 2503–2511. https://doi.org/10.1016/j.csbj.2022.05.018
Erian, A. M., Egermeier, M., Marx, H., & Sauer, M. (2022). Insights into the glycerol transport of Yarrowia lipolytica. Yeast, 39(5), 323–336. https://doi.org/10.1002/yea.3702
Hauvermale, A., Kuner, J., Rosenzweig, B., Guerra, D., Diltz, S., & Metz, J. G. (2006). Fatty acid production in Schizochytrium sp.: Involvement of a polyunsaturated fatty acid synthase and a type I fatty acid synthase. Lipids, 41(8), 739–747. https://doi.org/10.1007/s11745-006-5025-6

DHA Yield Model

Overview

Kinetic modelling is a common computational and mathematical approach used to investigate mechanisms of chemical reactions. To gain further insight into our DHA producing PUFA synthase, we did comprehensive research on the mechanism of the enzyme complex and derived simplified pathway reactions. Through building this model, we hope to provide a base model for future research on PUFA synthase, enabling in silico experiments to investigate deeper into the mechanisms and the behavior of the pathway.

Pathway Deviation

To model the production of DHA, we must first understand the mechanism of the PUFA Synthase. The PUFA synthase consists of several functional subdomains, and they work together in iterative steps to build the polyunsaturated fatty acids, which in our case are DHA and omega-6 DPA. We first investigated the subdomains and their respective functions.

Functional Domains

    Acyl-Carrier Protein (ACP) Domains:
    • a. There are nine copies of ACP domains in the PKS synthase. These domains are responsible for carrying the growing fatty acid chain during its synthesis. The number of ACP domains can influence the productivity of PUFAs, with more ACP units correlating with higher PUFA production.
    2. Malonyl-CoA: ACP Transacylase (MAT) Domain:
    • a. This domain is responsible for loading malonyl groups onto the ACP domains.
    3. Ketoacyl-ACP Synthase (KS) Domains:
    • a. These domains catalyze the condensation reaction between acyl-ACP and malonyl-ACP, which creates Œ≤-ketoacyl-ACP as an intermediate. This adds two carbons to the carbon chain. In the PKS synthase, there are two KS domains:
    4. KSA Domain:
  • a. It is located at the N-terminus of subunit A. Besides the KS activity, it may also involve in the decarboxylation of malonyl-ACP to form acetyl-ACP.
  • 5. KSB Domain:
    • a. It is located at the N-terminus of subunit B, adjacent to a Chain Length Factor (CLF) domain. Besides the KS activity, I also work with CLF domain to control the substrate specificity, which is critical for producing specific bioactive compounds.
    6. Ketoacyl-ACP Reductase (KR) Domain:
    • a. It is responsible for the reduction of the keto group to a hydroxyl group in the β-ketoacyl-ACP intermediates, after which the β-hydroxyacyl-ACP intermediates will be formed.
    7. Dehydratase (DH) Domains:
    • a. These domains catalyze the dehydration of the β-hydroxyacyl-ACP intermediates to form trans double bonds. There are two DH domains in the PKS synthase:
      • i. DHPKS Domain:
        • 1. Catalyzes the dehydration reactions in the cycles that will not add double bond to the fatty acid chains.
      • ii. DHFabA Domain:
        • 1. Catalyzes the dehydration reactions and the isomerization reactions in the cycles that will add double bond to the fatty acid chains.
    8. Enoyl-ACP Reductase (ER) Domains:
    • a. These domains are involved in the of trans-enoyl-ACP to form saturated acyl-ACP. There are two ER domains in the PKS synthase: ERB domain and ERC domain. ERB domain is more important for the synthesis of PUFAs while ERC domain is more related to synthesis of short-chain fatty acids.
    9. Acyltransferase (AT) Domain AT:
    • a. Responsible for the release of the completed PUFA from the ACP domain as a free fatty acid.

We then organized the order of steps conducted by specific subdomains per iterative cycle according to the following figure:

Figure 1. Diagram depicting the iterative cycles and steps specific subdomains catalyze of the PUFA synthase producing DHA and omega-6 DPA (Guo et al., 2022).

Figure 1. Diagram depicting the iterative cycles and steps specific subdomains catalyze of the PUFA synthase producing DHA and omega-6 DPA (Guo et al., 2022).

Steps Function Involved Domains Specific Notes
1 Loading of acetyl and malonyl groups MAT Loading of the malonyl group from Coenzyme A
ACP Carrier
KSA Decarboxylation of malonyl group
2 Condensation of malonyl-ACP and acetyl-ACP KSA/KSB
3 Reduction of ketone group KR
4 Dehydration to form double bond DHPKS Cycles that do not add double bonds
DHFabA Cycles that add double bonds
5 Adjustment of the structure ERB/C Reduction of the double bond
DHFabA 2,2 or 2,3-isomerization

Table 1. The functions of each subdomain, and the order in which they work together to complete a full iterative cycle. For steps 4 and 5, only one of the two listed domains would be involved each cycle, and if DHFabA was involved in step 4, it would also be the one involved in step 5.

Most PUFA synthases produce more than one major product. However, not much is known about the key within the iterative process of the enzyme complex that leads to the production of different molecules. The major products of our PUFA synthase, taken from Schizochytrium sp. 31 (ATCC88888) are DHA and omega-6 DPA. The two polyunsaturated fatty acids have the same carbon number but differ in the number and position of double bonds.

Based on such structural difference, we go through each iterative cycle to explore possible steps that would allow the formation of either omega 6 DPA and DHA.

Based on such structural difference, we go through each iterative cycle to explore possible steps that would allow the formation of either omega 6 DPA and DHA.

Types of Iterative cycles:

After formation of acetyl-ACP and malonyl-ACP, under the catalyzation of PKS synthase, the substrates will enter a cyclic process, where each cycle will add two carbon groups to the fatty acid chain. There are mainly three types of cycles.

First Type: Adds two carbon atoms without adding double bond, catalyzed by the KS, KR, DHPKS and EH domains. This type requires 2 NADPH, while the other two types only requires 1 NADPH.

Figure 2.  Diagram depicting the first type of cycle.

Figure 2. Diagram depicting the first type of cycle.

Second Type: Adds two carbon atoms and a cis-double bond at the nearest carbon to the carbonyl group. However, instead of DHPKS and EH domains, the last two steps are catalyzed by the DHFabA domain through 2,2-isomerization after dehydration.

Figure 3.  Diagram depicting the second type of cycle.

Figure 3. Diagram depicting the second type of cycle.

Third Type: Adds two carbon atoms and a cis-double bond at the second-nearest carbon to the carbonyl group, The last two steps are also catalyzed by the DHFabA domain but via 2,3-isomerization.

Figure 4.  Diagram depicting the third type of cycle.

Figure 4. Diagram depicting the third type of cycle.

In DHA and omega-6 DPA, the spacing between double bonds is three carbons. Therefore, to ensure the double bond are in their right place, when synthesizing DHA, the cycle sequence for DHA is the third type -> the second type -> the first type, repeated from the second cycle until completion. Omega-6 DPA lacks one double bond at the fourth carbon in comparison to DHA, and by evaluating the possible sequences of cycles that would lead to this result, we identify the branch point to be at the second cycle. This indicates that while the first cycle is the first type for both DHA and omega-6 DPA synthesis, the second cycle is the third type for DHA while omega-6 DPA’s is the first type.

When the synthesis of DHA and omega-6 DPA completes, CLF and AT domains catalyzes the release of the fatty acid chain from the ACP domain.

Pathway Simplification & Model Assumption

There are many steps within each cycle, and it would be a very complicated model if we were to write ordinary differential equations to model every reaction involved. Though discussions with our advisors and consulting Professor Henry, we decided to simplify the reactions into one single balanced reaction per cycle:

Reactions (DHA):

Reactions (omega-6 DPA):

Note:

    Cx:y-E
    • o x = no. of carbons
    • o y = no. of double bonds
    • o E = Enzyme Complex (including all subdomains)
    Acetyl-CoA binding and cycle 1 reaction is shared between DHA and omega-6 DPA pathways.

Simbiology And Matlab Model Construction

A. Model Assumptions:

  • Assume first order reactions
  • Assume that H+ and CO2 does not change on the timescale of other substrates
    d/dt = 0
  • Assume that reactions are irreversible
  • Assume only one activator for regulated genes
  • Assume that the molecule will not detach from the complex until the last step
  • Assume the enzyme can only handle one molecule at once
  • Assume that the enhanced production of NADPH from specific gene expression is only used for the PUFA synthesis
  • All the Acetyl-CoA produced in the model are consumed for the DHA production
  • Assume the enzyme only uses Acetyl-CoA from cytosol

B. Model Construction:

  • We constructed a model in MATLAB and in SimBiology:
  • Gene Expression Equations:
  • Pathway ODEs:
  • We set initial parameters and ran the model to ensure its function:

Figure 5: Model results with initial approximations; scattered points are experimental data from Chang et al., 2013. The predicted curve is of similar shape as the experimental data, suggesting that our model is functional.

PARAMETER PREDICTION

A. Parameters for our pathway have yet to be tested, and it is difficult to estimate the value based on currently limited knowledge of the mechanism. It is most ideal to test the values experimentally, which we may consider in the future. For an initial guess, we opted for machine learning methods to predict kcat values, which depicts the turnover rate of the enzymes when there is no substrate deficiency, for each simplified equation.

B. Justification for Using kcat is as Follows:

  • i. From our assumption, the intermediates in the catalytic cycles will not unbound from the enzyme complex until the final product is formed, and the product from one catalytic cycle will pass on directly and be the substrate for the next catalytic cycle. Thus, it is assumed that the enzyme complex is always saturated with substrates.

C. Additional Assumptions:

  • i. Since we predict kcat values for each step in each cycle, but simplified the steps into one reaction per cycle in the model, we took the lowest kcat value per cycle and assumed it to be the limiting reaction of the cycle.

D. Kcat Parameters

  • We used CATPRED, a deep learning framework for prediction of in vitro enzyme kinetic parameters, including kcat (Boorla et al, 2024).
  • See “predicted_kcat.pdf”

    Sensitivity Analysis

    A. After the construction of the model, we performed sensitivity analysis via the Simbiology model to investigate the effect of variation in input variables on outputs.

    Figure x: Sensitivity analysis of PUFA synthase pathway inputs and outputs

    B. From the results, NADPH is identified as an important metabolite to the DHA production pathway. This suggests that increasing the availability of NADPH may be a crucial strategy to optimizing the yield of DHA.

    References

    Guo, P., Dong, L., Wang, F., Chen, L., & Zhang, W. (2022). Deciphering and engineering the polyunsaturated fatty acid synthase pathway from eukaryotic microorganisms. Frontiers in Bioengineering and Biotechnology, 10. https://doi.org/10.3389/fbioe.2022.1052785

    Yeast Growth Model

    Overview

    The development of a batch reactor model aims to optimize the production of lipids in yeast by implementing a controlled expression strategy for a large enzyme. From Genome Scale Metabolic Model simulations and discussions with Prof. David Banfield, we recognized that constitutive expression of this enzyme may not be beneficial for yeast health and productivity. Instead, a regulated expression approach is preferred to separate growth and lipid accumulation to two independent stages. Our model divides the entire bioprocess into two distinct stages: a growth stage, where the yeast cells are cultivated under optimal conditions, and a production stage, where lipid production is maximized. This staged approach allows for better management of the yeast's metabolic pathways, ultimately enhancing the yield of lipids.

    Background

    In engineering cycle 0, we investigate how different concentrations of glucose and glycerol influence the growth of Yarrowia Lipolytica. Building upon this foundational understanding, our model aims to simulate how various substrate types and their concentrations affect cell growth and lipid accumulation. The development of this model is essential not only for refining our bioprocess but also for providing a predictive tool that can guide future experiments and scale-up processes.

    Impact of different carbon sources on growth:

    Glucose and glycerol are universal substrates commonly used for the growth of Y. Lipolytica. Research indicates that the yeast exhibits varying growth rates and biomass yields depending on the carbon source employed suggested by Workman, M., Holt, P., & Thykaer, J. in 2013. For instance, glucose has been identified as a preferred substrate, resulting in higher specific growth rates compared to glycerol. However, glucose also leads to increased biomass production. Therefore, we aim to balance these two parameters and investigate which substrate should be used to optimize growth.

    Factors affecting Lipid Accumulation:

    Y. lipolytica has been known for its remarkable ability to accumulate lipids founded by Robles-Rodríguez et al in 2018. According to Kuttiraja, M., et al. (2023), the minimum carbon-nitrogen (C/N) ratio for Y. lipolytica to initiate lipid accumulation is generally reported to be 20:1. Therefore, in our model design, the lipid accumulation would be initiated when the ratio of concentration of C to concentration of N reaches 20:1.

    Design

    Our initial approach involves developing a batch growth model based on ordinary differential equations, as described in (Robles-Rodríguez et al, 2018) which models a continuously stirred tank reactor system. We have removed the parameters related to inflow and outflow in order to simulate a batch model. The fundamental concept is that glucose and nitrogen are consumed and converted into functional biomass, which supports growth. Once the carbon-to-nitrogen (C/N) ratio reaches 20:1, lipid accumulation begins. Concurrently, the production of citric acid, a metabolic byproduct, is also initiated under these conditions. Citric acid is derived from the substrate and may compete with or inhibit lipid accumulation.

    Assumptions:

    1. It is an ideal batch reactor system.
    2. Lipid accumulation is triggered when the C/N ratio reaches 20:1.
    3. pH, temperature, agitation, dissolved oxygen remains constant; thus, it won't be the growth parameters
    4. The capacity of lipid accumulation would be 0.4 * functional biomass+ lipid biomass,
    5. Due to insufficient literature support, it is assumed that k3, YCIT/S and ki,3 are same even when using different carbon source.
    6. The model would stop operating when the lipid accumulation has reached the maximum capacity.
    7. The initial lipid biomass and citric acid concentration is assumed to be zero, as the yeast would not be in a metabolic overflow during inoculation.

    Parameters and Variables:

    Our initial input will be:

    Symbol Meaning Values (Glucose/Glycerol) Reference
    𝑆𝑖 Initial concentration of carbon source [g.L−1] 20, 40 N/A
    𝑆 Concentration of carbon source [g.L−1] N/A N/A
    𝑁𝑖 Initial concentration of nitrogen [g.L−1] 3.26 Our own media recipe
    𝑁 Concentration of nitrogen [g.L−1] N/A N/A
    𝑋𝐿𝑖 Initial concentration of lipid biomass [g.L−1] 0 Assumption
    𝑋𝐿 Concentration of lipid biomass [g.L−1] N/A N/A
    𝑋𝑓𝑖 Initial concentration of functional biomass [g.L−1] 0.001 Our culture
    𝑋𝑓 Concentration of functional biomass [g.L−1] N/A N/A
    𝐶𝐼𝑇𝑖 Initial concentration of citric acid [g.L−1] 0 Assumption
    𝐶𝐼𝑇 Concentration of citric acid [g.L−1] N/A N/A
    𝜋𝑚𝑎𝑥𝐶𝐼𝑇 Specific citric acid production rate [gCIT.(gXf.h)−1] 0.013 (Robles-Rodríguez et al, 2018)
    𝜋𝐶𝐼𝑇 Specific citric acid production rate [gCIT.(gXf.h)−1] N/A N/A
    𝜋𝑚𝑎𝑥𝐿𝑖𝑝𝑖𝑑 Maximum specific lipid production rate [gLIP.(gXfh)−1] 0.030 (Robles-Rodríguez et al, 2018)
    𝜋𝐿𝑖𝑝𝑖𝑑 Specific lipid production rate [gLIP.(gXfh)−1] N/A N/A
    𝜇 Growth rate [h−1] N/A (Robles-Rodríguez et al, 2018)
    𝜇𝑚𝑎𝑥 Maximum growth rate [h−1] 20g/L glucose: 0.2245
    40g/L glucose: 0.1840
    20g/L glycerol: 0.2410
    40g/L glucose: 0.2725
    Our culture
    𝑌𝑋/𝑆 Yield coefficient for functional biomass production with respect to carbon source [gXf.gS−1] Glucose: 0.603 / Glycerol: 0.166 (Robles-Rodríguez et al, 2018) / (Dobrowolski et al., 2023)
    𝑌𝐿𝐼𝑃/𝑆 Yield coefficient for lipid biomass production with respect to carbon source [gLIP.gS−1] Glucose: 0.276 / Glycerol: 0.3 (Robles-Rodríguez et al, 2018) / (Dobrowolski et al., 2023)
    𝑌𝐶𝐼𝑇/𝑆 Yield coefficient for citric acid production with respect to substrate [gCIT.gS−1] 0.88 (Robles-Rodríguez et al, 2018) / Assumption
    𝑌𝑋/𝑁 Yield coefficient for functional biomass production with respect to nitrogen [gXf.gN−1] 23.23 (Robles-Rodríguez et al, 2018)
    𝛾 Constant fraction of structural lipids to functional biomass [gXL.gXf−1] 0.06 (Robles-Rodríguez et al, 2018)
    𝑘1 Inhibition constant for nitrogen due to lipid accumulation due to growth [g.L−1] 9.24 (Robles-Rodríguez et al, 2018)
    𝑘2 Inhibition constant for lipid accumulation due to citric acid concentration [g.L−1] 8.3 (Robles-Rodríguez et al, 2018)
    𝑘3 Inhibition constant for citric acid concentration due to citric acid concentration [g.L−1] 40.44 (Robles-Rodríguez et al, 2018)
    𝐾𝑁 Saturation constants for nitrogen due to functional biomass growth [gN.L−1] 0.074 (Robles-Rodríguez et al, 2018)
    𝐾𝑆,1 Saturation constants for carbon source due to functional biomass growth [gS.L−1] 1.517 / 0.5 (Robles-Rodríguez et al, 2018) / Gao et al. (2016)
    𝐾𝑆,2 Saturation constants for carbon source due to lipid accumulation [gS.L−1] 0.978 / 0.7 (Robles-Rodríguez et al, 2018) / Li et al. (2019)
    𝐾𝑆,3 Saturation constants for carbon source due to citric acid production [gS.L−1] 11.26 (Robles-Rodríguez et al, 2018)
    𝑘𝑖,1 Carbon source inhibition constant for functional biomass growth [gS.L−1] 83.14 / 80 (Robles-Rodríguez et al, 2018) / (Kuttiraja et al., 2023; Dobrowolski et al., 2023)
    𝑘𝑖,2 Carbon source inhibition constant for lipid accumulation [gS.L−1] 35.2 / 112.5 (Robles-Rodríguez et al, 2018) / (Kuttiraja et al., 2023; Dobrowolski et al., 2023)
    𝑘𝑖,3 Carbon source inhibition constant for citric acid accumulation [gS.L−1] 28.8 (Robles-Rodríguez et al, 2018)
    𝑙𝑛𝑑𝑄 Boolean value to indicate if citric acid production is initiated 0: if S/N>20
    1: if S/N>20:1
    (Robles-Rodríguez et al, 2018)
    𝑙𝑛𝑑𝐿 Boolean value to indicate if lipid accumulation is initiated 0: if S/N>20
    1: if S/N>20:1
    (Robles-Rodríguez et al, 2018)

    The Ordinary Differential Equations (ODEs):

    Result:

    Here is the growth curve generated by our model

    In our bioreactor growth experiment, we achieved a final dry mass concentration of 2.646 g/L with 40 g/L glycerol, further validating the feasibility of our model.

    The results of our growth model, which align closely with the observed growth curves, indicate that glycerol serves as a superior carbon source for yeast compared to glucose. Specifically, glycerol enables yeast cells to enter the mid-exponential growth phase in a shorter timeframe, facilitating an earlier transition into the lipid accumulation stage. This rapid progression into the exponential phase is particularly beneficial for applications aimed at enhancing lipid production, as it allows for a more efficient use of resources and potentially increases overall yield. It allows us to further investigate the optimization of target product yield.

    Future Direction Improvement

    • Integrate our genetic circuit into the lipid accumulation part, seeing how the inducible promoter will affect the cell growth and the lipid accumulation
    • Obtain more experimental data to validate our model
    • Include more parameters in the model, e.g. the dissolved oxygen level, distribution of the substrate during agitation
    • Explore other possible metabolic mechanism we can use to manipulate our circuit to give the desired output

    References

    Workman, M., Holt, P., & Thykaer, J. (2013). Comparing cellular performance of Yarrowia lipolytica during growth on glycerol, glucose or a mixture of the two. AMB Express, 3(1), 58. https://doi.org/10.1186/2191-0855-3-58
    Robles-Rodríguez, C. E., Muñoz-Tamayo, R., Bideaux, C., Gorret, N., Guillouet, S. E., & Molina-Jouve, C. (2018). Modeling and optimization of lipid accumulation by Yarrowia lipolytica from glucose under nitrogen depletion conditions. Biotechnology and Bioengineering, 115(5), 1137-1151. https://doi.org/10.1002/bit.26537
    Kuttiraja, M., Workman, M., & Thykaer, J. (2023). Engineering strategies for efficient bioconversion of glycerol to value-added products by Yarrowia lipolytica. Catalysts, 13(4), 657. https://doi.org/10.3390/catal13040657
    Robles-Rodríguez, C. E., Muñoz-Tamayo, R., Bideaux, C., Gorret, N., Guillouet, S. E., & Molina-Jouve, C. (2018). Modeling and optimization of lipid accumulation by Yarrowia lipolytica from glucose under nitrogen depletion conditions. Biotechnology and Bioengineering, 115(5), 1137-1151. https://doi.org/10.1002/bit.26537
    Dobrowolski, A., Nicaud, J. M., & Kuttiraja, M. (2023). Novel evolved Yarrowia lipolytica strains for enhanced growth and lipid content. Microbial Cell Factories, 22(1), 62. https://doi.org/10.1186/s12934-023-02072-8
    Li, X., Yang, S., & Liu, J. (2019). Lipid production from glycerol by Yarrowia lipolytica: A review of its metabolic pathways and engineering strategies. Biotechnology Advances, 37(5), 107462. https://doi.org/10.1016/j.biotechadv.2019.107462
    Li, X., Yang, S., & Liu, J. (2019). Lipid production from glycerol by Yarrowia lipolytica: A review of its metabolic pathways and engineering strategies. Biotechnology Advances, 37(5), 107462. https://doi.org/10.1016/j.biotechadv.2019.107462