diff --git a/1-Model_construction.R b/1-Model_construction.R new file mode 100644 index 0000000..0e929b3 --- /dev/null +++ b/1-Model_construction.R @@ -0,0 +1,469 @@ + + +################### +# Set environment # +################### + +# load libraries and initialize environment +source("set_env.R") + + +###################### +# Model construction # +###################### + +# Test whether the different models (1-4) fit the data accurately, based on chi2 statistics + +# nb_points: number of experimental data points (152) +# nb_par: number of free parameters + +# Go to models directory +setwd(model_dir) + +# Model 1: no inhibition by acetate +loadModel("Millard2020_Ecoli_glc_ace_model_1.cps") +res_PE_model1 <- runParameterEstimation(method = "Statistics", model = getCurrentModel()) +res_PE_model1$chi2 <- test_khi2(nb_points=152, k_val=res_PE_model1$main$objective_value, nb_par=8) + +print("chi2-test for model 1:") +print(paste(" p-value = ", res_PE_model1$chi2$`p-value, i.e. P(X^2<=value)`, sep="")) +print(paste(" conclusion: ", res_PE_model1$chi2$conclusion, sep="")) + +# Model 2: inhibition of glycolysis +loadModel("Millard2020_Ecoli_glc_ace_model_2.cps") +res_PE_model2 <- runParameterEstimation(method = "Statistics", model = getCurrentModel()) +res_PE_model2$chi2 <- test_khi2(nb_points=152, k_val=res_PE_model2$main$objective_value, nb_par=9) + +print("chi2-test for model 2:") +print(paste(" p-value = ", res_PE_model2$chi2$`p-value, i.e. P(X^2<=value)`, sep="")) +print(paste(" conclusion: ", res_PE_model2$chi2$conclusion, sep="")) + +# Model 3: inhibition of TCA cycle +loadModel("Millard2020_Ecoli_glc_ace_model_3.cps") +res_PE_model3 <- runParameterEstimation(method = "Statistics", model = getCurrentModel()) +res_PE_model3$chi2 <- test_khi2(nb_points=152, k_val=res_PE_model3$main$objective_value, nb_par=9) + +print("chi2-test for model 3:") +print(paste(" p-value = ", res_PE_model3$chi2$`p-value, i.e. P(X^2<=value)`, sep="")) +print(paste(" conclusion: ", res_PE_model3$chi2$conclusion, sep="")) + +# Model 4: inhibition of TCA and glycolysis +loadModel("Millard2020_Ecoli_glc_ace_model_4.cps") +res_PE_model4 <- runParameterEstimation(method = "Statistics", model = getCurrentModel()) +res_PE_model4$chi2 <- test_khi2(nb_points=152, k_val=res_PE_model4$main$objective_value, nb_par=10) + +print("chi2-test for model 4:") +print(paste(" p-value = ", res_PE_model4$chi2$`p-value, i.e. P(X^2<=value)`, sep="")) +print(paste(" conclusion: ", res_PE_model4$chi2$conclusion, sep="")) + +# plot model evaluation results +setwd(results_dir) +pdf(file="Figure 1B.pdf", width = 5, height = 5) + +barplot(c(res_PE_model1$main$objective_value, + res_PE_model2$main$objective_value, + res_PE_model3$main$objective_value, + res_PE_model4$main$objective_value), + las=1, + ylab="variance-weighted SSR", + names.arg=c("model1", "model2", "model3", "model4")) +chi2_threshold <- 171 +abline(h=chi2_threshold) + +dev.off() + + +########### +# PLOT FIT vs MEASUREMENTS # +########### + +# In current version of CoRc package, custom weights cannot be provided when defining an experiment, +# hence we cannot define weight according to the experimental standard deviations. +# To overcome this limitation, we use a copasi model file with weights pre-defined in the +# Parameter estimation task. + +# MODEL 1 # +########### + +setwd(model_dir) +loadModel("Millard2020_Ecoli_glc_ace_model_1.cps") + +# get parameters for the best fit +res_PE_model <- runParameterEstimation(method = "Statistics", model = getCurrentModel()) + +# list used to set up global sensitivity analysis and store calculation results +fit_results_m1 <- list("1" = list("initial_concentrations"= c("Ace0_out"=0.92, "Glc"=12.9, "X"=0.0671), + "species_idx" = c("[Ace0_out]_0"=1, "[Glc]_0"=1, "[X]_0"=1), + "data_exp" = read.table("data_1mM.txt", header=TRUE, sep="\t"), + "sd" = c("Ace"=0.2, "Glc"=0.5, "X"=0.045, "Ace_enr"=0.03), + "mapping" = c("Ace"="Acet_out", "Glc"="Glc", "X"="X", "Ace_enr"="Values[Ace_enr]")), + "10" = list("initial_concentrations"= c("Ace0_out"=9, "Glc"=14.192, "X"=0.07), + "species_idx" = c("[Ace0_out]_0"=2, "[Glc]_0"=2, "[X]_0"=2), + "data_exp" = read.table("data_10mM.txt", header=TRUE, sep="\t"), + "sd" = c("Ace"=0.2, "Glc"=0.5, "X"=0.045, "Ace_enr"=0.03), + "mapping" = c("Ace"="Acet_out", "Glc"="Glc", "X"="X", "Ace_enr"="Values[Ace_enr]")), + "30" = list("initial_concentrations"= c("Ace0_out"=30.2448, "Glc"=12.5, "X"=0.07227), + "species_idx" = c("[Ace0_out]_0"=3, "[Glc]_0"=3, "[X]_0"=3), + "data_exp" = read.table("data_30mM.txt", header=TRUE, sep="\t"), + "sd" = c("Ace"=0.2, "Glc"=0.5, "X"=0.045, "Ace_enr"=0.03), + "mapping" = c("Ace"="Acet_out", "Glc"="Glc", "X"="X", "Ace_enr"="Values[Ace_enr]")), + "res_cost" = c(res_PE_model$main$objective_value), + "res_converged" = c(test_khi2(nb_points=152, k_val=res_PE_model$main$objective_value, nb_par=10)$`p-value, i.e. P(X^2<=value)` < 0.95)) + +# simulate data for the best fit +for (i in c("1", "10", "30")){ + # update initial conditions (Glc, Ace and biomass concentrations) + for (j in names(fit_results_m1[[i]][["initial_concentrations"]])){ + setSpecies(key=j, initial_concentration=fit_results_m1[[i]][["initial_concentrations"]][j]) + } + # apply initial state + applyInitialState() + # run time course simulation + tmp <- as.matrix(runTimeCourse()$result) + # get complete simulation results + idx_t_max <- (tmp[,"Time"] <= max(fit_results_m1[[i]][["data_exp"]]$time)) + fit_results_m1[[i]]$simulations <- tmp[idx_t_max,] + # get index of simulated data at experimental time points + fit_results_m1[[i]]$t_idx <- match(fit_results_m1[[i]][["data_exp"]]$time, fit_results_m1[[i]]$simulations[,"Time"]) + # get simulated data that match experimental time points + sim_data <- fit_results_m1[[i]][["data_exp"]] + for (k in names(fit_results_m1[[i]][["mapping"]])){ + sim_data[,k] <- fit_results_m1[[i]]$simulations[fit_results_m1[[i]]$t_idx,fit_results_m1[[i]][["mapping"]][[k]]] + } + fit_results_m1[[i]]$sim <- sim_data +} + +# Plot fitting results +setwd(results_dir) +pdf(file = "Figure 1-figure supplement 1.pdf", width = 7, height = 9) +par(mfrow=c(4,3)) + +plot_no_ci(fit_results_m1, "1", "X", col="#1F4E79", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,1), xaxs="i", yaxs="i", ylab="[biomass] (gDW/L)", xlab="", main="[acetate]=1mM") +plot_no_ci(fit_results_m1, "10", "X", col="#2E75B6", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,1), xaxs="i", yaxs="i", ylab="", xlab="", main="[acetate]=10mM") +plot_no_ci(fit_results_m1, "30", "X", col="#9DC3E6", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,1), xaxs="i", yaxs="i", ylab="", xlab="", main="[acetate]=30mM") + +plot_no_ci(fit_results_m1, "1", "Glc", col="#548235", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,15), xaxs="i", yaxs="i", ylab="[Glc] (mM)", xlab="", main="") +plot_no_ci(fit_results_m1, "10", "Glc", col="#70AD47", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,15), xaxs="i", yaxs="i", ylab="", xlab="", main="") +plot_no_ci(fit_results_m1, "30", "Glc", col="#A9D18E", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,15), xaxs="i", yaxs="i", ylab="", xlab="", main="") + +plot_no_ci(fit_results_m1, "1", "Ace", col="#B63A2D", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,5), xaxs="i", yaxs="i", ylab="[Ace] (mM)", xlab="", main="") +plot_no_ci(fit_results_m1, "10", "Ace", col="#D6685C", lwd=1.0, las=1, xlim=c(0,5), ylim=c(8,11), xaxs="i", yaxs="i", ylab="", xlab="", main="") +plot_no_ci(fit_results_m1, "30", "Ace", col="#E49890", lwd=1.0, las=1, xlim=c(0,5), ylim=c(25,32), xaxs="i", yaxs="i", ylab="", xlab="", main="") + +plot_no_ci(fit_results_m1, "1", "Ace_enr", col="#BF9000", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,1), xaxs="i", yaxs="i", ylab="[Ace] (mM)", xlab="time (h)", main="") +plot_no_ci(fit_results_m1, "10", "Ace_enr", col="#FFD966", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,0.7), xaxs="i", yaxs="i", ylab="", xlab="time (h)", main="") +plot_no_ci(fit_results_m1, "30", "Ace_enr", col="#FFE699", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,0.3), xaxs="i", yaxs="i", ylab="", xlab="time (h)", main="") + +dev.off() + +# MODEL 2 # +########### + +setwd(model_dir) +loadModel("Millard2020_Ecoli_glc_ace_model_2.cps") + +# get parameters for the best fit +res_PE_model <- runParameterEstimation(method = "Statistics", model = getCurrentModel()) + +# list used to set up global sensitivity analysis and store calculation results +fit_results_m2 <- list("1" = list("initial_concentrations"= c("Ace0_out"=0.92, "Glc"=12.9, "X"=0.0671), + "species_idx" = c("[Ace0_out]_0"=1, "[Glc]_0"=1, "[X]_0"=1), + "data_exp" = read.table("data_1mM.txt", header=TRUE, sep="\t"), + "sd" = c("Ace"=0.2, "Glc"=0.5, "X"=0.045, "Ace_enr"=0.03), + "mapping" = c("Ace"="Acet_out", "Glc"="Glc", "X"="X", "Ace_enr"="Values[Ace_enr]")), + "10" = list("initial_concentrations"= c("Ace0_out"=9, "Glc"=14.192, "X"=0.07), + "species_idx" = c("[Ace0_out]_0"=2, "[Glc]_0"=2, "[X]_0"=2), + "data_exp" = read.table("data_10mM.txt", header=TRUE, sep="\t"), + "sd" = c("Ace"=0.2, "Glc"=0.5, "X"=0.045, "Ace_enr"=0.03), + "mapping" = c("Ace"="Acet_out", "Glc"="Glc", "X"="X", "Ace_enr"="Values[Ace_enr]")), + "30" = list("initial_concentrations"= c("Ace0_out"=30.2448, "Glc"=12.5, "X"=0.07227), + "species_idx" = c("[Ace0_out]_0"=3, "[Glc]_0"=3, "[X]_0"=3), + "data_exp" = read.table("data_30mM.txt", header=TRUE, sep="\t"), + "sd" = c("Ace"=0.2, "Glc"=0.5, "X"=0.045, "Ace_enr"=0.03), + "mapping" = c("Ace"="Acet_out", "Glc"="Glc", "X"="X", "Ace_enr"="Values[Ace_enr]")), + "res_cost" = c(res_PE_model$main$objective_value), + "res_converged" = c(test_khi2(nb_points=152, k_val=res_PE_model$main$objective_value, nb_par=10)$`p-value, i.e. P(X^2<=value)` < 0.95)) + +# simulate data for the best fit +for (i in c("1", "10", "30")){ + # update initial conditions (Glc, Ace and biomass concentrations) + for (j in names(fit_results_m2[[i]][["initial_concentrations"]])){ + setSpecies(key=j, initial_concentration=fit_results_m2[[i]][["initial_concentrations"]][j]) + } + # apply initial state + applyInitialState() + # run time course simulation + tmp <- as.matrix(runTimeCourse()$result) + # get complete simulation results + idx_t_max <- (tmp[,"Time"] <= max(fit_results_m2[[i]][["data_exp"]]$time)) + fit_results_m2[[i]]$simulations <- tmp[idx_t_max,] + # get index of simulated data at experimental time points + fit_results_m2[[i]]$t_idx <- match(fit_results_m2[[i]][["data_exp"]]$time, fit_results_m2[[i]]$simulations[,"Time"]) + # get simulated data that match experimental time points + sim_data <- fit_results_m2[[i]][["data_exp"]] + for (k in names(fit_results_m2[[i]][["mapping"]])){ + sim_data[,k] <- fit_results_m2[[i]]$simulations[fit_results_m2[[i]]$t_idx,fit_results_m2[[i]][["mapping"]][[k]]] + } + fit_results_m2[[i]]$sim <- sim_data +} + +# Plot fitting results +setwd(results_dir) +pdf(file = "Figure 1-figure supplement 2.pdf", width = 7, height = 9) +par(mfrow=c(4,3)) + +plot_no_ci(fit_results_m2, "1", "X", col="#1F4E79", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,1), xaxs="i", yaxs="i", ylab="[biomass] (gDW/L)", xlab="", main="[acetate]=1mM") +plot_no_ci(fit_results_m2, "10", "X", col="#2E75B6", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,1), xaxs="i", yaxs="i", ylab="", xlab="", main="[acetate]=10mM") +plot_no_ci(fit_results_m2, "30", "X", col="#9DC3E6", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,1), xaxs="i", yaxs="i", ylab="", xlab="", main="[acetate]=30mM") + +plot_no_ci(fit_results_m2, "1", "Glc", col="#548235", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,15), xaxs="i", yaxs="i", ylab="[Glc] (mM)", xlab="", main="") +plot_no_ci(fit_results_m2, "10", "Glc", col="#70AD47", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,15), xaxs="i", yaxs="i", ylab="", xlab="", main="") +plot_no_ci(fit_results_m2, "30", "Glc", col="#A9D18E", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,15), xaxs="i", yaxs="i", ylab="", xlab="", main="") + +plot_no_ci(fit_results_m2, "1", "Ace", col="#B63A2D", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,5), xaxs="i", yaxs="i", ylab="[Ace] (mM)", xlab="", main="") +plot_no_ci(fit_results_m2, "10", "Ace", col="#D6685C", lwd=1.0, las=1, xlim=c(0,5), ylim=c(8,11), xaxs="i", yaxs="i", ylab="", xlab="", main="") +plot_no_ci(fit_results_m2, "30", "Ace", col="#E49890", lwd=1.0, las=1, xlim=c(0,5), ylim=c(25,32), xaxs="i", yaxs="i", ylab="", xlab="", main="") + +plot_no_ci(fit_results_m2, "1", "Ace_enr", col="#BF9000", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,1), xaxs="i", yaxs="i", ylab="[Ace] (mM)", xlab="time (h)", main="") +plot_no_ci(fit_results_m2, "10", "Ace_enr", col="#FFD966", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,0.7), xaxs="i", yaxs="i", ylab="", xlab="time (h)", main="") +plot_no_ci(fit_results_m2, "30", "Ace_enr", col="#FFE699", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,0.3), xaxs="i", yaxs="i", ylab="", xlab="time (h)", main="") + +dev.off() + +########### +# MODEL 3 # +########### + +setwd(model_dir) +loadModel("Millard2020_Ecoli_glc_ace_model_3.cps") + +# get parameters for the best fit +res_PE_model <- runParameterEstimation(method = "Statistics", model = getCurrentModel()) + +# list used to set up global sensitivity analysis and store calculation results +fit_results_m3 <- list("1" = list("initial_concentrations"= c("Ace0_out"=0.92, "Glc"=12.9, "X"=0.0671), + "species_idx" = c("[Ace0_out]_0"=1, "[Glc]_0"=1, "[X]_0"=1), + "data_exp" = read.table("data_1mM.txt", header=TRUE, sep="\t"), + "sd" = c("Ace"=0.2, "Glc"=0.5, "X"=0.045, "Ace_enr"=0.03), + "mapping" = c("Ace"="Acet_out", "Glc"="Glc", "X"="X", "Ace_enr"="Values[Ace_enr]")), + "10" = list("initial_concentrations"= c("Ace0_out"=9, "Glc"=14.192, "X"=0.07), + "species_idx" = c("[Ace0_out]_0"=2, "[Glc]_0"=2, "[X]_0"=2), + "data_exp" = read.table("data_10mM.txt", header=TRUE, sep="\t"), + "sd" = c("Ace"=0.2, "Glc"=0.5, "X"=0.045, "Ace_enr"=0.03), + "mapping" = c("Ace"="Acet_out", "Glc"="Glc", "X"="X", "Ace_enr"="Values[Ace_enr]")), + "30" = list("initial_concentrations"= c("Ace0_out"=30.2448, "Glc"=12.5, "X"=0.07227), + "species_idx" = c("[Ace0_out]_0"=3, "[Glc]_0"=3, "[X]_0"=3), + "data_exp" = read.table("data_30mM.txt", header=TRUE, sep="\t"), + "sd" = c("Ace"=0.2, "Glc"=0.5, "X"=0.045, "Ace_enr"=0.03), + "mapping" = c("Ace"="Acet_out", "Glc"="Glc", "X"="X", "Ace_enr"="Values[Ace_enr]")), + "res_cost" = c(res_PE_model$main$objective_value), + "res_converged" = c(test_khi2(nb_points=152, k_val=res_PE_model$main$objective_value, nb_par=10)$`p-value, i.e. P(X^2<=value)` < 0.95)) + +# simulate data for the best fit +for (i in c("1", "10", "30")){ + # update initial conditions (Glc, Ace and biomass concentrations) + for (j in names(fit_results_m3[[i]][["initial_concentrations"]])){ + setSpecies(key=j, initial_concentration=fit_results_m3[[i]][["initial_concentrations"]][j]) + } + # apply initial state + applyInitialState() + # run time course simulation + tmp <- as.matrix(runTimeCourse()$result) + # get complete simulation results + idx_t_max <- (tmp[,"Time"] <= max(fit_results_m3[[i]][["data_exp"]]$time)) + fit_results_m3[[i]]$simulations <- tmp[idx_t_max,] + # get index of simulated data at experimental time points + fit_results_m3[[i]]$t_idx <- match(fit_results_m3[[i]][["data_exp"]]$time, fit_results_m3[[i]]$simulations[,"Time"]) + # get simulated data that match experimental time points + sim_data <- fit_results_m3[[i]][["data_exp"]] + for (k in names(fit_results_m3[[i]][["mapping"]])){ + sim_data[,k] <- fit_results_m3[[i]]$simulations[fit_results_m3[[i]]$t_idx,fit_results_m3[[i]][["mapping"]][[k]]] + } + fit_results_m3[[i]]$sim <- sim_data +} + +# Plot fitting results +setwd(results_dir) +pdf(file = "Figure 1-figure supplement 3.pdf", width = 7, height = 9) +par(mfrow=c(4,3)) + +plot_no_ci(fit_results_m3, "1", "X", col="#1F4E79", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,1), xaxs="i", yaxs="i", ylab="[biomass] (gDW/L)", xlab="", main="[acetate]=1mM") +plot_no_ci(fit_results_m3, "10", "X", col="#2E75B6", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,1), xaxs="i", yaxs="i", ylab="", xlab="", main="[acetate]=10mM") +plot_no_ci(fit_results_m3, "30", "X", col="#9DC3E6", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,1), xaxs="i", yaxs="i", ylab="", xlab="", main="[acetate]=30mM") + +plot_no_ci(fit_results_m3, "1", "Glc", col="#548235", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,15), xaxs="i", yaxs="i", ylab="[Glc] (mM)", xlab="", main="") +plot_no_ci(fit_results_m3, "10", "Glc", col="#70AD47", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,15), xaxs="i", yaxs="i", ylab="", xlab="", main="") +plot_no_ci(fit_results_m3, "30", "Glc", col="#A9D18E", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,15), xaxs="i", yaxs="i", ylab="", xlab="", main="") + +plot_no_ci(fit_results_m3, "1", "Ace", col="#B63A2D", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,5), xaxs="i", yaxs="i", ylab="[Ace] (mM)", xlab="", main="") +plot_no_ci(fit_results_m3, "10", "Ace", col="#D6685C", lwd=1.0, las=1, xlim=c(0,5), ylim=c(8,11), xaxs="i", yaxs="i", ylab="", xlab="", main="") +plot_no_ci(fit_results_m3, "30", "Ace", col="#E49890", lwd=1.0, las=1, xlim=c(0,5), ylim=c(25,32), xaxs="i", yaxs="i", ylab="", xlab="", main="") + +plot_no_ci(fit_results_m3, "1", "Ace_enr", col="#BF9000", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,1), xaxs="i", yaxs="i", ylab="[Ace] (mM)", xlab="time (h)", main="") +plot_no_ci(fit_results_m3, "10", "Ace_enr", col="#FFD966", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,0.7), xaxs="i", yaxs="i", ylab="", xlab="time (h)", main="") +plot_no_ci(fit_results_m3, "30", "Ace_enr", col="#FFE699", lwd=1.0, las=1, xlim=c(0,5), ylim=c(0,0.3), xaxs="i", yaxs="i", ylab="", xlab="time (h)", main="") + +dev.off() + +############################################## +# MODEL 4 - with global sensitivity analysis # +############################################## + +setwd(model_dir) +loadModel("Millard2020_Ecoli_glc_ace_model_4_MC.cps") + +# create calibration data files for monte carlo simulations +file.copy("data_1mM.txt", "data_1mM_MC.txt", overwrite=TRUE) +file.copy("data_10mM.txt", "data_10mM_MC.txt", overwrite=TRUE) +file.copy("data_30mM.txt", "data_30mM_MC.txt", overwrite=TRUE) + +# number of Monte Carlo iterations +n_iter_mc <- 10 + +# get parameters for the best fit +res_PE_model4 <- runParameterEstimation(method = "Statistics", model = getCurrentModel()) + +# list used to set up global sensitivity analysis and store calculation results +fit_results <- list("1" = list("initial_concentrations"= c("Ace0_out"=0.92, "Glc"=12.9, "X"=0.0671), + "species_idx" = c("[Ace0_out]_0"=1, "[Glc]_0"=1, "[X]_0"=1), + "data_exp" = read.table("data_1mM.txt", header=TRUE, sep="\t"), + "sd" = c("Ace"=0.2, "Glc"=0.5, "X"=0.045, "Ace_enr"=0.03), + "mapping" = c("Ace"="Acet_out", "Glc"="Glc", "X"="X", "Ace_enr"="Values[Ace_enr]")), + "10" = list("initial_concentrations"= c("Ace0_out"=9, "Glc"=14.192, "X"=0.07), + "species_idx" = c("[Ace0_out]_0"=2, "[Glc]_0"=2, "[X]_0"=2), + "data_exp" = read.table("data_10mM.txt", header=TRUE, sep="\t"), + "sd" = c("Ace"=0.2, "Glc"=0.5, "X"=0.045, "Ace_enr"=0.03), + "mapping" = c("Ace"="Acet_out", "Glc"="Glc", "X"="X", "Ace_enr"="Values[Ace_enr]")), + "30" = list("initial_concentrations"= c("Ace0_out"=30.2448, "Glc"=12.5, "X"=0.07227), + "species_idx" = c("[Ace0_out]_0"=3, "[Glc]_0"=3, "[X]_0"=3), + "data_exp" = read.table("data_30mM.txt", header=TRUE, sep="\t"), + "sd" = c("Ace"=0.2, "Glc"=0.5, "X"=0.045, "Ace_enr"=0.03), + "mapping" = c("Ace"="Acet_out", "Glc"="Glc", "X"="X", "Ace_enr"="Values[Ace_enr]")), + "res_par" = res_PE_model4$parameters[, c("parameter", "value")], + "res_cost" = c(res_PE_model4$main$objective_value), + "res_converged" = c(test_khi2(nb_points=152, k_val=res_PE_model4$main$objective_value, nb_par=10)$`p-value, i.e. P(X^2<=value)` < 0.95)) + +# simulate data for the best fit +for (i in c("1", "10", "30")){ + # update initial conditions (Glc, Ace and biomass concentrations) + for (j in names(fit_results[[i]][["initial_concentrations"]])){ + setSpecies(key=j, initial_concentration=fit_results[[i]][["initial_concentrations"]][j]) + } + # apply initial state + applyInitialState() + # run time course simulation + tmp <- as.matrix(runTimeCourse()$result) + # get complete simulation results + idx_t_max <- (tmp[,"Time"] <= max(fit_results[[i]][["data_exp"]]$time)) + fit_results[[i]]$simulations <- array(NA, dim=c(n_iter_mc+1, sum(idx_t_max), ncol(tmp)), dimnames=list(iter=NULL, time=NULL, specie=colnames(tmp))) + fit_results[[i]]$simulations[1,,] <- tmp[idx_t_max,] + # get index of simulated data at experimental time points + fit_results[[i]]$t_idx <- match(fit_results[[i]][["data_exp"]]$time, fit_results[[i]]$simulations[1,,"Time"]) + # get simulated data that match experimental time points + sim_data <- fit_results[[i]][["data_exp"]] + for (k in names(fit_results[[i]][["mapping"]])){ + sim_data[,k] <- fit_results[[i]]$simulations[1,fit_results[[i]]$t_idx,fit_results[[i]][["mapping"]][[k]]] + } + fit_results[[i]]$sim_mc <- sim_data +} + +# create progress bar +pb <- txtProgressBar(min=0, max=n_iter_mc, style=3) + +# run global sensitivity analysis +j <- 0 +while (j < n_iter_mc){ + + # generate noisy datasets according to the experimental standard deviation + for (i in c("1", "10", "30")){ + # initialize dataset with simulations of the best fit + data_noise <- fit_results[[i]]$sim_mc + nr <- nrow(data_noise) + # add noise + for (k in names(fit_results[[i]]$sd)){ + data_noise[,k] <- rnorm(nr, mean=data_noise[,k], sd=fit_results[[i]]$sd[[k]]) + } + # set negative concentrations to 0 + data_noise[data_noise < 0] <- 0 + # save the noisy dataset + write.table(data_noise, paste("data_", i, "mM_MC.txt", sep=""), sep="\t", row.names = FALSE, quote = FALSE) + } + + # run parameter estimation (SRES followed by Hooke and Jeeves and Levenberg-Marquadt methods to refine convergence) + # note: this protocol ensures faster convergence in comparison to PSO, which is necessary for monte carlo analysis + res_PE_tmp <- runParameterEstimation(method = "SRES", randomize_start_values=TRUE, update_model=TRUE) + res_PE_tmp <- runParameterEstimation(method = "HookeJeeves", randomize_start_values=FALSE, update_model=TRUE) + res_PE_tmp <- runParameterEstimation(method = "LevenbergMarquardt", randomize_start_values=FALSE, update_model=TRUE) + + # check if optimization has converged, reiterate if not + p_val_chi2 <- test_khi2(nb_points=152, k_val=res_PE_tmp$main$objective_value, nb_par=10)$`p-value, i.e. P(X^2<=value)` + if (p_val_chi2 < 0.95){ + j = j+1 + }else{ + next + } + + # update the progress bar + setTxtProgressBar(pb, j) + + # save parameter estimation results + fit_results$res_par <- cbind(fit_results$res_par, res_PE_tmp$parameters[, "value"]) + fit_results$res_cost <- c(fit_results$res_cost, res_PE_tmp$main$objective_value) + fit_results$res_converged <- c(fit_results$res_converged, p_val_chi2 < 0.95) + + # simulate all experiments + for (i in c("1", "10", "30")){ + # update experiment-dependent parameters + for (l in names(fit_results[[i]]$species_idx)){ + idx <- which(res_PE_tmp$parameters$parameter==l)[fit_results[[i]]$species_idx[l]] + k <- sub(".*\\[(.*)\\].*", "\\1", l, perl=TRUE) + setSpecies(key=k, initial_concentration=as.numeric(res_PE_tmp$parameters[idx,"value"])) + } + # apply initial state + applyInitialState() + # run time course simulations + tmp <- as.matrix(runTimeCourse()$result) + # get complete simulation results + fit_results[[i]]$simulations[j+1,,] <- tmp[tmp[,"Time"] <= max(fit_results[[i]][["data_exp"]]$time),] + } +} + +# close progress bar +close(pb) + +# save global sensitivity analysis results for further use +setwd(results_dir) +save(fit_results, file = "mc_results_10.RData") + +# load global sensitivity analysis results +# (here we use results obtained with 100 iterations) +load("mc_results_100.RData") + +# Calculate statistics on parameters (mean, sd, rsd) +res_stats_params <- get_parameters_stats(fit_results) +print(res_stats_params) +write.table(res_stats_params, "parameters_stats.txt", sep="\t", quote = FALSE) + +# Plot fitting results +pdf(file = "Figure 1C.pdf", width = 7, height = 9) +par(mfrow=c(4,3)) + +plot_with_ci(fit_results, "1", "X", col="#1F4E79", lwd=1.2, las=1, xlim=c(0,5), ylim=c(0,1), xaxs="i", yaxs="i", ylab="[biomass] (gDW/L)", xlab="", main="[acetate]=1mM") +plot_with_ci(fit_results, "10", "X", col="#2E75B6", lwd=1.2, las=1, xlim=c(0,5), ylim=c(0,1), xaxs="i", yaxs="i", ylab="", xlab="", main="[acetate]=10mM") +plot_with_ci(fit_results, "30", "X", col="#9DC3E6", lwd=1.2, las=1, xlim=c(0,5), ylim=c(0,1), xaxs="i", yaxs="i", ylab="", xlab="", main="[acetate]=30mM") + +plot_with_ci(fit_results, "1", "Glc", col="#548235", lwd=1.2, las=1, xlim=c(0,5), ylim=c(0,15), xaxs="i", yaxs="i", ylab="[Glc] (mM)", xlab="", main="") +plot_with_ci(fit_results, "10", "Glc", col="#70AD47", lwd=1.2, las=1, xlim=c(0,5), ylim=c(0,15), xaxs="i", yaxs="i", ylab="", xlab="", main="") +plot_with_ci(fit_results, "30", "Glc", col="#A9D18E", lwd=1.2, las=1, xlim=c(0,5), ylim=c(0,15), xaxs="i", yaxs="i", ylab="", xlab="", main="") + +plot_with_ci(fit_results, "1", "Ace", col="#B63A2D", lwd=1.2, las=1, xlim=c(0,5), ylim=c(0,5), xaxs="i", yaxs="i", ylab="[Ace] (mM)", xlab="", main="") +plot_with_ci(fit_results, "10", "Ace", col="#D6685C", lwd=1.2, las=1, xlim=c(0,5), ylim=c(8,11), xaxs="i", yaxs="i", ylab="", xlab="", main="") +plot_with_ci(fit_results, "30", "Ace", col="#E49890", lwd=1.2, las=1, xlim=c(0,5), ylim=c(27,32), xaxs="i", yaxs="i", ylab="", xlab="", main="") + +plot_with_ci(fit_results, "1", "Ace_enr", col="#BF9000", lwd=1.2, las=1, xlim=c(0,5), ylim=c(0,1), xaxs="i", yaxs="i", ylab="[Ace] (mM)", xlab="time (h)", main="") +plot_with_ci(fit_results, "10", "Ace_enr", col="#FFD966", lwd=1.2, las=1, xlim=c(0,5), ylim=c(0,0.7), xaxs="i", yaxs="i", ylab="", xlab="time (h)", main="") +plot_with_ci(fit_results, "30", "Ace_enr", col="#FFE699", lwd=1.2, las=1, xlim=c(0,5), ylim=c(0,0.3), xaxs="i", yaxs="i", ylab="", xlab="time (h)", main="") + +dev.off() + + diff --git a/2-Model_validation.R b/2-Model_validation.R new file mode 100644 index 0000000..ef6b6f4 --- /dev/null +++ b/2-Model_validation.R @@ -0,0 +1,1049 @@ + + +################### +# Set environment # +################### + +# load libraries and initialize environment +source("set_env.R") + + +############################################ +# Load global sensitivity analysis results # +############################################ + +# This file is generated by script "1-Model_construction.R". +setwd(results_dir) +load("mc_results_100.RData") + + +################### +# Testing model 4 # +################### + +setwd(results_dir) +pdf(file="Figure 3.pdf", width = 7, height = 9) +par(mfrow=c(4,3)) + +# Pulse experiment (Enjalbert et al., 2017) +########################################### + +setwd(model_dir) +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model.cps") +setwd(results_dir) + +# set biomass concentration at which the pulse is performed +setGlobalQuantities(key = "_X_conc_pulse", initial_value = 0.9) +# simulate response to acetate pulse +res_pulse <- runTimeCourse() +# get simulation results +id_start <- which(res_pulse$result$Ace_out >= 30)[1]+1 +id_end <- which(res_pulse$result$Time >= (res_pulse$result$Time[id_start]+8/60))[1]+1 +t_pulse <- (res_pulse$result$Time[seq(id_start, id_end)] - res_pulse$result$Time[id_start]) +ace_conc_pulse <- res_pulse$result$Ace_out[seq(id_start, id_end)] - res_pulse$result$Ace_out[id_start] +glc_conc_pulse <- res_pulse$result$Glc[seq(id_start, id_end)] - res_pulse$result$Glc[id_start] + +# set biomass concentration at which the pulse is performed +setGlobalQuantities(key = "_X_conc_pulse", initial_value = 1000) +res_nopulse <- runTimeCourse() +# get simulation results +ace_conc_nopulse <- res_nopulse$result$Ace_out[seq(id_start, id_end)] - res_nopulse$result$Ace_out[id_start] +glc_conc_nopulse <- res_nopulse$result$Glc[seq(id_start, id_end)] - res_nopulse$result$Glc[id_start] + +run_sa <- function(){ + # set biomass concentration at which the pulse is performed + setGlobalQuantities(key = "_X_conc_pulse", initial_value = 0.9) + + # simulate response to acetate pulse + res_pulse <- runTimeCourse() + + # get simulation results + id_start <- which(res_pulse$result$Ace_out >= 30)[1]+1 + id_end <- which(res_pulse$result$Time >= (res_pulse$result$Time[id_start]+8/60))[1]+1 + + t_pulse <- (res_pulse$result$Time[seq(id_start, id_end)] - res_pulse$result$Time[id_start]) + ace_conc_pulse <- res_pulse$result$Ace_out[seq(id_start, id_end)] - res_pulse$result$Ace_out[id_start] + glc_conc_pulse <- res_pulse$result$Glc[seq(id_start, id_end)] - res_pulse$result$Glc[id_start] + + # set biomass concentration at which the pulse is performed + setGlobalQuantities(key = "_X_conc_pulse", initial_value = 1000) + + # simulate response to water pulse (control experiment) + res_nopulse <- runTimeCourse() + + # get simulation results + ace_conc_nopulse <- res_nopulse$result$Ace_out[seq(id_start, id_end)] - res_nopulse$result$Ace_out[id_start] + glc_conc_nopulse <- res_nopulse$result$Glc[seq(id_start, id_end)] - res_nopulse$result$Glc[id_start] + + return(list(ace_conc_pulse=ace_conc_pulse, glc_conc_pulse=glc_conc_pulse, ace_conc_nopulse=ace_conc_nopulse, glc_conc_nopulse=glc_conc_nopulse, t_pulse=t_pulse)) +} + +# create progress bar +pb <- txtProgressBar(min=0, max=ncol(fit_results$res_par)-1, style=3) + +# sensitivity analysis +# run simulations for all sets of parameters +res <- array(NA, dim=c(ncol(fit_results$res_par)-1, 4, length(t_pulse)), dimnames=list(iter=NULL, var=c("ace_conc_pulse", "glc_conc_pulse", "ace_conc_nopulse", "glc_conc_nopulse"), time=NULL)) +for (i in seq(ncol(fit_results$res_par)-1)){ + rp <- c(fit_results$res_par[,i+1]) + names(rp) <- fit_results$res_par[,"parameter"] + model <- update_params(getCurrentModel(), rp) + rsa <- run_sa() + for (j in names(rsa)){ + if (j != "t_pulse"){ + res[i,j,] <- rsa[[j]] + } + } + # update the progress bar + setTxtProgressBar(pb, i) +} + +# close progress bar +close(pb) + +# measurements +# sampling time +time_meas <- seq(0,8)/60 +# acetate pulse experiment +glc <- c(0, -0.06931183, -0.151415145, -0.189227994, -0.269451057, -0.290764495, -0.230785281, -0.464084162, -0.551878527) +sd_glc <- c(0, 0.032476344, 0.073133915, 0.113018846, 0.049485284, 0.005325541, 0.163377704, 0.034786419, 0.048477157) +ace <- c(0, -0.027907926, -0.078000853, -0.155334163, -0.165031608, -0.111098424, -0.182877548, -0.237262298, -0.276903255) +sd_ace <- c(0, 0.002740145, 0.025693594, 0.053641876, 0.089975321, 0.005875669, 0.085604161, 0.061930626, 0.099140975) +# control experiment (pulse without acetate) +glc_nopulse <- c(0, -0.195637774, -0.325512845, -0.433785334, -0.628518958, -0.726913548, -0.892908748, -1.071230579, -1.16285575) +sd_glc_nopulse <- c(0, 0.058023617, 0.109115511, 0.047090371, 0.052331837, 0.065619906, 0.133896355, 0.16828754, 0.03515465) +ace_nopulse <- c(0, 0.01012067, 0.05974009, 0.086787283, 0.086690257, 0.104329693, 0.128507087, 0.130026354, 0.187336965) +sd_ace_nopulse <- c(0, 0.005117257, 0.022507218, 0.040319045, 0.037548873, 0.037235285, 0.044494365, 0.045029023, 0.023982374) + +# plot simulations vs measurements +plot_with_ci_2(t_pulse, ace_conc_pulse, res[,"ace_conc_pulse",], time_meas, ace, sd_ace, col="#D6685C", h=0, ylim=c(-0.4,0.3), las=1, main="Enjalbert_2017 (pulse 30mM ace)", ylab="change of [Ace]", xlab="time (h)", lwd=1.2) +plot_with_ci_2(t_pulse, glc_conc_pulse, res[,"glc_conc_pulse",], time_meas, glc, sd_glc, col="#70AD47", h=NULL, ylim=c(-1.4,0), las=1, main="Enjalbert_2017 (pulse 30mM ace)", ylab="change of [Glc]", xlab="time (h)", lwd=1.2) +plot_with_ci_2(t_pulse, ace_conc_nopulse, res[,"ace_conc_nopulse",], time_meas, ace_nopulse, sd_ace_nopulse, col="#D6685C", h=0, ylim=c(-0.4,0.3), las=1, main="Enjalbert_2017 (control)", ylab="change of [Ace]", xlab="time (h)", lwd=1.2) +plot_with_ci_2(t_pulse, glc_conc_nopulse, res[,"glc_conc_nopulse",], time_meas, glc_nopulse, sd_glc_nopulse, col="#70AD47", h=NULL, ylim=c(-1.4,0), las=1, main="Enjalbert_2017 (control)", ylab="change of [Glc]", xlab="time (h)", lwd=1.2) + +# Chemostat experiment (Renilla et al., 2012) +############################################# + +setwd(model_dir) +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model.cps") +setwd(results_dir) + +# delete events and set glc feed to 1 +deleteEvent(getEvents()$key) +setGlobalQuantities(key = "_feed", initial_value = 1) + +n_step <- 50 +dilution_rates <- seq(0.1, 0.5, length.out = n_step) +fluxes <- c("Values[v_growth_rate]", "Values[v_glc_uptake]", "Values[v_ace_net]") + +# sensitivity analysis + +# create progress bar +pb <- txtProgressBar(min=0, max=ncol(fit_results$res_par)-1, style=3) + +# run simulations for all sets of parameters +res <- array(NA, dim=c(ncol(fit_results$res_par)-1, n_step, length(fluxes)+1), dimnames=list(iter=NULL, dilution_rate=NULL, specie=c("dilution_rate", "Values[v_growth_rate]", "Values[v_glc_uptake]", "Values[v_ace_net]"))) +for (i in seq(ncol(fit_results$res_par)-1)){ + rp <- c(fit_results$res_par[,i+1]) + names(rp) <- fit_results$res_par[,"parameter"] + model <- update_params(getCurrentModel(), rp) + res_chemostat <- matrix(NA, nrow=n_step, ncol=length(fluxes)+1, dimnames=list(r=NULL, c=c("dilution_rate", fluxes))) + for (j in seq(n_step)){ + setGlobalQuantities(key = "_dilution_rate", initial_value = dilution_rates[j], model=model) + res_ss <- runSteadyState(model=model) + res_chemostat[j,] <- c(dilution_rates[j], unlist(res_ss$global_quantities[res_ss$global_quantities$key %in% fluxes, "value"])) + } + res[i,,] <- res_chemostat + # update the progress bar + setTxtProgressBar(pb, i) +} + +# close progress bar +close(pb) + +# acetate flux as function of dilution rate +dilution_rate <- c(0.09586056644880175, 0.20043572984749464, 0.2997821350762528, 0.3468409586056645, 0.39912854030501094, 0.44618736383442276, 0.5002178649237474) +q_ace <- c(0.40472342596168076, 0.7396156614846294, 1.477019302736899, 1.2154669794005626, 1.96934635755591, 2.0929143843289824, 2.006569318707304) +q_ace_sd <- c(0.120440687823878, 0.15316333195082, 0.04595662724122, 0.0274214232252601, 0.28885933172199, 0.0777639526513201, 0.28664731148963) + +plot_with_ci_3(res, "dilution_rate", "Values[v_ace_net]", "#D6685C", las=1, xlim=c(0,0.6), ylim=c(0,4), main="Renilla_2012", lwd=1.2, xlab="dilution rate (h-1)", ylab="ace flux (mmol/gDW/h)") +plot_points(dilution_rate, q_ace, q_ace_sd, offset=0.01, col="#D6685C") + +# acetate flux as function of glc uptake +q_glc <- c(1.43478260869565, 2.7391304347826, 4.59130434782608, 4.69565217391304, 5.7391304347826, 5.92173913043478, 6.10434782608695, 7.12173913043477, 8.34782608695652) +q_ace <- c(0.381358340437624, 0.762906128635029, 1.40509614473808, 1.17836506583309, 1.99602159704461, 2.08290233967983, 2.01292033721701, 2.44717249218528, 2.5848252344416) +q_ace_sd <- c(0.139547219854126, 0.139471440750213, 0.15697641375391, 0.1045751633987, 0.19179691200152, 0.0174291938997801, 0.24397082504499, 0.52295159609737, 0.33130624230368) +q_glc_sd <- c(0.20234760540444, 0.58899555210581, 0.36827270271402, 0.14607356337583, 0.630432119199, 0.13069807968161, 0.47630364392495, 0.99361847036576, 0.66185122193432) + +plot_with_ci_3(res, "Values[v_glc_uptake]", "Values[v_ace_net]", col="#D6685C", las=1, xlim=c(0,10), ylim=c(0,4), main="Renilla_2012", lwd=1.2, xlab="glc uptake (mmol/gDW/h)", ylab="ace flux (mmol/gDW/h)") +plot_points(q_glc, q_ace, q_ace_sd, offset=0.2, col="#D6685C") +plot_points(q_glc, q_ace, q_glc_sd, offset=0.08, mode="h", col="#D6685C") + +# Steady-state fluxes for acetate concentration between 0.1 and 100mM (Enjalbert et al., 2017; Pinhal et al., 2019) +################################################################################################################### + +setwd(model_dir) +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model.cps") +setwd(results_dir) + +# remove events and fix concentrations of actate, glucose and biomass +deleteEvent(getEvents()$key) +setSpecies(key="Ace_out", type="fixed") +setSpecies(key="Glc", type="fixed") +setSpecies(key="X", type="fixed") + +# run simulations +n_step <- 50 +ace_concentration <- 10**seq(-1, 2, length.out = n_step) +fluxes <- c("Values[v_growth_rate]", "Values[v_glc_uptake]", "Values[v_ace_net]") + +# sensitivity analysis +# run simulations for all sets of parameters +res_e2 <- array(NA, dim=c(ncol(fit_results$res_par)-1, n_step, length(fluxes)+1), dimnames=list(iter=NULL, ace_concentration=NULL, specie=c("ace_concentration", fluxes))) + +# create progress bar +pb <- txtProgressBar(min=0, max=ncol(fit_results$res_par)-1, style=3) + +for (i in seq(ncol(fit_results$res_par)-1)){ + rp <- c(fit_results$res_par[,i+1]) + names(rp) <- fit_results$res_par[,"parameter"] + model <- update_params(getCurrentModel(), rp) + + res_ace_range <- matrix(NA, nrow=n_step, ncol=length(fluxes)+1, dimnames=list(r=NULL, c=c("ace_concentration", fluxes))) + for (j in seq(n_step)){ + setSpecies(key="Ace_out", initial_concentration=ace_concentration[j]) + applyInitialState() + res_ss <- runSteadyState() + res_ace_range[j,] <- c(ace_concentration[j], unlist(res_ss$global_quantities[res_ss$global_quantities$key %in% fluxes, "value"])) + } + res_e2[i,,] <- res_ace_range + # update the progress bar + setTxtProgressBar(pb, i) +} + +# close progress bar +close(pb) + +xlab_main <- c(0.1, 1, 10, 100) +xlab_sec <- c(seq(0.2, 0.9, by=0.1), seq(2, 9, by=1), seq(20, 90, by=10)) + +# growth rate as function of acetate concentration +growth_rates <- c(0.521128511, 0.611148842625582, 0.613161998174498, 0.502533817, 0.496290415, 0.488201506, 0.547635665, 0.499830448, 0.474554197, 0.425356578, 0.377534684, 0.645724326, 0.618475601, 0.554887936, 0.564811523, 0.527571192, 0.434972836, 0.3824734, 0.583623355, 0.620905534, 0.564259247, 0.532148135, 0.483885215, 0.557074418, 0.630654409249223) +sd_growth_rates <- c(0.001793104, 0.00204807928657914, 0.00219396182484705, 0.001709207, 0.001846205, 0.001757403, 0.001821375, 0.001025702, 0.001940912, 0.001204707, 0.001999188, 0.001418374, 0.001932601, 0.001455791, 0.001574234, 0.001206265, 0.001292476, 0.001068259, 0.001804648, 0.001643459, 0.001598405, 0.001121218, 0.000912408, 0.00194896, 0.00203369597574686) +ace_conc <- c(0.451964624, 1.11600286648471, 2.04718732931708, 2.86252145, 5.285907977, 9.106164204, 16.67476528, 16.76626787, 17.00218707, 30.8667961, 57.92292091, 0.478352574, 4.55673229, 8.10163028, 8.22100734, 7.829591756, 33.53244905, 66.20361403, 0.436956014, 1.28468189, 1.555875222, 12.02564968, 30.24566673, 0.548011282, 2.29562227069566) +sd_ace_conc <- c(0.262318006, 0.281361538208953, 0.289527601555302, 0.302351163, 0.327705513, 0.330782201, 0.277259011, 0.233956798, 0.31883929, 0.300578057, 0.313371784, 0.231202155, 0.264687437, 0.243480317, 0.289821733, 0.263808862, 0.289478134, 0.264861034, 0.22461248, 0.229031308, 0.241718918, 0.254111384, 0.187394292, 0.011107606, 0.290519090995892) + +plot_with_ci_3(res_e2, "ace_concentration", "Values[v_growth_rate]", col="#2E75B6", xaxt='n', las=1, lwd=1.2, xlim=c(0.1,100), ylim=c(0,0.8), log="x", main="Enjalbert_2017", xlab="[acetate] (mM)", ylab="growth rate (h-1)") +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +plot_points(ace_conc, growth_rates, sd_growth_rates, offset=0.01, col="#2E75B6") +plot_points(ace_conc, growth_rates, sd_ace_conc, offset=0.02, mode="h", col="#2E75B6") + +# glc uptake as function of acetate concentration +glc_upt <- c(8.654860167, 8.36127542981722, 7.98010111285252, 9.236935826, 8.274418986, 7.560431219, 7.339194455, 5.775312502, 6.423391263, 5.1544758, 3.938631573, 8.115447647, 9.28067031, 6.737153424, 7.172748804, 5.884186033, 5.684201497, 4.811576974, 9.632702365, 8.055042777, 9.708342814, 7.100081588, 5.505759496, 9.242859752, 8.18621623190759) +sd_glc_upt <- c(0.337812425, 0.38531328268303, 0.373770045721031, 0.356787032, 0.334672954, 0.317509322, 0.288025925, 0.16053276, 0.375934255, 0.293148172, 0.359225607, 0.197331684, 0.360984112, 0.229372278, 0.241820396, 0.20450532, 0.260869273, 0.216134352, 0.34289286, 0.350305744, 0.293144783, 0.220135755, 0.153471508, 0.25245346, 0.396815184905029) +ace_conc <- c(0.451964624, 1.11600286648471, 2.04718732931708, 2.86252145, 5.285907977, 9.106164204, 16.67476528, 16.76626787, 17.00218707, 30.8667961, 57.92292091, 0.478352574, 4.55673229, 8.10163028, 8.22100734, 7.829591756, 33.53244905, 66.20361403, 0.436956014, 1.28468189, 1.555875222, 12.02564968, 30.24566673, 0.548011282, 2.29562227069566) +sd_ace_conc <- c(0.262318006, 0.281361538208953, 0.289527601555302, 0.302351163, 0.327705513, 0.330782201, 0.277259011, 0.233956798, 0.31883929, 0.300578057, 0.313371784, 0.231202155, 0.264687437, 0.243480317, 0.289821733, 0.263808862, 0.289478134, 0.264861034, 0.22461248, 0.229031308, 0.241718918, 0.254111384, 0.187394292, 0.011107606, 0.290519090995892) + +plot_with_ci_3(res_e2, "ace_concentration", "Values[v_glc_uptake]", col="#70AD47", xaxt='n', las=1, lwd=1.2, xlim=c(0.1,100), ylim=c(0,10), log="x", main="Enjalbert_2017", xlab="[acetate] (mM)", ylab="glc uptake (mmol/gDW/h)") +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +plot_points(ace_conc, glc_upt, sd_glc_upt, offset=0.02, col="#70AD47") +plot_points(ace_conc, glc_upt, sd_ace_conc, offset=0.2, mode="h", col="#70AD47") + +# ace flux as function of acetate concentration +ace_flx <- c(3.5, -2.7, 1.516999356, 1.26845123082679, 0.775821380507016, 0.678877137, 0.017366464, -0.991478151, -1.286687213, -2.078474994, -1.530841439, -1.525342269, -1.253581266, 1.984679487, 0.546462624, -0.136780389, -0.393883917, -0.610240984, -1.120767885, -1.277455315, 2.574285211, 2.051935093, 1.828415596, -1.262442483, -1.317987733, 2.333568565, 1.85234639824858) +sd_ace_flx <- c(0.35, 0.27, 0.316118066, 0.388752117258161, 0.40715278851436, 0.33357638, 0.37333751, 0.347029894, 0.280501612, 0.195031303, 0.376252463, 0.226182385, 0.303661317, 0.253610517, 0.385450715, 0.243880325, 0.30665695, 0.257983739, 0.23844407, 0.198458448, 0.299832036, 0.334956504, 0.317134334, 0.263807154, 0.195219648, 0.016120887, 0.386174129654754) +ace_conc <- c(0.2, 100, 0.451964624, 1.11600286648471, 2.04718732931708, 2.86252145, 5.285907977, 9.106164204, 16.67476528, 16.76626787, 17.00218707, 30.8667961, 57.92292091, 0.478352574, 4.55673229, 8.10163028, 8.22100734, 7.829591756, 33.53244905, 66.20361403, 0.436956014, 1.28468189, 1.555875222, 12.02564968, 30.24566673, 0.548011282, 2.29562227069566) +sd_ace_conc <- c(0.05, 10, 0.262318006, 0.281361538208953, 0.289527601555302, 0.302351163, 0.327705513, 0.330782201, 0.277259011, 0.233956798, 0.31883929, 0.300578057, 0.313371784, 0.231202155, 0.264687437, 0.243480317, 0.289821733, 0.263808862, 0.289478134, 0.264861034, 0.22461248, 0.229031308, 0.241718918, 0.254111384, 0.187394292, 0.011107606, 0.290519090995892) +plot_with_ci_3(res_e2, "ace_concentration", "Values[v_ace_net]", col="#D6685C", xaxt='n', las=1, lwd=1.2, xlim=c(0.1,100), ylim=c(-4,4), log="x", main="Enjalbert_2017", xlab="[acetate] (mM)", ylab="ace flux (mmol/gDW/h)") +abline(h=0) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +plot_points(ace_conc, ace_flx, sd_ace_flx, offset=0.04, col="#D6685C") +plot_points(ace_conc, ace_flx, sd_ace_conc, offset=0.2, mode="h", col="#D6685C") + +dev.off() + + +############################## +# Testing alternative models # +############################## + +########### +# Model 1 # +########### + +setwd(results_dir) +pdf(file="Figure 3-figure supplement 1.pdf", width = 7, height = 9) +par(mfrow=c(4,3)) + +# Pulse experiment (Enjalbert et al., 2017) +########################################### + +setwd(model_dir) +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_1.cps") + +# set biomass concentration at which the pulse is performed +setGlobalQuantities(key = "_X_conc_pulse", initial_value = 0.9) + +# simulate response to acetate pulse +res_pulse <- runTimeCourse() + +# get simulation results +id_start <- which(res_pulse$result$Ace_out >= 30)[1]+1 +id_end <- which(res_pulse$result$Time >= (res_pulse$result$Time[id_start]+8/60))[1]+1 + +t_pulse <- (res_pulse$result$Time[seq(id_start, id_end)] - res_pulse$result$Time[id_start]) +ace_conc_pulse <- res_pulse$result$Ace_out[seq(id_start, id_end)] - res_pulse$result$Ace_out[id_start] +glc_conc_pulse <- res_pulse$result$Glc[seq(id_start, id_end)] - res_pulse$result$Glc[id_start] + +# measurements +time_meas <- seq(0,8)/60 +glc <- c(0, -0.06931183, -0.151415145, -0.189227994, -0.269451057, -0.290764495, -0.230785281, -0.464084162, -0.551878527) +sd_glc <- c(0, 0.032476344, 0.073133915, 0.113018846, 0.049485284, 0.005325541, 0.163377704, 0.034786419, 0.048477157) +ace <- c(0, -0.027907926, -0.078000853, -0.155334163, -0.165031608, -0.111098424, -0.182877548, -0.237262298, -0.276903255) +sd_ace <- c(0, 0.002740145, 0.025693594, 0.053641876, 0.089975321, 0.005875669, 0.085604161, 0.061930626, 0.099140975) + +# plot simulations vs measurements +plot(t_pulse, ace_conc_pulse, type="l", ylim=c(-0.4,0), las=1, main="Enjalbert_2017 (pulse 30mM ace)", ylab="change of [Ace]", xlab="time (h)", col="#D6685C", lwd=2) +plot_points(time_meas, ace, sd_ace, offset=0.002, col="#D6685C") +plot(t_pulse, glc_conc_pulse, type="l", main="Enjalbert_2017 (pulse 30mM ace)", las=1, ylab="change of [Glc]", xlab="time (h)", col="#70AD47", lwd=2) +plot_points(time_meas, glc, sd_glc, offset=0.002, col="#70AD47") + +# simulate control experiment (i.e. no acetate pulse) +deleteEvent(getEvents()$key) + +res_nopulse <- runTimeCourse() + +# get simulation results +ace_conc_nopulse <- res_nopulse$result$Ace_out[seq(id_start, id_end)] - res_nopulse$result$Ace_out[id_start] +glc_conc_nopulse <- res_nopulse$result$Glc[seq(id_start, id_end)] - res_nopulse$result$Glc[id_start] + +# measurements +glc_nopulse <- c(0, -0.195637774, -0.325512845, -0.433785334, -0.628518958, -0.726913548, -0.892908748, -1.071230579, -1.16285575) +sd_glc_nopulse <- c(0, 0.058023617, 0.109115511, 0.047090371, 0.052331837, 0.065619906, 0.133896355, 0.16828754, 0.03515465) +ace_nopulse <- c(0, 0.01012067, 0.05974009, 0.086787283, 0.086690257, 0.104329693, 0.128507087, 0.130026354, 0.187336965) +sd_ace_nopulse <- c(0, 0.005117257, 0.022507218, 0.040319045, 0.037548873, 0.037235285, 0.044494365, 0.045029023, 0.023982374) + +# plot simulations vs measurements +plot(t_pulse, ace_conc_nopulse, type="l", ylim=c(0,0.3), main="Enjalbert_2017 (control)", las=1, ylab="change of [Ace]", xlab="time (h)", col="#D6685C", lwd=2) +plot_points(time_meas, ace_nopulse, sd_ace_nopulse, offset=0.002, col="#D6685C") +plot(t_pulse, glc_conc_nopulse, type="l", ylim=c(-1.4,0), main="Enjalbert_2017 (control)", las=1, ylab="change of [Glc]", xlab="time (h)", col="#70AD47", lwd=2) +plot_points(time_meas, glc_nopulse, sd_glc_nopulse, offset=0.002, col="#70AD47") + +# Chemostat experiment (Renilla et al., 2012) +############################################# + +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_1.cps") + +# delete events and set glc feed to 1 +deleteEvent(getEvents()$key) +setGlobalQuantities(key = "_feed", initial_value = 1) + +n_step <- 50 +dilution_rates <- seq(0.1, 0.5, length.out = n_step) +fluxes <- c("Values[v_growth_rate]", "Values[v_glc_uptake]", "Values[v_ace_net]") +res_chemostat <- matrix(NA, nrow=n_step, ncol=length(fluxes)+1, dimnames=list(r=NULL, c=c("dilution_rate", fluxes))) +for (i in seq(n_step)){ + setGlobalQuantities(key = "_dilution_rate", initial_value = dilution_rates[i]) + res_ss <- runSteadyState() + res_chemostat[i,] <- c(dilution_rates[i], unlist(res_ss$global_quantities[res_ss$global_quantities$key %in% fluxes, "value"])) +} + +# acetate flux as function of dilution rate +dilution_rate <- c(0.09586056644880175, 0.20043572984749464, 0.2997821350762528, 0.3468409586056645, 0.39912854030501094, 0.44618736383442276, 0.5002178649237474) +q_ace <- c(0.40472342596168076, 0.7396156614846294, 1.477019302736899, 1.2154669794005626, 1.96934635755591, 2.0929143843289824, 2.006569318707304) +q_ace_sd <- c(0.120440687823878, 0.15316333195082, 0.04595662724122, 0.0274214232252601, 0.28885933172199, 0.0777639526513201, 0.28664731148963) + +plot(res_chemostat[,"dilution_rate"], res_chemostat[,"Values[v_ace_net]"], type="l", las=1, xlim=c(0,0.6), ylim=c(0,4), main="Renilla_2012", lwd=2, xlab="dilution rate (h-1)", ylab="ace flux (mmol/gDW/h)", col="#D6685C") +plot_points(dilution_rate, q_ace, q_ace_sd, offset=0.01, col="#D6685C") + +# acetate flux as function of glc uptake +q_glc <- c(1.43478260869565, 2.7391304347826, 4.59130434782608, 4.69565217391304, 5.7391304347826, 5.92173913043478, 6.10434782608695, 7.12173913043477, 8.34782608695652) +q_ace <- c(0.381358340437624, 0.762906128635029, 1.40509614473808, 1.17836506583309, 1.99602159704461, 2.08290233967983, 2.01292033721701, 2.44717249218528, 2.5848252344416) +q_ace_sd <- c(0.139547219854126, 0.139471440750213, 0.15697641375391, 0.1045751633987, 0.19179691200152, 0.0174291938997801, 0.24397082504499, 0.52295159609737, 0.33130624230368) +q_glc_sd <- c(0.20234760540444, 0.58899555210581, 0.36827270271402, 0.14607356337583, 0.630432119199, 0.13069807968161, 0.47630364392495, 0.99361847036576, 0.66185122193432) + +plot(res_chemostat[,"Values[v_glc_uptake]"], res_chemostat[,"Values[v_ace_net]"], type="l", las=1, xlim=c(0,10), ylim=c(0,4), main="Renilla_2012", lwd=2, xlab="glc uptake (mmol/gDW/h)", ylab="ace flux (mmol/gDW/h)", col="#D6685C") +plot_points(q_glc, q_ace, q_ace_sd, offset=0.2, col="#D6685C") +plot_points(q_glc, q_ace, q_glc_sd, offset=0.08, mode="h", col="#D6685C") + +# Steady-state fluxes for acetate concentration between 0.1 and 100mM (Enjalbert et al., 2017; Pinhal et al., 2019) +############################################################################################## + +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_1.cps") + +# remove events and fix concentrations of actate, glucose and biomass +deleteEvent(getEvents()$key) +setSpecies(key="Ace_out", type="fixed") +setSpecies(key="Glc", type="fixed") +setSpecies(key="X", type="fixed") + +# run simulations +n_step <- 50 +ace_concentration <- 10**seq(-1, 2, length.out = n_step) +fluxes <- c("Values[v_growth_rate]", "Values[v_glc_uptake]", "Values[v_ace_net]") +res_ace_range <- matrix(NA, nrow=n_step, ncol=length(fluxes)+1, dimnames=list(r=NULL, c=c("ace_concentration", fluxes))) +for (i in seq(n_step)){ + setSpecies(key="Ace_out", initial_concentration=ace_concentration[i]) + applyInitialState() + res_ss <- runSteadyState() + res_ace_range[i,] <- c(ace_concentration[i], unlist(res_ss$global_quantities[res_ss$global_quantities$key %in% fluxes, "value"])) +} + +xlab_main <- c(0.1, 1, 10, 100) +xlab_sec <- c(seq(0.2, 0.9, by=0.1), seq(2, 9, by=1), seq(20, 90, by=10)) + +# growth rate as function of acetate concentration +growth_rates <- c(0.521128511, 0.611148842625582, 0.613161998174498, 0.502533817, 0.496290415, 0.488201506, 0.547635665, 0.499830448, 0.474554197, 0.425356578, 0.377534684, 0.645724326, 0.618475601, 0.554887936, 0.564811523, 0.527571192, 0.434972836, 0.3824734, 0.583623355, 0.620905534, 0.564259247, 0.532148135, 0.483885215, 0.557074418, 0.630654409249223) +sd_growth_rates <- c(0.001793104, 0.00204807928657914, 0.00219396182484705, 0.001709207, 0.001846205, 0.001757403, 0.001821375, 0.001025702, 0.001940912, 0.001204707, 0.001999188, 0.001418374, 0.001932601, 0.001455791, 0.001574234, 0.001206265, 0.001292476, 0.001068259, 0.001804648, 0.001643459, 0.001598405, 0.001121218, 0.000912408, 0.00194896, 0.00203369597574686) +ace_conc <- c(0.451964624, 1.11600286648471, 2.04718732931708, 2.86252145, 5.285907977, 9.106164204, 16.67476528, 16.76626787, 17.00218707, 30.8667961, 57.92292091, 0.478352574, 4.55673229, 8.10163028, 8.22100734, 7.829591756, 33.53244905, 66.20361403, 0.436956014, 1.28468189, 1.555875222, 12.02564968, 30.24566673, 0.548011282, 2.29562227069566) +sd_ace_conc <- c(0.262318006, 0.281361538208953, 0.289527601555302, 0.302351163, 0.327705513, 0.330782201, 0.277259011, 0.233956798, 0.31883929, 0.300578057, 0.313371784, 0.231202155, 0.264687437, 0.243480317, 0.289821733, 0.263808862, 0.289478134, 0.264861034, 0.22461248, 0.229031308, 0.241718918, 0.254111384, 0.187394292, 0.011107606, 0.290519090995892) +plot(res_ace_range[,"ace_concentration"], xaxt='n', res_ace_range[,"Values[v_growth_rate]"], las=1, col="#2E75B6", lwd=2, type="l", xlim=c(0.1,100), ylim=c(0,0.8), log="x", main="Enjalbert_2017", xlab="[acetate] (mM)", ylab="growth rate (h-1)") +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +plot_points(ace_conc, growth_rates, sd_growth_rates, offset=0.01, col="#2E75B6") +plot_points(ace_conc, growth_rates, sd_ace_conc, offset=0.02, mode="h", col="#2E75B6") + +# glc uptake as function of acetate concentration +glc_upt <- c(8.654860167, 8.36127542981722, 7.98010111285252, 9.236935826, 8.274418986, 7.560431219, 7.339194455, 5.775312502, 6.423391263, 5.1544758, 3.938631573, 8.115447647, 9.28067031, 6.737153424, 7.172748804, 5.884186033, 5.684201497, 4.811576974, 9.632702365, 8.055042777, 9.708342814, 7.100081588, 5.505759496, 9.242859752, 8.18621623190759) +sd_glc_upt <- c(0.337812425, 0.38531328268303, 0.373770045721031, 0.356787032, 0.334672954, 0.317509322, 0.288025925, 0.16053276, 0.375934255, 0.293148172, 0.359225607, 0.197331684, 0.360984112, 0.229372278, 0.241820396, 0.20450532, 0.260869273, 0.216134352, 0.34289286, 0.350305744, 0.293144783, 0.220135755, 0.153471508, 0.25245346, 0.396815184905029) +ace_conc <- c(0.451964624, 1.11600286648471, 2.04718732931708, 2.86252145, 5.285907977, 9.106164204, 16.67476528, 16.76626787, 17.00218707, 30.8667961, 57.92292091, 0.478352574, 4.55673229, 8.10163028, 8.22100734, 7.829591756, 33.53244905, 66.20361403, 0.436956014, 1.28468189, 1.555875222, 12.02564968, 30.24566673, 0.548011282, 2.29562227069566) +sd_ace_conc <- c(0.262318006, 0.281361538208953, 0.289527601555302, 0.302351163, 0.327705513, 0.330782201, 0.277259011, 0.233956798, 0.31883929, 0.300578057, 0.313371784, 0.231202155, 0.264687437, 0.243480317, 0.289821733, 0.263808862, 0.289478134, 0.264861034, 0.22461248, 0.229031308, 0.241718918, 0.254111384, 0.187394292, 0.011107606, 0.290519090995892) +plot(res_ace_range[,"ace_concentration"], res_ace_range[,"Values[v_glc_uptake]"], xaxt='n', las=1, type="l", lwd=2, xlim=c(0.1,100), col="#70AD47", ylim=c(0,10), log="x", main="Enjalbert_2017", xlab="[acetate] (mM)", ylab="glc uptake (mmol/gDW/h)") +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +plot_points(ace_conc, glc_upt, sd_glc_upt, offset=0.02, col="#70AD47") +plot_points(ace_conc, glc_upt, sd_ace_conc, offset=0.2, mode="h", col="#70AD47") + +# ace flux as function of acetate concentration +ace_flx <- c(3.5, -2.7, 1.516999356, 1.26845123082679, 0.775821380507016, 0.678877137, 0.017366464, -0.991478151, -1.286687213, -2.078474994, -1.530841439, -1.525342269, -1.253581266, 1.984679487, 0.546462624, -0.136780389, -0.393883917, -0.610240984, -1.120767885, -1.277455315, 2.574285211, 2.051935093, 1.828415596, -1.262442483, -1.317987733, 2.333568565, 1.85234639824858) +sd_ace_flx <- c(0.35, 0.27, 0.316118066, 0.388752117258161, 0.40715278851436, 0.33357638, 0.37333751, 0.347029894, 0.280501612, 0.195031303, 0.376252463, 0.226182385, 0.303661317, 0.253610517, 0.385450715, 0.243880325, 0.30665695, 0.257983739, 0.23844407, 0.198458448, 0.299832036, 0.334956504, 0.317134334, 0.263807154, 0.195219648, 0.016120887, 0.386174129654754) +ace_conc <- c(0.2, 100, 0.451964624, 1.11600286648471, 2.04718732931708, 2.86252145, 5.285907977, 9.106164204, 16.67476528, 16.76626787, 17.00218707, 30.8667961, 57.92292091, 0.478352574, 4.55673229, 8.10163028, 8.22100734, 7.829591756, 33.53244905, 66.20361403, 0.436956014, 1.28468189, 1.555875222, 12.02564968, 30.24566673, 0.548011282, 2.29562227069566) +sd_ace_conc <- c(0.05, 10, 0.262318006, 0.281361538208953, 0.289527601555302, 0.302351163, 0.327705513, 0.330782201, 0.277259011, 0.233956798, 0.31883929, 0.300578057, 0.313371784, 0.231202155, 0.264687437, 0.243480317, 0.289821733, 0.263808862, 0.289478134, 0.264861034, 0.22461248, 0.229031308, 0.241718918, 0.254111384, 0.187394292, 0.011107606, 0.290519090995892) +plot(res_ace_range[,"ace_concentration"], res_ace_range[,"Values[v_ace_net]"], xaxt='n', las=1, type="l", lwd=2, xlim=c(0.1,100), ylim=c(-4,4), log="x", main="Enjalbert_2017, Pinhal_2019", col="#D6685C", xlab="[acetate] (mM)", ylab="ace flux (mmol/gDW/h)") +abline(h=0) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +plot_points(ace_conc, ace_flx, sd_ace_flx, offset=0.04, col="#D6685C") +plot_points(ace_conc, ace_flx, sd_ace_conc, offset=0.2, mode="h", col="#D6685C") + +dev.off() + +########### +# Model 2 # +########### + +setwd(results_dir) +pdf(file="Figure 3-figure supplement 2.pdf", width = 7, height = 9) +par(mfrow=c(4,3)) + +# Pulse experiment (Enjalbert et al., 2017) +########################################### + +setwd(model_dir) +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_2.cps") + +# set biomass concentration at which the pulse is performed +setGlobalQuantities(key = "_X_conc_pulse", initial_value = 0.9) + +# simulate response to acetate pulse +res_pulse <- runTimeCourse() + +# get simulation results +id_start <- which(res_pulse$result$Ace_out >= 30)[1]+1 +id_end <- which(res_pulse$result$Time >= (res_pulse$result$Time[id_start]+8/60))[1]+1 + +t_pulse <- (res_pulse$result$Time[seq(id_start, id_end)] - res_pulse$result$Time[id_start]) +ace_conc_pulse <- res_pulse$result$Ace_out[seq(id_start, id_end)] - res_pulse$result$Ace_out[id_start] +glc_conc_pulse <- res_pulse$result$Glc[seq(id_start, id_end)] - res_pulse$result$Glc[id_start] + +# measurements +time_meas <- seq(0,8)/60 +glc <- c(0, -0.06931183, -0.151415145, -0.189227994, -0.269451057, -0.290764495, -0.230785281, -0.464084162, -0.551878527) +sd_glc <- c(0, 0.032476344, 0.073133915, 0.113018846, 0.049485284, 0.005325541, 0.163377704, 0.034786419, 0.048477157) +ace <- c(0, -0.027907926, -0.078000853, -0.155334163, -0.165031608, -0.111098424, -0.182877548, -0.237262298, -0.276903255) +sd_ace <- c(0, 0.002740145, 0.025693594, 0.053641876, 0.089975321, 0.005875669, 0.085604161, 0.061930626, 0.099140975) + +# plot simulations vs measurements +plot(t_pulse, ace_conc_pulse, type="l", ylim=c(-0.4,0), las=1, main="Enjalbert_2017 (pulse 30mM ace)", ylab="change of [Ace]", xlab="time (h)", col="#D6685C", lwd=2) +plot_points(time_meas, ace, sd_ace, offset=0.002, col="#D6685C") +plot(t_pulse, glc_conc_pulse, type="l", main="Enjalbert_2017 (pulse 30mM ace)", las=1, ylab="change of [Glc]", xlab="time (h)", col="#70AD47", lwd=2) +plot_points(time_meas, glc, sd_glc, offset=0.002, col="#70AD47") + +# simulate control experiment (i.e. no acetate pulse) +deleteEvent(getEvents()$key) + +res_nopulse <- runTimeCourse() + +# get simulation results +ace_conc_nopulse <- res_nopulse$result$Ace_out[seq(id_start, id_end)] - res_nopulse$result$Ace_out[id_start] +glc_conc_nopulse <- res_nopulse$result$Glc[seq(id_start, id_end)] - res_nopulse$result$Glc[id_start] + +# measurements +glc_nopulse <- c(0, -0.195637774, -0.325512845, -0.433785334, -0.628518958, -0.726913548, -0.892908748, -1.071230579, -1.16285575) +sd_glc_nopulse <- c(0, 0.058023617, 0.109115511, 0.047090371, 0.052331837, 0.065619906, 0.133896355, 0.16828754, 0.03515465) +ace_nopulse <- c(0, 0.01012067, 0.05974009, 0.086787283, 0.086690257, 0.104329693, 0.128507087, 0.130026354, 0.187336965) +sd_ace_nopulse <- c(0, 0.005117257, 0.022507218, 0.040319045, 0.037548873, 0.037235285, 0.044494365, 0.045029023, 0.023982374) + +# plot simulations vs measurements +plot(t_pulse, ace_conc_nopulse, type="l", ylim=c(0,0.3), main="Enjalbert_2017 (control)", las=1, ylab="change of [Ace]", xlab="time (h)", col="#D6685C", lwd=2) +plot_points(time_meas, ace_nopulse, sd_ace_nopulse, offset=0.002, col="#D6685C") +plot(t_pulse, glc_conc_nopulse, type="l", ylim=c(-1.4,0), main="Enjalbert_2017 (control)", las=1, ylab="change of [Glc]", xlab="time (h)", col="#70AD47", lwd=2) +plot_points(time_meas, glc_nopulse, sd_glc_nopulse, offset=0.002, col="#70AD47") + +# Chemostat experiment (Renilla et al., 2012) +############################################# + +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_2.cps") + +# delete events and set glc feed to 1 +deleteEvent(getEvents()$key) +setGlobalQuantities(key = "_feed", initial_value = 1) + +n_step <- 50 +dilution_rates <- seq(0.1, 0.5, length.out = n_step) +fluxes <- c("Values[v_growth_rate]", "Values[v_glc_uptake]", "Values[v_ace_net]") +res_chemostat <- matrix(NA, nrow=n_step, ncol=length(fluxes)+1, dimnames=list(r=NULL, c=c("dilution_rate", fluxes))) +for (i in seq(n_step)){ + setGlobalQuantities(key = "_dilution_rate", initial_value = dilution_rates[i]) + res_ss <- runSteadyState() + res_chemostat[i,] <- c(dilution_rates[i], unlist(res_ss$global_quantities[res_ss$global_quantities$key %in% fluxes, "value"])) +} + +# acetate flux as function of dilution rate +dilution_rate <- c(0.09586056644880175, 0.20043572984749464, 0.2997821350762528, 0.3468409586056645, 0.39912854030501094, 0.44618736383442276, 0.5002178649237474) +q_ace <- c(0.40472342596168076, 0.7396156614846294, 1.477019302736899, 1.2154669794005626, 1.96934635755591, 2.0929143843289824, 2.006569318707304) +q_ace_sd <- c(0.120440687823878, 0.15316333195082, 0.04595662724122, 0.0274214232252601, 0.28885933172199, 0.0777639526513201, 0.28664731148963) + +plot(res_chemostat[,"dilution_rate"], res_chemostat[,"Values[v_ace_net]"], type="l", las=1, xlim=c(0,0.6), ylim=c(0,4), main="Renilla_2012", lwd=2, xlab="dilution rate (h-1)", ylab="ace flux (mmol/gDW/h)", col="#D6685C") +plot_points(dilution_rate, q_ace, q_ace_sd, offset=0.01, col="#D6685C") + +# acetate flux as function of glc uptake +q_glc <- c(1.43478260869565, 2.7391304347826, 4.59130434782608, 4.69565217391304, 5.7391304347826, 5.92173913043478, 6.10434782608695, 7.12173913043477, 8.34782608695652) +q_ace <- c(0.381358340437624, 0.762906128635029, 1.40509614473808, 1.17836506583309, 1.99602159704461, 2.08290233967983, 2.01292033721701, 2.44717249218528, 2.5848252344416) +q_ace_sd <- c(0.139547219854126, 0.139471440750213, 0.15697641375391, 0.1045751633987, 0.19179691200152, 0.0174291938997801, 0.24397082504499, 0.52295159609737, 0.33130624230368) +q_glc_sd <- c(0.20234760540444, 0.58899555210581, 0.36827270271402, 0.14607356337583, 0.630432119199, 0.13069807968161, 0.47630364392495, 0.99361847036576, 0.66185122193432) + +plot(res_chemostat[,"Values[v_glc_uptake]"], res_chemostat[,"Values[v_ace_net]"], type="l", las=1, xlim=c(0,10), ylim=c(0,4), main="Renilla_2012", lwd=2, xlab="glc uptake (mmol/gDW/h)", ylab="ace flux (mmol/gDW/h)", col="#D6685C") +plot_points(q_glc, q_ace, q_ace_sd, offset=0.2, col="#D6685C") +plot_points(q_glc, q_ace, q_glc_sd, offset=0.08, mode="h", col="#D6685C") + +# Steady-state fluxes for acetate concentration between 0.1 and 100mM (Enjalbert et al., 2017; Pinhal et al., 2019) +############################################################################################## + +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_2.cps") + +# remove events and fix concentrations of actate, glucose and biomass +deleteEvent(getEvents()$key) +setSpecies(key="Ace_out", type="fixed") +setSpecies(key="Glc", type="fixed") +setSpecies(key="X", type="fixed") + +# run simulations +n_step <- 50 +ace_concentration <- 10**seq(-1, 2, length.out = n_step) +fluxes <- c("Values[v_growth_rate]", "Values[v_glc_uptake]", "Values[v_ace_net]") +res_ace_range <- matrix(NA, nrow=n_step, ncol=length(fluxes)+1, dimnames=list(r=NULL, c=c("ace_concentration", fluxes))) +for (i in seq(n_step)){ + setSpecies(key="Ace_out", initial_concentration=ace_concentration[i]) + applyInitialState() + res_ss <- runSteadyState() + res_ace_range[i,] <- c(ace_concentration[i], unlist(res_ss$global_quantities[res_ss$global_quantities$key %in% fluxes, "value"])) +} + +xlab_main <- c(0.1, 1, 10, 100) +xlab_sec <- c(seq(0.2, 0.9, by=0.1), seq(2, 9, by=1), seq(20, 90, by=10)) + +# growth rate as function of acetate concentration +growth_rates <- c(0.521128511, 0.611148842625582, 0.613161998174498, 0.502533817, 0.496290415, 0.488201506, 0.547635665, 0.499830448, 0.474554197, 0.425356578, 0.377534684, 0.645724326, 0.618475601, 0.554887936, 0.564811523, 0.527571192, 0.434972836, 0.3824734, 0.583623355, 0.620905534, 0.564259247, 0.532148135, 0.483885215, 0.557074418, 0.630654409249223) +sd_growth_rates <- c(0.001793104, 0.00204807928657914, 0.00219396182484705, 0.001709207, 0.001846205, 0.001757403, 0.001821375, 0.001025702, 0.001940912, 0.001204707, 0.001999188, 0.001418374, 0.001932601, 0.001455791, 0.001574234, 0.001206265, 0.001292476, 0.001068259, 0.001804648, 0.001643459, 0.001598405, 0.001121218, 0.000912408, 0.00194896, 0.00203369597574686) +ace_conc <- c(0.451964624, 1.11600286648471, 2.04718732931708, 2.86252145, 5.285907977, 9.106164204, 16.67476528, 16.76626787, 17.00218707, 30.8667961, 57.92292091, 0.478352574, 4.55673229, 8.10163028, 8.22100734, 7.829591756, 33.53244905, 66.20361403, 0.436956014, 1.28468189, 1.555875222, 12.02564968, 30.24566673, 0.548011282, 2.29562227069566) +sd_ace_conc <- c(0.262318006, 0.281361538208953, 0.289527601555302, 0.302351163, 0.327705513, 0.330782201, 0.277259011, 0.233956798, 0.31883929, 0.300578057, 0.313371784, 0.231202155, 0.264687437, 0.243480317, 0.289821733, 0.263808862, 0.289478134, 0.264861034, 0.22461248, 0.229031308, 0.241718918, 0.254111384, 0.187394292, 0.011107606, 0.290519090995892) +plot(res_ace_range[,"ace_concentration"], xaxt='n', res_ace_range[,"Values[v_growth_rate]"], las=1, col="#2E75B6", lwd=2, type="l", xlim=c(0.1,100), ylim=c(0,0.8), log="x", main="Enjalbert_2017", xlab="[acetate] (mM)", ylab="growth rate (h-1)") +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +plot_points(ace_conc, growth_rates, sd_growth_rates, offset=0.01, col="#2E75B6") +plot_points(ace_conc, growth_rates, sd_ace_conc, offset=0.02, mode="h", col="#2E75B6") + +# glc uptake as function of acetate concentration +glc_upt <- c(8.654860167, 8.36127542981722, 7.98010111285252, 9.236935826, 8.274418986, 7.560431219, 7.339194455, 5.775312502, 6.423391263, 5.1544758, 3.938631573, 8.115447647, 9.28067031, 6.737153424, 7.172748804, 5.884186033, 5.684201497, 4.811576974, 9.632702365, 8.055042777, 9.708342814, 7.100081588, 5.505759496, 9.242859752, 8.18621623190759) +sd_glc_upt <- c(0.337812425, 0.38531328268303, 0.373770045721031, 0.356787032, 0.334672954, 0.317509322, 0.288025925, 0.16053276, 0.375934255, 0.293148172, 0.359225607, 0.197331684, 0.360984112, 0.229372278, 0.241820396, 0.20450532, 0.260869273, 0.216134352, 0.34289286, 0.350305744, 0.293144783, 0.220135755, 0.153471508, 0.25245346, 0.396815184905029) +ace_conc <- c(0.451964624, 1.11600286648471, 2.04718732931708, 2.86252145, 5.285907977, 9.106164204, 16.67476528, 16.76626787, 17.00218707, 30.8667961, 57.92292091, 0.478352574, 4.55673229, 8.10163028, 8.22100734, 7.829591756, 33.53244905, 66.20361403, 0.436956014, 1.28468189, 1.555875222, 12.02564968, 30.24566673, 0.548011282, 2.29562227069566) +sd_ace_conc <- c(0.262318006, 0.281361538208953, 0.289527601555302, 0.302351163, 0.327705513, 0.330782201, 0.277259011, 0.233956798, 0.31883929, 0.300578057, 0.313371784, 0.231202155, 0.264687437, 0.243480317, 0.289821733, 0.263808862, 0.289478134, 0.264861034, 0.22461248, 0.229031308, 0.241718918, 0.254111384, 0.187394292, 0.011107606, 0.290519090995892) +plot(res_ace_range[,"ace_concentration"], res_ace_range[,"Values[v_glc_uptake]"], xaxt='n', las=1, type="l", lwd=2, xlim=c(0.1,100), col="#70AD47", ylim=c(0,10), log="x", main="Enjalbert_2017", xlab="[acetate] (mM)", ylab="glc uptake (mmol/gDW/h)") +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +plot_points(ace_conc, glc_upt, sd_glc_upt, offset=0.02, col="#70AD47") +plot_points(ace_conc, glc_upt, sd_ace_conc, offset=0.2, mode="h", col="#70AD47") + +# ace flux as function of acetate concentration +ace_flx <- c(3.5, -2.7, 1.516999356, 1.26845123082679, 0.775821380507016, 0.678877137, 0.017366464, -0.991478151, -1.286687213, -2.078474994, -1.530841439, -1.525342269, -1.253581266, 1.984679487, 0.546462624, -0.136780389, -0.393883917, -0.610240984, -1.120767885, -1.277455315, 2.574285211, 2.051935093, 1.828415596, -1.262442483, -1.317987733, 2.333568565, 1.85234639824858) +sd_ace_flx <- c(0.35, 0.27, 0.316118066, 0.388752117258161, 0.40715278851436, 0.33357638, 0.37333751, 0.347029894, 0.280501612, 0.195031303, 0.376252463, 0.226182385, 0.303661317, 0.253610517, 0.385450715, 0.243880325, 0.30665695, 0.257983739, 0.23844407, 0.198458448, 0.299832036, 0.334956504, 0.317134334, 0.263807154, 0.195219648, 0.016120887, 0.386174129654754) +ace_conc <- c(0.2, 100, 0.451964624, 1.11600286648471, 2.04718732931708, 2.86252145, 5.285907977, 9.106164204, 16.67476528, 16.76626787, 17.00218707, 30.8667961, 57.92292091, 0.478352574, 4.55673229, 8.10163028, 8.22100734, 7.829591756, 33.53244905, 66.20361403, 0.436956014, 1.28468189, 1.555875222, 12.02564968, 30.24566673, 0.548011282, 2.29562227069566) +sd_ace_conc <- c(0.05, 10, 0.262318006, 0.281361538208953, 0.289527601555302, 0.302351163, 0.327705513, 0.330782201, 0.277259011, 0.233956798, 0.31883929, 0.300578057, 0.313371784, 0.231202155, 0.264687437, 0.243480317, 0.289821733, 0.263808862, 0.289478134, 0.264861034, 0.22461248, 0.229031308, 0.241718918, 0.254111384, 0.187394292, 0.011107606, 0.290519090995892) +plot(res_ace_range[,"ace_concentration"], res_ace_range[,"Values[v_ace_net]"], xaxt='n', las=1, type="l", lwd=2, xlim=c(0.1,100), ylim=c(-4,4), log="x", main="Enjalbert_2017, Pinhal_2019", col="#D6685C", xlab="[acetate] (mM)", ylab="ace flux (mmol/gDW/h)") +abline(h=0) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +plot_points(ace_conc, ace_flx, sd_ace_flx, offset=0.04, col="#D6685C") +plot_points(ace_conc, ace_flx, sd_ace_conc, offset=0.2, mode="h", col="#D6685C") + +dev.off() + +########### +# Model 3 # +########### + +setwd(results_dir) +pdf(file="Figure 3-figure supplement 3.pdf", width = 7, height = 9) +par(mfrow=c(4,3)) + +# Pulse experiment (Enjalbert et al., 2017) +########################################### + +setwd(model_dir) +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_3.cps") + +# set biomass concentration at which the pulse is performed +setGlobalQuantities(key = "_X_conc_pulse", initial_value = 0.9) + +# simulate response to acetate pulse +res_pulse <- runTimeCourse() + +# get simulation results +id_start <- which(res_pulse$result$Ace_out >= 30)[1]+1 +id_end <- which(res_pulse$result$Time >= (res_pulse$result$Time[id_start]+8/60))[1]+1 + +t_pulse <- (res_pulse$result$Time[seq(id_start, id_end)] - res_pulse$result$Time[id_start]) +ace_conc_pulse <- res_pulse$result$Ace_out[seq(id_start, id_end)] - res_pulse$result$Ace_out[id_start] +glc_conc_pulse <- res_pulse$result$Glc[seq(id_start, id_end)] - res_pulse$result$Glc[id_start] + +# measurements +time_meas <- seq(0,8)/60 +glc <- c(0, -0.06931183, -0.151415145, -0.189227994, -0.269451057, -0.290764495, -0.230785281, -0.464084162, -0.551878527) +sd_glc <- c(0, 0.032476344, 0.073133915, 0.113018846, 0.049485284, 0.005325541, 0.163377704, 0.034786419, 0.048477157) +ace <- c(0, -0.027907926, -0.078000853, -0.155334163, -0.165031608, -0.111098424, -0.182877548, -0.237262298, -0.276903255) +sd_ace <- c(0, 0.002740145, 0.025693594, 0.053641876, 0.089975321, 0.005875669, 0.085604161, 0.061930626, 0.099140975) + +# plot simulations vs measurements +plot(t_pulse, ace_conc_pulse, type="l", ylim=c(-0.4,0), las=1, main="Enjalbert_2017 (pulse 30mM ace)", ylab="change of [Ace]", xlab="time (h)", col="#D6685C", lwd=2) +plot_points(time_meas, ace, sd_ace, offset=0.002, col="#D6685C") +plot(t_pulse, glc_conc_pulse, type="l", main="Enjalbert_2017 (pulse 30mM ace)", las=1, ylab="change of [Glc]", xlab="time (h)", col="#70AD47", lwd=2) +plot_points(time_meas, glc, sd_glc, offset=0.002, col="#70AD47") + +# simulate control experiment (i.e. no acetate pulse) +deleteEvent(getEvents()$key) + +res_nopulse <- runTimeCourse() + +# get simulation results +ace_conc_nopulse <- res_nopulse$result$Ace_out[seq(id_start, id_end)] - res_nopulse$result$Ace_out[id_start] +glc_conc_nopulse <- res_nopulse$result$Glc[seq(id_start, id_end)] - res_nopulse$result$Glc[id_start] + +# measurements +glc_nopulse <- c(0, -0.195637774, -0.325512845, -0.433785334, -0.628518958, -0.726913548, -0.892908748, -1.071230579, -1.16285575) +sd_glc_nopulse <- c(0, 0.058023617, 0.109115511, 0.047090371, 0.052331837, 0.065619906, 0.133896355, 0.16828754, 0.03515465) +ace_nopulse <- c(0, 0.01012067, 0.05974009, 0.086787283, 0.086690257, 0.104329693, 0.128507087, 0.130026354, 0.187336965) +sd_ace_nopulse <- c(0, 0.005117257, 0.022507218, 0.040319045, 0.037548873, 0.037235285, 0.044494365, 0.045029023, 0.023982374) + +# plot simulations vs measurements +plot(t_pulse, ace_conc_nopulse, type="l", ylim=c(0,0.3), main="Enjalbert_2017 (control)", las=1, ylab="change of [Ace]", xlab="time (h)", col="#D6685C", lwd=2) +plot_points(time_meas, ace_nopulse, sd_ace_nopulse, offset=0.002, col="#D6685C") +plot(t_pulse, glc_conc_nopulse, type="l", ylim=c(-1.4,0), main="Enjalbert_2017 (control)", las=1, ylab="change of [Glc]", xlab="time (h)", col="#70AD47", lwd=2) +plot_points(time_meas, glc_nopulse, sd_glc_nopulse, offset=0.002, col="#70AD47") + +# Chemostat experiment (Renilla et al., 2012) +############################################# + +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_3.cps") + +# delete events and set glc feed to 1 +deleteEvent(getEvents()$key) +setGlobalQuantities(key = "_feed", initial_value = 1) + +n_step <- 50 +dilution_rates <- seq(0.1, 0.5, length.out = n_step) +fluxes <- c("Values[v_growth_rate]", "Values[v_glc_uptake]", "Values[v_ace_net]") +res_chemostat <- matrix(NA, nrow=n_step, ncol=length(fluxes)+1, dimnames=list(r=NULL, c=c("dilution_rate", fluxes))) +for (i in seq(n_step)){ + setGlobalQuantities(key = "_dilution_rate", initial_value = dilution_rates[i]) + res_ss <- runSteadyState() + res_chemostat[i,] <- c(dilution_rates[i], unlist(res_ss$global_quantities[res_ss$global_quantities$key %in% fluxes, "value"])) +} + +# acetate flux as function of dilution rate +dilution_rate <- c(0.09586056644880175, 0.20043572984749464, 0.2997821350762528, 0.3468409586056645, 0.39912854030501094, 0.44618736383442276, 0.5002178649237474) +q_ace <- c(0.40472342596168076, 0.7396156614846294, 1.477019302736899, 1.2154669794005626, 1.96934635755591, 2.0929143843289824, 2.006569318707304) +q_ace_sd <- c(0.120440687823878, 0.15316333195082, 0.04595662724122, 0.0274214232252601, 0.28885933172199, 0.0777639526513201, 0.28664731148963) + +plot(res_chemostat[,"dilution_rate"], res_chemostat[,"Values[v_ace_net]"], type="l", las=1, xlim=c(0,0.6), ylim=c(0,4), main="Renilla_2012", lwd=2, xlab="dilution rate (h-1)", ylab="ace flux (mmol/gDW/h)", col="#D6685C") +plot_points(dilution_rate, q_ace, q_ace_sd, offset=0.01, col="#D6685C") + +# acetate flux as function of glc uptake +q_glc <- c(1.43478260869565, 2.7391304347826, 4.59130434782608, 4.69565217391304, 5.7391304347826, 5.92173913043478, 6.10434782608695, 7.12173913043477, 8.34782608695652) +q_ace <- c(0.381358340437624, 0.762906128635029, 1.40509614473808, 1.17836506583309, 1.99602159704461, 2.08290233967983, 2.01292033721701, 2.44717249218528, 2.5848252344416) +q_ace_sd <- c(0.139547219854126, 0.139471440750213, 0.15697641375391, 0.1045751633987, 0.19179691200152, 0.0174291938997801, 0.24397082504499, 0.52295159609737, 0.33130624230368) +q_glc_sd <- c(0.20234760540444, 0.58899555210581, 0.36827270271402, 0.14607356337583, 0.630432119199, 0.13069807968161, 0.47630364392495, 0.99361847036576, 0.66185122193432) + +plot(res_chemostat[,"Values[v_glc_uptake]"], res_chemostat[,"Values[v_ace_net]"], type="l", las=1, xlim=c(0,10), ylim=c(0,4), main="Renilla_2012", lwd=2, xlab="glc uptake (mmol/gDW/h)", ylab="ace flux (mmol/gDW/h)", col="#D6685C") +plot_points(q_glc, q_ace, q_ace_sd, offset=0.2, col="#D6685C") +plot_points(q_glc, q_ace, q_glc_sd, offset=0.08, mode="h", col="#D6685C") + +# Steady-state fluxes for acetate concentration between 0.1 and 100mM (Enjalbert et al., 2017; Pinhal et al., 2019) +############################################################################################## + +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_3.cps") + +# remove events and fix concentrations of actate, glucose and biomass +deleteEvent(getEvents()$key) +setSpecies(key="Ace_out", type="fixed") +setSpecies(key="Glc", type="fixed") +setSpecies(key="X", type="fixed") + +# run simulations +n_step <- 50 +ace_concentration <- 10**seq(-1, 2, length.out = n_step) +fluxes <- c("Values[v_growth_rate]", "Values[v_glc_uptake]", "Values[v_ace_net]") +res_ace_range <- matrix(NA, nrow=n_step, ncol=length(fluxes)+1, dimnames=list(r=NULL, c=c("ace_concentration", fluxes))) +for (i in seq(n_step)){ + setSpecies(key="Ace_out", initial_concentration=ace_concentration[i]) + applyInitialState() + res_ss <- runSteadyState() + res_ace_range[i,] <- c(ace_concentration[i], unlist(res_ss$global_quantities[res_ss$global_quantities$key %in% fluxes, "value"])) +} + +xlab_main <- c(0.1, 1, 10, 100) +xlab_sec <- c(seq(0.2, 0.9, by=0.1), seq(2, 9, by=1), seq(20, 90, by=10)) + +# growth rate as function of acetate concentration +growth_rates <- c(0.521128511, 0.611148842625582, 0.613161998174498, 0.502533817, 0.496290415, 0.488201506, 0.547635665, 0.499830448, 0.474554197, 0.425356578, 0.377534684, 0.645724326, 0.618475601, 0.554887936, 0.564811523, 0.527571192, 0.434972836, 0.3824734, 0.583623355, 0.620905534, 0.564259247, 0.532148135, 0.483885215, 0.557074418, 0.630654409249223) +sd_growth_rates <- c(0.001793104, 0.00204807928657914, 0.00219396182484705, 0.001709207, 0.001846205, 0.001757403, 0.001821375, 0.001025702, 0.001940912, 0.001204707, 0.001999188, 0.001418374, 0.001932601, 0.001455791, 0.001574234, 0.001206265, 0.001292476, 0.001068259, 0.001804648, 0.001643459, 0.001598405, 0.001121218, 0.000912408, 0.00194896, 0.00203369597574686) +ace_conc <- c(0.451964624, 1.11600286648471, 2.04718732931708, 2.86252145, 5.285907977, 9.106164204, 16.67476528, 16.76626787, 17.00218707, 30.8667961, 57.92292091, 0.478352574, 4.55673229, 8.10163028, 8.22100734, 7.829591756, 33.53244905, 66.20361403, 0.436956014, 1.28468189, 1.555875222, 12.02564968, 30.24566673, 0.548011282, 2.29562227069566) +sd_ace_conc <- c(0.262318006, 0.281361538208953, 0.289527601555302, 0.302351163, 0.327705513, 0.330782201, 0.277259011, 0.233956798, 0.31883929, 0.300578057, 0.313371784, 0.231202155, 0.264687437, 0.243480317, 0.289821733, 0.263808862, 0.289478134, 0.264861034, 0.22461248, 0.229031308, 0.241718918, 0.254111384, 0.187394292, 0.011107606, 0.290519090995892) +plot(res_ace_range[,"ace_concentration"], xaxt='n', res_ace_range[,"Values[v_growth_rate]"], las=1, col="#2E75B6", lwd=2, type="l", xlim=c(0.1,100), ylim=c(0,0.8), log="x", main="Enjalbert_2017", xlab="[acetate] (mM)", ylab="growth rate (h-1)") +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +plot_points(ace_conc, growth_rates, sd_growth_rates, offset=0.01, col="#2E75B6") +plot_points(ace_conc, growth_rates, sd_ace_conc, offset=0.02, mode="h", col="#2E75B6") + +# glc uptake as function of acetate concentration +glc_upt <- c(8.654860167, 8.36127542981722, 7.98010111285252, 9.236935826, 8.274418986, 7.560431219, 7.339194455, 5.775312502, 6.423391263, 5.1544758, 3.938631573, 8.115447647, 9.28067031, 6.737153424, 7.172748804, 5.884186033, 5.684201497, 4.811576974, 9.632702365, 8.055042777, 9.708342814, 7.100081588, 5.505759496, 9.242859752, 8.18621623190759) +sd_glc_upt <- c(0.337812425, 0.38531328268303, 0.373770045721031, 0.356787032, 0.334672954, 0.317509322, 0.288025925, 0.16053276, 0.375934255, 0.293148172, 0.359225607, 0.197331684, 0.360984112, 0.229372278, 0.241820396, 0.20450532, 0.260869273, 0.216134352, 0.34289286, 0.350305744, 0.293144783, 0.220135755, 0.153471508, 0.25245346, 0.396815184905029) +ace_conc <- c(0.451964624, 1.11600286648471, 2.04718732931708, 2.86252145, 5.285907977, 9.106164204, 16.67476528, 16.76626787, 17.00218707, 30.8667961, 57.92292091, 0.478352574, 4.55673229, 8.10163028, 8.22100734, 7.829591756, 33.53244905, 66.20361403, 0.436956014, 1.28468189, 1.555875222, 12.02564968, 30.24566673, 0.548011282, 2.29562227069566) +sd_ace_conc <- c(0.262318006, 0.281361538208953, 0.289527601555302, 0.302351163, 0.327705513, 0.330782201, 0.277259011, 0.233956798, 0.31883929, 0.300578057, 0.313371784, 0.231202155, 0.264687437, 0.243480317, 0.289821733, 0.263808862, 0.289478134, 0.264861034, 0.22461248, 0.229031308, 0.241718918, 0.254111384, 0.187394292, 0.011107606, 0.290519090995892) +plot(res_ace_range[,"ace_concentration"], res_ace_range[,"Values[v_glc_uptake]"], xaxt='n', las=1, type="l", lwd=2, xlim=c(0.1,100), col="#70AD47", ylim=c(0,10), log="x", main="Enjalbert_2017", xlab="[acetate] (mM)", ylab="glc uptake (mmol/gDW/h)") +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +plot_points(ace_conc, glc_upt, sd_glc_upt, offset=0.02, col="#70AD47") +plot_points(ace_conc, glc_upt, sd_ace_conc, offset=0.2, mode="h", col="#70AD47") + +# ace flux as function of acetate concentration +ace_flx <- c(3.5, -2.7, 1.516999356, 1.26845123082679, 0.775821380507016, 0.678877137, 0.017366464, -0.991478151, -1.286687213, -2.078474994, -1.530841439, -1.525342269, -1.253581266, 1.984679487, 0.546462624, -0.136780389, -0.393883917, -0.610240984, -1.120767885, -1.277455315, 2.574285211, 2.051935093, 1.828415596, -1.262442483, -1.317987733, 2.333568565, 1.85234639824858) +sd_ace_flx <- c(0.35, 0.27, 0.316118066, 0.388752117258161, 0.40715278851436, 0.33357638, 0.37333751, 0.347029894, 0.280501612, 0.195031303, 0.376252463, 0.226182385, 0.303661317, 0.253610517, 0.385450715, 0.243880325, 0.30665695, 0.257983739, 0.23844407, 0.198458448, 0.299832036, 0.334956504, 0.317134334, 0.263807154, 0.195219648, 0.016120887, 0.386174129654754) +ace_conc <- c(0.2, 100, 0.451964624, 1.11600286648471, 2.04718732931708, 2.86252145, 5.285907977, 9.106164204, 16.67476528, 16.76626787, 17.00218707, 30.8667961, 57.92292091, 0.478352574, 4.55673229, 8.10163028, 8.22100734, 7.829591756, 33.53244905, 66.20361403, 0.436956014, 1.28468189, 1.555875222, 12.02564968, 30.24566673, 0.548011282, 2.29562227069566) +sd_ace_conc <- c(0.05, 10, 0.262318006, 0.281361538208953, 0.289527601555302, 0.302351163, 0.327705513, 0.330782201, 0.277259011, 0.233956798, 0.31883929, 0.300578057, 0.313371784, 0.231202155, 0.264687437, 0.243480317, 0.289821733, 0.263808862, 0.289478134, 0.264861034, 0.22461248, 0.229031308, 0.241718918, 0.254111384, 0.187394292, 0.011107606, 0.290519090995892) +plot(res_ace_range[,"ace_concentration"], res_ace_range[,"Values[v_ace_net]"], xaxt='n', las=1, type="l", lwd=2, xlim=c(0.1,100), ylim=c(-4,4), log="x", main="Enjalbert_2017, Pinhal_2019", col="#D6685C", xlab="[acetate] (mM)", ylab="ace flux (mmol/gDW/h)") +abline(h=0) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +plot_points(ace_conc, ace_flx, sd_ace_flx, offset=0.04, col="#D6685C") +plot_points(ace_conc, ace_flx, sd_ace_conc, offset=0.2, mode="h", col="#D6685C") + +dev.off() + + +################################################## +# Calculate variance weighted SSR for all models # +################################################## + +# Pulse experiment (Enjalbert et al., 2017) +########################################### + +# measurements +# sampling time +time_meas <- seq(0,8)/60 +# acetate pulse experiment +glc <- c(0, -0.06931183, -0.151415145, -0.189227994, -0.269451057, -0.290764495, -0.230785281, -0.464084162, -0.551878527) +sd_glc <- c(0.1, 0.032476344, 0.073133915, 0.113018846, 0.049485284, 0.005325541, 0.163377704, 0.034786419, 0.048477157) +ace <- c(0, -0.027907926, -0.078000853, -0.155334163, -0.165031608, -0.111098424, -0.182877548, -0.237262298, -0.276903255) +sd_ace <- c(0.01, 0.002740145, 0.025693594, 0.053641876, 0.089975321, 0.005875669, 0.085604161, 0.061930626, 0.099140975) +# control experiment (pulse without acetate) +glc_nopulse <- c(0, -0.195637774, -0.325512845, -0.433785334, -0.628518958, -0.726913548, -0.892908748, -1.071230579, -1.16285575) +sd_glc_nopulse <- c(0.01, 0.058023617, 0.109115511, 0.047090371, 0.052331837, 0.065619906, 0.133896355, 0.16828754, 0.03515465) +ace_nopulse <- c(0, 0.01012067, 0.05974009, 0.086787283, 0.086690257, 0.104329693, 0.128507087, 0.130026354, 0.187336965) +sd_ace_nopulse <- c(0.01, 0.005117257, 0.022507218, 0.040319045, 0.037548873, 0.037235285, 0.044494365, 0.045029023, 0.023982374) + +setwd(model_dir) +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_4.cps") +# set biomass concentration at which the pulse is performed +setGlobalQuantities(key = "_X_conc_pulse", initial_value = 0.9) +# simulate response to acetate pulse +res_pulse <- runTimeCourse(intervals=10000) +# get simulation results +id_start <- which(res_pulse$result$Ace_out >= 30)[1]+1 +id_end <- which(res_pulse$result$Time >= (res_pulse$result$Time[id_start]+8/60))[1]+1 +t_pulse <- (res_pulse$result$Time[seq(id_start, id_end)] - res_pulse$result$Time[id_start]) +ace_conc_pulse <- res_pulse$result$Ace_out[seq(id_start, id_end)] - res_pulse$result$Ace_out[id_start] +glc_conc_pulse <- res_pulse$result$Glc[seq(id_start, id_end)] - res_pulse$result$Glc[id_start] +# set biomass concentration at which the pulse is performed +setGlobalQuantities(key = "_X_conc_pulse", initial_value = 1000) +res_nopulse <- runTimeCourse(intervals=10000) +# get simulation results +ace_conc_nopulse <- res_nopulse$result$Ace_out[seq(id_start, id_end)] - res_nopulse$result$Ace_out[id_start] +glc_conc_nopulse <- res_nopulse$result$Glc[seq(id_start, id_end)] - res_nopulse$result$Glc[id_start] +# get simulation times +time_sim <- res_nopulse$result$Time[seq(id_start, id_end)] - res_nopulse$result$Time[id_start] +idx_sim <- get_index_closest(time_meas, time_sim) +# calculate variance weighted SSR +SSR_model_4 <- sum((glc-glc_conc_pulse[idx_sim])**2/sd_glc**2, (ace-ace_conc_pulse[idx_sim])**2/sd_ace**2, (glc_nopulse-glc_conc_nopulse[idx_sim])**2/sd_glc_nopulse**2, (ace_nopulse-ace_conc_nopulse[idx_sim])**2/sd_ace_nopulse**2 ) + +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_1.cps") +# set biomass concentration at which the pulse is performed +setGlobalQuantities(key = "_X_conc_pulse", initial_value = 0.9) +# simulate response to acetate pulse +res_pulse <- runTimeCourse(intervals=10000) +# get simulation results +id_start <- which(res_pulse$result$Ace_out >= 30)[1]+1 +id_end <- which(res_pulse$result$Time >= (res_pulse$result$Time[id_start]+8/60))[1]+1 +t_pulse <- (res_pulse$result$Time[seq(id_start, id_end)] - res_pulse$result$Time[id_start]) +ace_conc_pulse <- res_pulse$result$Ace_out[seq(id_start, id_end)] - res_pulse$result$Ace_out[id_start] +glc_conc_pulse <- res_pulse$result$Glc[seq(id_start, id_end)] - res_pulse$result$Glc[id_start] +# set biomass concentration at which the pulse is performed +setGlobalQuantities(key = "_X_conc_pulse", initial_value = 1000) +res_nopulse <- runTimeCourse(intervals=10000) +# get simulation results +ace_conc_nopulse <- res_nopulse$result$Ace_out[seq(id_start, id_end)] - res_nopulse$result$Ace_out[id_start] +glc_conc_nopulse <- res_nopulse$result$Glc[seq(id_start, id_end)] - res_nopulse$result$Glc[id_start] +# get simulation times +time_sim <- res_nopulse$result$Time[seq(id_start, id_end)] - res_nopulse$result$Time[id_start] +idx_sim <- get_index_closest(time_meas, time_sim) +# calculate variance weighted SSR +SSR_model_1 <- sum((glc-glc_conc_pulse[idx_sim])**2/sd_glc**2, (ace-ace_conc_pulse[idx_sim])**2/sd_ace**2, (glc_nopulse-glc_conc_nopulse[idx_sim])**2/sd_glc_nopulse**2, (ace_nopulse-ace_conc_nopulse[idx_sim])**2/sd_ace_nopulse**2 ) + +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_2.cps") +# set biomass concentration at which the pulse is performed +setGlobalQuantities(key = "_X_conc_pulse", initial_value = 0.9) +# simulate response to acetate pulse +res_pulse <- runTimeCourse(intervals=10000) +# get simulation results +id_start <- which(res_pulse$result$Ace_out >= 30)[1]+1 +id_end <- which(res_pulse$result$Time >= (res_pulse$result$Time[id_start]+8/60))[1]+1 +t_pulse <- (res_pulse$result$Time[seq(id_start, id_end)] - res_pulse$result$Time[id_start]) +ace_conc_pulse <- res_pulse$result$Ace_out[seq(id_start, id_end)] - res_pulse$result$Ace_out[id_start] +glc_conc_pulse <- res_pulse$result$Glc[seq(id_start, id_end)] - res_pulse$result$Glc[id_start] +# set biomass concentration at which the pulse is performed +setGlobalQuantities(key = "_X_conc_pulse", initial_value = 1000) +res_nopulse <- runTimeCourse(intervals=10000) +# get simulation results +ace_conc_nopulse <- res_nopulse$result$Ace_out[seq(id_start, id_end)] - res_nopulse$result$Ace_out[id_start] +glc_conc_nopulse <- res_nopulse$result$Glc[seq(id_start, id_end)] - res_nopulse$result$Glc[id_start] +# get simulation times +time_sim <- res_nopulse$result$Time[seq(id_start, id_end)] - res_nopulse$result$Time[id_start] +idx_sim <- get_index_closest(time_meas, time_sim) +# calculate variance weighted SSR +SSR_model_2 <- sum((glc-glc_conc_pulse[idx_sim])**2/sd_glc**2, (ace-ace_conc_pulse[idx_sim])**2/sd_ace**2, (glc_nopulse-glc_conc_nopulse[idx_sim])**2/sd_glc_nopulse**2, (ace_nopulse-ace_conc_nopulse[idx_sim])**2/sd_ace_nopulse**2 ) + +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_3.cps") +# set biomass concentration at which the pulse is performed +setGlobalQuantities(key = "_X_conc_pulse", initial_value = 0.9) +# simulate response to acetate pulse +res_pulse <- runTimeCourse(intervals=10000) +# get simulation results +id_start <- which(res_pulse$result$Ace_out >= 30)[1]+1 +id_end <- which(res_pulse$result$Time >= (res_pulse$result$Time[id_start]+8/60))[1]+1 +t_pulse <- (res_pulse$result$Time[seq(id_start, id_end)] - res_pulse$result$Time[id_start]) +ace_conc_pulse <- res_pulse$result$Ace_out[seq(id_start, id_end)] - res_pulse$result$Ace_out[id_start] +glc_conc_pulse <- res_pulse$result$Glc[seq(id_start, id_end)] - res_pulse$result$Glc[id_start] +# set biomass concentration at which the pulse is performed +setGlobalQuantities(key = "_X_conc_pulse", initial_value = 1000) +res_nopulse <- runTimeCourse(intervals=10000) +# get simulation results +ace_conc_nopulse <- res_nopulse$result$Ace_out[seq(id_start, id_end)] - res_nopulse$result$Ace_out[id_start] +glc_conc_nopulse <- res_nopulse$result$Glc[seq(id_start, id_end)] - res_nopulse$result$Glc[id_start] +# get simulation times +time_sim <- res_nopulse$result$Time[seq(id_start, id_end)] - res_nopulse$result$Time[id_start] +idx_sim <- get_index_closest(time_meas, time_sim) +# calculate variance weighted SSR +SSR_model_3 <- sum((glc-glc_conc_pulse[idx_sim])**2/sd_glc**2, (ace-ace_conc_pulse[idx_sim])**2/sd_ace**2, (glc_nopulse-glc_conc_nopulse[idx_sim])**2/sd_glc_nopulse**2, (ace_nopulse-ace_conc_nopulse[idx_sim])**2/sd_ace_nopulse**2 ) + + +# Steady-state fluxes for acetate concentration between 0.1 and 100mM (Enjalbert et al., 2017; Pinhal et al., 2019) +############################################################################################## + +growth_rates <- c(0.521128511, 0.611148842625582, 0.613161998174498, 0.502533817, 0.496290415, 0.488201506, 0.547635665, 0.499830448, 0.474554197, 0.425356578, 0.377534684, 0.645724326, 0.618475601, 0.554887936, 0.564811523, 0.527571192, 0.434972836, 0.3824734, 0.583623355, 0.620905534, 0.564259247, 0.532148135, 0.483885215, 0.557074418, 0.630654409249223) +sd_growth_rates <- c(0.001793104, 0.00204807928657914, 0.00219396182484705, 0.001709207, 0.001846205, 0.001757403, 0.001821375, 0.001025702, 0.001940912, 0.001204707, 0.001999188, 0.001418374, 0.001932601, 0.001455791, 0.001574234, 0.001206265, 0.001292476, 0.001068259, 0.001804648, 0.001643459, 0.001598405, 0.001121218, 0.000912408, 0.00194896, 0.00203369597574686) +sd_growth_rates <- sapply(sd_growth_rates, FUN=function(x) max(x,0.01)) +glc_upt <- c(8.654860167, 8.36127542981722, 7.98010111285252, 9.236935826, 8.274418986, 7.560431219, 7.339194455, 5.775312502, 6.423391263, 5.1544758, 3.938631573, 8.115447647, 9.28067031, 6.737153424, 7.172748804, 5.884186033, 5.684201497, 4.811576974, 9.632702365, 8.055042777, 9.708342814, 7.100081588, 5.505759496, 9.242859752, 8.18621623190759) +sd_glc_upt <- c(0.337812425, 0.38531328268303, 0.373770045721031, 0.356787032, 0.334672954, 0.317509322, 0.288025925, 0.16053276, 0.375934255, 0.293148172, 0.359225607, 0.197331684, 0.360984112, 0.229372278, 0.241820396, 0.20450532, 0.260869273, 0.216134352, 0.34289286, 0.350305744, 0.293144783, 0.220135755, 0.153471508, 0.25245346, 0.396815184905029) +sd_glc_upt <- sapply(sd_glc_upt, FUN=function(x) max(x,0.1)) +ace_flx <- c(3.5, -2.7, 1.516999356, 1.26845123082679, 0.775821380507016, 0.678877137, 0.017366464, -0.991478151, -1.286687213, -2.078474994, -1.530841439, -1.525342269, -1.253581266, 1.984679487, 0.546462624, -0.136780389, -0.393883917, -0.610240984, -1.120767885, -1.277455315, 2.574285211, 2.051935093, 1.828415596, -1.262442483, -1.317987733, 2.333568565, 1.85234639824858) +sd_ace_flx <- c(0.35, 0.27, 0.316118066, 0.388752117258161, 0.40715278851436, 0.33357638, 0.37333751, 0.347029894, 0.280501612, 0.195031303, 0.376252463, 0.226182385, 0.303661317, 0.253610517, 0.385450715, 0.243880325, 0.30665695, 0.257983739, 0.23844407, 0.198458448, 0.299832036, 0.334956504, 0.317134334, 0.263807154, 0.195219648, 0.016120887, 0.386174129654754) +sd_ace_flx <- sapply(sd_ace_flx, FUN=function(x) max(x,0.1)) + +ace_conc <- c(0.451964624, 1.11600286648471, 2.04718732931708, 2.86252145, 5.285907977, 9.106164204, 16.67476528, 16.76626787, 17.00218707, 30.8667961, 57.92292091, 0.478352574, 4.55673229, 8.10163028, 8.22100734, 7.829591756, 33.53244905, 66.20361403, 0.436956014, 1.28468189, 1.555875222, 12.02564968, 30.24566673, 0.548011282, 2.29562227069566) + +fluxes <- c("Values[v_growth_rate]", "Values[v_glc_uptake]", "Values[v_ace_net]") + +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_1.cps") +# remove events and fix concentrations of actate, glucose and biomass +deleteEvent(getEvents()$key) +setSpecies(key="Ace_out", type="fixed") +setSpecies(key="Glc", type="fixed") +setSpecies(key="X", type="fixed") +Enj_SS_SSR_model_1 <- c() +for (j in seq(length(ace_conc))){ + setSpecies(key="Ace_out", initial_concentration=ace_conc[j]) + applyInitialState() + res_ss <- runSteadyState() + res_flx <- unlist(res_ss$global_quantities[res_ss$global_quantities$key %in% fluxes, "value"]) + Enj_SS_SSR_model_1 <- c(Enj_SS_SSR_model_1, (res_flx[1]-growth_rates[j])**2/sd_growth_rates[j]**2 + (res_flx[2]-glc_upt[j])**2/sd_glc_upt[j]**2 + (res_flx[3]-ace_flx[j])**2/sd_ace_flx[j]**2) +} +Enj_SS_SSR_model_1 <- sum(Enj_SS_SSR_model_1) + +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_2.cps") +# remove events and fix concentrations of actate, glucose and biomass +deleteEvent(getEvents()$key) +setSpecies(key="Ace_out", type="fixed") +setSpecies(key="Glc", type="fixed") +setSpecies(key="X", type="fixed") +Enj_SS_SSR_model_2 <- c() +for (j in seq(length(ace_conc))){ + setSpecies(key="Ace_out", initial_concentration=ace_conc[j]) + applyInitialState() + res_ss <- runSteadyState() + res_flx <- unlist(res_ss$global_quantities[res_ss$global_quantities$key %in% fluxes, "value"]) + Enj_SS_SSR_model_2 <- c(Enj_SS_SSR_model_2, (res_flx[1]-growth_rates[j])**2/sd_growth_rates[j]**2 + (res_flx[2]-glc_upt[j])**2/sd_glc_upt[j]**2 + (res_flx[3]-ace_flx[j])**2/sd_ace_flx[j]**2) +} +Enj_SS_SSR_model_2 <- sum(Enj_SS_SSR_model_2) + +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_3.cps") +# remove events and fix concentrations of actate, glucose and biomass +deleteEvent(getEvents()$key) +setSpecies(key="Ace_out", type="fixed") +setSpecies(key="Glc", type="fixed") +setSpecies(key="X", type="fixed") +Enj_SS_SSR_model_3 <- c() +for (j in seq(length(ace_conc))){ + setSpecies(key="Ace_out", initial_concentration=ace_conc[j]) + applyInitialState() + res_ss <- runSteadyState() + res_flx <- unlist(res_ss$global_quantities[res_ss$global_quantities$key %in% fluxes, "value"]) + Enj_SS_SSR_model_3 <- c(Enj_SS_SSR_model_3, (res_flx[1]-growth_rates[j])**2/sd_growth_rates[j]**2 + (res_flx[2]-glc_upt[j])**2/sd_glc_upt[j]**2 + (res_flx[3]-ace_flx[j])**2/sd_ace_flx[j]**2) +} +Enj_SS_SSR_model_3 <- sum(Enj_SS_SSR_model_3) + +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_4.cps") +# remove events and fix concentrations of actate, glucose and biomass +deleteEvent(getEvents()$key) +setSpecies(key="Ace_out", type="fixed") +setSpecies(key="Glc", type="fixed") +setSpecies(key="X", type="fixed") +Enj_SS_SSR_model_4 <- c() +for (j in seq(length(ace_conc))){ + setSpecies(key="Ace_out", initial_concentration=ace_conc[j]) + applyInitialState() + res_ss <- runSteadyState() + res_flx <- unlist(res_ss$global_quantities[res_ss$global_quantities$key %in% fluxes, "value"]) + Enj_SS_SSR_model_4 <- c(Enj_SS_SSR_model_4, (res_flx[1]-growth_rates[j])**2/sd_growth_rates[j]**2 + (res_flx[2]-glc_upt[j])**2/sd_glc_upt[j]**2 + (res_flx[3]-ace_flx[j])**2/sd_ace_flx[j]**2) +} +Enj_SS_SSR_model_4 <- sum(Enj_SS_SSR_model_4) + + +# Chemostat experiment (Renilla et al., 2012) +############################################# + +fluxes <- c("Values[v_growth_rate]", "Values[v_glc_uptake]", "Values[v_ace_net]") + +# acetate flux as function of dilution rate +dilution_rate <- c(0.09586056644880175, 0.20043572984749464, 0.2997821350762528, 0.3468409586056645, 0.39912854030501094, 0.44618736383442276, 0.5002178649237474) +# acetate flux as function of glc uptake +q_glc <- c(1.43478260869565, 2.7391304347826, 4.59130434782608, 4.69565217391304, 5.7391304347826, 5.92173913043478, 6.10434782608695, 7.12173913043477, 8.34782608695652) +q_ace <- c(0.381358340437624, 0.762906128635029, 1.40509614473808, 1.17836506583309, 1.99602159704461, 2.08290233967983, 2.01292033721701, 2.44717249218528, 2.5848252344416) +q_ace_sd <- c(0.139547219854126, 0.139471440750213, 0.15697641375391, 0.1045751633987, 0.19179691200152, 0.0174291938997801, 0.24397082504499, 0.52295159609737, 0.33130624230368) +q_glc_sd <- c(0.20234760540444, 0.58899555210581, 0.36827270271402, 0.14607356337583, 0.630432119199, 0.13069807968161, 0.47630364392495, 0.99361847036576, 0.66185122193432) +q_ace_sd <- sapply(q_ace_sd, FUN=function(x) max(x,0.1)) +q_glc_sd <- sapply(q_glc_sd, FUN=function(x) max(x,0.1)) + + +setwd(model_dir) +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_1.cps") +# delete events and set glc feed to 1 +deleteEvent(getEvents()$key) +setGlobalQuantities(key = "_feed", initial_value = 1) +Ren_SSR_model_1 <- c() +for (j in seq(length(dilution_rate))){ + setGlobalQuantities(key = "_dilution_rate", initial_value = dilution_rate[j], model=getCurrentModel()) + applyInitialState(model=getCurrentModel()) + res_ss <- runSteadyState(model=getCurrentModel()) + val <- unlist(res_ss$global_quantities[res_ss$global_quantities$key %in% fluxes, "value"]) + Ren_SSR_model_1 <- c(Ren_SSR_model_1, + (val[2]-q_glc[j])**2/q_glc_sd[j]**2 + (val[3]-q_ace[j])**2/q_ace_sd[j]**2) +} +Ren_SSR_model_1 <- sum(Ren_SSR_model_1) + +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_2.cps") +# delete events and set glc feed to 1 +deleteEvent(getEvents()$key) +setGlobalQuantities(key = "_feed", initial_value = 1) +Ren_SSR_model_2 <- c() +for (j in seq(length(dilution_rate))){ + setGlobalQuantities(key = "_dilution_rate", initial_value = dilution_rate[j]) + applyInitialState() + res_ss <- runSteadyState() + val <- unlist(res_ss$global_quantities[res_ss$global_quantities$key %in% fluxes, "value"]) + Ren_SSR_model_2 <- c(Ren_SSR_model_2, + (val[2]-q_glc[j])**2/q_glc_sd[j]**2 + (val[3]-q_ace[j])**2/q_ace_sd[j]**2) +} +Ren_SSR_model_2 <- sum(Ren_SSR_model_2) + +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_3.cps") +# delete events and set glc feed to 1 +deleteEvent(getEvents()$key) +setGlobalQuantities(key = "_feed", initial_value = 1) +Ren_SSR_model_3 <- c() +for (j in seq(length(dilution_rate))){ + setGlobalQuantities(key = "_dilution_rate", initial_value = dilution_rate[j], model=getCurrentModel()) + applyInitialState(model=getCurrentModel()) + res_ss <- runSteadyState(model=getCurrentModel()) + val <- unlist(res_ss$global_quantities[res_ss$global_quantities$key %in% fluxes, "value"]) + Ren_SSR_model_3 <- c(Ren_SSR_model_3, + (val[2]-q_glc[j])**2/q_glc_sd[j]**2 + (val[3]-q_ace[j])**2/q_ace_sd[j]**2) +} +Ren_SSR_model_3 <- sum(Ren_SSR_model_3) + + +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_4.cps") +# delete events and set glc feed to 1 +deleteEvent(getEvents()$key) +setGlobalQuantities(key = "_feed", initial_value = 1) +Ren_SSR_model_4 <- c() +for (j in seq(length(dilution_rate))){ + setGlobalQuantities(key = "_dilution_rate", initial_value = dilution_rate[j], model=getCurrentModel()) + applyInitialState(model=getCurrentModel()) + res_ss <- runSteadyState(model=getCurrentModel()) + val <- unlist(res_ss$global_quantities[res_ss$global_quantities$key %in% fluxes, "value"]) + Ren_SSR_model_4 <- c(Ren_SSR_model_4, + (val[2]-q_glc[j])**2/q_glc_sd[j]**2 + (val[3]-q_ace[j])**2/q_ace_sd[j]**2) +} +Ren_SSR_model_4 <- sum(Ren_SSR_model_4) + + +# Plot results +setwd(results_dir) +pdf(file="Figure 3-figure supplement 4.pdf", width = 5, height = 5) + +data <- c(SSR_model_1, SSR_model_2, SSR_model_3, SSR_model_4) +xx <- barplot(data, las=1, ylab="variance-weighted SSR", log="y", names.arg=c("model_1", "model_2", "model_3", "model_4")) +text(x = xx, y = round(data), label = round(data), pos = c(1,3,1,3)) + +data <- c(Enj_SS_SSR_model_1, Enj_SS_SSR_model_2, Enj_SS_SSR_model_3, Enj_SS_SSR_model_4) +xx <- barplot(data, las=1, ylab="variance-weighted SSR (steady-state data from Enjalbert 2017)", log="y", names.arg=c("model_1", "model_2", "model_3", "model_4")) +text(x = xx, y = round(data), label = round(data), pos = c(1,3,1,3)) + +data <- c(Ren_SSR_model_1, Ren_SSR_model_2, Ren_SSR_model_3, Ren_SSR_model_4) +xx <- barplot(data, las=1, ylab="variance-weighted SSR", names.arg=c("model_1", "model_2", "model_3", "model_4")) +text(x = xx, y = round(data), label = round(data), pos = 1) + +dev.off() diff --git a/3-Metabolic_control_analyses.R b/3-Metabolic_control_analyses.R new file mode 100644 index 0000000..b26a03b --- /dev/null +++ b/3-Metabolic_control_analyses.R @@ -0,0 +1,379 @@ + + +################### +# Set environment # +################### + +# load libraries and initialize environment +source("set_env.R") + + +############################################ +# Load global sensitivity analysis results # +############################################ + +# This file is generated by script "1-Model_construction.R". +setwd(results_dir) +load("mc_results_100.RData") + + +######################################## +# Metabolic control analyses - Model 4 # +######################################## + +setwd(model_dir) +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model.cps") +setwd(results_dir) + +# delete events and fix concentrations of biomass and extracellular glc and acetate +deleteEvent(getEvents()$key) +setSpecies(key="Ace_out", type="fixed") +setSpecies(key="Glc", type="fixed") +setSpecies(key="X", type="fixed") + +# calculate flux control coefficients +n_step <- 50 +ace_range <- 10**(seq(-1,2,length.out = n_step)) +res_ace_MCA <- matrix(NA, nrow=length(ace_range), ncol=7) +res_ace_MCA[,1] <- ace_range +for (i in seq(n_step)){ + setSpecies(key="Ace_out{cell}", initial_concentration = ace_range[i], model = getCurrentModel()) + applyInitialState(model = getCurrentModel()) + res_MCA <- runMetabolicControlAnalysis(perform_steady_state_analysis=TRUE, model = getCurrentModel()) + if (res_MCA$result_ss != "found"){ + print("error") + }else{ + fcc <- res_MCA$flux_control_coefficients_scaled + res_ace_MCA[i,seq(2,7)] <- fcc["(ackA)", colnames(fcc) != "'Summation Error'"] + } +} +colnames(res_ace_MCA) <- c("ace_conc", colnames(fcc)[colnames(fcc) != "'Summation Error'"]) + +# plot results +pdf(file = "Figure 4.pdf", height = 4, width = 5) + +setSpecies(key="Ace_out{cell}", initial_concentration = 0.1, model = getCurrentModel()) +applyInitialState(model = getCurrentModel()) +res_MCA <- runMetabolicControlAnalysis(model = getCurrentModel()) + +hmc <- res_MCA$flux_control_coefficients_scaled[c("(glc_upt)", "(ackA)", "(sink)"), c("(glc_upt)", "(sink)", "(pta)", "(ackA)", "(ace_xch)")] +colnames(hmc) <- c("glycolysis", "TCA", "Pta", "AckA", "Ace_xch") +rownames(hmc) <- c("glycolysis", "acetate", "TCA") +rgb.palette <- colorRampPalette(brewer.pal(n = 11, name = "RdBu")) +breaks <- c(-1, -0.5, -0.2, -0.1, -0.05, -0.01, 0.01, 0.05, 0.1, 0.2, 0.5, 1, 1.5) +col_scale <- rgb.palette(length(breaks)-1) +col_scale[6] <- "#DDDDDD" +heatmap.2(hmc, Rowv=FALSE, Colv=FALSE, dendrogram="none", col=col_scale, breaks=breaks, + scale="none", trace="none", colsep=seq(0,ncol(hmc)), rowsep=seq(0,nrow(hmc)), sepwidth=c(0.04,0.02)) + +dev.off() + +# sensitivity analysis +run_mca_sa <- function(){ + n_step <- 50 + ace_range <- 10**(seq(-1,2,length.out = n_step)) + res_ace_MCA <- matrix(NA, nrow=length(ace_range), ncol=7) + res_ace_MCA[,1] <- ace_range + for (i in seq(n_step)){ + setSpecies(key="Ace_out{cell}", initial_concentration = ace_range[i], model = getCurrentModel()) + applyInitialState(model = getCurrentModel()) + res_MCA <- runMetabolicControlAnalysis(perform_steady_state_analysis=TRUE, model = getCurrentModel()) + if (res_MCA$result_ss != "found"){ + print("error") + }else{ + fcc <- res_MCA$flux_control_coefficients_scaled + res_ace_MCA[i,seq(2,7)] <- fcc["(ackA)", colnames(fcc) != "'Summation Error'"] + } + } + colnames(res_ace_MCA) <- c("ace_conc", colnames(fcc)[colnames(fcc) != "'Summation Error'"]) + return(res_ace_MCA) +} + +# create progress bar +pb <- txtProgressBar(min=0, max=ncol(fit_results$res_par)-1, style=3) + +# sensitivity analysis +res <- array(NA, dim=c(ncol(fit_results$res_par)-1, length(ace_range), 7), dimnames=list(iter=NULL, ace_conc=NULL, fcc=c("ace_conc", colnames(fcc)[colnames(fcc) != "'Summation Error'"]))) +for (i in seq(ncol(fit_results$res_par)-1)){ + rp <- c(fit_results$res_par[,i+1]) + names(rp) <- fit_results$res_par[,"parameter"] + model <- update_params(getCurrentModel(), rp) + res[i,,] <- run_mca_sa() + setTxtProgressBar(pb, i) +} + +# close progress bar +close(pb) + +# plot results +pdf(file="Figure 5.pdf", width = 7, height = 9) +par(mfrow=c(4,3)) + +xlab_main <- c(0.1, 1, 10, 100) +xlab_sec <- c(seq(0.2, 0.9, by=0.1), seq(2, 9, by=1), seq(20, 90, by=10)) +conc_threshold <- 13.27174 + +ctrl_by_acetate_pathway <- apply(apply(res[,,c("(ackA)", "(pta)", "(ace_xch)")], c(1,2), sum), 2, mean) +ctrl_by_glc_upt <- apply(res[,,"(glc_upt)"], 2, mean) +ctrl_by_TCA <- apply(res[,,"(sink)"], 2, mean) + +plot(x=ace_range, y=ctrl_by_acetate_pathway, type="l", xaxt="n", las=1, ylim=c(0,1), log="x", xaxs="i", yaxs="i", xlab="[acetate] (mM)", ylab="control by acetate pathway", lwd=1.2) +polygon(x=c(ace_range, rev(ace_range)), + y=c(apply(apply(res[,,c("(ackA)", "(pta)", "(ace_xch)")], c(1,2), sum), 2, max), rev(apply(apply(res[,,c("(ackA)", "(pta)", "(ace_xch)")], c(1,2), sum), 2, min))), + col="#00000055", border=NA) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) + +fconc <- 0.8 +lines_threshold(x=ace_range, y=ctrl_by_glc_upt, threshold=conc_threshold, new=TRUE, xaxt="n", xaxs="i", yaxs="i", las=1, xlim=c(0.1,100), type="l", ylim=c(-3,3), log="x", xlab="[acetate] (mM)", ylab="control by glycolysis", lwd=1.2) +polygon(x=c(ace_range[ace_range < conc_threshold*fconc], rev(ace_range[ace_range < conc_threshold*fconc])), + y=c(apply(res[,ace_range < conc_threshold*fconc,"(glc_upt)"], 2, max), rev(apply(res[,ace_range < conc_threshold*fconc,"(glc_upt)"], 2, min))), + col="#00000055", border=NA) + +polygon(x=c(ace_range[ace_range > conc_threshold/fconc], rev(ace_range[ace_range > conc_threshold/fconc])), + y=c(apply(res[,ace_range > conc_threshold/fconc,"(glc_upt)"], 2, max), rev(apply(res[,ace_range > conc_threshold/fconc,"(glc_upt)"], 2, min))), + col="#00000055", border=NA) + +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +abline(h=0) + +lines_threshold(x=ace_range, y=ctrl_by_TCA, threshold=conc_threshold, new=TRUE, xaxt="n", xaxs="i", yaxs="i", las=1, xlim=c(0.1,100), type="l", ylim=c(-3,3), log="x", xlab="[acetate] (mM)", ylab="control by TCA", lwd=1.2) +polygon(x=c(ace_range[ace_range < conc_threshold*fconc], rev(ace_range[ace_range < conc_threshold*fconc])), + y=c(apply(res[,ace_range < conc_threshold*fconc,"(sink)"], 2, max), rev(apply(res[,ace_range < conc_threshold*fconc,"(sink)"], 2, min))), + col="#00000055", border=NA) + +polygon(x=c(ace_range[ace_range > conc_threshold/fconc], rev(ace_range[ace_range > conc_threshold/fconc])), + y=c(apply(res[,ace_range > conc_threshold/fconc,"(sink)"], 2, max), rev(apply(res[,ace_range > conc_threshold/fconc,"(sink)"], 2, min))), + col="#00000055", border=NA) + +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +abline(h=0) + +plot(x=ace_range, y=ctrl_by_TCA+ctrl_by_glc_upt, type="l", ylim=c(0,1), xaxt="n", xaxs="i", yaxs="i", las=1, log="x", xlab="[acetate] (mM)", ylab="control by TCA+glycolysis", lwd=1.2) +polygon(x=c(ace_range, rev(ace_range)), + y=c(apply(res[,,"(sink)"]+res[,,"(glc_upt)"], 2, max), rev(apply(res[,,"(sink)"]+res[,,"(glc_upt)"], 2, min))), + col="#00000055", border=NA) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) + +plot(x=ace_range, y=-ctrl_by_TCA/ctrl_by_glc_upt, type="l", ylim=c(0,2), xaxt="n", xaxs="i", yaxs="i", las=1, log="x", xlab="[acetate] (mM)", ylab="control by TCA/glycolysis", lwd=1.2) +polygon(x=c(ace_range, rev(ace_range)), + y=c(apply(-res[,,"(sink)"]/res[,,"(glc_upt)"], 2, max), rev(apply(-res[,,"(sink)"]/res[,,"(glc_upt)"], 2, min))), + col="#00000055", border=NA) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +abline(h=1) + +dev.off() + + +##################################################### +# Metabolic control analyses for alternative models # +##################################################### + +########### +# Model 1 # +########### + +setwd(model_dir) +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_1.cps") +setwd(results_dir) + +# delete events and fix concentrations of biomass and extracellular glc and acetate +deleteEvent(getEvents()$key) +setSpecies(key="Ace_out", type="fixed") +setSpecies(key="Glc", type="fixed") +setSpecies(key="X", type="fixed") + +# calculate flux control coefficients +n_step <- 50 +ace_range <- 10**(seq(-1,2,length.out = n_step)) +res_ace_MCA <- matrix(NA, nrow=length(ace_range), ncol=7) +res_ace_MCA[,1] <- ace_range +for (i in seq(n_step)){ + setSpecies(key="Ace_out{cell}", initial_concentration = ace_range[i], model = getCurrentModel()) + applyInitialState(model = getCurrentModel()) + res_MCA <- runMetabolicControlAnalysis(model = getCurrentModel()) + if (res_MCA$result_ss != "found"){ + print("error") + }else{ + fcc <- res_MCA$flux_control_coefficients_scaled + res_ace_MCA[i,seq(2,7)] <- fcc["(ackA)", colnames(fcc) != "'Summation Error'"] + } +} +colnames(res_ace_MCA) <- c("ace_conc", colnames(fcc)[colnames(fcc) != "'Summation Error'"]) + +pdf(file="Figure 5-figure supplement 1.pdf", width = 7, height = 9) +par(mfrow=c(4,3)) + +xlab_main <- c(0.1, 1, 10, 100) +xlab_sec <- c(seq(0.2, 0.9, by=0.1), seq(2, 9, by=1), seq(20, 90, by=10)) +conc_threshold <- 16 + +ctrl_by_acetate_pathway <- apply(res_ace_MCA[,c("(ackA)", "(pta)", "(ace_xch)")], 1, sum) +ctrl_by_glc_upt <- res_ace_MCA[,"(glc_upt)"] +ctrl_by_TCA <- res_ace_MCA[,"(sink)"] + +plot(x=ace_range, y=ctrl_by_acetate_pathway, type="l", xaxt="n", las=1, ylim=c(0,1), log="x", xaxs="i", yaxs="i", xlab="[acetate] (mM)", ylab="control by acetate pathway", lwd=2) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) + +lines_threshold(x=ace_range, y=ctrl_by_glc_upt, threshold=conc_threshold, new=TRUE, xaxt="n", xaxs="i", yaxs="i", las=1, xlim=c(0.1,100), type="l", ylim=c(-3,3), log="x", xlab="[acetate] (mM)", ylab="control by glycolysis", lwd=2) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +abline(h=0) + +lines_threshold(x=ace_range, y=ctrl_by_TCA, threshold=conc_threshold, new=TRUE, xaxt="n", xaxs="i", yaxs="i", las=1, xlim=c(0.1,100), type="l", ylim=c(-3,3), log="x", xlab="[acetate] (mM)", ylab="control by TCA", lwd=2) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +abline(h=0) + +plot(x=ace_range, y=ctrl_by_TCA+ctrl_by_glc_upt, type="l", ylim=c(0,1), xaxt="n", xaxs="i", yaxs="i", las=1, log="x", xlab="[acetate] (mM)", ylab="control by TCA+glycolysis", lwd=2) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) + +plot(x=ace_range, y=-ctrl_by_TCA/ctrl_by_glc_upt, type="l", ylim=c(0,2), xaxt="n", xaxs="i", yaxs="i", las=1, log="x", xlab="[acetate] (mM)", ylab="control by TCA/glycolysis", lwd=2) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +abline(h=1) + +dev.off() + +########### +# Model 2 # +########### + +setwd(model_dir) +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_2.cps") +setwd(results_dir) + +# delete events and fix concentrations of biomass and extracellular glc and acetate +deleteEvent(getEvents()$key) +setSpecies(key="Ace_out", type="fixed") +setSpecies(key="Glc", type="fixed") +setSpecies(key="X", type="fixed") + +# calculate flux control coefficients +n_step <- 50 +ace_range <- 10**(seq(-1,2,length.out = n_step)) +res_ace_MCA <- matrix(NA, nrow=length(ace_range), ncol=7) +res_ace_MCA[,1] <- ace_range +for (i in seq(n_step)){ + setSpecies(key="Ace_out{cell}", initial_concentration = ace_range[i], model = getCurrentModel()) + applyInitialState(model = getCurrentModel()) + res_MCA <- runMetabolicControlAnalysis(model = getCurrentModel()) + if (res_MCA$result_ss != "found"){ + print("error") + }else{ + fcc <- res_MCA$flux_control_coefficients_scaled + res_ace_MCA[i,seq(2,7)] <- fcc["(ackA)", colnames(fcc) != "'Summation Error'"] + } +} +colnames(res_ace_MCA) <- c("ace_conc", colnames(fcc)[colnames(fcc) != "'Summation Error'"]) + +pdf(file="Figure 5-figure supplement 2.pdf", width = 7, height = 9) +par(mfrow=c(4,3)) + +xlab_main <- c(0.1, 1, 10, 100) +xlab_sec <- c(seq(0.2, 0.9, by=0.1), seq(2, 9, by=1), seq(20, 90, by=10)) +conc_threshold <- 14 + +ctrl_by_acetate_pathway <- apply(res_ace_MCA[,c("(ackA)", "(pta)", "(ace_xch)")], 1, sum) +ctrl_by_glc_upt <- res_ace_MCA[,"(glc_upt)"] +ctrl_by_TCA <- res_ace_MCA[,"(sink)"] + +plot(x=ace_range, y=ctrl_by_acetate_pathway, type="l", xaxt="n", las=1, ylim=c(0,1), log="x", xaxs="i", yaxs="i", xlab="[acetate] (mM)", ylab="control by acetate pathway", lwd=2) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) + +lines_threshold(x=ace_range, y=ctrl_by_glc_upt, threshold=conc_threshold, new=TRUE, xaxt="n", xaxs="i", yaxs="i", las=1, xlim=c(0.1,100), type="l", ylim=c(-3,3), log="x", xlab="[acetate] (mM)", ylab="control by glycolysis", lwd=2) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +abline(h=0) + +lines_threshold(x=ace_range, y=ctrl_by_TCA, threshold=conc_threshold, new=TRUE, xaxt="n", xaxs="i", yaxs="i", las=1, xlim=c(0.1,100), type="l", ylim=c(-3,3), log="x", xlab="[acetate] (mM)", ylab="control by TCA", lwd=2) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +abline(h=0) + +plot(x=ace_range, y=ctrl_by_TCA+ctrl_by_glc_upt, type="l", ylim=c(0,1), xaxt="n", xaxs="i", yaxs="i", las=1, log="x", xlab="[acetate] (mM)", ylab="control by TCA+glycolysis", lwd=2) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) + +plot(x=ace_range, y=-ctrl_by_TCA/ctrl_by_glc_upt, type="l", ylim=c(0,2), xaxt="n", xaxs="i", yaxs="i", las=1, log="x", xlab="[acetate] (mM)", ylab="control by TCA/glycolysis", lwd=2) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +abline(h=1) + +dev.off() + +########### +# Model 3 # +########### + +setwd(model_dir) +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model_3.cps") +setwd(results_dir) + +# delete events and fix concentrations of biomass and extracellular glc and acetate +deleteEvent(getEvents()$key) +setSpecies(key="Ace_out", type="fixed") +setSpecies(key="Glc", type="fixed") +setSpecies(key="X", type="fixed") + +# calculate flux control coefficients +n_step <- 50 +ace_range <- 10**(seq(-1,2,length.out = n_step)) +res_ace_MCA <- matrix(NA, nrow=length(ace_range), ncol=7) +res_ace_MCA[,1] <- ace_range +for (i in seq(n_step)){ + setSpecies(key="Ace_out{cell}", initial_concentration = ace_range[i], model = getCurrentModel()) + applyInitialState(model = getCurrentModel()) + res_MCA <- runMetabolicControlAnalysis(model = getCurrentModel()) + if (res_MCA$result_ss != "found"){ + print("error") + }else{ + fcc <- res_MCA$flux_control_coefficients_scaled + res_ace_MCA[i,seq(2,7)] <- fcc["(ackA)", colnames(fcc) != "'Summation Error'"] + } +} +colnames(res_ace_MCA) <- c("ace_conc", colnames(fcc)[colnames(fcc) != "'Summation Error'"]) + +pdf(file="Figure 5-figure supplement 3.pdf", width = 7, height = 9) +par(mfrow=c(4,3)) + +xlab_main <- c(0.1, 1, 10, 100) +xlab_sec <- c(seq(0.2, 0.9, by=0.1), seq(2, 9, by=1), seq(20, 90, by=10)) +conc_threshold <- 2.3 + +ctrl_by_acetate_pathway <- apply(res_ace_MCA[,c("(ackA)", "(pta)", "(ace_xch)")], 1, sum) +ctrl_by_glc_upt <- res_ace_MCA[,"(glc_upt)"] +ctrl_by_TCA <- res_ace_MCA[,"(sink)"] + +plot(x=ace_range, y=ctrl_by_acetate_pathway, type="l", xaxt="n", las=1, ylim=c(0,1), log="x", xaxs="i", yaxs="i", xlab="[acetate] (mM)", ylab="control by acetate pathway", lwd=2) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) + +lines_threshold(x=ace_range, y=ctrl_by_glc_upt, threshold=conc_threshold, new=TRUE, xaxt="n", xaxs="i", yaxs="i", las=1, xlim=c(0.1,100), type="l", ylim=c(-3,3), log="x", xlab="[acetate] (mM)", ylab="control by glycolysis", lwd=2) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +abline(h=0) + +lines_threshold(x=ace_range, y=ctrl_by_TCA, threshold=conc_threshold, new=TRUE, xaxt="n", xaxs="i", yaxs="i", las=1, xlim=c(0.1,100), type="l", ylim=c(-3,3), log="x", xlab="[acetate] (mM)", ylab="control by TCA", lwd=2) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +abline(h=0) + +plot(x=ace_range, y=ctrl_by_TCA+ctrl_by_glc_upt, type="l", ylim=c(0,1), xaxt="n", xaxs="i", yaxs="i", las=1, log="x", xlab="[acetate] (mM)", ylab="control by TCA+glycolysis", lwd=2) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) + +plot(x=ace_range, y=-ctrl_by_TCA/ctrl_by_glc_upt, type="l", ylim=c(0,2), xaxt="n", xaxs="i", yaxs="i", las=1, log="x", xlab="[acetate] (mM)", ylab="control by TCA/glycolysis", lwd=2) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +abline(h=1) + +dev.off() + diff --git a/4-Regulation_analyses.R b/4-Regulation_analyses.R new file mode 100644 index 0000000..4555c50 --- /dev/null +++ b/4-Regulation_analyses.R @@ -0,0 +1,186 @@ + + +################### +# Set environment # +################### + +# load libraries and initialize environment +source("set_env.R") + + +############################################ +# Load global sensitivity analysis results # +############################################ + +# This file is generated by script "1-Model_construction.R". +setwd(results_dir) +load("mc_results_100.RData") + + +################################# +# Metabolic regulation analyses # +################################# + +setwd(model_dir) +loadModel("Millard2020_Ecoli_glc_ace_kinetic_model.cps") +setwd(results_dir) +# delete events and fix concentrations of biomass and extracellular glc and acetate +deleteEvent(getEvents()$key) +setSpecies(key="Ace_out", type="fixed") +setSpecies(key="Glc", type="fixed") +setSpecies(key="X", type="fixed") + +n_step <- 300 +delta_p <- 0.001 +conc_threshold <- 14.27174 +ace_range <- 10**(seq(-1, 2, length.out = n_step)) +ace_range <- ace_range[abs(ace_range - conc_threshold) > 0.1] + +res_reg <- array(NA, dim=c(ncol(fit_results$res_par)-1, length(ace_range), 4), dimnames=list(iter=NULL, r=NULL, c=c("ace_conc", "via_acetate_pathway", "via_glc_upt", "via_tca"))) + +# create progress bar +pb <- txtProgressBar(min=0, max=ncol(fit_results$res_par)-1, style=3) + +for (j in seq(ncol(fit_results$res_par)-1)){ + + res_ace_regulation <- matrix(NA, nrow=length(ace_range), ncol=4, dimnames=list(r=NULL, c=c("ace_conc", "via_acetate_pathway", "via_glc_upt", "via_tca"))) + + for (i in seq(length(ace_range))){ + + rp <- c(fit_results$res_par[,j+1]) + names(rp) <- fit_results$res_par[,"parameter"] + model <- update_params(getCurrentModel(), rp) + + # set ace concentration + setSpecies(key="Ace_out{cell}", initial_concentration = ace_range[i], model=model) + applyInitialState(model=model) + + # get steady-state + res_ss_i <- runSteadyState(update_model=TRUE, model=model)$global_quantities + + # calculate control coefficients + res_MCA_R <- runMCA(model=model)$flux_control_coefficients_scaled + + # calculate elasticities + + # fix acetylCoA concentration to calculate elasticity of each pathway wrt acetate + setSpecies(key="AcCoA", type="fixed", model=model) + + # change acetate concentration + setSpecies(key="Ace_out{cell}", initial_concentration = ace_range[i]*(1+delta_p), model=model) + applyInitialState(model=model) + + # get steady-state + res_ss_i_eps <- runSteadyState(model=model)$global_quantities + + # calculate elasticities (using the more stable numerical method, both being equivalent) + #elasticities <- (log(abs(res_ss_i_eps$value)) - log(abs(res_ss_i$value))) / log(1+delta_p) + #print(elasticities) + elasticities <- (res_ss_i_eps$value - res_ss_i$value) / delta_p / res_ss_i_eps$value + #print(elasticities) + + # reset balance on accoa + setSpecies(key="AcCoA", type="reactions", model=model) + + # calculate response coefficient + # acetate via Pta-AckA + res_reg_ace_ace <- sum(res_MCA_R["(ackA)", c("(ackA)", "(pta)", "(ace_xch)")]) * elasticities[res_ss_i_eps$key == "Values[v_ace_net]"] + # acetate via glc uptake + res_reg_ace_glc <- res_MCA_R["(ackA)", "(glc_upt)"] * elasticities[res_ss_i_eps$key == "Values[v_glc_uptake]"] + # acetate via sink + res_reg_ace_tca <- res_MCA_R["(ackA)", "(sink)"] * elasticities[res_ss_i_eps$key == "Values[v_growth_rate]"] + + res_ace_regulation[i,] <- c(ace_range[i], res_reg_ace_ace, res_reg_ace_glc, res_reg_ace_tca) + } + + # save results + res_reg[j,,] <- res_ace_regulation + + # update the progress bar + setTxtProgressBar(pb, j) +} + +# close progress bar +close(pb) + +# plot regulation results +pdf(file="Figure 6.pdf", width = 7, height = 9) +par(mfrow=c(4,3)) + +xlab_main <- c(0.1, 1, 10, 100) +xlab_sec <- c(seq(0.2, 0.9, by=0.1), seq(2, 9, by=1), seq(20, 90, by=10)) +conc_threshold <- 14.5 + + +# plot partitioned response coefficients +fconc_max <- 1.12 +fconc_min <- 0.82 +lines_threshold(ace_range, apply(res_reg[,,"via_acetate_pathway"], 2, median), threshold=conc_threshold, new=TRUE, xaxt="n", las=1, xaxs="i", yaxs="i", col="#2E75B6", xlim=c(0.1,100), type="l", log="x", ylim=c(-5, 5), xlab="[acetate] (mM)", ylab="R_ace_pathway", lwd=1.2) +polygon(x=c(ace_range[ace_range < conc_threshold*fconc_min], rev(ace_range[ace_range < conc_threshold*fconc_min])), + y=c(apply(res_reg[,ace_range < conc_threshold*fconc_min,"via_acetate_pathway"], 2, max), rev(apply(res_reg[,ace_range < conc_threshold*fconc_min,"via_acetate_pathway"], 2, min))), + col="#2E75B655", border=NA) +polygon(x=c(ace_range[ace_range > conc_threshold*fconc_max], rev(ace_range[ace_range > conc_threshold*fconc_max])), + y=c(apply(res_reg[,ace_range > conc_threshold*fconc_max,"via_acetate_pathway"], 2, max), rev(apply(res_reg[,ace_range > conc_threshold*fconc_max,"via_acetate_pathway"], 2, min))), + col="#2E75B655", border=NA) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +abline(h=0) + +fconc_max <- 1.12 +fconc_min <- 0.84 +lines_threshold(ace_range, apply(res_reg[,,"via_glc_upt"], 2, median), threshold=conc_threshold, new=TRUE, xaxt="n", las=1, xaxs="i", yaxs="i", col="#D6685C", xlim=c(0.1,100), type="l", log="x", ylim=c(-5, 5), xlab="[acetate] (mM)", ylab="R_Glc_uptake", lwd=1.2) +polygon(x=c(ace_range[ace_range < conc_threshold*fconc_min], rev(ace_range[ace_range < conc_threshold*fconc_min])), + y=c(apply(res_reg[,ace_range < conc_threshold*fconc_min,"via_glc_upt"], 2, max), rev(apply(res_reg[,ace_range < conc_threshold*fconc_min,"via_glc_upt"], 2, min))), + col="#D6685C55", border=NA) +polygon(x=c(ace_range[ace_range > conc_threshold*fconc_max], rev(ace_range[ace_range > conc_threshold*fconc_max])), + y=c(apply(res_reg[,ace_range > conc_threshold*fconc_max,"via_glc_upt"], 2, max), rev(apply(res_reg[,ace_range > conc_threshold*fconc_max,"via_glc_upt"], 2, min))), + col="#D6685C55", border=NA) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +abline(h=0) + +lines_threshold(ace_range, apply(res_reg[,,"via_tca"], 2, median), threshold=conc_threshold, new=TRUE, xaxt="n", las=1, xaxs="i", yaxs="i", col="#70AD47", xlim=c(0.1,100), type="l", log="x", ylim=c(-5, 5), xlab="[acetate] (mM)", ylab="R_TCA", lwd=1.2) +polygon(x=c(ace_range[ace_range < conc_threshold*fconc_min], rev(ace_range[ace_range < conc_threshold*fconc_min])), + y=c(apply(res_reg[,ace_range < conc_threshold*fconc_min,"via_tca"], 2, max), rev(apply(res_reg[,ace_range < conc_threshold*fconc_min,"via_tca"], 2, min))), + col="#70AD4755", border=NA) +polygon(x=c(ace_range[ace_range > conc_threshold*fconc_max], rev(ace_range[ace_range > conc_threshold*fconc_max])), + y=c(apply(res_reg[,ace_range > conc_threshold*fconc_max,"via_tca"], 2, max), rev(apply(res_reg[,ace_range > conc_threshold*fconc_max,"via_tca"], 2, min))), + col="#70AD4755", border=NA) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +abline(h=0) + +# plot contribution of each pathway +fconc_max <- 1.14 +fconc_min <- 0.84 +contributio_ace <- res_reg[,,"via_acetate_pathway"]/apply(res_reg[,,2:4], 1:2, FUN=function(x) sum(abs(x))) +lines_threshold(ace_range, apply(contributio_ace, 2, median), threshold=conc_threshold, new=TRUE, xaxt="n", las=1, xaxs="i", yaxs="i", xlim=c(0.1,100), type="l", log="x", xlab="[acetate] (mM)", ylab="relative_R", ylim=c(-0.7, 0.7), lwd=1.2, col="#2E75B6") +polygon(x=c(ace_range[ace_range < conc_threshold*fconc_min], rev(ace_range[ace_range < conc_threshold*fconc_min])), + y=c(apply(contributio_ace[,ace_range < conc_threshold*fconc_min], 2, max), rev(apply(contributio_ace[,ace_range < conc_threshold*fconc_min], 2, min))), + col="#2E75B655", border=NA) +polygon(x=c(ace_range[ace_range > conc_threshold*fconc_max], rev(ace_range[ace_range > conc_threshold*fconc_max])), + y=c(apply(contributio_ace[,ace_range > conc_threshold*fconc_max], 2, max), rev(apply(contributio_ace[,ace_range > conc_threshold*fconc_max], 2, min))), + col="#2E75B655", border=NA) +axis(side = 1, at = xlab_main, labels = TRUE) +axis(side = 1, at = xlab_sec, labels = FALSE, tcl=-0.3) +abline(h=0) +contributio_glc_upt <- res_reg[,,"via_glc_upt"]/apply(res_reg[,,2:4], 1:2, FUN=function(x) sum(abs(x))) +lines_threshold(ace_range, apply(contributio_glc_upt, 2, median), threshold=conc_threshold, new=FALSE, type="l", col="#D6685C", lwd=1.2) +polygon(x=c(ace_range[ace_range < conc_threshold*fconc_min], rev(ace_range[ace_range < conc_threshold*fconc_min])), + y=c(apply(contributio_glc_upt[,ace_range < conc_threshold*fconc_min], 2, max), rev(apply(contributio_glc_upt[,ace_range < conc_threshold*fconc_min], 2, min))), + col="#D6685C55", border=NA) +polygon(x=c(ace_range[ace_range > conc_threshold*fconc_max], rev(ace_range[ace_range > conc_threshold*fconc_max])), + y=c(apply(contributio_glc_upt[,ace_range > conc_threshold*fconc_max], 2, max), rev(apply(contributio_glc_upt[,ace_range > conc_threshold*fconc_max], 2, min))), + col="#D6685C55", border=NA) +contributio_tca <- res_reg[,,"via_tca"]/apply(res_reg[,,2:4], 1:2, FUN=function(x) sum(abs(x))) +lines_threshold(ace_range, apply(contributio_tca, 2, median), threshold=conc_threshold, new=FALSE, type="l", col="#70AD47", lwd=1.2) +polygon(x=c(ace_range[ace_range < conc_threshold*fconc_min], rev(ace_range[ace_range < conc_threshold*fconc_min])), + y=c(apply(contributio_tca[,ace_range < conc_threshold*fconc_min], 2, max), rev(apply(contributio_tca[,ace_range < conc_threshold*fconc_min], 2, min))), + col="#70AD4755", border=NA) +polygon(x=c(ace_range[ace_range > conc_threshold*fconc_max], rev(ace_range[ace_range > conc_threshold*fconc_max])), + y=c(apply(contributio_tca[,ace_range > conc_threshold*fconc_max], 2, max), rev(apply(contributio_tca[,ace_range > conc_threshold*fconc_max], 2, min))), + col="#70AD4755", border=NA) + +dev.off() + + diff --git a/README.md b/README.md index 0edb188..d95af71 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,8 @@ -# Kinetic modeling of glucose and acetate metabolisms in *E. coli* - Millard et al., 2020. +# Kinetic modeling of glucose and acetate metabolisms in *E. coli* - Millard et al., 2021. ## Overview -This R script performs all analyses detailed in the following publication: +These R scripts perform all analyses detailed in the following publication: > Control and regulation of acetate overflow in *Escherichia coli* > @@ -11,17 +11,17 @@ This R script performs all analyses detailed in the following publication: All models are available in COPASI format in the directory `/model/cps/`, with the experimental data used for model calibration. The final kinetic model is also available in SBML format in the directory `/model/sbml/` and from the Biomodels database (http://www.ebi.ac.uk/biomodels/) under identifier MODEL2005050001. -Details on the calculations can be found in the [publication](https://doi.org/10.1101/2020.08.18.255356) and in the script `run_analysis.R`. +Details on the calculations can be found in the [original publication](https://doi.org/10.1101/2020.08.18.255356) and in the R scripts. ## Dependencies Some R packages are required. -`RColorBrewer` and `gplots` can be installed +`RColorBrewer`, `stringr` and `gplots` can be installed by running the following command in an R console: ```bash -install.packages(c("RColorBrewer", "gplots")) +install.packages(c("RColorBrewer", "gplots", "stringr")) ``` `CoRC` can be installed @@ -53,10 +53,22 @@ cd /home/usr/data/acetate_regulation/ R ``` -- run calculations: +- run calculations, starting from model construction to regulation analyses: ```bash -source("run_analysis.R") +source("1-Model_construction.R") +``` + +```bash +source("2-Model_validation.R") +``` + +```bash +source("3-Metabolic_control_analyses.R") +``` + +```bash +source("4-Regulation_analyses.R") ``` The code is open-source and available under GPLv3 license. diff --git a/misc.R b/misc.R new file mode 100644 index 0000000..9cb1723 --- /dev/null +++ b/misc.R @@ -0,0 +1,168 @@ +library(RColorBrewer) +library(gplots) +library(CoRC) +library(stringr) + +test_khi2 <- function(nb_points, k_val, nb_par){ + # Perform Khi2 statistical test. + # + # Args: + # nb_points (int): number of data points + # k_val (float): khi2 value (cost) + # nb_par (int): number of free parameters + # + # Returns (list): + # $'khi2 value' (float): khi2 value (cost) + # $'data points' (int): number of data points + # $'fitted parameters' (int): number of free parameters + # $'degrees of freedom' (int): degrees of freedom + # $'khi2 reduced value' (float): chi2 reduced value + # $'p-value, i.e. P(X^2<=value)' (float): p value + # $conclusion (str): message indicating whether the models fits (or not) the data at 95% confidence interval + # + df <- nb_points - nb_par + p_val <- pchisq(k_val, df=df) + khi2test <- list("khi2 value" = k_val, + "data points" = nb_points, + "fitted parameters" = nb_par, + "degrees of freedom" = df, + "khi2 reduced value" = k_val/df, + "p-value, i.e. P(X^2<=value)" = p_val) + if (p_val > 0.95){ + khi2test$conclusion <- "At level of 95% confidence, the model does not fit the data good enough with respect to the provided measurement SD." + }else{ + khi2test$conclusion <- "At level of 95% confidence, the model fits the data good enough with respect to the provided measurement SD." + } + return(khi2test) +} + +plot_points <- function(x, y, sd, col="black", offset=1.5, mode="v", cex=1){ + # Scatterplot with error bars. + # + # Args: + # x (vector): x coordinates + # y (vector): y coordinates + # sd (vector): error bars + # col (color code): color of points + # offset (float): width (or height if mode = 'v') of error bars + # mode ('v' or 'h'): errors of y (if mode='v') or x (if mode='h') + # + if (mode == "v"){ + segments(x0=x, y0=y-sd, x1=x, y1=y+sd) + segments(x0=x-offset, y0=y+sd, x1=x+offset, y1=y+sd) + segments(x0=x-offset, y0=y-sd, x1=x+offset, y1=y-sd) + }else if (mode == "h"){ + segments(x0=x-sd, y0=y, x1=x+sd, y1=y) + segments(x0=x+sd, y0=y-offset, x1=x+sd, y1=y+offset) + segments(x0=x-sd, y0=y-offset, x1=x-sd, y1=y+offset) + } + points(x, y, pch=21, bg=col, col="black", cex=cex) +} + +lines_threshold <- function(x, y, threshold, new, ...){ + # Split data according to a given threshold, and plot lines for + # each set. + # + # Args: + # x (vector): x coordinates + # y (vector): y coordinates + # threshold (float): value of x at which lines should not + # be connected + # new (bool): create a new plot if TRUE, otherwise add + # lines to an existing plot + # + id_1 <- (x < threshold) + id_2 <- (x > threshold) + if (new){ + plot(x[id_1], y[id_1], ...) + }else{ + suppressWarnings(lines(x[id_1], y[id_1], ...)) + } + suppressWarnings(lines(x[id_2], y[id_2], ...)) +} + +get_parameters_stats <- function(fit_results){ + li <- grep("]_0", fit_results$res_par$parameter, fixed=TRUE, invert=TRUE) + tmp <- matrix(NA, nrow=length(li), ncol=6, dimnames=list(par=fit_results$res_par$parameter[li], stats=c("mean", "median", "ci95_lb", "ci95_up", "sd", "rsd"))) + for (i in li){ + data <- as.numeric(unlist(fit_results$res_par[i,-1])) + tmp[fit_results$res_par[i,1],] <- c(mean(data), median(data), quantile(data, probs = c(0.025, 0.975)), sd(data), sd(data)/mean(data)) + } + return(tmp) +} + +update_params <- function(model, rp){ + setCurrentModel(model) + for (i in names(rp)){ + if (grepl("_0", i, fixed = TRUE)){ + next + }else if (grepl(".InitialValue", i, fixed = TRUE)){ + k <- str_remove(i, ".InitialValue") + #print(getGlobalQuantities(key=k)) + setGlobalQuantities(key=k, initial_value=rp[i]) + #print(getGlobalQuantities(key=k)) + }else{ + #print(getParameters(key=i)) + setParameters(key=i, value=rp[i]) + #print(getParameters(key=i)) + } + } + applyInitialState() + return(model) +} + +plot_with_ci <- function(fit_results, cond, specie, col, ...){ + if (specie %in% dimnames(fit_results[[cond]]$simulations)$specie){ + specie_id <- specie + }else{ + specie_id <- fit_results[[cond]]$mapping[specie] + } + plot(fit_results[[cond]]$simulations[1,,"Time"], apply(fit_results[[cond]]$simulations[,,specie_id], 2, mean), col=col, type="l", ...) + #polygon(x=c(fit_results[[cond]]$simulations[1,,"Time"], rev(fit_results[[cond]]$simulations[1,,"Time"])), + # y=c(apply(fit_results[[cond]]$simulations[,,specie_id], 2, max), rev(apply(fit_results[[cond]]$simulations[,,specie_id], 2, min))), + # col=paste(col, "33", sep=""), border=NA) + polygon(x=c(fit_results[[cond]]$simulations[1,,"Time"], rev(fit_results[[cond]]$simulations[1,,"Time"])), + y=c(apply(fit_results[[cond]]$simulations[,,specie_id], 2, max), rev(apply(fit_results[[cond]]$simulations[,,specie_id], 2, min))), + col=paste(col, "55", sep=""), border=NA) + plot_points(fit_results[[cond]]$data_exp$time, + fit_results[[cond]]$data_exp[, specie], + fit_results[[cond]]$sd[specie], offset=0.03, col=col, cex=1.2) +} + +plot_with_ci_2 <- function(x1, y1, y2, x2, y3, sd_y3, col, h=NULL, ...){ + plot(x1, y1, type="l", col=col, ...) + if (!is.null(h)){ + abline(h=h) + } + polygon(x=c(x1, rev(x1)), + y=c(apply(y2, 2, max), rev(apply(y2, 2, min))), + col=paste(col, "55", sep=""), border=NA) + plot_points(x2, y3, sd_y3, offset=0.002, col=col) +} + +plot_with_ci_3 <- function(sim_results, x, specie, col, ...){ + plot(sim_results[1,,x], apply(sim_results[,,specie], 2, mean), col=col, type="l", ...) + polygon(x=c(sim_results[1,,x], rev(sim_results[1,,x])), + y=c(apply(sim_results[,,specie], 2, max), rev(apply(sim_results[,,specie], 2, min))), + col=paste(col, "55", sep=""), border=NA) +} + +plot_no_ci <- function(fit_results, cond, specie, col, ...){ + if (specie %in% dimnames(fit_results[[cond]]$simulations)$specie){ + specie_id <- specie + }else{ + specie_id <- fit_results[[cond]]$mapping[specie] + } + plot(fit_results[[cond]]$simulations[,"Time"], fit_results[[cond]]$simulations[,specie_id], col=col, type="l", ...) + plot_points(fit_results[[cond]]$data_exp$time, + fit_results[[cond]]$data_exp[, specie], + fit_results[[cond]]$sd[specie], offset=0.03, col=col, cex=1.2) +} + +get_index_closest <- function(x, v){ + idx <- c() + for (i in x){ + idx <- c(idx, which.min(abs(v - i))) + } + return(idx) +} diff --git a/model/cps/Millard2020_Ecoli_glc_ace_model_1.cps b/model/cps/Millard2020_Ecoli_glc_ace_isotopic_model_1.cps similarity index 100% rename from model/cps/Millard2020_Ecoli_glc_ace_model_1.cps rename to model/cps/Millard2020_Ecoli_glc_ace_isotopic_model_1.cps diff --git a/model/cps/Millard2020_Ecoli_glc_ace_model_2.cps b/model/cps/Millard2020_Ecoli_glc_ace_isotopic_model_2.cps similarity index 100% rename from model/cps/Millard2020_Ecoli_glc_ace_model_2.cps rename to model/cps/Millard2020_Ecoli_glc_ace_isotopic_model_2.cps diff --git a/model/cps/Millard2020_Ecoli_glc_ace_model_3.cps b/model/cps/Millard2020_Ecoli_glc_ace_isotopic_model_3.cps similarity index 100% rename from model/cps/Millard2020_Ecoli_glc_ace_model_3.cps rename to model/cps/Millard2020_Ecoli_glc_ace_isotopic_model_3.cps diff --git a/model/cps/Millard2020_Ecoli_glc_ace_model_4.cps b/model/cps/Millard2020_Ecoli_glc_ace_isotopic_model_4.cps similarity index 100% rename from model/cps/Millard2020_Ecoli_glc_ace_model_4.cps rename to model/cps/Millard2020_Ecoli_glc_ace_isotopic_model_4.cps diff --git a/model/cps/Millard2020_Ecoli_glc_ace_isotopic_model_4_MC.cps b/model/cps/Millard2020_Ecoli_glc_ace_isotopic_model_4_MC.cps new file mode 100644 index 0000000..432883f --- /dev/null +++ b/model/cps/Millard2020_Ecoli_glc_ace_isotopic_model_4_MC.cps @@ -0,0 +1,3941 @@ + + + + + + + + + + + + 2019-09-24T12:10:06Z + + + + + + + Y*X*sink + + + + + + + + + + + + + + 2019-09-24T12:15:09Z + + + + + + + V*S/((KM+S)*(1+I/Ki)) + + + + + + + + + + + + + + + + 2019-09-24T11:32:52Z + + + + + + + AcP_0/AcP*v_ACK_Vmax*(ACP*ADP)/ (v_ACK_KmACP*v_ACK_KmADP)/((1+ACP/ v_ACK_KmACP+ACE/v_ACK_KmACE)*(1+ADP/v_ACK_KmADP+ATP/v_ACK_KmATP)) + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T11:32:52Z + + + + + + + AcP_1/AcP*v_ACK_Vmax*(ACP*ADP)/ (v_ACK_KmACP*v_ACK_KmADP)/((1+ACP/ v_ACK_KmACP+ACE/v_ACK_KmACE)*(1+ADP/v_ACK_KmADP+ATP/v_ACK_KmATP)) + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T11:33:45Z + + + + + + + Ace_0/Ace*v_ACK_Vmax*(ACE*ATP/v_ACK_Keq)/ (v_ACK_KmACP*v_ACK_KmADP)/((1+ACP/ v_ACK_KmACP+ACE/v_ACK_KmACE)*(1+ADP/v_ACK_KmADP+ATP/v_ACK_KmATP)) + + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T11:33:45Z + + + + + + + Ace_1/Ace*v_ACK_Vmax*(ACE*ATP/v_ACK_Keq)/ (v_ACK_KmACP*v_ACK_KmADP)/((1+ACP/ v_ACK_KmACP+ACE/v_ACK_KmACE)*(1+ADP/v_ACK_KmADP+ATP/v_ACK_KmATP)) + + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T11:31:22Z + + + + + + + (ACCOA_0/(ACCOA))*v_PTA_Vmax*(ACCOA*P)/ (v_PTA_KiACCOA*v_PTA_KmP)/(1+ACCOA/v_PTA_KiACCOA+P/ v_PTA_KiP+ACP/v_PTA_KiACP+COA/v_PTA_KiCOA+ACCOA*P/ (v_PTA_KiACCOA*v_PTA_KmP)+ACP*COA/(v_PTA_KmACP*v_PTA_KiCOA)) + + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T11:40:53Z + + + + + + + (ACCOA_1/(ACCOA))*v_PTA_Vmax*(ACCOA*P)/ (v_PTA_KiACCOA*v_PTA_KmP)/(1+ACCOA/v_PTA_KiACCOA+P/ v_PTA_KiP+ACP/v_PTA_KiACP+COA/v_PTA_KiCOA+ACCOA*P/ (v_PTA_KiACCOA*v_PTA_KmP)+ACP*COA/(v_PTA_KmACP*v_PTA_KiCOA)) + + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T11:32:27Z + + + + + + + (AcP_0/AcP)*v_PTA_Vmax*(ACP*COA/v_PTA_Keq)/ (v_PTA_KiACCOA*v_PTA_KmP)/(1+ACCOA/v_PTA_KiACCOA+P/ v_PTA_KiP+ACP/v_PTA_KiACP+COA/v_PTA_KiCOA+ACCOA*P/ (v_PTA_KiACCOA*v_PTA_KmP)+ACP*COA/(v_PTA_KmACP*v_PTA_KiCOA)) + + + + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T11:32:27Z + + + + + + + (AcP_1/AcP)*v_PTA_Vmax*(ACP*COA/v_PTA_Keq)/ (v_PTA_KiACCOA*v_PTA_KmP)/(1+ACCOA/v_PTA_KiACCOA+P/ v_PTA_KiP+ACP/v_PTA_KiACP+COA/v_PTA_KiCOA+ACCOA*P/ (v_PTA_KiACCOA*v_PTA_KmP)+ACP*COA/(v_PTA_KmACP*v_PTA_KiCOA)) + + + + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T14:21:54Z + + + + + + + A/At*V*At^n/((At^n+Km)*(1+I/Ki)) + + + + + + + + + + + + + + + + + + 2020-02-25T11:27:00Z + + + + + + + Sx/S*(Vmax*(S/KM/((1+S/KM+P/KM)))) + + + + + + + + + + + + + + + + + 2019-09-24T11:27:30Z + + + + + + + + Model of E. coli glucose and acetate metabolism + +Kinetic and isotopic model, inhibition of glycolysis and TCA by acetate. + +model #4 of Millard et al., 2020 + + + + + + + + + 2019-09-24T11:27:54Z + + + + + + + + + + + + + + + 2019-09-24T11:27:54Z + + + + + + + if(<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[_STEADY_STATE],Reference=InitialValue> == 0,-<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[X],Reference=Concentration>*<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Reactions[glc_upt],Reference=Flux>*<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[volume],Reference=InitialValue>,0) + + + + + + + + + 2019-09-24T11:43:28Z + + + + + + + + + + + + + 2019-09-24T11:43:28Z + + + + + + + + + + + + + 2019-09-24T11:46:12Z + + + + + + + + + + + + + 2019-09-24T11:46:12Z + + + + + + + + + + + + + 2019-09-24T11:47:36Z + + + + + + + + + + + + + 2019-09-24T11:47:47Z + + + + + + + + + + + + + 2019-09-24T12:08:24Z + + + + + + + if(<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[_STEADY_STATE],Reference=InitialValue> == 0,<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Reactions[growth],Reference=Flux>,0) + + + + + + + + + 2019-09-24T14:21:04Z + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[AcCoA0],Reference=Concentration>+<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[AcCoA1],Reference=Concentration> + + + + + + + + + 2019-09-24T14:55:11Z + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[Ace0],Reference=Concentration>+<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[Ace1],Reference=Concentration> + + + + + + + + + 2019-09-24T14:55:41Z + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[AcP0],Reference=Concentration>+<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[AcP1],Reference=Concentration> + + + + + + + + + 2020-02-26T10:12:37Z + + + + + + + if(<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[_STEADY_STATE],Reference=InitialValue> == 0,<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[X],Reference=Concentration>*<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[volume],Reference=InitialValue>*(<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Reactions[trp_ace_f0],Reference=Flux>-<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Reactions[trp_ace_r0],Reference=Flux>),0) + + + + + + + + + 2020-02-26T10:12:50Z + + + + + + + if(<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[_STEADY_STATE],Reference=InitialValue> == 0,<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[X],Reference=Concentration>*<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[volume],Reference=InitialValue>*(<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Reactions[trp_ace_f1],Reference=Flux>-<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Reactions[trp_ace_r1],Reference=Flux>),0) + + + + + + + + + 2020-02-26T10:12:59Z + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[Ace0_out],Reference=Concentration>+<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[Ace1_out],Reference=Concentration> + + + + + + + + + + + 2019-09-24T11:36:27Z + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[AcCoAt],Reference=Concentration> + + + + + + + + + 2019-09-24T11:50:33Z + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[AcP0],Reference=Concentration>+<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[AcP1],Reference=Concentration> + + + + + + + + + 2019-09-24T11:51:40Z + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[Ace0],Reference=Concentration>+<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[Ace1],Reference=Concentration> + + + + + + + + + 2019-09-24T11:53:52Z + + + + + + + + + + + + + 2019-09-24T11:54:01Z + + + + + + + + + + + + + 2019-09-24T11:54:06Z + + + + + + + + + + + + + 2019-09-24T11:54:10Z + + + + + + + + + + + + + 2019-09-24T11:54:14Z + + + + + + + + + + + + + 2019-09-24T11:54:18Z + + + + + + + + + + + + + 2019-09-24T11:54:21Z + + + + + + + + + + + + + 2019-09-24T11:54:41Z + + + + + + + + + + + + + 2019-09-24T11:54:47Z + + + + + + + + + + + + + 2019-09-24T11:56:34Z + + + + + + + + + + + + + 2019-09-24T11:59:26Z + + + + + + + + + + + + + 2019-09-24T11:59:31Z + + + + + + + + + + + + + 2019-09-24T11:59:35Z + + + + + + + + + + + + + 2019-09-24T11:59:38Z + + + + + + + + + + + + + 2019-09-24T11:59:42Z + + + + + + + + + + + + + 2019-09-24T11:59:45Z + + + + + + + + + + + + + 2019-09-24T11:59:50Z + + + + + + + + + + + + + 2019-09-24T11:59:53Z + + + + + + + + + + + + + 2019-09-24T12:03:16Z + + + + + + + + + + + + + 2019-09-24T12:03:26Z + + + + + + + + + + + + + 2019-09-24T12:08:35Z + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Reactions[sink0],Reference=Flux>+<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Reactions[sink1],Reference=Flux> + + + + + + + + + 2019-09-24T12:10:39Z + + + + + + + + + + + + + 2019-09-24T12:53:57Z + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[Ace1],Reference=Concentration>/(<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[Ace0],Reference=Concentration>+<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Compartments[cell],Vector=Metabolites[Ace1],Reference=Concentration>) + + + + + + + + + 2019-09-24T15:58:00Z + + + + + + + + + + + + + 2019-09-24T19:22:58Z + + + + + + + + + + + + + 2019-09-24T19:23:05Z + + + + + + + + + + + + + 2020-02-25T11:25:57Z + + + + + + + + + + + + + 2020-02-26T10:41:17Z + + + + + + + + + + + + + 2020-02-26T10:41:34Z + + + + + + + + + + + + + 2020-02-26T13:37:18Z + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[v_sink],Reference=Value>*<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[Y],Reference=InitialValue> + + + + + + + + + 2020-02-26T14:22:55Z + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Reactions[trp_ace_f0],Reference=Flux>+<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Reactions[trp_ace_f1],Reference=Flux>-<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Reactions[trp_ace_r0],Reference=Flux>-<CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Reactions[trp_ace_r1],Reference=Flux> + + + + + + + + + 2020-04-29T19:12:37Z + + + + + + + + + + + + + + + 2019-09-24T11:27:37Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T11:28:06Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T11:28:35Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T11:28:46Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T11:34:27Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T11:35:06Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T11:28:35Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T11:34:27Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T11:28:06Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T11:35:06Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T11:28:46Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2019-09-24T12:07:54Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2020-02-26T10:33:46Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2020-02-26T10:33:46Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2020-02-26T10:33:46Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2020-02-26T10:33:46Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +2020-12-21T11:38:18Z + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_Vmax],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[P],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[CoA],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KiACCOA],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KmP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KiP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KiACP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KiCOA],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KmACP],Reference=InitialValue> + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_Vmax],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[ADP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_KmACP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_KmADP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_KmACE],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[ATP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_KmATP],Reference=InitialValue> + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_sink_Vmax],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[sink_n],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_sink_Km],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_sink_Ki],Reference=InitialValue> + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_Vmax],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[ATP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_Keq],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_KmACP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_KmADP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_KmACE],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[ADP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_KmATP],Reference=InitialValue> + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_Vmax],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[CoA],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_Keq],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KiACCOA],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KmP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[P],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KiP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KiACP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KiCOA],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KmACP],Reference=InitialValue> + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_Vmax],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[ADP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_KmACP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_KmADP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_KmACE],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[ATP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_KmATP],Reference=InitialValue> + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_Vmax],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[ATP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_Keq],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_KmACP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_KmADP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_KmACE],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[ADP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_ACK_KmATP],Reference=InitialValue> + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_Vmax],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[P],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[CoA],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KiACCOA],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KmP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KiP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KiACP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KiCOA],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KmACP],Reference=InitialValue> + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_Vmax],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[CoA],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_Keq],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KiACCOA],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KmP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[P],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KiP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KiACP],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KiCOA],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_PTA_KmACP],Reference=InitialValue> + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[sink_n],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_sink_Km],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_sink_Ki],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_sink_Vmax],Reference=InitialValue> + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[Y],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[v_sink],Reference=InitialValue> + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_TRP_Vmax],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_TRP_Km],Reference=InitialValue> + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_TRP_Vmax],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_TRP_Km],Reference=InitialValue> + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_TRP_Vmax],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_TRP_Km],Reference=InitialValue> + + + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_TRP_Vmax],Reference=InitialValue> + + + + + <CN=Root,Model=E. coli model - Glc-Ace metabolism,Vector=Values[r_TRP_Km],Reference=InitialValue> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 7.7685596306705115e+21 4.2154949148271837e+19 5.5403817796619515e+20 6.0153929090582774e+18 6.0221408570000005e+18 5.7661288689345823e+20 0 8.7923256512199996e+20 0 5.9326521219652067e+21 8.8525470597899996e+20 5.7661288689345823e+20 5.9326521219652067e+21 5.6005357087525346e+20 1.47 9.8514004617962829 0.95748820990000005 29662.247991208762 0 2.9604186319476407 376.69231041352646 1 9565521.7634551357 0.20000000000000001 2.6000000000000001 2.6000000000000001 0.20000000000000001 0.029000000000000001 0.69999999999999996 10 1.22 0.0281 336940.01763224584 0.60599999999999998 0.16 0.5 7 2.3999999999999999 0.070000000000000007 174 740859.87349901034 24.759179607035055 9.980425734506176e-05 0.0017700000000000001 1 2.3261654310710522 0 33.153580532659561 480034.55363364267 60 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Automatically generated report. + +