Generation method and application of multi-constraint metabolic network model toolkit

Abstract

The invention relates to the technical field of network model toolkits, in particular to a generation method and application of a multi-constraint metabolic network model toolkit. The invention discloses a method for generating a multi-constraint metabolic network model toolkit. The method comprises the following steps: constructing a multi-constraint metabolic network model; performing parameter correction on the multi-constraint metabolic network model to obtain a model after parameter correction; performing model analysis by using the model after parameter correction to obtain a model analysis result; the multi-constraint metabolic network model comprises an enzyme constraint model, an enzyme heat constraint model and a thermodynamic constraint model; the parameter correction comprises enzyme kinetic parameter correction and thermodynamic parameter correction. According to the generation method disclosed by the invention, the influence of the metabolite concentration on reversibility can be effectively reflected by integrating thermodynamic factors, so that the ubiquitous problems of constraint relaxation, prediction distortion, single research angle and the like caused by excessive simplification of a current enzyme constraint model are solved.

Description

Technical Field

[0001] The present invention relates to the technical field of network model toolkits, and in particular to a generation method and application of a multi-constraint metabolic network model toolkit. Background Art

[0002] Genome-scale metabolic models (GEMs / GSMM) are mathematical models that cover the biochemical reaction network in cells and can quantitatively describe the relationship between genes and phenotypes. They are widely used to simulate and predict metabolic fluxes at different system levels. At present, metabolic network model analysis algorithms can be divided into two categories: one is the flux balance analysis algorithm, which is based on the flux balance analysis (FBA), which is mainly used to predict the flux distribution of metabolic pathways.

[0003] In order to improve the accuracy of predicting flux distribution and expand the scope of application of the model, a series of algorithms that add additional constraints to metabolic flux have emerged, such as Regulatory FBA (rFBA), Integrated FBA (iFBA), Dynamic FBA (dFBA) and Community FBA (cFBA). Another type of algorithm is related to the prediction of metabolic engineering targets, which can be divided into gene knockout algorithm (OptKnock), gene addition algorithm (OptSwap), gene expression regulation algorithm (OptForce and k-OptForce), heterologous pathway expression algorithm (OptStrain and SimOptStrain) and flux scanning algorithm (FESOF). These algorithms identify overexpressed targets, weakened targets and knockout targets in metabolic pathways that can be used to improve the flux of target products through simulation.

[0004] In the flux balance analysis of metabolic networks based on metrology, most reaction rates are subject to broad upper and lower limits. Therefore, when the substrate input rate is set to a high value, the internal reaction flux value may be too large and deviate from the experimental value. In actual cells, the rate of each reaction is limited by the content and specific enzyme activity of the enzyme that catalyzes the reaction. In addition, the total amount of enzyme per unit mass, unit volume or single cell also needs to be limited to a certain range. The enzyme constraint model can be obtained by introducing these system-level enzyme resource constraints into the metrology model in a specific mathematical form.

[0005] At present, enzyme-constrained metabolic network model analysis algorithms are mainly divided into two categories: one type of algorithm is used to simulate and predict biological phenomena, such as overflow metabolism, substrate hierarchical utilization, energy efficiency and other biological phenomena. The other type of algorithm is metabolic engineering target simulation, which can be used to identify key enzymes in the pathway, such as enzyme cost-based ranking method, FESOF-based enzyme amount scanning algorithm, and enzyme cost difference method based on different experimental conditions (such as high growth and low product production (HGLP) and low growth and high product production (LGHP) conditions).

[0006] In addition, due to the lack of thermodynamically relevant parameters in the above models, it is often difficult to accurately evaluate the reversibility of the reaction, which may lead to the lack of biological feasibility of the pathway predicted by the model, such as the infinite cycle problem of energy and reducing power. At present, there are two main methods to introduce the thermodynamic parameters of the reaction into the metabolic network model: one is to add thermodynamic constraints to each metabolic reaction in the GEM, and to achieve the introduction of pathway thermodynamic constraints through mixed integer linear programming. The other method is to maximize the thermodynamic driving force as the objective function (MDF method). The core principle of this method follows the barrel principle, that is, in the solution space that satisfies the metabolite concentration constraint, find the concentration distribution that can make the thermodynamic bottleneck reaction in the pathway the highest. Through this method, the thermodynamic bottleneck (referring to the reaction with the worst thermodynamic feasibility in the pathway) and the limiting metabolite (referring to the metabolite that plays a limiting role in the thermodynamic level of the pathway) in the pathway can be predicted, and the corresponding optimization strategy can be formulated. However, the evaluation method of thermodynamic bottlenecks and limiting metabolites has been proposed a long time ago, but there has been no algorithm that can be used directly.

[0007] In summary, as the number of constraints increases, the simulation calculation will be very slow. To this end, the present application provides a method for generating and applying a multi-constraint metabolic network model toolkit. Summary of the invention

[0008] In order to overcome the shortcomings of the prior art, the present invention provides a method for generating and applying a multi-constraint metabolic network model toolkit. The generation method can effectively reflect the effect of metabolite concentration on reversibility by integrating thermodynamic factors, thereby compensating for the problems of constraint relaxation, prediction distortion and single research angle caused by over-simplification that are common in current enzyme constraint models.

[0009] The technical solution adopted by the present invention to solve its technical problem is:

[0010] The first purpose of this application is a method for generating a multi-constraint metabolic network model toolkit, comprising the following steps:

[0011] Construct multi-constraint metabolic network models;

[0012] Performing parameter correction on the multi-constraint metabolic network model to obtain a model with corrected parameters;

