A dynamic enzyme-constrained genome-scale metabolic network model construction method

Abstract

The application discloses a kind of dynamic enzyme constraint genome scale metabolic network model construction methods, comprising the following steps: step S1, determine kinetic model, obtain the specific numerical value of reaction flux and the change of metabolite with time;Step S2, according to the same reaction in ecGEMs upper and lower bound constraint of different reaction flux, the flux change condition of all reactions of organism is obtained by FBA solution with time change;Step S3, the model is verified.The ecGEMs can be made dynamic, the problem that ecGEMs cannot obtain metabolite concentration and intracellular reaction flux change with time is solved.And greatly improve the accuracy of ecGEMs model, it has important significance to the research system biology and synthetic biology research.

Description

Technical Field

[0001] This invention relates to the interdisciplinary field of biology and computational science, and in particular to a method for constructing a dynamic enzyme-constrained genome-scale metabolic network model. Background Technology

[0002] Genome-scale metabolic models (GEMs) have been widely used in systems biology and synthetic biology over the past few decades. GEMs link genes encoding enzymes to the reactions catalyzed by these enzymes, allowing for the prediction of phenotypes from cellular genotypes and thus providing a systematic characterization of various metabolic activities in organisms. GEMs have played a significant role in the study of cellular metabolism over the past decade.

[0003] GEMs use the Constraint-Based Reconstruction and Analysis Toolbox (COBRA) to simulate the metabolic states of organisms in different environments. Flux balance analysis (FBA) is the most widely used algorithm in COBRA because it can accurately quantify cellular metabolism. FBA simplifies complex metabolic networks into a linear programming problem, thus significantly reducing the number of parameters and model complexity. This simplification greatly aids in the analysis of cellular metabolism. The mathematical representation of FBA is shown below:

[0004] max c T v

[0005] stSv=0

[0006] lb i ≤v i ≤ub i

[0007] Where S∈R m×n This is a reaction matrix, where m represents the number of metabolites, n represents the number of reactions, and v∈R. n×1 For reaction flux, lb i ub is the lower bound of the i-th reaction. i c is the upper bound of the i-th reaction. T ∈R 1×n It is the target vector.

[0008] FBA (Fluidized Absorption Regulator) introduces various constraints to characterize the properties and behavior of a system, enabling accurate model descriptions without heavily relying on biological information such as kinetics and thermodynamics. FBA assumes the system remains in steady state and imposes constraints on external fluxes to obtain objective function values ​​and a set of flux distributions.

[0009] While GEMs provide a convenient framework for large-scale metabolic studies, they still have some limitations. For example, GEMs cannot accurately simulate metabolic changes that occur at high growth rates, such as spillover metabolism in bacteria and the Crabtree effect in Saccharomyces cerevisiae. Enzyme-constrained genome-scale metabolic network models (ecGEMs) can solve this problem by integrating enzyme kinetic parameters and enzyme concentration, and ecGEMs are therefore widely used. ecGEMs are described below:

[0010] Sv=0

[0011] lb≤v≤ub

[0012]

[0013] Where S∈R m×n This is a reaction matrix, where m represents the number of metabolites, n represents the number of reactions, and v∈R. n×1 Let lb be the reaction flux, ub be the lower bound of the reaction, and E be the upper bound of the reaction. i k is the enzyme concentration. cat Enzyme turnover number is a kinetic parameter of the enzyme.

[0014] However, ecGEMs still have certain limitations. While ecGEMs can provide the distribution of metabolic fluxes, they cannot obtain information on the changes in metabolite concentrations and intracellular responses over time, which is essential for studying the dynamic metabolic behavior of organisms. Secondly, ecGEMs struggle to quantitatively analyze the interactions between various components in somatic cells, hindering further in-depth research on cellular metabolism. To address these limitations, kinetic models have emerged. Kinetic models use ordinary differential equations (ODEs) to capture the dynamic changes of all components included in the model (metabolites, proteins, mRNA, etc.). The general expression is shown below:

[0015]

[0016] Where S and v represent the stoichiometric matrix and the vector of metabolic reaction flux, respectively. v is determined by the metabolite concentration (x), kinetic parameter (k), and enzyme concentration (E). x0 is the initial metabolite concentration, where the kinetic parameter mainly involves the maximum reaction rate constant v. max enzyme catalytic reaction constant k cat Michaelis constant k m equilibrium constant k eq etc.

