Model

Overview

The intracellular metabolic activity of Yarrowia lipolytica is complex, and the synthesis path of intracellular polyunsaturated fatty acids is long, which makes it difficult to accumulate EPA directionally. Traditional metabolic engineering strategies are mostly based on a certain metabolic pathway or specific genes, without considering the dynamic impact of genetic disturbance on cell growth and product synthesis. As a mathematical model, genome scale metabolic network model is widely used in analyzing network characteristics, predicting cell phenotypes, guiding metabolic engineering, driving model discovery, studying evolutionary processes and analyzing interactions. We tried to use the genome scale metabolic network model of Yarrowia lipolytica to systematically analyze the distribution of its intracellular metabolic flow, and combined with different algorithms to predict the transformation target of efficient EPA synthesis, providing guidance for the subsequent metabolic engineering transformation.

1.Selection and Improvement of the Starting Model.

Through literature mining, by 2022, people have constructed seven genome scale metabolic network models (Table 1) to understand Yarrowia lipolytica.The newly published genome scale metabolic network model iYli21 of Yarrowia lipolytica is selected as the starting model. (https://ars.els-cdn.com/content/image/1-s2.0-S200103702200174X-mmc6.xml)

Table 1 Statistics of published Yarrowia lipolytica models

Use the readCbModel instruction of COBRA toolbox to read the yeast model in SBML format into MATLAB, and get the basic information of iYli21. (Figure1)

Figure1 MATLAB read model iYli21

In order to view the model structure and specific content intuitively, use the writeCbModel instruction to convert the SBML format model to EXCEL format. The specific instruction is writeCbModel (model, 'format', 'xls',' fileName ',' iYli21. xlsx '), and finally get the model file named iYli21. xlsx.

The iYli21 model in EXCEL format has some problems, which need to be improved, specifically involving.

(1)Improvement of Cell Interval Information:

Improve the metabolite range so that it can be read by COBRA Toolbox.

Table2 Modification of metabolite partition information

(2)Improvement of Metabolite Information

As the metabolites in the original iYli21 model are all represented by m1, m2 and other strings, it is difficult to view the metabolite information conveniently, and its metabolites need to be replaced, such as replacing m32 representing H2O with h2o. Combined with BIGG Model, KEGG, MetaCyC and other databases, 1868 metabolite formats in iYli21 model are replaced.In the process of replacement, it is found that m200 [c] and m1855 [c] also represent malonyl CoA. Therefore, these two metabolites are replaced with malcoa [c], and the information related to m1855 [c] in the metabolite list is deleted. Similar situations also occur in m2017 [p] (h2o2 [p]), m2038 [c] (ahcys [c]), m2037 [c] (amet [c]), m2043 [c] (trp__L [c]), m2206 [c] (pac [c]), m1959 [c] (34droxpeg [c]), and m1963 [c] (acgam1p [c]). Use ultrareplace software to replace the metabolite format in iYli21 model. (Figure2)

Figure2 Model iYli21 after metabolite replacement

(3)Determination of Biomass Equation

In the iYli21 model, there are five reactions related to the biomass equation, three of which are the biomass equation of Saccharomyces cerevisiae (yeast 5 biomass pseudoaction, yeast 6 biomass pseudoaction, yeast 8 biomass pseudoaction), which have no direct relationship with Yarrowia lipolytica. Therefore, these three reactions are deleted first. In addition, a general biomass equation of microorganisms is also deleted (xBIOMASS). Finally, the replaced model is obtained, including 2281 reactions, 1860 metabolites and 1056 genes. (Figure3)

Figure3 Model after deleting redundant biomass equation

(4)Improvement of EPA Synthesis Pathway of Yarrowia Lipolytica

Since there is no EPA synthesis pathway in Yarrowia lipolytica (Figure 4), it is necessary to introduce relevant metabolic pathways manually in iYli21 model.

Figure4 Synthetic pathway of polyunsaturated fatty acids in Yarrowia lipolytica

Synthetic pathway of EPA 1(c184 -> c204 -> c205):

① lincoa[em] + nadph[em] + o2[em] + 2 h[c] <=> c183coa[em] + nadh[em] + 2 h2o[em]

② c183 coa [em] + nadph[em] + o2[em] + 2 h[c] <=> c184 coa [em] + nadh[em] + 2 h2o[em]

③ h[em] + malcoa[em] + c184coa[em] -> coa[em] + co2[em] + c184_ocoa[em]

④h[em] + c184_ocoa[em] + nadph[em] -> nadp[em] + c184_hcoa[em]

⑤ c184_hcoa[em] <=> h2o[em] + c184_dcoa[em]

⑥ h[em] + nadph[em] + c184_dcoa[em] -> c204coa[em] + nadp[em]

⑦ c204coa[em] + nadph[em] + o2[em] + 2 h[c] <=> c205coa[em] + nadh[em] + 2 h2o[em]

⑧ c205coa[c] <=> c205coa[em]

⑨ c205coa[c] <=> c205coa[p]

⑩ c205coa[p] + h2o[p] -> h[p] + coa[p] + c205[p]

Synthetic pathway of EPA 2(c183 -> c203 -> c204 -> c205):

① h[em] + malcoa[em] + c183coa[em] -> coa[em] + co2[em] + c183_ocoa[em]

②h[em] + c183_ocoa[em] + nadph[em] -> nadp[em] + c183_hcoa[em]

③ c183_hcoa[em] <=> h2o[em] + c183_dcoa[em]

④ h[em] + nadph[em] + c183_dcoa[em] -> c203coa[em] + nadp[em]

⑤ c203coa[em] + nadph[em] + o2[em] + 2 h[c] <=> c204coa[em] + nadh[em] + 2 h2o[em]

⑥ c204coa[em] + nadph[em] + o2[em] + 2 h[c] <=> c205coa[em] + nadh[em] + 2 h2o[em]

Transport reaction of EPA:

c205[p] <=> c205[c]

c205[c] <=> c205[e]

Exchange reaction of EPA:

c205[c] <=> c205[e]

After the introduction of EPA related synthetic pathway, the model iYli21 contains 2299 reactions, 1876 metabolites and 1058 genes. (Figure 5)

Figure5 Model after importing EPA synthesis pathway iYli21

2.Verification of iYli21 Model's Ability to Synthesize EPA

After completing the revision of model iYli21, we verify the prediction ability of the model. First, the biomass equation is taken as the target equation, and the flow balance analysis (FBA) algorithm is used to solve it. The specific growth rate of Yarrowia lipolytica yeast is 2.9231h-1, which is consistent with the starting model, indicating that the improvement work does not affect the prediction ability of the model. (Figure 6)

Figure6 Simulation of cell growth rate of Yarrowia lipolytica

In order to verify the ability of the model to predict EPA synthesis, by setting the upper and lower bounds of the biomass equation as 0.1 mmol/gDW/h, the exchange reaction EX of EPA_ C205 is the target equation, which is solved by FBA, and the resultant rate of EPA is 7.6744 mmol/gDW/h, which indicates that the model can synthesize EPA from scratch. (Figure7)

Figure7 Simulation of EPA synthesis rate of Yarrowia lipolytica

3.Using iYLi21 Model to Predict Targets for Improving EPA Synthesis

OptForce algorithm (Ranganathan et al. 2010) is widely used to predict targets of up regulation, down regulation and gene knockout, and has been successfully used to predict targets of product metabolic engineering transformation, such as fatty acids (Xu et al. 2011), amino acids (Ye et al. 2020), nucleotides (Deng et al. 2022). With reference to the COBRA toolbox 3.0 (Heirendt et al. 2019) tutorial, we try to use OptForce algorithm to predict potential targets for improving EPA synthesis capability.

Using OptForce algorithm, the upper and lower boundary flow of biomass equation is set as 0.1 mmol/gDW/h, so as to meet the basic growth demand of Yarrowia lipolytica. At the same time, the exchange reaction of EPA is selected as the target equation to solve. Finally, 8 over expression modification targets were obtained, of which 5 were EPA synthesis related genes and 3 were cell growth related (Table 3). Among these targets, some increase the EPA synthesis rate by increasing the supply of EPA synthesis cofactors (acetyl coenzyme A and NADPH), some directly strengthen the key enzymes in the EPA synthesis pathway, such as ∆ - 9 elongase, and some indirectly increase the EPA synthesis rate by increasing the cell growth rate.

Table 3 Modification Targets for Improving EPA Content Predicted by OptForce Algorithm

4.Materials and Methods

Table 4 Server Configuration

Table 5 List of software used for model analysis

Reference

Deng AH, Qiu QD, Sun QY, Chen ZX, Wang JY, Zhang Y, Liu SW, Wen TY (2022) In silico-guided metabolic engineering of Bacillus subtilis for efficient biosynthesis of purine nucleosides by blocking the key backflow nodes. Biotechnol Biof Biop 15(1) doi:ARTN 82 10.1186/s13068-022-02179-x

Guo YF, Su LQ, Liu Q, Zhu Y, Dai ZJ, Wang QH (2022) Dissecting carbon metabolism of Yarrowia lipolytica type strain W29 using genome-scale metabolic modelling. Comput Struct Biotec 20:2503-2511 doi:10.1016/j.csbj.2022.05.018

Heirendt L, Arreckx S, Pfau T, Mendoza SN, Richelle A, Heinken A, Haraldsdottir HS, Wachowiak J, Keating SM, Vlasov V, Magnusdottir S, Ng CY, Preciat G, Zagare A, Chan SHJ, Aurich MK, Clancy CM, Modamio J, Sauls JT, Noronha A, Bordbar A, Cousins B, El Assal DC, Valcarcel LV, Apaolaza I, Ghaderi S, Ahookhosh M, Ben Guebila M, Kostromins A, Sompairac N, Le HM, Ma D, Sun Y, Wang L, Yurkovich JT, Oliveira MAP, Vuong PT, El Assal LP, Kuperstein I, Zinovyev A, Hinton HS, Bryant WA, Aragon Artacho FJ, Planes FJ, Stalidzans E, Maass A, Vempala S, Hucka M, Saunders MA, Maranas CD, Lewis NE, Sauter T, Palsson BØ, Thiele I, Fleming RMT (2019) Creation and analysis of biochemical constraint-based models using the COBRA Toolbox v.3.0. Nat Protoc 14(3):639-702 doi:10.1038/s41596-018-0098-2

Kavscek M, Bhutada G, Madl T, Natter K (2015) Optimization of lipid production with a genome-scale model of Yarrowia lipolytica. BMC Syst Biol 9:72 doi:10.1186/s12918-015-0217-4

Kerkhoven EJ, Pomraning KR, Baker SE, Nielsen J (2016) Regulation of amino-acid metabolism controls flux to lipid accumulation in Yarrowia lipolytica. NPJ Syst Biol Appl 2:16005 doi:10.1038/npjsba.2016.5

Loira N, Dulermo T, Nicaud JM, Sherman DJ (2012) A genome-scale metabolic model of the lipid-accumulating yeast Yarrowia lipolytica. BMC Syst Biol 6:35

Mishra P, Lee NR, Lakshmanan M, Kim M, Kim BG, Lee DY (2018) Genome-scale model-driven strain design for dicarboxylic acid production in Yarrowia lipolytica. BMC Syst Biol 12:12 doi:10.1186/s12918-018-0542-5

Pan P, Hua Q (2012) Reconstruction and in silico analysis of metabolic network for an oleaginous yeast, Yarrowia lipolytica. PLoS ONE 7(12):e51535

Ranganathan S, Suthers PF, Maranas CD (2010) OptForce: an optimization procedure for identifying all genetic manipulations leading to targeted overproductions. PLoS Comput Biol 6(4):e1000744 doi:10.1371/journal.pcbi.1000744

Wei SS, Jian XX, Chen J, Zhang C, Hua Q (2017) Reconstruction of genome-scale metabolic model of Yarrowia lipolytica and its application in overproduction of triacylglycerol. Bioresour Bioprocess 4:51 doi:10.1186/s40643-017-0180-6

Xu P, Ranganathan S, Fowler ZL, Maranas CD, Koffas MAG (2011) Genome-scale metabolic network modeling results in minimal interventions that cooperatively force carbon flux towards malonyl-CoA. Metab Eng 13(5):578-587

Ye C, Luo Q, Guo L, Gao C, Xu N, Zhang L, Liu L, Chen X (2020) Improving lysine production through construction of an Escherichia coli enzyme-constrained model. Biotechnol Bioeng 117(11):3533-3544 doi:10.1002/bit.27485

Overview

In the era of sustainable development, the use of microorganisms to produce various compounds through fermentation has attracted extensive attention. For a commercialized fermentation project, efficient production strains are the prerequisite, and suitable fermentation condition as well as large-scale production are also essential. In this project, an engineered strain with high yield of PUFA was successfully obtained through synthetic biology. However, in actual production, we also found that extracellular conditions (such as stirring and ventilation) would have a great impact on the yield. To obtain suitable fermentation conditions, the traditional method is to repeatedly test the effect of different conditions on fermentation results through single factor experiment. There is no doubt that this will be a time-consuming and laborious work, and the interrelationship between different parameters as well as the size of different sizes of bioreactors will be more difficult to experiment. To solve this problem, computational fluid dynamics (CFD) technology was used in this project. The flow field characteristics of bioreactor can be quantitatively characterized by mathematical modeling and computer numerical solution. In addition, CFD can also be used to investigate the oxygen mass transfer in the reactor under different conditions, which is crucial for aerobic PUFA fermentation. Finally, through CFD numerical simulation, we compared different fermentation conditions in advance and determined the optimal fermentation conditions.

1. Introduction of Computational Fluid Dynamics (CFD)

Computational Fluid Dynamics (CFD) is a simulation technology integrating mathematics, fluid mechanics and computer, which can simulate and analyze various fluid problems. The main research content is to solve the control equations of fluid mechanics through computer and numerical methods, and simulate and analyze the fluid mechanics problems.

The basic steps of CFD are: (1) abstract the actual fermentation tank into a 3D model, using the software SpaceClaim; (2) The 3D model is divided into grids one by one, because the CFD software calculation is based on a single grid. All grids are connected in series to form the flow field data of the entire fermentation tank. The number of grids is 106 orders of magnitude. The software used is ICEM and Fluent Mesh; (3) Import the grid into the solver, and select an appropriate solution model for solution (different models use different equations, and different equations have different solution accuracy and time). The software used is Fluent. (Note: All software is included in the ANSYS 2020 R2 software package)

2、Build model

（1）3D modeling: In order to abstract the actual fermentation tank into a 3D model, we need to use SpaceClaim software for 3D modeling, and draw all parts and details of the fermentation tank. The fermentation tank is equipped with three layers of agitators and four baffles, and uses an annular air distributor. The vent holes with a diameter of 1mm are evenly distributed in the vent ring, four vertical baffles are evenly distributed on the reactor wall, and the blades are six flat bladed propellers.

（2）Grid generation: Although the computer is intelligent, it is not what you ask it to do. So is computational fluid dynamics. We can't tell the computer to operate on this fermenter. We can only give the computer instructions to operate on a more regular geometry. This is the role of grid. Through ICEM software or Fluent Meshing software, We can divide the drawn 3D CAD geometry into millions of grids with regular shapes for computer calculation. In the mesh generation step, we need to first divide the surface mesh, describe the details of each face of the geometry to the computer, and then perform the 3D mesh on the drawn surface mesh, and describe the entire 3D geometry to the computer, so that we can continue to calculate.

The number of grids in fermentation tanks is generally between 1.5 and 2 million, including tetrahedron and hexahedron grids.

3、Calculation preparation

In order to achieve our goal (EPA fermentation by Yarrowia lipolytica), we need to select appropriate models (hydrodynamic equations and computer simulation equations) for CFD simulation, that is, tell the computer what kind of fluid and what kind of environment we want to simulate through corresponding instructions. Considering the calculation force of the workstation, we will use MRF multiple reference system method (this setting is because our model has a rotating area: agitator) to solve.

（1）First, determine the relationship between air and liquid. In the solver, we regard air and liquid as two mobile phases and treat them with Euler Euler multiphase flow model. The material and energy exchange between the two phases is expressed by the following equation:

Mass transport equation

Where: t-time, s; ρ-density, kg/m3; u-speed, m/s; α is phase holdup, %; subscript α and β are gas phase and liquid phase respectively

Momentum transfer equation

Where: P-pressure, Pa; g-gravitational acceleration, m/s2; μeff-effective ni, Pa·s; Mα、 Mβ are respectively gas-liquid interface forces,N.

（2）The equation of motion of each phase is called turbulence equation. In the past research, the k-e turbulence equation was generally selected for gas-liquid two-phase flow in fermentation tanks. This double equation model is more consistent with the actual situation in fermentation tanks, and has low requirements for computing resources:

For k equation:

For ε equation:

Where:Gk is the turbulent kinetic energy generated by the average velocity gradient; Gb is the generation term of turbulent kinetic energy brought by buoyancy; Ym is the contribution of excessive dissipation to the total dissipation rate in compressible turbulence; Rε is an add-on; αk and αε are reverse effective Prandtl number of k and ε respectively; μeff is effective viscosity; C1ε、C2ε and C3ε are constant, C1ε=1.42, C2ε=1.68, C3ε=1.3。

（3）Mass transfer equation between the air phase and the liquid phase: For fermentation, it is important to simulate the mass transfer in the tank, so KLA equation is indispensable. KLA equation is the characterization of oxygen mass transfer coefficient. The equation is composed of two parts. For different multiphase fluid simulation, different KLA equations are often used. Here we choose the equation that meets the conditions in the fermentation tank, and the relevant coefficients include local gas holdup and average bubble diameter.

Where Dl is the oxygen diffusivity in the liquid phase, with a value of 2.54×10-9m2/s; k, n and ρ respectively denote the consistency coefficient, flow index and density of the fluid, ε is the local energy dissipation rate, αg is the local gas holdup, and d is the mean diameter of the bubbles.

（4）For the drawn grid, the most important area where computing resources are used is the blade part. The fluid motion near the blade is the most complex. We divide the flow area into two parts, the internal moving area and the external static area. The two areas transmit data through the grid interface, which is why the MRF method is used. For the grid interface, we also need to set the boundary conditions at the inlet and outlet. In order to get the most simulated simulation, we set the phase holdup settings at the inlet and outlet, the speed selection, and the inlet bubble diameter. After all the preliminary preparations for the calculation are made, the grid is initialized, that is, the simulation will inject the corresponding gas and liquid into our fermentation tank model, and then the calculation can begin. The operation generally requires more than 1000 steps of iteration. When the residuals of all variables are less than e-4 and the result meets the change rate of less than 1%, the result converges and can be used as a reference.

4、Result processing

（1）After the iterative calculation results are obtained, the data can be imported into the CFD-POST software to analyze and present the results.

Select the vector map output to observe the fluid velocity simulation results in the fermentation tank.

Select cloud image output to analyze gas holdup distribution and velocity distribution of fluid.

(2)concrete analysis:

The mixed flow in the tank is obviously separated, the corresponding mixing ability is poor, and it is difficult to obtain a uniform flow field by simply using the radial flow paddle. Therefore, the blade combination is improved, including the combination of axial flow paddle and increasing the number of blades to enhance the mixing of fluid in each rotating area.

Figure Fluid velocity diagram

A and B flow patterns are seriously separated.

The velocity distribution of C is significantly higher and more uniform.

The velocity distribution of D is relatively uniform, but the effect of axial flow propeller is too obvious, which makes the velocity at the bottom faster.

Gas holdup distribution diagram

The bubbles in A are more concentrated in the lower part.

Bubbles in B are mostly concentrated in the lower propeller area.

The gas holdup distribution of C is obviously more uniform.

The gas holdup of D is relatively uniform, but the bubbles in the lower propeller area are relatively concentrated.

Although the average kla of C is lower than that of D, its velocity and gas holdup distribution are more uniform, so C is the best combination.