[0013] Using the model after parameter correction, performing model analysis to obtain model analysis results;

[0014] The multi-constraint metabolic network model includes an enzyme constraint model, an enzyme thermal constraint model, and a thermodynamic constraint model;

[0015] The parameter correction includes enzyme kinetic parameter correction and thermodynamic parameter correction.

[0016] Preferably, the model analysis includes phenotypic simulation analysis and metabolic engineering target prediction analysis;

[0017] The metabolic engineering target prediction analysis includes bottleneck reaction mining, isozyme coupling mining, continuous reaction coupling mining, energy coupling mining, limiting metabolite mining, and key enzyme mining.

[0018] Preferably, the model analysis result is to visualize the result;

[0019] The result visualization includes MDF line graph, cumulative distribution graph, phenotypic phase plane graph, metabolic overflow graph, tradeoff graph, and pathway visualization.

[0020] Preferably, the thermodynamic parameter correction includes kinetic parameter correction and thermodynamic parameter correction.

[0021] Preferably, the kinetic parameter correction includes using an enzyme kinetic parameter automated correction algorithm, the process is:

[0022] First, the enzyme usage is calculated according to the current model, and the reaction that uses the most enzyme is selected;

[0023] Then, according to the EC code of the reaction, find the corresponding kcat maximum value in the kcat_max dictionary and apply it to the model;

[0024] Finally, update the loop counter, record the result of each operation, and print the result.

[0025] Preferably, according to the generation method according to claim 4, it is characterized in that the correction of the kinetic parameters includes using an automated correction algorithm for thermodynamic parameters to perform correction.

[0026] Preferably, the phenotypic simulation analysis includes growth rate simulation and product simulation.

[0027] Preferably, the algorithm for bottleneck reaction mining has the following process:

[0028] Set a lower limit on the flux of the biomass synthesis reaction and then determine the MDF level that the pathway can achieve;

[0029] This MDF value was then used as a constraint to solve for the maximum growth rate level;

[0030] Then, the maximum biomass synthesis rate and MDF value were used as constraints to calculate the upper and lower limits of the thermodynamic driving force of each reaction in the model;

[0031] Finally, a reaction was identified as a bottleneck reaction if both the upper and lower limits were equal to the MDF value; where the variability of the reaction was analyzed within the pathway output by the pFBA method and the thermodynamic driving force value was equal to the MDF.

[0032] More preferably, isozyme coupling is sought for the bottleneck reaction, and the process is as follows: based on the identified bottleneck reaction, a reaction associated with the same gene-protein-reaction (GPR) and having a standard Gibbs free energy less than 0 is searched in a multi-constraint model.

[0033] The second purpose of this application is to provide an application of a multi-constraint metabolic network model toolkit, and to apply the multi-constraint metabolic network model toolkit to predict the bottlenecks and key regulatory transformation targets of cellular metabolism.

[0034] Compared with the prior art, the present invention has the following beneficial effects:

[0035] 1. The generation method described in this application, the generated toolkit, contains a variety of simulation analysis methods and visualization methods, such as thermodynamic bottleneck reaction analysis algorithm, restriction enzyme analysis algorithm, and restriction metabolite analysis algorithm. And also includes isozyme coupling mining algorithm, continuous reaction coupling mining algorithm and energy coupling mining algorithm; the toolkit can mine the pathway bottleneck reaction that has not been widely concerned from the perspective of pathway thermodynamic feasibility. In addition, the simulation result visualization method of the toolkit enables researchers without programming ability to intuitively see the relevant results of model analysis.

[0036] 2. The generation method described in this application and the toolkit generated can help researchers analyze multi-constraint network models from multiple perspectives, thereby obtaining various types of targets such as thermodynamic bottleneck reactions, restrictive enzymes and restrictive metabolites in growth and chemical production pathways, thereby improving the coverage and selectivity of targets and the efficiency of strain transformation.

[0037] 3. The generation method described in the present application, when applied, can discover various types of targets for assisting metabolic engineering transformation and developing industrial strains that can produce high value-added products. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] The present invention is further described below in conjunction with the accompanying drawings and embodiments.

[0039] Figure 1 A workflow diagram of a method for generating a multi-constraint metabolic network model toolkit described in this application;

[0040] Figure 2-9 A workflow diagram of corresponding modules of a method for generating a multi-constraint metabolic network model toolkit described in this application;

[0041] Fig.10 MDF line graph calculated for the inflection point of this application;

[0042] Fig.11 Example diagram for visualizing the calculation results of this application;

[0043] Fig.12 The reaction screening results in the metabolic engineering target calculation for this application;

[0044] Fig.13 After coupling for this application, the simulation results of the model are compared;

[0045] Fig.14 For this application case, the MDF level changes in the serine synthesis process;

[0046] Fig.15 The application case of this application is a strategy to solve the thermodynamic bottleneck in the L-serine synthesis pathway. DETAILED DESCRIPTION