[0017] Since the introduction of kinetic models, determining the ODE (Optical Derivative Expression) to describe a reaction has been a major challenge. One solution is to determine different ODEs for different reaction mechanisms based on the reaction principle; another is to use the same ODE for all reactions. Both methods require extensive experimental data for parameter fitting, consuming significant computational resources and time, which is a key difficulty in kinetic modeling. Furthermore, the structural and quantitative uncertainties in the kinetic representation further complicate kinetic modeling. Therefore, current kinetic modeling is limited to small, well-studied metabolic networks and has not yet achieved full-scale metabolic dynamics modeling at the cellular level. In contrast, ecGEMs (eco-cellular gene sequencing and analysis) model and analyze at the cellular level. Therefore, integrating kinetic models with ecGEMs to leverage their respective strengths is crucial. Summary of the Invention

[0018] The purpose of this invention is to provide a method for constructing a dynamic enzyme-constrained genome-scale metabolic network model. This method uses the time-varying reaction flux obtained from a biological dynamics model as an input variable, and uses this as a reference to constrain the same reactions in ecGEMs. The constrained ecGEMs are then subjected to FBA operations to obtain the distribution of all intracellular reaction fluxes of the organism over time, thereby solving the problems mentioned in the background art.

[0019] To achieve the above-mentioned objectives, this invention provides a method for constructing a dynamic enzyme-constrained genome-scale metabolic network model, comprising the following steps:

[0020] Step S1: Determine the kinetic model and obtain specific values ​​for reaction flux and metabolite changes over time;

[0021] Step S2: Apply upper and lower bound constraints to the same reaction in ecGEMs based on different reaction fluxes, and use FBA to solve for the flux changes of all reactions in the organism over time.

[0022] Step S3: Validate the model.

[0023] Furthermore, in step S1, organisms with well-established ecGEMs are collected, including Escherichia coli, Saccharomyces cerevisiae, Bacillus subtilis, Corynebacterium glutamicum, Agropyron lipolytica, and Aspergillus niger. The kinetic models involve pathways that have an important impact on the metabolism of organisms, including glycolysis, the tricarboxylic acid cycle, the pentose phosphate pathway, and lipid synthesis.

[0024] Furthermore, in step S1, the input constraints of ecGEMs are obtained by running the biodynamic model, as expressed by the following formula:

[0025] Y i =f i (xi ,t),i=1,2...,m

[0026] Where Y i The reaction flux is calculated using the kinetic model of the i-th pathway, and x is... i The i-th dynamic model input variable is t, which is the time variable, and m is the number of dynamic models that have been constructed.

[0027] Furthermore, in step S2, the reaction fluxes in the kinetic models obtained at different dilution rates are used as reference values. These fluxes are multiplied by a1 (a1∈[0.5,0.9]) to obtain the lower limit of the flux for the same reaction in ecGEMs, and multiplied by a2 (a2∈[1.1,1.5]) to obtain the upper limit of the flux for the same reaction in ecGEMs. FBA is then run to obtain the distribution of the fluxes of all reactions in Saccharomyces cerevisiae over time.

[0028] Furthermore, the FBA formula is expressed as follows:

[0029] max c T v

[0030] stSv=0

[0031] a1v di ≤v i ≤a2v di

[0032] lb j ≤v j ≤ub j

[0033]

[0034] Where S∈R m×n It is a reaction matrix, v∈R n×1 c is the reaction flux. T ∈R 1×n It is the target vector, v i For the reactions involved in the dynamic model, v di Let lb be the reaction flux value obtained from the kinetic model at the current moment. j ub is the lower bound of the j-th reaction. j v is the upper bound of the j-th reaction. j E represents the remaining reaction flux. i k is the enzyme concentration. cat Enzyme turnover is a kinetic parameter of the enzyme.

[0035] Furthermore, in step S3, model validation includes three methods: flux variability analysis, model self-validation, and experimental data validation.

[0036] Furthermore, the formula for the flux variability analysis method is expressed as follows:

[0037] max / min v i

