Simulations
For each simulation step, the function receives a metabolic model and performs:
update fluxes by metabolites concentrations
update fluxes by coregnet and influence value
update fluxes by gene state from the GRN simulator
The simulation result is a list containing:
objective_history: time series of objective function value for the linear program
metabolites: metabolites concentrations over time
fluxes_history: time series of the fluxes values for all the time series
metabolites_concentration_history: time series of metabolite concentrations
metabolites_fluxes_history: time series of the metabolites fluxes during the simulation
rate_history: time series of the growth rate values for all simulation
time: vector containing the simulation times
gene_state_history: list containing the values for the gene state during the simulation
The fluxes for the simulation time are stored in a matrix which row names are the fluxes reaction id.
data("aliases_SC")
data("iMM904")
metabolites<-data.frame("names" = c("D-Glucose","Ethanol"),
"concentrations" = c(16.6,0))
Simulation1<-Simulation(model = iMM904,
time = seq(1,20,by = 1),
metabolites = metabolites,
initial_biomass = 0.45,
aliases = aliases_SC)
# Default biomass flux index use is 1577 corresponding to Biomass SC5 notrace
# Joining by: metabolites_id
# simulation step
# simulation step
# simulation step
# simulation step
# simulation step
# simulation step
# simulation step
# simulation step
# simulation step
# simulation step
# simulation step
# simulation step
# simulation step
# simulation step
# simulation step
# simulation step
# simulation step
# simulation step
# simulation step
Simulation1$fluxes_history[1:10,1:5]
# [,1] [,2] [,3] [,4] [,5]
# 13BGH 0.00000 0.00000 0.0000000 0.00000000 0.00000000
# 13BGHe 0.00000 0.00000 0.0000000 0.00000000 0.00000000
# 13GS 0.32667 0.32667 0.2063497 0.03309391 0.03309391
# 16GS 0.00000 0.00000 0.0000000 0.00000000 0.00000000
# 23CAPPD 0.00000 0.00000 0.0000000 0.00000000 0.00000000
# 2DDA7Ptm 0.00000 0.00000 0.0000000 0.00000000 0.00000000
# 2DHPtm 0.00000 0.00000 0.0000000 0.00000000 0.00000000
# 2DOXG6PP 0.00000 0.00000 0.0000000 0.00000000 0.00000000
# 2HBO 0.00000 0.00000 0.0000000 0.00000000 0.00000000
# 2HBt2 0.00000 0.00000 0.0000000 0.00000000 0.00000000
To have access to the gprRules users can use the sybil package, which returns a vector of size equal to the number of fluxes and the associated genes.
If you only wish to know which gene affects which reaction; the sybil objects have a slot for obtaining the flux-gene matrix.
rxnGeneMat(iMM904)[1:10,1:10]
# 10 x 10 sparse Matrix of class "lgCMatrix"
#
# [1,] | . . . . . . . . .
# [2,] . | | | . . . . . .
# [3,] . . . . | | | | | |
# [4,] . . . . . . . . . .
# [5,] . . . . . . . . . .
# [6,] . . . . . . . . . .
# [7,] . . . . . . . . . .
# [8,] . . . . . . . . . .
# [9,] . . . . . . . . . .
# [10,] . . . . . . . . . .
Simulate a dFBA over time (here 20h) without constraint
metabolites<-data.frame("names" = c("D-Glucose","Ethanol"),
"concentrations" = c(16.6,0))
Simulation1<-Simulation(model = iMM904,
time = seq(1,20,by = 1),
metabolites = metabolites,
initial_biomass = 0.45,
aliases = aliases_SC)
Simulation1$biomass_history
# [1] 0.4500000 0.6001102 0.8002939 0.9598884 0.9882935 1.0175393 1.0476505
# [8] 1.0786527 1.1105724 1.1434366 1.1772734 1.2121114 1.2479804 1.2849109
# [15] 1.3229342 1.3620827 1.4023896 1.4175428 1.4175428 1.4175428
Simulation1$met_concentration_history
# [,1] [,2] [,3] [,4] [,5] [,6] [,7]
# D-Glucose 16.6 11.385409 4.431344 0.00000 0.00000 0.00000 0.00000
# Ethanol 0.0 8.247123 19.245307 26.42949 24.77683 23.07527 21.32335
# [,8] [,9] [,10] [,11] [,12] [,13] [,14]
# D-Glucose 0.00000 0.00000 0.00000 0.00000 0.00000 0.000000 0.000000
# Ethanol 19.51959 17.66245 15.75036 13.78168 11.75474 9.667829 7.519157
# [,15] [,16] [,17] [,18] [,19] [,20]
# D-Glucose 0.000000 0.000000 0.0000000 0 0 0
# Ethanol 5.306901 3.029179 0.6840547 0 0 0
Constraining the model
When provided with different kind of constraints, CoRegFlux process the given information in the following order:
- Gene expression is first integrated
- TF KO or OV is carried out, starting with the first line in the regulator table and going along the rows.
- Gene KO or OV is carried out. starting with the first line in the gene table and going along the rows. Thus order in the regulator table and gene table might play a role and potentially give different results.
The different functions used by CoRegFlux to constraint the model are individually accessible to allow the combination of CoRegFlux’s models with other algorithms and parameters, provided by sybil for instance. A model can be constrain iteratively through the different function. In that case, the recommended order is as follows: uptake constraint, gene expression, TF KO or OV, gene KO or OV.
Simulate a dFBA with gene expression as a constraint
Simulation2<-Simulation(model = iMM904,
time = seq(1,20,by = 1),
metabolites = metabolites,
initial_biomass = 0.45,
aliases = aliases_SC,
gene_state_function = function(a,b){GeneState})
Simulation2$biomass_history
# [1] 0.4500000 0.5248670 0.6121898 0.7140405 0.8328362 0.9713961 1.0679966
# [8] 1.0996010 1.1321406 1.1656430 1.2001369 1.2356516 1.2722172 1.3098649
# [15] 1.3486266 1.3866631 1.3866631 1.3866631 1.3866631 1.3866631
Simulation2$met_concentration_history
# [,1] [,2] [,3] [,4] [,5] [,6] [,7]
# D-Glucose 16.6 14.477830 12.002591 9.115544 5.748176 1.820576 0.00000
# Ethanol 0.0 3.005114 6.510192 10.598414 15.366798 20.928504 22.88805
# [,8] [,9] [,10] [,11] [,12] [,13] [,14]
# D-Glucose 0.00000 0.00000 0.00000 0.00000 0.00000 0.000000 0.000000
# Ethanol 21.04926 19.15605 17.20682 15.19991 13.13361 7.204882 5.014481
# [,15] [,16] [,17] [,18] [,19] [,20]
# D-Glucose 0.000000 0 0 0 0 0
# Ethanol 1.416063 0 0 0 0 0
Simulate a dFBA with TF knock-out (KO) while constraining the model with gene expression
If the simulated mutant have several TFs KO or OV, CoRegFlux will constrain the model according to the order of the TFs in the regulator table. While this example also constraint the model with gene expression, it is possible to run the simulation without such constraints.
regulator_table <- data.frame("regulator" = "MET32",
"influence" = -1.20322,
"expression" = 0,
stringsAsFactors = FALSE)
SimulationTFKO<-Simulation(model = iMM904,
time = seq(1,20,by = 1),
metabolites = metabolites,
initial_biomass = 0.45,
aliases = aliases_SC,
coregnet = SC_GRN_1,
regulator_table = regulator_table ,
gene_state_function = function(a,b){GeneState})
SimulationTFKO$biomass_history ## This KO is predicted as non-lethal
# [1] 0.4500000 0.5248670 0.6121898 0.7140405 0.8328362 0.9713961 1.0679966
# [8] 1.0996010 1.1321406 1.1656430 1.2001369 1.2356516 1.2722172 1.3098649
# [15] 1.3332605 1.3332605 1.3332605 1.3332605 1.3332605 1.3332605
Simulate a dFBA with TF over-expression (OV) while constraining the model with gene expression
If the simulated mutant have several TFs KO or OV, CoRegFlux will constrain the model according to the order of the TFs in the regulator table. While this example also constraint the model with gene expression, it is possible to run the simulation without such constraints.
regulator_table <- data.frame("regulator" = "MET32",
"influence" = -1.20322 ,
"expression" = 3,
stringsAsFactors = FALSE)
SimulationTFOV<-Simulation(model = iMM904,
time = seq(1,20,by = 1),
metabolites = metabolites,
initial_biomass = 0.45,
aliases = aliases_SC,
coregnet = SC_GRN_1,
regulator_table = regulator_table,
gene_state_function = function(a,b){GeneState})
SimulationTFOV$biomass_history ## This OV is predicted as non-lethal
# [1] 0.4500000 0.5248670 0.6121898 0.7140405 0.8198573 0.8441187 0.8690980
# [8] 0.8948165 0.9212961 0.9485592 0.9766291 0.9766291 0.9766291 0.9766291
# [15] 0.9766291 0.9766291 0.9766291 0.9766291 0.9766291 0.9766291
Simulate a dFBA with gene(s) knock-out or over-expression simulation while constraining the model with gene expression
If the simulated mutant have several gene KO or gene OV, CoRegFlux will constrain the model according to the order of the genes in the gene table. While this example also constraint the model with gene expression, it is possible to run the simulation without such constraints.
gene_table <- data.frame("gene" = c("YJL026W","YIL162W"),
"expression" =c(2,0),
stringsAsFactors = FALSE)
SimulationGeneKO_OV<-Simulation(model = iMM904,
time = seq(1,20,by = 1),
metabolites = metabolites,
initial_biomass = 0.45,
aliases = aliases_SC,
coregnet = SC_GRN_1,
gene_table = gene_table,
gene_state_function = function(a,b){GeneState})
SimulationGeneKO_OV$biomass_history ## This OV is predicted as non-lethal
# [1] 0.4500000 0.5248670 0.6121898 0.7140405 0.8075046 0.8314004 0.8560033
# [8] 0.8813343 0.9074149 0.9342673 0.9619142 0.9903794 1.0196868 1.0488524
# [15] 1.0488524 1.0488524 1.0488524 1.0488524 1.0488524 1.0488524
Constraining the model according to gene expression, TF KO or OV, gene KO or OV to run various FBA using sybil
The different functions used by CoRegFlux to constraint the model are individually accessible to allow the combination of CoRegFlux’s models with other algorithms and parameters, provided by sybil for instance. A model can be constrain iteratively through the different function. In that case, the recommended order is as follows: uptake constraint, gene expression, TF KO or OV, gene KO or OV.
metabolites_rates <- data.frame("name"=c("D-Glucose"),
"concentrations"=c(16.6),
"rates"=c(-2.81))
model_uptake_constraints <- adjust_constraints_to_observed_rates(model = iMM904,
metabolites_with_rates = metabolites_rates)
model_gene_constraints <- coregflux_static(model= iMM904,
predicted_gene_expression =
PredictedGeneState,
aliases = aliases_SC)$model
model_TF_KO_OV_constraints <- update_fluxes_constraints_influence(model= iMM904,
coregnet = SC_GRN_1,
regulator_table = regulator_table,
aliases = aliases_SC )
model_gene_KO_OV_constraints <- update_fluxes_constraints_geneKOOV(
model= iMM904,
gene_table = gene_table,
aliases = aliases_SC)
sol <- sybil::optimizeProb(model_TF_KO_OV_constraints)
sol
# solver: glpkAPI
# method: simplex
# algorithm: fba
# number of variables: 1577
# number of constraints: 1226
# return value of solver: solution process was successful
# solution status: solution is optimal
# value of objective function (fba): 0.287866
# value of objective function (model): 0.287866
From observations to fluxes
Here we will compute the fluxes from the observed growth rates (which can be obtained directly from the growth curves)
Assuming we have an observed growth rate of 0.3
Given that the fba solution is not unique, if you wish to see the intervals of maximum and minimum allowed fluxes for a reaction, flux variability analysis should be used
It worth noting that none of the two methods guarantee that the observed growth rate will be reached.
This could mean that the uptake rates for the limiting substrates (most commonly glucose uptake rate) does not allow for higher growth.
To constraint the model using the substrate uptake rate, the user must also provide the metabolites_rates argument
Calibration: identifying the softplus parameter using bayesian optimization
To translate the gene expression to fluxes in the GEM, CoRegFlux use the softplus function
where \(\theta\) is the softplus parameter applied to all fluxes, \(gpr_{i}\left(X\right)\) is the result of evaluating the gene-protein-reaction rules for a set of gene expression levels of the metabolic genes \(X\). These rules relate genes to reactions and are logical form. CoRegflux transform these rules as follows
- AND are substitued by MIN()
- OR are substituted by SUM().
Given a known growth rate and predicted gene expressions obtained through the function predict_linear_model_influence, the users have the possibility to adjust the softplus parameter \(\theta\) to calibrate the integration of the gene expression in the GEM. This step requires the installation of the package rBayesianOptimization.
library(rBayesianOptimization)
gRates <- 0.1
opF<-function(p){
CoRegFlux_model<-coregflux_static(model = model_uptake_constraints,
gene_parameter = p,
predicted_gene_expression =
PredictedGeneState)
ts<-optimizeProb(CoRegFlux_model$model)
list(Score=-1*log(abs(lp_obj(ts)-gRates)/gRates),Pred=0)
}
result<-BayesianOptimization(FUN = opF,
bounds = list(p = c(-10,10)),
data.frame(p = seq(-10,10,by = 0.5)),
n_iter = 10,
verbose = TRUE)