[0047] For the convenience of understanding by those skilled in the art, the present invention is further described below in conjunction with the embodiments, and the contents mentioned in the embodiments are not limitations of the present invention. As used herein, "and / or" includes terms of all combinations of any one or more associated listed items. The terms used here are only used to describe specific embodiments, and are not intended to limit the present invention. As used herein, the singular forms "one", "a", "a kind" and "the" are also intended to include plural forms, unless the context clearly indicates otherwise. It is further understood that "including" when used in this specification specifies the stated features, integers, steps, operations, elements and / or compositions, but does not exclude the presence or addition of one or more other features, integers, steps, operations, elements, compositions and / or combinations thereof. Unless otherwise defined, all terms (including technical and scientific terms) used herein have the same meanings generally understood by ordinary technicians in the technical field to which the present invention belongs. It is further understood that terms, such as defined in common dictionaries, are interpreted in accordance with their meanings in the context of the relevant field, and are not idealized or overly formal meanings, unless explicitly defined here. The exemplary invention described herein may appropriately lack any one or more element limitations, which are not particularly disclosed here. Therefore, the terms "comprising", "including", "containing", etc. should be interpreted broadly and non-restrictively. In addition, the terminology used herein is used as a description, not a limitation, and the use of these terminology expressions that do not include any equivalent characteristics is not intended to describe only part of their characteristics, but various modifications are possible within the scope of the present invention according to the rights. Therefore, although the present invention has been specifically disclosed through preferred embodiments and optional features, the modifications disclosed herein to embody the changes of the present invention may be recorded by those skilled in the art, and such modifications and changes will be considered to be within the scope of the present invention.

[0048] Initially, the inventors proposed a semi-automated modeling framework (ETGEMs) that integrates thermodynamic constraints with enzyme resource constraints, and implemented analysis algorithms for multiple pathway targets for this framework (Table 1). The ETGEMs analysis process provides a variety of information, including the distribution of enzyme costs in the pathway and the distribution of the optimal thermodynamic driving force level of the reaction in addition to the flux distribution. In addition, the limiting metabolites in the pathway can be determined by the changes in metabolite concentrations, and the bottleneck reactions and key enzymes in the pathway (the reactions that occupy the majority of enzyme costs in the pathway, the increase in their enzyme activity has a very obvious effect on the flux of the pathway, and are usually the implementation targets of strategies such as overexpression and directed evolution of enzymes) can be determined based on the thermodynamic driving force and enzyme cost distribution information.

[0049] Table 1 Algorithms implemented by ETGEMs

[0050]

[0051] Recently, the inventors have restructured the ETGEMs framework and defined the model structure (including five objects: reactions, genes, metabolites, models and environments) in an object-oriented way to facilitate future model modification and expansion. We call this toolkit for analyzing multi-constraint model systems MCM_Analyzer.

[0052] Specifically, the multi-constraint metabolic network model is constructed based on modular thinking. In FlexModeler, a multi-constraint template function Template_Concretemodel is built. By controlling the parameters passed in, different types of models can be constructed, such as GEM (metabolic network model), ECM (enzyme constraint model), TCM (thermodynamic constraint model), ETM (enzyme thermal constraint model), etc. The inputs are all metabolic network models, but the parameters passed are different.

[0053] For example, metabolic network models

[0054] GEM=Template_Concretemodel(model, obj_target=obj_target, set_obj_value=True, set_substrate_ini=True, set_stoi_matrix=True, set_bound=True);

[0055] Enzyme thermal constraint model

[0056] ETM=Template_Concretemodel(model,obj_target=obj_target,set_obj_value=Tru e, set_substrate_ini=True, set_stoi_matrix=True, set_bound=True, set_metabolite=True, set_Df=True, set_integer=True, set_metabolite_ratio=True, set_thermodynamics=True, B_value=B_value, set_enzyme_constraint=True)

[0057] like Figure 1 As shown, the present application provides a method for generating a multi-constraint metabolic network model toolkit (MCM_Analyzer), comprising the following steps:

[0058] Construct multi-constraint metabolic network models;

[0059] Performing parameter correction on the multi-constraint metabolic network model to obtain a model with corrected parameters;

[0060] Using the model after parameter correction, performing model analysis to obtain model analysis results;

[0061] The multi-constraint metabolic network model includes an enzyme constraint model, an enzyme thermal constraint model, and a thermodynamic constraint model;

[0062] The parameter correction includes enzyme kinetic parameter correction and thermodynamic parameter correction.

[0063] Specifically, the parameter correction is to run the kinetic parameter correction process, perform parameter correction, and obtain a model after the kinetic parameter correction. Then, the thermodynamic parameter correction process is performed to perform parameter correction and obtain a model after the thermodynamic parameter correction; more specifically, Figure 2 As shown, when performing kinetic parameter correction, the model, correction times, threshold, and correction file are input, and the model after kinetic parameter correction is output; when performing thermodynamic parameter correction, the model, optimization target, minimum MDF value, and concurrent number are input, and the model after thermodynamic parameter correction is output.

[0064] Specifically, model analysis is performed using a model corrected by thermodynamic parameters. The model analysis includes two parts: phenotypic simulation and metabolic engineering target prediction. Phenotypic simulation involves growth rate simulation and product simulation. The simulation process requires the input of model files, substrates, substrate concentrations, objective functions and other parameter information, and the output result is the flux distribution result. Metabolic engineering target prediction includes bottleneck reaction mining, isozyme coupling mining, continuous reaction coupling mining, energy coupling mining and limiting metabolite mining. The common input parameters of these target mining algorithms are models and post-selection translation lists, and the output results are various metabolic engineering strategies. When mining bottleneck reactions, limiting metabolites and key enzymes, there are specific parameters such as B value, objective function, maximum objective function and optimization target. More specifically, such as Figure 3 As shown in , when performing phenotypic simulation, the model after thermodynamic parameter correction, substrate concentration, objective function, reaction Gibbs free energy, metabolite concentration, B value, and total enzyme amount are input, and the pathway flux distribution results are output. Figure 4 As shown in , when conducting bottleneck reaction mining, the model, candidate reaction list, and bottleneck reaction are input; the model, B value, maximum value of objective function, candidate reaction list, objective function, and optimization target are input, and the bottleneck reaction information is output; isozyme coupling, the model and candidate reaction list are input, and the reaction coupling strategy is output; as shown in Figure 5 As shown in the continuous reaction coupling, the input model and candidate reaction list are used to output the reaction coupling strategy; the energy coupling, the input model and candidate reaction list are used to output the reaction coupling strategy; Figure 6For the limiting metabolite shown, input the model, B value, maximum value of the objective function, candidate reaction list, objective function, and optimization target, and output the limiting metabolite information; for the key enzyme, input the model, B value, maximum value of the objective function, candidate reaction list, objective function, and optimization target, and output the key enzyme information.