[0038] stSv=0

[0039] Z0≤γw T

[0040] lb j ≤v j ≤ub j

[0041]

[0042] Where Z0 = w T It is the optimal solution for the following formula:

[0043] max c T v

[0044] stSv=0

[0045] lb i ≤v i ≤ub i

[0046] γ is a parameter that takes a value between 0 and 1, S∈R m×n This is a reaction matrix, where m represents the number of metabolites, n represents the number of reactions, and v∈R. n×1 For reaction flux, lb j ub is the lower bound of the j-th reaction. j c is the upper bound of the j-th reaction. T ∈R 1×n It is the target vector, v j E represents the remaining reaction flux. i k is the enzyme concentration. cat Enzyme turnover is a kinetic parameter of the enzyme.

[0047] Furthermore, the model self-validation is based on actual throughput experimental data from organisms. Sensitivity analysis is performed on n common reactions of the kinetic model and the organism's ecGEMs. The m reactions with the highest and lowest sensitivity, totaling 2m reactions, are selected. The remaining n-2m reactions obtained from the dynamic EMP (Embden-Meyerhof-Parnas pathway) model are compared with the remaining n-2m reactions obtained from FBA running on the 2m fixed reactions selected from ecGEMs for model self-validation.

[0048] Furthermore, the experimental data verification involves comparing the reaction flux of the organism's ecGEMs with the reaction flux obtained from the integrated dynamic model and the organism's ecGEMs, respectively, with the actual experimental data of the organism to obtain the error before and after the integration.

[0049] Compared with existing technologies, this system and method have the following advantages:

[0050] (1) Using flux as the point of integration between dynamic and steady-state models provides a way to improve the accuracy of ecGEMs of organisms.

[0051] (2) By integrating the dynamic model with the organism's ecGEMs, the problem that steady-state ecGEMs cannot obtain the metabolite concentration and the change of intracellular reaction flux over time is solved. Attached Figure Description

[0052] Figure 1 This is a flowchart of a method for constructing a dynamic enzyme-constrained genome-scale metabolic network model.

[0053] Figure 2 A schematic diagram illustrating the methodology for integrating a kinetic model of *Saccharomyces cerevisiae* with an enzyme-constrained genome-scale metabolic network model.

[0054] Figure 3 The EMP kinetic model includes the main reactions and metabolite concentrations, as well as a graph showing the change in reaction flux.

[0055] Figure 4 This is a comparison chart of model self-validation errors.

[0056] Figure 5 This is a comparison chart of the errors between the model before and after integration and the experimental values. Detailed Implementation

[0057] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0058] like Figure 1 The diagram shown is a flowchart of the method of the present invention, and the specific steps are as follows:

[0059] Step S1: Determine the kinetic model and obtain specific values ​​for reaction flux and metabolite changes over time.

[0060] First, kinetic models are collected, involving organisms with well-established ecGEMs, such as *Escherichia coli*, *Saccharomyces cerevisiae*, *Bacillus subtilis*, *Corynebacterium glutamicum*, *Agrobacterium lipolyticum*, and *Aspergillus niger*. These kinetic models can cover pathways that significantly impact the organism's metabolism, including glycolysis, the tricarboxylic acid cycle, the pentose phosphate pathway, and lipid synthesis. Then, the kinetic models are appropriately modified according to the corresponding organism's ecGEMs to ensure consistency between the two models. The input constraints of the ecGEMs are obtained by running the kinetic models. The mathematical description is as follows:

[0061] Y i =f i (x i ,t),i=1,2...,m

[0062] Where Y i The reaction flux is calculated using the kinetic model of the i-th pathway, and x is... i The i-th dynamic model input variable is t, which is the time variable, and m is the number of dynamic models that have been constructed.

[0063] Combination method such as Figure 2 The diagram shown illustrates the methodology for integrating the kinetic model of Saccharomyces cerevisiae with an enzyme-constrained genome-scale metabolic network model.

[0064] Step S2: Apply upper and lower bound constraints to the same reaction in ecGEMs based on different reaction fluxes, and use FBA to solve for the flux changes of all reactions in the organism over time.