[0065] Specifically, the model analysis results are obtained, and the result visualization function is set accordingly, which can draw MDF line graphs, metabolic overflow graphs, cumulative distribution graphs, trade off graphs, phenotypic phase plane graphs, and pathway visualization graphs. These visualization functions require different parameters, involving models, input files, output files, coordinate axis settings, and color settings. More specifically, such as Figure 7 As shown, pathway visualization, input flux results, custom coordinates, flux value display, reaction settings, cofactor settings, node settings, output pathway diagram; MDF line graph, input model, MDF file, output MDF line graph; as shown Figure 8 For the cumulative distribution diagram shown in the figure, input the data file, coordinate axis settings, output file name, scale settings, and output the cumulative distribution diagram; for the phenotypic phase plane diagram, input the flux value, maximum oxygen concentration, maximum substrate value, and maximum biomass rate, and output the phenotypic phase plane diagram; Fig. 9 The metabolic overflow diagram shown includes input flux value, axis settings, substrate name, overflow maximum value, displayed metabolite list, color settings, biomass maximum value, overflow output file, and output phenotypic phase plane diagram.

[0066] Specifically, a method for generating a multi-constraint metabolic network model toolkit (MCM_Analyzer) is described as follows:

[0067] 1. Model modification

[0068] (1) Automated calibration algorithm for enzyme kinetic parameters

[0069] Code:

[0070] optimize_enzyme_parameters(model,expected_value,rounds_num,kcat_max_json_file)

[0071] Parameter Description:

[0072] Model: Models containing kinetic parameters, such as the enzyme constraint model (ecGEM) and the enzyme thermal constraint model (etcGEM)

[0073] expected_value: corrected expected value

[0074] rounds_num: cycle number threshold

[0075] kcat_max_json: JSON file containing the maximum value of kcat

[0076] Main function: In the while loop, the enzyme parameters are continuously optimized according to the expected value and the loop number threshold. First, the enzyme usage is calculated according to the current model, and the reaction that uses the most enzyme is selected. Then, according to the ECcode of the reaction, the corresponding kcat maximum value is found in the kcat_max dictionary and applied to the model. Finally, the loop counter is updated, the result of each operation is recorded, and the result is printed.

[0077] (2) Automatic correction algorithm for thermodynamic parameters

[0078] Code:

[0079] find_bottleneck_reaction(model,obj_target,B_value,max_workers)

[0080] Parameter Description:

[0081] Model: Models containing thermodynamic parameters, such as the thermodynamic constraint model (etGEM) and the enzyme thermal constraint model (etcGEM)

[0082] obj_target: generally maximize, that is, maximize the objective function

[0083] B_value: minimum thermodynamic driving force in the pathway

[0084] max_workers: specifies the maximum number of worker processes in the process pool, that is, the number of tasks that can be executed concurrently at the same time.

[0085] Main function: Check the model for thermodynamically overconstrained reactions and adjust their standard Gibbs free energies.

[0086] 2. Model Analysis and Visualization

[0087] (1) Metabolic network model analysis

[0088] from script.analysis.flux_balance_analysis import*

[0089] #Model parameter settings

[0090] model.add_medium("EX_glc__D_e_reverse",6)

[0091] model.obj_name="BIOMASS_Ec_iML1515_core_75p37M"

[0092] #Model calculation——The algorithm is repackaged on the original ETGEMs framework

[0093] from script.analysis.flux_balance_analysis import*

[0094] opt_gem=calculate_GEM(model,"maximize")

[0095] #Objective function flux result

[0096] opt_gem.obj()

[0097] #Get the flux results for a specific reaction

[0098] opt_gem.reaction["PDH"].value

[0099] #Calculation result formatting

[0100] from script.io.convert import*

[0101] flux_solution_df=result2flux_solution_df(model,opt_gem)

[0102] (2) Analysis of enzyme-constrained metabolic network model

[0103] The parameter settings and calculation results are the same as those of the metabolic network model calculation, with the following differences:

[0104] #Modify the total enzyme amount constraint

[0105] model.E_total = 0.3

[0106] #Model calculation

[0107] from script.analysis.flux_balance_analysis import*

[0108] opt_gem=calculate_ECM(model,"maximize")

[0109] (3) Analysis of the thermodynamic (enzyme heat) constrained metabolic network model

[0110] The parameter settings and calculation results are the same as those of the metabolic network model calculation, with the following differences:

[0111] #Other parameter settings, such as total enzyme amount, reaction Gibbs free energy, reaction metabolite concentration model.get_reaction_by_id('AIRC3_reverse').set_g0(0) #Reaction Gibbs free energy setting model.get_reaction_by_id('ATPS4rpp_num1').g0=0

[0112] opt_etm.metabolite["atp_c"].value=-3.91#Metabolite concentration setting

[0113] #When calculating the thermodynamic constraint model or enzyme thermal constraint model, it is necessary to calculate the B value to determine the minimum thermodynamic drive of the model - the new algorithm process for determining the inflection point is reduced from 500 fixed operations to more than a dozen times (one more than the actual inflection point)

[0114] Inflection point calculation process:

[0115] 1) Determine the maximum product yield without thermodynamic constraints using enzyme-constrained models and metabolic network models;

[0116] 2) Start with a minimum product output of 0.1 and calculate the MDF at that point. Then, calculate the maximum output at this MDF, which represents the first inflection point. Continue to add a very small value (0.02) to this output and repeat the above process until it exceeds the maximum value obtained in 1);

[0117] 3) Get the output and MDF of each point and draw a line chart, such as Fig.10 shown.

[0118] from script.analysis.max_driving_force import*

[0119] calculate_MDF_list(model,"ETM",MDF_result_file)

[0120] #Determine the inflection point value based on the MDF value

[0121] from script.plotting.plotting import*

[0122] Draw_MDF_By_Product_rate(MDF_dict,xmin,xlab,xmax,step,ytick,save_file)

[0123] #choose Fig.10 Calculate the B value for a certain value in the curve (Product rate in the figure).

[0124] B_value=MDF_Calculation(model,obj_value,obj_name)

[0125] #Model calculation

[0126] opt_tcm=calculate_TCM(model,"maximize",B_value)#Thermodynamic constraint model

[0127] opt_etm=calculate_ETM(model,"maximize",B_value)#enzyme thermal constraint model

[0128] (4) Visualization of calculation results

[0129] Code:

[0130] from script.cave.spav_function import*

[0131] flux2map(flux_solution_df,rxn_switch,coordinates,cofactor_switch,id_or_name,flow,col_equ,col_flux,html_file)

[0132] Parameter Description:

[0133] flux_solution_df: Pathway flux results

[0134] rxn_switch: If the value is 'T', the reaction is shown in the diagram

[0135] coordinates: The coordinate values ​​of the nodes in the graph. This can be empty (the backend code automatically calculates the node coordinates)

[0136] cofactor_switch: If the value is 'T', cofactors are shown in the plot

[0137] id_or_name: Whether the node is displayed as ID or Name

[0138] flow: Whether to display the flux value in the graph

[0139] col_equ: the column name corresponding to the reaction equation in flux_solution_df

[0140] col_flux: the column name corresponding to the flux value in flux_solution_df

[0141] html_file: Output HTML file format

[0142] The results are as follows Fig.11 shown.

[0143] (5) COBRA method rewrite

[0144] #pFBA——Cobrapy cobra.flux_analysis.pfba(model)

[0145] from script.analysis.flux_balance_analysis import*

[0146] Min_Flux_Sum_Calculation(model,model_type,obj_max_value,obj_id,B_value)

[0147] Parameter Description:

[0148] Model: Model of the new structure

[0149] model_type: can be metabolic network model (GEM), enzyme constraint model (ECM), thermodynamic constraint model (TCM) and enzyme thermal constraint model (ETM), etc.

[0150] obj_max_value: the maximum value of the objective function

[0151] obj_id: ID of the objective function in the model

[0152] B_value: For GEM and ECM, B_value can be empty

[0153] #FVA——Cobrapy: cobra.flux_analysis.flux_variability_analysis(model,fraction_of_optimum=0.9)

[0154] from script.analysis.flux_variability_analysis import*

[0155] flux_variability_analysis(model,model_type,obj_max_value,fraction_of_optimum,B_value,obj_id,max_workers)

[0156] Parameter Description:

[0157] Model: Model of the new structure

[0158] model_type: can be metabolic network model (GEM), enzyme constraint model (ECM), thermodynamic constraint model (TCM) and enzyme thermal constraint model (ETM), etc.

[0159] obj_max_value: the maximum value of the objective function

[0160] fraction_of_optimum: Setting the parameter fraction_of_optimum = 0.90 will give the flux range where the reaction is 90% optimal.

[0161] B_value: For GEM and ECM, B_value can be empty

[0162] obj_id: ID of the objective function in the model

[0163] max_workers: specifies the maximum number of worker processes in the process pool, that is, the number of tasks that can be executed concurrently at the same time.

[0164] 3. Calculation of metabolic engineering targets

[0165] (1) Bottleneck reaction mining: rewriting the BRA algorithm implemented in the original ETGEM framework

[0166] Algorithm flow:

[0167] 1) Specifically, a flux lower limit is set for the biomass synthesis reaction, and then the MDF level that the pathway can reach is determined. Subsequently, this MDF value is used as a constraint to solve for the maximum growth rate level. Then, the upper and lower limits of the thermodynamic driving force for each reaction in the model are calculated using both the maximum biomass synthesis rate and the MDF value as constraints. Finally, if both the upper and lower limits are equal to the MDF value, the reaction is identified as a bottleneck reaction.

[0168] 2) Given that bottleneck reactions are always involved in a pathway and their driving force (Dfi) is immutable, only reactions within the pathway output by the pFBA method need to be analyzed for variability, where the thermodynamic driving force value is equal to the MDF. This significantly reduces the size of the reaction list used for variability analysis, thereby shortening the time required to identify bottleneck reactions.

[0169] use_result=Get_Results_Thermodynamics(model,Concretemodel)

[0170] use_result_tmp=use_result[use_result['f']>-model.K_value]