[0065] Using the reaction fluxes from the EMP kinetic model at different dilution rates as reference values, multiply them by a1 (a1∈[0.5,0.9]) to obtain the lower limit of flux for the same reaction in ecGEMs, and multiply them by a2 (a2∈[1.1,1.5]) to obtain the upper limit of flux for the same reaction in ecGEMs. Then, run FBA to obtain the distribution of fluxes for all reactions in Saccharomyces cerevisiae over time. Here, FBA becomes:

[0066] max c T v

[0067] stSv=0

[0068] a1v di ≤v i ≤a2v di

[0069] lb j ≤v j ≤ub j

[0070]

[0071] Where S∈R m×n It is a reaction matrix, v∈R n×1 c is the reaction flux. T ∈R 1×n It is the target vector, v i For the reactions involved in the dynamic model, v di Let lb be the reaction flux value obtained from the kinetic model at the current moment. j ub is the lower bound of the j-th reaction. j v is the upper bound of the j-th reaction. j E represents the remaining reaction flux. i k is the enzyme concentration. cat Enzyme turnover is a kinetic parameter of the enzyme.

[0072] Step S3, Model Validation.

[0073] Model validation is mainly carried out through three methods:

[0074] Method 1, Flux variability analysis (FVA).

[0075] Flux variability analysis is commonly used to determine the robustness of metabolic models under various simulation conditions. The mathematical expression for FVA is:

[0076] max / min v i

[0077] stSv=0

[0078] Z0≤γw T

[0079] lb j ≤v j ≤ub j

[0080]

[0081] Where Z0 = w T It is the optimal solution for the following formula:

[0082] max c T v

[0083] stSv=0

[0084] lb i ≤v i ≤ub i

[0085] γ is a parameter that takes a value between 0 and 1, S∈Rm×n This is a reaction matrix, where m represents the number of metabolites, n represents the number of reactions, and v∈R. n×1 For reaction flux, lb j ub is the lower bound of the j-th reaction. j c is the upper bound of the j-th reaction. T ∈R 1×n It is the target vector, v j E represents the remaining reaction flux. i k is the enzyme concentration. cat Enzyme turnover is a kinetic parameter of the enzyme. If FVA analysis shows a decrease in model variability after binding, it indicates that the model performance has been improved.

[0086] Method 2, model self-validation.

[0087] Using actual throughput experimental data from organisms as a reference, sensitivity analysis was performed on n reactions common to both the kinetic model and the organism's ecGEMs. The m reactions with the highest and lowest sensitivity, totaling 2m reactions, were selected. The remaining n-2m reactions obtained from the dynamic EMP model were compared with the remaining n-2m reactions obtained from FBA running on the fixed 2m reactions selected from the ecGEMs for model self-validation. The errors before and after the combination can be obtained, thus indicating whether the combination is reasonable.

[0088] Method 3: Validation using experimental data.

[0089] The reaction fluxes of biological ecGEMs and the reaction fluxes obtained from the integrated dynamic model and biological ecGEMs of this invention are compared with real experimental data of organisms to obtain the errors before and after integration. This allows us to determine whether the integration method proposed in this invention is beneficial to improving the accuracy of biological ecGEMs.

[0090] To clearly describe the integration approach, the kinetic model of the Embden-Meyerhof-Parnas pathway (EMP) of Saccharomyces cerevisiae is integrated with Saccharomyces cerevisiae ecGEMs to describe the specific implementation method. The integration method of the kinetic models of other pathways is the same as that of ecGEMs.