[0171] reaction_result_select=use_result_tmp[abs(use_result_tmp['f']-B_value)<=1e-05]

[0172] # Reaction screening results such as Fig.12 shown.

[0173] #Determine bottleneck response

[0174] reaction_Df_list_fixed=Get_Max_Min_Df_Complete(model,obj_target,reac_list,B_value,max_value_obj,obj_id)

[0175] Bottleneck_reaction=reaction_Df_list_fixed[(reaction_Df_list_fixed['maximize']-B_value)<=0.01]

[0176] Parameter Description:

[0177] Model: Model of the new structure

[0178] obj_target: optimization target, such as 'maximize', 'minimize'

[0179] reac_list: list of reactions to be calculated

[0180] B_value: Maximum minimum thermodynamic driving force

[0181] max_value_obj: The maximum value of the objective function determined under the current MDF value

[0182] obj_id: The ID corresponding to the objective function in the model

[0183] (2) Searching for isozyme coupling strategies for bottleneck reactions - newly developed algorithms

[0184] Principle: Based on the identified bottleneck reactions, reactions associated with the same gene-protein-reaction (GPR) and with a standard Gibbs free energy less than 0 are searched in the multi-constraint model. These reactions may be potential targets of the same enzyme-coupled reaction.

[0185] #Determine isoenzymes

[0186] from script.analysis.coupling_strategy import*

[0187] Get_Fusion_Protein(Bottleneck_reaction,model)

[0188] Parameter Description:

[0189] Bottleneck_reaction: List of bottleneck reactions

[0190] Model: Model of the new structure

[0191] #Perform isozyme coupling analysis, taking the reactions catalyzed by the two isozymes "PGCD" and "AHGDx_reverse" as examples

[0192] model.merge_reaction("PGCD","AHGDx_reverse")

[0193] calculate_MDF_list(model,"ETM",MDF_result_file)

[0194] Draw_MDF_By_Product_rate(MDF_dict,0,xlab,xlab_max+1,1,ytick,png_file)

[0195] #After coupling, the simulation results of the model are compared Fig.13 Shown

[0196] (3) Finding a continuous reaction coupling strategy for bottleneck reactions - a newly developed algorithm

[0197] Principle: A continuous reaction is a situation where the product of one reaction is consumed by another reaction. Based on the identified bottleneck reaction, search for continuous reactions with a standard Gibbs free energy less than 0.

[0198] from script.analysis.continuous_reaction import*

[0199] Get_Continuous_Reaction(Bottleneck_reaction_list,model)

[0200] Parameter Description:

[0201] Bottleneck_reaction: List of bottleneck reactions

[0202] Model: Model of the new structure

[0203] #Reaction coupling and the results after coupling show the same isozyme coupling strategy

[0204] (4) Finding energy coupling strategies for bottleneck reactions - newly developed algorithms

[0205] Principle: The energy coupling strategy refers to the situation where energy-related metabolites (such as ATP, NADPH, etc.) produced in one reaction are consumed by another reaction. Based on the identified bottleneck reaction, search for energy coupling reactions with a standard Gibbs free energy less than 0.

[0206] from script.analysis.continuous_reaction import*

[0207] Find_Energy_Coupling_Strategy(Bottleneck_reaction_list,model)

[0208] Parameter Description:

[0209] Bottleneck_reaction: List of bottleneck reactions

[0210] Model: Model of the new structure

[0211] #Reaction coupling and the results after coupling show the same isozyme coupling strategy

[0212] (5) Restricted metabolite mining: rewriting the LMA algorithm implemented in the original ETGEM framework

[0213] Principle: Limiting metabolites are those whose concentrations reach a set boundary or have a dual role as substrate and product in different bottleneck reactions. This dual role results in their concentration being limited, eliminating variability. Therefore, there is an optimal concentration level. Since limiting metabolites are bound to participate in the bottleneck reaction, a list of metabolites participating in the bottleneck reaction can be compiled. Subsequently, the concentration variation range of these metabolites is calculated, and metabolites with unchanged concentrations are identified as limiting metabolites.

[0214] from script.analysis.metabolite_concentration_variability_analysisimport*

[0215] Get_Max_Min_Met_Concentration_for_list(model,B_value,max_value_obj,obj_id,Bottleneck_reaction,Bottleneck_reaction_met,max_workers)

[0216] Parameter Description:

[0217] Model: Model of the new structure

[0218] B_value: Maximum minimum thermodynamic driving force

[0219] max_value_obj: The maximum value of the objective function determined under the current MDF value

[0220] obj_id: The ID corresponding to the objective function in the model

[0221] Bottleneck_reaction: List of bottleneck reactions

[0222] Bottleneck_reaction_met: List of metabolites involved in the bottleneck reaction

[0223] max_workers: specifies the maximum number of worker processes in the process pool, that is, the number of tasks that can be executed concurrently at the same time.

[0224] (6) Key enzyme mining: rewriting the KEA algorithm implemented in the original ETGEM framework

[0225] Principle: In the comprehensive enzyme constraint model, for reactions with higher enzyme costs, they can be considered to exert stronger control on the pathway, making them the key enzymes in the pathway. Therefore, by determining the growth and MDF levels, the minimum enzyme cost required for each reaction can be calculated. Higher values ​​indicate a more critical role. Generally, in enzyme constraint models, when enzyme amount constraints are in effect, the variation in enzyme levels in pathway reactions becomes smaller or even non-existent. In contrast, when enzyme amount constraints are not met, the enzyme costs of reactions in the pathway will be more variable. Therefore, minimizing enzyme costs makes sense regardless of whether the enzyme amount constraints are met, while maximizing enzyme costs may be less meaningful.

[0226] from script.analysis.enzyme_cost_variability_analysis import*

[0227] max_Max_Min_E_for_list(model,B_value,obj_id,max_value_obj,enz_use_reaction_list,max_workers)

[0228] Parameter Description:

[0229] Model: Model of the new structure

[0230] B_value: Maximum minimum thermodynamic driving force

[0231] max_value_obj: The maximum value of the objective function determined under the current MDF value

[0232] obj_id: The ID corresponding to the objective function in the model

[0233] enz_use_reaction_list: List of reactions that need to be evaluated

[0234] max_workers: specifies the maximum number of worker processes in the process pool, that is, the number of tasks that can be executed concurrently at the same time.

[0235] 4. Application Cases:

[0236] Prediction of serine biosynthesis pathway using EcoETM - bottleneck reaction mining and isozyme coupling mining

[0237] Using EcoETM (an enzyme thermal constraint model for Escherichia coli), we analyzed the tryptophan synthesis pathway using the toolkit and identified five distributed bottleneck reactions (PGCD, PGK_reverse, GAPD, FBA, and TPI), whose joint participation made the pathway thermodynamically infeasible. Among them, PGCD is a well-known reaction in the L-serine synthesis pathway, and the other four reactions are also core to the glycolysis pathway. This subset is suitable as an example to illustrate the concept of distributed bottleneck reactions, but the conclusion that the classical EMP pathway cannot coexist with the conventional L-serine synthesis process is inconsistent with the observations. Therefore, we analyzed the L-serine synthesis pathway in detail.

[0238] We found that the maximum thermodynamic driving force (MDF) of the L-serine biosynthesis pathway is variable, with 13 turning points, such as Fig.14 As shown in the figure, the serine synthesis process is divided into 13 stages according to the MDF level. The MDF level of the first seven stages can reach more than 4 kJ / mol, and the thermodynamic feasibility is ideal. The MDF level of the middle four stages is close to zero, and the feasibility is quite low. The MDF of the last two stages is negative, so they can be considered thermodynamically infeasible. Further analysis shows that there are a total of 2 local bottleneck reactions and 11 distributed bottleneck reactions, representing a combination of reactions located in central metabolic pathways (mainly in EMP, PPP and TCA). For the yield index, the highest flux level that can be achieved at the end of the seventh stage is 8.99 mmol / gDW / h, which is only 42.7% compared with the highest flux of 21.03 mmol / gDW / h, which is not satisfactory because it represents a 57% carbon loss. Starting from the eighth stage, the MDF dropped to 1.57 kJ / mol, which was caused by the local bottleneck reaction of PGCD alone. Although the participation of this reaction reduces the MDF level of the pathway, it makes a significant contribution to the increase in yield, thus bringing the synthetic flux of the pathway to 17.175mmol / gDW / h, equivalent to 81.7% of the maximum yield. From stage 9 to stage 12, the thermodynamic feasibility of the pathway gradually decreases, which is caused by several distributed bottleneck reactions. The common point of these distributed bottleneck reaction combinations is that the PGCD reaction always shares the limiting metabolite with other bottleneck reactions, and the overall thermodynamic feasibility of the pathway is lower than that of stage 8 (when PGCD is a local bottleneck reaction). According to the second law of thermodynamics (MDF ≥ 0), the flux limit of the last thermodynamically feasible pathway at the end of stage 11 is 20.66mmol / gDW / h, equivalent to 98.2% of the maximum yield. In stage 13, there are still three distributed bottleneck reactions, but PGCD is not among them. By investigating the individual thermodynamic driving force levels of these three reactions, it is found that they are all thermodynamically feasible in principle. However, due to the close connection in sharing limiting metabolites, their synchronous activities can make the pathway thermodynamically infeasible.

[0239] According to these results, it can be found that the PGCD reaction has a great influence on the L-serine yield. With its participation, the maximum yield of the pathway can be increased from 42.7% of the theoretical yield predicted by the initial iML1515 to 98.2%. According to the MetaCyc database, the PGCD reaction should indeed appear in the L-serine synthesis pathway of Escherichia coli, which has been confirmed by experimental results. In the combination of bottleneck reactions in stages 9-13, the distributed bottleneck reactions with or without PGCD are always located outside the core L-serine synthesis pathway. This suggests that when specific pathways are placed in a cellular environment, it may be inappropriate to perform thermodynamic analysis and optimization on these pathways alone, which reflects the difference between the thermodynamic evaluation of the pathway MDF solution and the specific pathway, indicating that the thermodynamic evaluation of the preset pathway is likely to misjudge the actual situation of the complex intracellular system.

[0240] The PGCD reaction is very important for producing high yields of L-serine, but its unfavorable thermodynamics is also a problem that must be solved. Fig.15 As shown, the PGCD reaction will release the reducing power of NADH. (In 2017, Zhang et al. demonstrated that when PGCD is coupled with a reducing power consumption reaction catalyzed by the same enzyme expressed by the same gene, its thermodynamic feasibility is significantly improved)