[0091] The kinetic model adopted is the well-established EMP kinetic model of *Saccharomyces cerevisiae*. To facilitate subsequent integration with its ecGEMs, we modified the original EMP kinetic model. Since integration requires kinetic models at different dilution rates, and the kinetic parameters involved in these models vary with different dilution rates, we calculated the kinetic parameters at different dilution rates using piecewise cubic polynomial interpolation, obtaining eight different dilution rates (0.027 h⁻¹). -1 0.051h -10.118h -1 0.154h -1 0.207h -1 0.249h -1 0.316h -1 0.355h -1 The EMP kinetic model under this model includes 29 reactions (vGLT, vGLK, vPGI, vPFK, vALD, vTPI, vGAPDH, vPGK, vPGM, vENO, vPYK, vPDC, vADH, vG3PDH, vHOR2, vVacPi, vADK1, vmitoNADH, vETOHt, vGLYCt, vATPase, vmito, vsinkG6P, vsinkF6P, vsinkGAP, vsinkP3G, vsinkPEP, vsinkPYR, vsinkACE) and 21 metabolites (G6P, GLCi, NAD, NADH, ATP, P2G, P3G, PEP, PYR, GLYCERAL3P, ADP, AMP, DHAP, GLYC3P, GL YCEROL, ETOH, G1P, TRE, GLYCec, BPG, PHOS), 45 kinetic parameters (G6P, GLCi, NAD, NADH, ATP, P2G, P3G, PEP, PYR, GLYCERAL3P, ADP, AMP, DHAP, GLYC3P, GLYCEROL, ETOH, G1P, TRE, GLYCec, p_GLT_VmGLT, p.HXK_ExprsCor, p.PGI_ExprsCor, p.PFK_ExprsCor, p.FBA_ExprsCor, p.TPI_ExprsCor, p.GAPDH_ExprsCor, p.PGK_ExprsCor, p.PG M_ExprsCor, p.ENO_ExprsCor, p.PYK_ExprsCor, p.PDC_ExprsCor, p.ADH_ExprsCor, f_GlCo, p.poly_sinkACE, p.poly_sinkF6P, p.poly_sin kG6P, p.poly_sinkGAP, p.poly_sinkP3G, p.poly_sinkPEP, p.poly_sinkPYR, p_mitoNADHVmax, RQ, p.gam_ngam, adp_atp_ratio, p_T6P_ss).

[0092] For ease of description, the following only considers 0.027h from 8 different dilution rates. -1The kinetic model at the dilution rate was integrated with e cGEMs, and the integration method for other dilution rates was the same. Piecewise cubic polynomial interpolation of the 45 kinetic parameters in the model yielded a value of 0.027h. -1The kinetic model at the dilution rate yielded 45 kinetic parameters: G6P = 1.09684093517648, GLCi = 0.118979261665238, NAD = 1.57950584410507, NADH = 0.0104941558949249, ATP = 2.79653743424957, P2G = 0.286343064496289, P3G = 2.38938534937561, PEP=1.15570412745331, PYR=0.266866352602992, GLYCERAL3P=0.00063497 9086903727, ADP=0.766679639293019, AMP=0.218513060483924, DHAP=0.0838634387000135, GL YC3P=0.0402357801663814, GLYCEROL=0.122214112919759, ETOH=0, G1P=0.1, TRE=125.347006856529, GLYCe c=0.100005858162466, p_GLT_VmGLT=1.72719667166439, p.HXK_ExprsCor=1.31243129845865, p.PGI_Exprs Cor=1.00172412081006,p.PFK_ExprsCor=1.00290408474836,p.FBA_ExprsCor=1.11908785294106,p.TPI_E xprsCor=1.21746764132189,p.GAPDH_ExprsCor=1.4727106219162,p.PGK_ExprsCor=1.2970865370221,p.PG M_ExprsCor=1.31190013040651, p.ENO_ExprsCor=1.58021536154459, p.PYK_ExprsCor= 1.15718799856897, p.PDC_ExprsCor=1.34472986275393, p.ADH_ExprsCor=1.284773775 24626, f_GlCo=0.132214112919759, p.poly_sinkACE=-0.00937597409185371, p.poly_s inkF6P=0.00608215002899966, p.poly_sinkG6P=-0.0222991353718, p.poly_sinkGAP=0.00545677667953178, p.poly_sinkP3G=-0.0041405749, p.poly_sinkPEP=-0.0025123373, p.poly_sinkPYR=-0.0484338781019264, p_mitoNADHVmax=0.141225290814443, RQ=1.0467349537037, p.gam_ngam=0.322231922115345, adp_atp_ratio=0.274153183112584, p_T6P_ss=0.153874697609578. Based on the calculation, 0.027h is obtained. -1 The dynamic EMP model at dilution rate can produce curves showing the changes in metabolites and reaction flux over time, such as... Figure 3 The figure shows the flux changes of the EMP kinetic model, including the main reactions and metabolite concentrations. The flux at each time step can be used as input to the ecGEMs model at that time step to obtain the flux distribution of all reactions at that time step.