[0241] In the iML1515 model, there are two potential natural coupling reactions associated with PGCD, AHGDx and ARHGDx ( Fig.15 Part A in the figure). The above three reactions have the same gene number (b2913) and therefore encode multifunctional enzymes. In the initial iML1515 model, the AHGDx and ARHGDx reactions are reversible, and the only difference between them is the configuration of the substrate. However, in the entire metabolic network, there is no configuration conversion process between their substrates, and there is no consumption or silencing reaction with the R-configuration substrate (R-2-hydroxyglutamate, r2hglut), so ARHGDx is redundant and cannot participate in the flux balance analysis. The S-configuration substrate (S-2-hydroxyglutamate, s2hglut) of the AHGDx reaction can only be used as a substrate for another irreversible reaction SHGO, that is, there is no other way to consume it. Therefore, it can be considered that the forward reaction direction of AHGDx in the model is also redundant, and only its reverse reaction is valid. The AHGDx_reverse reaction consumes reducing power, so when it is coupled with the SHGO reaction, the generation and consumption of carbon compounds can also be offset. This means that when the combined reaction continues to couple with PGCD, the thermodynamic feasibility of the overall reaction can be promoted by balancing the reducing force, and the metabolite α-ketoglutarate (α-KG) can also be supplemented, such as Fig.15 As shown in part A.

[0242] We added the combined reaction catalyzed by SerA with a ΔrGi'0 of 2.1 kJ / mol to the model and set the equation to "3pg_c+akg_c-->3php_c+S2hglut_c" to replace the original two partial reactions PGCD and AHGDx_reverse. Since the catalytic rate of the enzyme is limited by the step with the weakest kinetic parameter, the lower kcat value of the AHGDx_reverse reaction is due to the combined reaction of SerA, such as Fig.15 After the reaction combination, we recalculated the change of MDF in the L-serine synthesis process (as shown in Fig.15 (as shown in part B of the figure), it was found that the thermodynamic driving force of the pathway was significantly improved. At the stage where the original thermodynamic feasibility reached a good level (left side of the red dashed line), the achievable flux increased from 8.99 to 18.77 mmol / gDW / h. In addition, after correction by merging reactions, the thermodynamic driving force in the entire yield solution space was optimized. Due to the dual adjustment of reaction stoichiometry and enzyme kinetic parameters, the yield solution space was reduced by 9.1% (relative to the original thermodynamically feasible yield solution space, blue shaded area) to 10.7% (relative to the original total yield solution space, the entire shaded area), as shown in Fig.15 This is shown in part B. In addition, the number of MDF levels was reduced and the pathway was divided into only five stages.

[0243] The above embodiments are preferred implementation schemes of the present invention. In addition, the present invention may also be implemented in other ways. Any obvious replacement without departing from the concept of the present invention is within the protection scope of the present invention.

Claims

1. A method for generating a multi-constraint metabolic network model toolkit, characterized in that: The following steps are involved: Construct multi-constraint metabolic network models; Performing parameter correction on the multi-constraint metabolic network model to obtain a model with corrected parameters; Using the model after parameter correction, performing model analysis to obtain model analysis results; The multi-constraint metabolic network model includes an enzyme constraint model, an enzyme thermal constraint model, and a thermodynamic constraint model; The parameter correction includes enzyme kinetic parameter correction and thermodynamic parameter correction.

2. The generation method according to claim 1, characterized in that: The model analysis includes phenotypic simulation analysis and metabolic engineering target prediction analysis; The metabolic engineering target prediction analysis includes bottleneck reaction mining, isozyme coupling mining, continuous reaction coupling mining, energy coupling mining, limiting metabolite mining, and key enzyme mining.

3. The generation method according to claim 1, characterized in that: The model analysis result is to visualize the result; The result visualization includes MDF line graph, cumulative distribution graph, phenotypic phase plane graph, metabolic overflow graph, trade off graph, and pathway visualization.

4. The generation method according to claim 1, characterized in that: The thermodynamic parameter correction includes kinetic parameter correction and thermodynamic parameter correction.

5. The generation method according to claim 4, characterized in that: The kinetic parameter correction includes using an enzyme kinetic parameter automatic correction algorithm, the process is: First, the enzyme usage is calculated according to the current model, and the reaction that uses the most enzyme is selected; Then, according to the EC code of the reaction, find the corresponding kcat maximum value in the kcat_max dictionary and apply it to the model; Finally, update the loop counter, record the result of each operation, and print the result.

6. The generation method according to claim 1, characterized in that: The generation method according to claim 4 is characterized in that the correction of the kinetic parameters includes using an automated correction algorithm for thermodynamic parameters to perform correction.

7. The generation method according to claim 2, characterized in that: The phenotypic simulation analysis includes growth rate simulation and product simulation.

8. The generation method according to claim 2, characterized in that: The algorithm for bottleneck reaction mining has the following process: Set a lower limit on the flux of the biomass synthesis reaction and then determine the MDF level that the pathway can achieve; This MDF value was then used as a constraint to solve for the maximum growth rate level; Then, the maximum biomass synthesis rate and MDF value were used as constraints to calculate the upper and lower limits of the thermodynamic driving force of each reaction in the model; Finally, a reaction was identified as a bottleneck reaction if both the upper and lower limits were equal to the MDF value; where the variability of the reaction was analyzed within the pathway output by the pFBA method and the thermodynamic driving force value was equal to the MDF.

9. The generation method according to claim 8, characterized in that: The isozyme coupling is sought for the bottleneck reaction, and the process is as follows: based on the identified bottleneck reaction, a reaction associated with the same gene-protein-reaction (GPR) and having a standard Gibbs free energy less than 0 is searched in a multi-constraint model.

10. An application of a multi-constraint metabolic network model toolkit, characterized in that: The toolkit generated by the generation method according to any one of claims 1 to 9 is applied to predict the bottlenecks of cell metabolism and key regulatory modification targets.

Images