[0093] Comparing the reaction fluxes calculated by the ecGEMs model with those obtained by the integrated method proposed in this invention, and then with experimental values, it can be seen that in the self-validation, only one of the seven common reactions has a large error, while the errors of the rest are very small. This also proves the rationality of the combined method. Specific errors are as follows: Figure 4 As shown, this is a comparison chart of model self-validation errors. The error calculated directly using the ecGEMs model is 58.12%, while the error of the integrated method proposed in this invention is 35.82%, representing a relative improvement of 38.26%. The errors are as follows... Figure 5 The figure shown is a comparison of the error between the model before and after integration and the experimental values.

[0094] In summary, the method of this invention can dynamically represent ecGEMs, solving the problems that ecGEMs cannot obtain information on metabolite concentrations and changes in intracellular reaction flux over time. Furthermore, it significantly improves the accuracy of ecGEM models, which is of great significance for research in systems biology and synthetic biology.

[0095] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for constructing a dynamic enzyme-constrained genome-scale metabolic network model, characterized in that, Includes the following steps: Step S1: Determine the kinetic model and obtain specific values ​​for reaction flux and metabolite changes over time; In step S1, organisms with well-established ecGEMs are collected, including Escherichia coli, Saccharomyces cerevisiae, Bacillus subtilis, Corynebacterium glutamicum, Agropyron lipolytica, and Aspergillus niger. The kinetic models involve pathways that have an important impact on the metabolism of organisms, including glycolysis, the tricarboxylic acid cycle, the pentose phosphate pathway, and lipid synthesis. In step S1, the input constraints of ecGEMs are obtained by running the biodynamic model, as expressed by the following formula: ， in It is the first The reaction flux was calculated using a kinetic model of one pathway. No. There are two input variables for the dynamic model, where t is the time variable. This represents the number of already constructed dynamic models; Step S2 involves applying upper and lower bound constraints to the same reactions in ecGEMs based on different reaction fluxes, and using FBA to obtain the flux changes of all reactions in the organism over time. In step S2, the reaction fluxes obtained in the kinetic models at different dilution rates are used as reference values ​​and multiplied by [the appropriate factor]. , As the lower limit of flux for the same reaction in ecGEMs, multiplied by , As the upper limit of the flux for the same reaction in ecGEMs, FBA was run to obtain the distribution of the flux of all reactions in Saccharomyces cerevisiae over time. The FBA formula is expressed as follows: ， in, It is a reaction matrix. For reaction flux, It is the target vector. For the reactions involved in the dynamic model, This represents the reaction flux value obtained from the kinetic model at the current moment. This is the lower bound of the j-th reaction. This is the upper bound of the j-th reaction. This represents the remaining reaction flux value. Enzyme concentration, Enzyme turnover number is a kinetic parameter of the enzyme. Step S3: Validate the model. Model validation includes three methods: flux variability analysis, model self-validation, and experimental data validation. The formula for the flux variability analysis method is expressed as follows: ， in It is the optimal solution for the following formula: ， It is a parameter that takes a value between 0 and 1. This is a reaction matrix, where m represents the number of metabolites and n represents the number of reactions. For reaction flux, This is the lower bound of the j-th reaction. This is the upper bound of the j-th reaction. It is the target vector. This represents the remaining reaction flux value. Enzyme concentration, Enzyme turnover number is a kinetic parameter of the enzyme. The model self-validation is based on actual flux experimental data of organisms. Sensitivity analysis is performed on n common reactions of the kinetic model and the organism's ecGEMs. The m reactions with the highest and lowest sensitivity, totaling 2m reactions, are selected. The remaining n-2m reactions obtained from the dynamic EMP model are compared with the remaining n-2m reactions obtained by running FBA on the 2m fixed reactions selected from the ecGEMs for model self-validation.

2. The method for constructing a dynamic enzyme-constrained genome-scale metabolic network model according to claim 1, characterized in that, The experimental data verification involves comparing the reaction flux of the organism's ecGEMs with the reaction flux obtained from the integrated dynamic model and the organism's ecGEMs, respectively, with the actual experimental data of the organism, to obtain the error before and after the integration.

Images