Sensitivities of OFMs in E. coli core

using ElementaryFluxModes

import AbstractFBCModels as A
import AbstractFBCModels.CanonicalModel as CM
using JSONFBCModels
import FastDifferentiation as F
const Ex = F.Node
import DifferentiableMetabolism as D
using COBREXA
import COBREXA as X
using JSON
import Tulip as T
using CairoMakie

It has been proven that enzyme constrained models with K constraints will use a maximum of K EFMs in their optimal solution (de Groot 2019).

Here, we take the E. coli core model with two enzyme constraints, find the optimal solution, and then find the EFMs in the optimal solution.

!isfile("e_coli_core.json") &&
    download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json")

include("../../test/data_static.jl")

model = convert(CM.Model, X.load_model("e_coli_core.json"))

AbstractFBCModels.CanonicalModel.Model(
  reactions = Dict{String, AbstractFBCModels.CanonicalModel.Reaction}("ACALD" => AbstractFBCModels.CanonicalModel.Reaction("Acetaldehyde dehydrogenas…
  metabolites = Dict{String, AbstractFBCModels.CanonicalModel.Metabolite}("glu__L_c" => AbstractFBCModels.CanonicalModel.Metabolite("L-Glutamate", "c…
  genes = Dict{String, AbstractFBCModels.CanonicalModel.Gene}("b4301" => AbstractFBCModels.CanonicalModel.Gene("sgcE", Dict("sbo" => ["SBO:0000243"],…
  couplings = Dict{String, AbstractFBCModels.CanonicalModel.Coupling}(),
)

unconstrain glucose

model.reactions["EX_glc__D_e"].lower_bound = -1000.0

-1000.0

constrain PFL to zero

model.reactions["PFL"].upper_bound = 0.0

rid_kcat = Dict(k => Ex(Symbol(k)) for (k, _) in ecoli_core_reaction_kcats)
kcats = Symbol.(keys(ecoli_core_reaction_kcats))

float_reaction_isozymes = Dict{String,Dict{String,X.Isozyme}}()
for (rid, rxn) in model.reactions
    grrs = rxn.gene_association_dnf
    isnothing(grrs) && continue # skip if no grr available
    haskey(ecoli_core_reaction_kcats, rid) || continue # skip if no kcat data available
    for (i, grr) in enumerate(grrs)

        kcat = ecoli_core_reaction_kcats[rid] * 3.6 # change unit to k/h

        d = get!(float_reaction_isozymes, rid, Dict{String,X.Isozyme}())
        d["isozyme_$i"] = X.Isozyme(
            gene_product_stoichiometry = Dict(grr .=> fill(1.0, size(grr))), # assume subunit stoichiometry of 1 for all isozymes
            kcat_forward = kcat, # assume forward and reverse have the same kcat
            kcat_reverse = kcat,
        )
    end
end

gene_product_molar_masses = Dict(k => v for (k, v) in ecoli_core_gene_product_masses)

capacity = [
    (
        "membrane", # enzyme group names must be of type Symbol
        [x for (x, y) in ecoli_core_subcellular_location if y == "membrane"],
        15.0,
    ),
    ("cytosol", [x for (x, y) in ecoli_core_subcellular_location if y == "cytosol"], 35.0),
]

ec_solution = X.enzyme_constrained_flux_balance_analysis(
    model;
    reaction_isozymes = float_reaction_isozymes,
    gene_product_molar_masses,
    capacity,
    optimizer = T.Optimizer,
)

ec_solution

ConstraintTrees.Tree{Float64} with 12 elements:
  :coupling                     => ConstraintTrees.Tree{Float64}(#= 0 elements =#)
  :flux_stoichiometry           => ConstraintTrees.Tree{Float64}(#= 72 elements =#)
  :fluxes                       => ConstraintTrees.Tree{Float64}(#= 95 elements =#)
  :fluxes_forward               => ConstraintTrees.Tree{Float64}(#= 95 elements =#)
  :fluxes_reverse               => ConstraintTrees.Tree{Float64}(#= 95 elements =#)
  :gene_product_amounts         => ConstraintTrees.Tree{Float64}(#= 130 elements =#)
  :gene_product_capacity        => ConstraintTrees.Tree{Float64}(#= 2 elements =#)
  :isozyme_flux_forward_balance => ConstraintTrees.Tree{Float64}(#= 62 elements =#)
  :isozyme_flux_reverse_balance => ConstraintTrees.Tree{Float64}(#= 62 elements =#)
  :isozyme_forward_amounts      => ConstraintTrees.Tree{Float64}(#= 62 elements =#)
  :isozyme_reverse_amounts      => ConstraintTrees.Tree{Float64}(#= 62 elements =#)
  :objective                    => 0.740724

This solution is both producing acetate and consuming oxygen, therefore it looks like overflow metabolism.

ec_solution.fluxes["EX_ac_e"]
ec_solution.fluxes["EX_o2_e"]

-11.119604376947802

Since we had two enzyme pools, the membrane and the cytosol, this optimal solution will be composed of at most two EFMs. In order to find out how small parameter changes will affect this optimal composition, we must first remove inactive reactions, so that the optimal solution is a unique superposition of the two EFMs, and we can differentiate it.

This solution contains many inactive reactions

sort(collect(ec_solution.fluxes), by = ComposedFunction(abs, last))

95-element Vector{Pair{Symbol, Union{Float64, ConstraintTrees.Tree{Float64}}}}:
    :EX_fru_e => 0.0
    :EX_fum_e => 0.0
 :EX_gln__L_e => 0.0
 :EX_mal__L_e => 0.0
     :FRUpts2 => 0.0
     :FUMt2_2 => 0.0
      :GLNabc => 0.0
     :MALt2_2 => 0.0
         :PFL => 0.0
    :EX_for_e => 1.1751401296095953e-14
              ⋮
    :SUCCt2_2 => 45.1017515034113
      :SUCCt3 => 45.101751508128025
         :PDH => 48.64218610958069
         :ENO => 53.247640002202544
         :PGM => -53.24764000220255
    :EX_co2_e => 53.669000943520395
        :CO2t => -53.66900094352041
        :GAPD => 54.355763660340095
         :PGK => -54.355763660340095

And also many inactive gene products.

sort(collect(ec_solution.gene_product_amounts), by = last)

130-element Vector{Pair{Symbol, Union{Float64, ConstraintTrees.Tree{Float64}}}}:
 :b0902 => 0.0
 :b0903 => 0.0
 :b2579 => 0.0
 :b3114 => 0.0
 :b3951 => 0.0
 :b3952 => 0.0
 :b0809 => 2.9237052619941317e-16
 :b0810 => 2.9237052619941317e-16
 :b0811 => 2.9237052619941317e-16
 :b1702 => 3.9759214856671975e-12
        ⋮
 :b2415 => 0.0347023816798353
 :b2416 => 0.0347023816798353
 :b0978 => 0.040328743170719945
 :b0979 => 0.040328743170719945
 :b3528 => 0.053498438408503594
 :b2779 => 0.07065207134510235
 :b1779 => 0.11726330582717176
 :b1478 => 0.15248718085537333
 :b2926 => 0.2619504379992874

Let us start by pruning the model.

flux_zero_tol = 1e-6 # these bounds make a real difference!
gene_zero_tol = 1e-6
pruned_model, pruned_reaction_isozymes = D.prune_model(
    model,
    ec_solution.fluxes,
    ec_solution.gene_product_amounts,
    float_reaction_isozymes,
    ec_solution.isozyme_forward_amounts,
    ec_solution.isozyme_reverse_amounts,
    flux_zero_tol,
    gene_zero_tol,
)

(AbstractFBCModels.CanonicalModel.Model(Dict{String, AbstractFBCModels.CanonicalModel.Reaction}("ACALD" => AbstractFBCModels.CanonicalModel.Reaction("Acetaldehyde dehydrogenase (acetylating)", 0.0, 1000.0, Dict("coa_c" => 1.0, "nad_c" => 1.0, "accoa_c" => -1.0, "acald_c" => 1.0, "h_c" => -1.0, "nadh_c" => -1.0), 0.0, [["b0351"], ["b1241"]], Dict("bigg.reaction" => ["ACALD"], "sabiork" => ["163"], "metanetx.reaction" => ["MNXR95210"], "rhea" => ["23290", "23291", "23289", "23288"], "sbo" => ["SBO:0000176"], "seed.reaction" => ["rxn00171"], "kegg.reaction" => ["R00228"], "biocyc" => ["META:ACETALD-DEHYDROG-RXN"], "ec-code" => ["1.2.1.10"]), Dict("original_bigg_ids" => ["ACALD"])), "PTAr" => AbstractFBCModels.CanonicalModel.Reaction("Phosphotransacetylase", 0.0, 1000.0, Dict("coa_c" => 1.0, "accoa_c" => -1.0, "actp_c" => 1.0, "pi_c" => -1.0), 0.0, [["b2297"], ["b2458"]], Dict("bigg.reaction" => ["PTAr"], "sabiork" => ["72"], "metanetx.reaction" => ["MNXR103319"], "rhea" => ["19521", "19523", "19522", "19524"], "sbo" => ["SBO:0000176"], "seed.reaction" => ["rxn00173"], "kegg.reaction" => ["R00230"], "biocyc" => ["META:PHOSACETYLTRANS-RXN"], "ec-code" => ["2.3.1.8"]), Dict("original_bigg_ids" => ["PTAr"])), "ALCD2x" => AbstractFBCModels.CanonicalModel.Reaction("Alcohol dehydrogenase (ethanol)", 0.0, 1000.0, Dict("nad_c" => 1.0, "acald_c" => -1.0, "etoh_c" => 1.0, "h_c" => -1.0, "nadh_c" => -1.0), 0.0, [["b0356"], ["b1241"], ["b1478"]], Dict("bigg.reaction" => ["ALCD2x"], "sabiork" => ["597"], "metanetx.reaction" => ["MNXR95725"], "rhea" => ["25292", "25290", "25291", "25293"], "sbo" => ["SBO:0000176"], "seed.reaction" => ["rxn00543"], "kegg.reaction" => ["R00754"], "biocyc" => ["META:ALCOHOL-DEHYDROG-RXN"], "reactome.reaction" => ["R-MMU-71707", "R-XTR-71707", "R-CFA-71707", "R-OSA-71707", "R-RNO-71707", "R-HSA-71707", "R-SSC-71707", "R-GGA-71707", "R-SPO-71707", "R-DME-71707", "R-DRE-71707", "R-BTA-71707", "R-ATH-71707", "R-SCE-71707", "R-TGU-71707", "R-CEL-71707", "R-DDI-71707"], "ec-code" => ["1.1.1.71", "1.1.1.1"]…), Dict("original_bigg_ids" => ["ALCD2x"])), "PDH" => AbstractFBCModels.CanonicalModel.Reaction("Pyruvate dehydrogenase", 0.0, 1000.0, Dict("coa_c" => -1.0, "nad_c" => -1.0, "accoa_c" => 1.0, "co2_c" => 1.0, "pyr_c" => -1.0, "nadh_c" => 1.0), 0.0, [["b0114", "b0115", "b0116"]], Dict("bigg.reaction" => ["PDH"], "sabiork" => ["523"], "metanetx.reaction" => ["MNXR102425"], "rhea" => ["28043", "28045", "28044", "28042"], "sbo" => ["SBO:0000176"], "seed.reaction" => ["rxn00154"], "kegg.reaction" => ["R00209"], "biocyc" => ["META:PYRUVDEH-RXN"], "reactome.reaction" => ["R-RNO-71397", "R-OSA-71397", "R-GGA-373177", "R-DME-71397", "R-TGU-71397", "R-XTR-71397", "R-CFA-71397", "R-BTA-71397", "R-HSA-71397", "R-GGA-71397", "R-SPO-71397", "R-CEL-71397", "R-DDI-71397", "R-DRE-71397", "R-SSC-71397", "R-SCE-71397", "R-ATH-71397", "R-MMU-71397"], "ec-code" => ["1.2.1", "1.8.1.4", "1.2.1.51", "1.2.4.1", "2.3.1.12"]…), Dict("original_bigg_ids" => ["PDH"])), "CO2t" => AbstractFBCModels.CanonicalModel.Reaction("CO2 transporter via diffusion", 0.0, 1000.0, Dict("co2_e" => 1.0, "co2_c" => -1.0), 0.0, [["s0001"]], Dict("metanetx.reaction" => ["MNXR96810"], "sbo" => ["SBO:0000185"], "seed.reaction" => ["rxn05467", "rxn08237", "rxn09706", "rxn09821", "rxn09775", "rxn09876", "rxn08238", "rxn09860"], "biocyc" => ["META:TRANS-RXN0-545"], "reactome.reaction" => ["R-DDI-1247645", "R-CFA-1237069", "R-TGU-1247649", "R-CFA-1237042", "R-CEL-1247645", "R-XTR-1237069", "R-HSA-1247645", "R-OSA-1237042", "R-GGA-1237042", "R-BTA-1247645"  …  "R-MMU-1247645", "R-DME-1237069", "R-RNO-1237042", "R-TGU-1247645", "R-BTA-1237042", "R-GGA-1237069", "R-DME-1247649", "R-SSC-1237069", "R-DRE-1237069", "R-ATH-1247649"], "bigg.reaction" => ["CO2t"]), Dict("original_bigg_ids" => ["CO2t"])), "PYK" => AbstractFBCModels.CanonicalModel.Reaction("Pyruvate kinase", 0.0, 1000.0, Dict("adp_c" => -1.0, "atp_c" => 1.0, "pep_c" => -1.0, "pyr_c" => 1.0, "h_c" => -1.0), 0.0, [["b1676"], ["b1854"]], Dict("bigg.reaction" => ["PYK"], "sabiork" => ["9"], "metanetx.reaction" => ["MNXR103371"], "rhea" => ["18160", "18159", "18158", "18157"], "sbo" => ["SBO:0000176"], "seed.reaction" => ["rxn00148"], "kegg.reaction" => ["R00200"], "biocyc" => ["META:PEPDEPHOS-RXN"], "reactome.reaction" => ["R-SSC-71670", "R-GGA-353056", "R-DRE-71670", "R-OSA-71670", "R-SPO-71670", "R-DDI-71670", "R-SCE-71670", "R-CEL-71670", "R-PFA-71670", "R-HSA-71670", "R-TGU-71670", "R-DME-71670", "R-ATH-71670", "R-BTA-71670", "R-MMU-71670", "R-RNO-71670", "R-GGA-71670", "R-CFA-71670", "R-XTR-71670"], "ec-code" => ["2.7.1.40"]…), Dict("original_bigg_ids" => ["PYK"])), "EX_nh4_e" => AbstractFBCModels.CanonicalModel.Reaction("Ammonia exchange", 0.0, 1000.0, Dict("nh4_e" => 1.0), 0.0, nothing, Dict("sabiork" => ["11683"], "metanetx.reaction" => ["MNXR101950"], "rhea" => ["28749", "28748", "28750", "28747"], "sbo" => ["SBO:0000627"], "seed.reaction" => ["rxn13364", "rxn09835", "rxn05466", "rxn08987", "rxn08986", "rxn09736"], "biocyc" => ["META:TRANS-RXN0-206", "META:RXN-9615", "META:TRANS-RXN0-544"], "reactome.reaction" => ["R-CEL-444416", "R-SSC-444416", "R-DDI-444416", "R-MMU-444416", "R-DRE-444416", "R-GGA-444416", "R-CFA-444416", "R-HSA-444416", "R-XTR-444416", "R-BTA-444416", "R-RNO-444416", "R-TGU-444416", "R-DME-444416"], "bigg.reaction" => ["EX_nh4_e"]), Dict("original_bigg_ids" => ["EX_nh4_e"])), "CS" => AbstractFBCModels.CanonicalModel.Reaction("Citrate synthase", 0.0, 1000.0, Dict("coa_c" => 1.0, "h2o_c" => -1.0, "accoa_c" => -1.0, "cit_c" => 1.0, "oaa_c" => -1.0, "h_c" => 1.0), 0.0, [["b0720"]], Dict("bigg.reaction" => ["CS"], "sabiork" => ["267"], "metanetx.reaction" => ["MNXR96920"], "rhea" => ["16847", "16846", "16845", "16848"], "sbo" => ["SBO:0000176"], "seed.reaction" => ["rxn00256"], "kegg.reaction" => ["R00351"], "biocyc" => ["META:CITSYN-RXN", "META:RXN-14905"], "reactome.reaction" => ["R-GGA-373006", "R-GGA-70975", "R-CEL-70975", "R-HSA-70975", "R-DRE-70975", "R-MMU-70975", "R-ATH-70975", "R-SPO-70975", "R-OSA-70975", "R-SCE-70975", "R-SSC-70975", "R-DDI-70975", "R-RNO-70975", "R-DME-70975", "R-PFA-70975", "R-XTR-70975", "R-CFA-70975", "R-BTA-70975"], "ec-code" => ["2.3.3.16", "2.3.3.1", "2.3.3.3"]…), Dict("original_bigg_ids" => ["CS"])), "PGM" => AbstractFBCModels.CanonicalModel.Reaction("Phosphoglycerate mutase", 0.0, 1000.0, Dict("3pg_c" => -1.0, "2pg_c" => 1.0), 0.0, [["b0755"], ["b3612"], ["b4395"]], Dict("bigg.reaction" => ["PGM"], "sabiork" => ["7641"], "metanetx.reaction" => ["MNXR102547"], "rhea" => ["15902", "15903", "15901", "15904"], "sbo" => ["SBO:0000176"], "seed.reaction" => ["rxn01106"], "kegg.reaction" => ["R01518"], "biocyc" => ["META:3PGAREARR-RXN", "META:RXN-15513"], "reactome.reaction" => ["R-XTR-71654", "R-SSC-71654", "R-DME-71445", "R-XTR-71445", "R-RNO-71445", "R-GGA-71654", "R-SSC-71445", "R-GGA-71445", "R-DDI-71654", "R-SPO-71654"  …  "R-DDI-71445", "R-RNO-71654", "R-DRE-71445", "R-OSA-71445", "R-PFA-71654", "R-DME-71654", "R-TGU-71445", "R-SCE-71445", "R-ATH-71445", "R-SCE-71654"], "ec-code" => ["5.4.2.1", "5.4.2.11", "5.4.2.12"]…), Dict("original_bigg_ids" => ["PGM"])), "TKT1" => AbstractFBCModels.CanonicalModel.Reaction("Transketolase", 0.0, 1000.0, Dict("g3p_c" => 1.0, "r5p_c" => -1.0, "s7p_c" => 1.0, "xu5p__D_c" => -1.0), 0.0, [["b2465"], ["b2935"]], Dict("bigg.reaction" => ["TKT1"], "metanetx.reaction" => ["MNXR104868"], "sbo" => ["SBO:0000176"], "ec-code" => ["2.2.1.1"]), Dict("original_bigg_ids" => ["TKT1"]))…), Dict{String, AbstractFBCModels.CanonicalModel.Metabolite}("glu__L_c" => AbstractFBCModels.CanonicalModel.Metabolite("L-Glutamate", "c", Dict("C" => 5, "N" => 1, "H" => 8, "O" => 4), -1, 0.0, Dict("kegg.drug" => ["D00007", "D04341"], "kegg.compound" => ["C00025", "C00302"], "sbo" => ["SBO:0000247"], "biocyc" => ["META:Glutamates", "META:GLT"], "chebi" => ["CHEBI:76051", "CHEBI:21301", "CHEBI:29985", "CHEBI:42825", "CHEBI:29987", "CHEBI:18237", "CHEBI:24314", "CHEBI:16015", "CHEBI:13107", "CHEBI:5431", "CHEBI:21304", "CHEBI:6224", "CHEBI:14321", "CHEBI:29988"], "bigg.metabolite" => ["glu__L"], "seed.compound" => ["cpd27177", "cpd00023", "cpd19002"], "sabiork" => ["73", "2010"], "metanetx.chemical" => ["MNXM89557"], "inchi_key" => ["WHUUTDBJXJRKMK-VKHMYHEASA-M"]…), Dict("original_bigg_ids" => ["glu_L_c"])), "icit_c" => AbstractFBCModels.CanonicalModel.Metabolite("Isocitrate", "c", Dict("C" => 6, "H" => 5, "O" => 7), -3, 0.0, Dict("chebi" => ["CHEBI:16087", "CHEBI:36454", "CHEBI:5998", "CHEBI:24886", "CHEBI:36453", "CHEBI:14465", "CHEBI:30887", "CHEBI:24884"], "metanetx.chemical" => ["MNXM89661"], "inchi_key" => ["ODBLHEXUDAPZAU-UHFFFAOYSA-K"], "kegg.compound" => ["C00311"], "hmdb" => ["HMDB00193"], "sabiork" => ["2013"], "sbo" => ["SBO:0000247"], "bigg.metabolite" => ["icit"], "seed.compound" => ["cpd00260"]), Dict("original_bigg_ids" => ["icit_c"])), "co2_c" => AbstractFBCModels.CanonicalModel.Metabolite("CO2 CO2", "c", Dict("C" => 1, "O" => 2), 0, 0.0, Dict("envipath" => ["650babc9-9d68-4b73-9332-11972ca26f7b/compound/2ec3da94-5f50-4525-81b1-5607c5c7a3d3", "32de3cf4-e3e6-4168-956e-32fa5ddb0ce1/compound/05f60af4-0a3f-4ead-9a29-33bb0f123379"], "kegg.drug" => ["D00004"], "kegg.compound" => ["C00011"], "sbo" => ["SBO:0000247"], "biocyc" => ["META:CARBON-DIOXIDE"], "chebi" => ["CHEBI:23011", "CHEBI:3283", "CHEBI:48829", "CHEBI:16526", "CHEBI:13283", "CHEBI:13285", "CHEBI:13284", "CHEBI:13282"], "bigg.metabolite" => ["co2"], "seed.compound" => ["cpd00011"], "sabiork" => ["1266"], "metanetx.chemical" => ["MNXM13"]…), Dict("original_bigg_ids" => ["co2_c"])), "6pgl_c" => AbstractFBCModels.CanonicalModel.Metabolite("6-phospho-D-glucono-1,5-lactone", "c", Dict("C" => 6, "P" => 1, "H" => 9, "O" => 9), -2, 0.0, Dict("kegg.compound" => ["C01236"], "sbo" => ["SBO:0000247"], "biocyc" => ["META:D-6-P-GLUCONO-DELTA-LACTONE"], "chebi" => ["CHEBI:12233", "CHEBI:4160", "CHEBI:57955", "CHEBI:16938", "CHEBI:20989", "CHEBI:12958"], "bigg.metabolite" => ["6pgl"], "seed.compound" => ["cpd00911"], "sabiork" => ["1366"], "metanetx.chemical" => ["MNXM429"], "inchi_key" => ["IJOJIVNDFQSGAB-SQOUGZDYSA-L"], "hmdb" => ["HMDB62628", "HMDB01127"]…), Dict("original_bigg_ids" => ["6pgl_c"])), "cit_c" => AbstractFBCModels.CanonicalModel.Metabolite("Citrate", "c", Dict("C" => 6, "H" => 5, "O" => 7), -3, 0.0, Dict("kegg.drug" => ["D00037"], "kegg.compound" => ["C00158"], "sbo" => ["SBO:0000247"], "biocyc" => ["META:CIT"], "chebi" => ["CHEBI:35802", "CHEBI:133748", "CHEBI:35809", "CHEBI:76049", "CHEBI:41523", "CHEBI:35810", "CHEBI:30769", "CHEBI:13999", "CHEBI:23322", "CHEBI:35808", "CHEBI:42563", "CHEBI:16947", "CHEBI:132362", "CHEBI:35804", "CHEBI:23321", "CHEBI:3727", "CHEBI:35806"], "bigg.metabolite" => ["cit"], "seed.compound" => ["cpd00137"], "sabiork" => ["1952"], "metanetx.chemical" => ["MNXM131"], "inchi_key" => ["KRKNYBCHXYNGOX-UHFFFAOYSA-K"]…), Dict("original_bigg_ids" => ["cit_c"])), "glc__D_e" => AbstractFBCModels.CanonicalModel.Metabolite("D-Glucose", "e", Dict("C" => 6, "H" => 12, "O" => 6), 0, 0.0, Dict("kegg.drug" => ["D00009"], "kegg.compound" => ["C00031"], "sbo" => ["SBO:0000247"], "biocyc" => ["META:Glucopyranose"], "chebi" => ["CHEBI:12965", "CHEBI:20999", "CHEBI:4167", "CHEBI:17634"], "bigg.metabolite" => ["glc__D"], "seed.compound" => ["cpd26821", "cpd00027"], "sabiork" => ["1406", "1407"], "metanetx.chemical" => ["MNXM41"], "inchi_key" => ["WQZGKKKJIJFFOK-GASJEMHNSA-N"]…), Dict("original_bigg_ids" => ["glc_D_e"])), "nadh_c" => AbstractFBCModels.CanonicalModel.Metabolite("Nicotinamide adenine dinucleotide - reduced", "c", Dict("C" => 21, "N" => 7, "P" => 2, "H" => 27, "O" => 14), -2, 0.0, Dict("kegg.compound" => ["C00004"], "sbo" => ["SBO:0000247"], "biocyc" => ["META:NADH"], "chebi" => ["CHEBI:13395", "CHEBI:21902", "CHEBI:7423", "CHEBI:44216", "CHEBI:57945", "CHEBI:16908", "CHEBI:13396"], "bigg.metabolite" => ["nadh"], "seed.compound" => ["cpd00004"], "sabiork" => ["38"], "metanetx.chemical" => ["MNXM10"], "inchi_key" => ["BOPGDPNILDQYTO-NNYOXOHSSA-L"], "hmdb" => ["HMDB01487"]…), Dict("original_bigg_ids" => ["nadh_c"])), "coa_c" => AbstractFBCModels.CanonicalModel.Metabolite("Coenzyme A", "c", Dict("S" => 1, "C" => 21, "N" => 7, "P" => 3, "H" => 32, "O" => 16), -4, 0.0, Dict("envipath" => ["32de3cf4-e3e6-4168-956e-32fa5ddb0ce1/compound/19310484-6aa5-4dcf-b1da-855a8c21ecfd"], "kegg.compound" => ["C00010"], "sbo" => ["SBO:0000247"], "biocyc" => ["META:COA-GROUP", "META:CO-A"], "chebi" => ["CHEBI:41631", "CHEBI:41597", "CHEBI:741566", "CHEBI:23355", "CHEBI:13295", "CHEBI:13298", "CHEBI:57287", "CHEBI:15346", "CHEBI:3771", "CHEBI:13294"], "bigg.metabolite" => ["coa"], "seed.compound" => ["cpd22528", "cpd00010"], "sabiork" => ["1265"], "metanetx.chemical" => ["MNXM12"], "inchi_key" => ["RGJOEKWQDUBAIZ-IBOSZNHHSA-J"]…), Dict("original_bigg_ids" => ["coa_c"])), "g3p_c" => AbstractFBCModels.CanonicalModel.Metabolite("Glyceraldehyde 3-phosphate", "c", Dict("C" => 3, "P" => 1, "H" => 5, "O" => 6), -2, 0.0, Dict("kegg.compound" => ["C00661", "C00118"], "sbo" => ["SBO:0000247"], "biocyc" => ["META:GAP"], "chebi" => ["CHEBI:5446", "CHEBI:14333", "CHEBI:12984", "CHEBI:181", "CHEBI:17138", "CHEBI:21026", "CHEBI:12983", "CHEBI:18324", "CHEBI:58027", "CHEBI:59776", "CHEBI:29052"], "bigg.metabolite" => ["g3p"], "seed.compound" => ["cpd00102", "cpd19005"], "sabiork" => ["27", "1687"], "metanetx.chemical" => ["MNXM74"], "inchi_key" => ["LXJXRIRHZLFYRP-VKHMYHEASA-L"], "hmdb" => ["HMDB01112"]…), Dict("original_bigg_ids" => ["g3p_c"])), "nad_c" => AbstractFBCModels.CanonicalModel.Metabolite("Nicotinamide adenine dinucleotide", "c", Dict("C" => 21, "N" => 7, "P" => 2, "H" => 26, "O" => 14), -1, 0.0, Dict("kegg.drug" => ["D00002"], "kegg.compound" => ["C00003"], "sbo" => ["SBO:0000247"], "biocyc" => ["META:NAD"], "chebi" => ["CHEBI:7422", "CHEBI:44215", "CHEBI:13394", "CHEBI:21901", "CHEBI:57540", "CHEBI:15846", "CHEBI:44281", "CHEBI:44214", "CHEBI:13393", "CHEBI:29867"], "bigg.metabolite" => ["nad"], "seed.compound" => ["cpd00003"], "sabiork" => ["37"], "metanetx.chemical" => ["MNXM8"], "inchi_key" => ["BAWFJGJZGIEFAR-NNYOXOHSSA-M"]…), Dict("original_bigg_ids" => ["nad_c"]))…), Dict{String, AbstractFBCModels.CanonicalModel.Gene}("b4301" => AbstractFBCModels.CanonicalModel.Gene("sgcE", Dict("sbo" => ["SBO:0000243"], "uniprot" => ["P39362"], "ecogene" => ["EG12553"], "ncbigene" => ["948829"], "ncbigi" => ["16132122"], "refseq_locus_tag" => ["b4301"], "refseq_name" => ["sgcE"], "asap" => ["ABE-0014097"], "refseq_synonym" => ["JW4263", "ECK4290", "yjhK"]), Dict("original_bigg_ids" => ["b4301"])), "b1779" => AbstractFBCModels.CanonicalModel.Gene("gapA", Dict("sbo" => ["SBO:0000243"], "uniprot" => ["P0A9B2"], "ecogene" => ["EG10367"], "ncbigene" => ["947679"], "ncbigi" => ["16129733"], "refseq_locus_tag" => ["b1779"], "refseq_name" => ["gapA"], "asap" => ["ABE-0005920"], "refseq_synonym" => ["ECK1777", "JW1768"]), Dict("original_bigg_ids" => ["b1779"])), "b2276" => AbstractFBCModels.CanonicalModel.Gene("nuoN", Dict("sbo" => ["SBO:0000243"], "uniprot" => ["P0AFF0"], "ecogene" => ["EG12093"], "ncbigene" => ["945136"], "ncbigi" => ["145698289"], "refseq_locus_tag" => ["b2276"], "refseq_name" => ["nuoN"], "asap" => ["ABE-0007526"], "refseq_synonym" => ["JW2271", "ECK2270"]), Dict("original_bigg_ids" => ["b2276"])), "b1761" => AbstractFBCModels.CanonicalModel.Gene("gdhA", Dict("sbo" => ["SBO:0000243"], "uniprot" => ["P00370"], "ecogene" => ["EG10372"], "ncbigene" => ["946802"], "ncbigi" => ["16129715"], "refseq_locus_tag" => ["b1761"], "refseq_name" => ["gdhA"], "asap" => ["ABE-0005865"], "refseq_synonym" => ["ECK1759", "JW1750"]), Dict("original_bigg_ids" => ["b1761"])), "b3493" => AbstractFBCModels.CanonicalModel.Gene("pitA", Dict("sbo" => ["SBO:0000243"], "uniprot" => ["P0AFJ7"], "ecogene" => ["EG12230"], "ncbigene" => ["948009"], "ncbigi" => ["16131365"], "refseq_locus_tag" => ["b3493"], "refseq_name" => ["pitA"], "asap" => ["ABE-0011407"], "refseq_synonym" => ["pit", "ECK3478", "JW3460"]), Dict("original_bigg_ids" => ["b3493"])), "b2926" => AbstractFBCModels.CanonicalModel.Gene("pgk", Dict("sbo" => ["SBO:0000243"], "uniprot" => ["P0A799"], "ecogene" => ["EG10703"], "ncbigene" => ["947414"], "ncbigi" => ["16130827"], "refseq_locus_tag" => ["b2926"], "refseq_name" => ["pgk"], "asap" => ["ABE-0009605"], "refseq_synonym" => ["ECK2922", "JW2893"]), Dict("original_bigg_ids" => ["b2926"])), "b0979" => AbstractFBCModels.CanonicalModel.Gene("cbdB", Dict("sbo" => ["SBO:0000243"], "uniprot" => ["P26458"], "ecogene" => ["EG11379"], "ncbigene" => ["947547"], "ncbigi" => ["16128945"], "refseq_locus_tag" => ["b0979"], "refseq_name" => ["cbdB"], "asap" => ["ABE-0003302"], "refseq_synonym" => ["JW0961", "appB", "cyxB", "ECK0970"]), Dict("original_bigg_ids" => ["b0979"])), "b2282" => AbstractFBCModels.CanonicalModel.Gene("nuoH", Dict("sbo" => ["SBO:0000243"], "uniprot" => ["P0AFD4"], "ecogene" => ["EG12088"], "ncbigene" => ["946761"], "ncbigi" => ["16130217"], "refseq_locus_tag" => ["b2282"], "refseq_name" => ["nuoH"], "asap" => ["ABE-0007541"], "refseq_synonym" => ["JW2277", "ECK2276"]), Dict("original_bigg_ids" => ["b2282"])), "b2283" => AbstractFBCModels.CanonicalModel.Gene("nuoG", Dict("sbo" => ["SBO:0000243"], "uniprot" => ["P33602"], "ecogene" => ["EG12087"], "ncbigene" => ["946762"], "ncbigi" => ["145698290"], "refseq_locus_tag" => ["b2283"], "refseq_name" => ["nuoG"], "asap" => ["ABE-0007543"], "refseq_synonym" => ["JW2278", "ECK2277"]), Dict("original_bigg_ids" => ["b2283"])), "b2416" => AbstractFBCModels.CanonicalModel.Gene("ptsI", Dict("sbo" => ["SBO:0000243"], "uniprot" => ["P08839"], "ecogene" => ["EG10789"], "ncbigene" => ["946879"], "ncbigi" => ["16130342"], "refseq_locus_tag" => ["b2416"], "refseq_name" => ["ptsI"], "asap" => ["ABE-0007967"], "refseq_synonym" => ["ECK2411", "ctr", "JW2409"]), Dict("original_bigg_ids" => ["b2416"]))…), Dict{String, AbstractFBCModels.CanonicalModel.Coupling}()), Dict{String, Dict{String, COBREXA.Isozyme}}("PFK" => Dict("isozyme_1" => COBREXA.Isozyme(Dict("b1723" => 1.0), 3601.6560000000004, nothing)), "PTAr" => Dict("isozyme_2" => COBREXA.Isozyme(Dict("b2458" => 1.0), 4219.092000000001, nothing)), "TKT2" => Dict("isozyme_2" => COBREXA.Isozyme(Dict("b2935" => 1.0), 1682.712, nothing)), "PGL" => Dict("isozyme_1" => COBREXA.Isozyme(Dict("b0767" => 1.0), 7633.512000000001, nothing)), "ACALD" => Dict("isozyme_1" => COBREXA.Isozyme(Dict("b0351" => 1.0), 2045.1960000000001, nothing)), "PGI" => Dict("isozyme_1" => COBREXA.Isozyme(Dict("b4025" => 1.0), 1685.1960000000001, nothing)), "ALCD2x" => Dict("isozyme_3" => COBREXA.Isozyme(Dict("b1478" => 1.0), 273.42, nothing)), "PGK" => Dict("isozyme_1" => COBREXA.Isozyme(Dict("b2926" => 1.0), 207.50400000000002, nothing)), "ICDHyr" => Dict("isozyme_1" => COBREXA.Isozyme(Dict("b1136" => 1.0), 142.632, nothing)), "PDH" => Dict("isozyme_1" => COBREXA.Isozyme(Dict("b0114" => 1.0, "b0115" => 1.0, "b0116" => 1.0), 1907.136, nothing))…))

Since we will later be differentiating our solution, we need to make parameter isozymes and a kinetic model

parameter_values = Dict(
    Symbol(x) => iso.kcat_forward for (x, y) in pruned_reaction_isozymes for (_, iso) in y
)

rid_kcat = Dict(k => Ex(Symbol(k)) for (k, _) in parameter_values)
parameter_isozymes = Dict(
    x => Dict(
        "isozyme" => X.IsozymeT{Ex}(
            iso.gene_product_stoichiometry,
            rid_kcat[Symbol(x)],
            nothing,
        ) for (_, iso) in y
    ) for (x, y) in pruned_reaction_isozymes
)

pkm = X.enzyme_constrained_flux_balance_constraints( # kinetic model
    pruned_model;
    reaction_isozymes = parameter_isozymes,
    gene_product_molar_masses,
    capacity,
)

pruned_solution = D.optimized_values(
    pkm,
    parameter_values;
    objective = pkm.objective.value,
    optimizer = T.Optimizer,
)

(tree = ConstraintTrees.Tree{Float64}(#= 12 elements =#), primal_values = [41.69304582306297, 3.3738864917780025, 0.7991675278074312, 0.7991675278074535, 3.373886491778005, 41.69304582306298, 8.390000014263562, 0.7407243746482269, 53.669001513849246, 0.7991675278074185  …  0.0029929638145053954, 0.019613928519562524, 0.070652071827256, 0.04032874991388847, 0.001959493158657936, 0.15248718390411423, 0.001389701890214975, 0.0011621347499833816, 0.001689819318361655, 0.0203858436174642], equality_dual_values = [0.02492930687014489, 0.014760868846430996, 0.014596725288290495, 0.014965230475163616, 0.015197498948580045, -0.00031185145643069606, 0.0, 0.004922650431641902, 0.00602434591003517, 0.03343676205069087  …  0.0011008155494533034, 0.0027816413336696046, 0.00014303200199665392, 6.0870518424461934e-5, 0.0014429753769459965, 0.0009355543548025853, 0.0015010372952537607, 0.0010786131569188657, 0.0007180292162103088, 0.00017951934251539064], inequality_dual_values = [-1.3851828240284269e-12, -5.011795376759092e-11, -3.709882083649073e-11, -3.874176899575392e-11, -5.072428333116834e-11, -1.4179922937351746e-12, -0.012547259306716914, -4.26844003327565e-11, -6.870164106469038e-13, -3.324292443370856e-11  …  -3.7987563045117077e-14, -3.790555906101654e-14, -9.777251415316644e-14, 0.0, -3.7867994935268234e-14, -3.785540430565886e-14, -3.781840397270636e-14, -4.6375268684640907e-14, -0.016731799175476156, -0.01735886145395185])

All reactions with zero flux have been removed in the pruned model

sort(collect(pruned_solution.tree.fluxes), by = ComposedFunction(abs, last))

52-element Vector{Pair{Symbol, Union{Float64, ConstraintTrees.Tree{Float64}}}}:
                     :GLNS => 0.18940322259736744
 :BIOMASS_Ecoli_core_w_GAM => 0.7407243746482269
                       :CS => 0.7991675278074185
                   :ACONTa => 0.7991675278074312
                   :ACONTb => 0.7991675278074535
                   :ICDHyr => 0.7991675278074873
                     :TKT2 => 1.716838610153485
                     :TALA => 1.9842401094015332
                     :TKT1 => 1.9842401094015338
                      :PPC => 2.1226197679914267
                           ⋮
                 :SUCCt2_2 => 45.10175150771673
                   :SUCCt3 => 45.101751507716735
                      :PDH => 48.642186653954965
                      :ENO => 53.247640453329765
                      :PGM => 53.247640453329765
                     :CO2t => 53.669001513849246
                 :EX_co2_e => 53.669001513849246
                      :PGK => 54.35576411780352
                     :GAPD => 54.35576411780354

Genes with zero concentration have also been removed.

sort(abs.(collect(values(pruned_solution.tree.gene_product_amounts))))

49-element Vector{Float64}:
 0.0005800203793505351
 0.0005856841395960682
 0.0007996712306205315
 0.0008318932491454554
 0.001689819318361655
 0.001959493158657936
 0.0025518366401983566
 0.0027916453170659383
 0.0029929638145053954
 0.0032366263256029797
 ⋮
 0.034702381966990506
 0.034702381966990506
 0.04032874991388847
 0.04032874991388847
 0.05349843841361214
 0.070652071827256
 0.1172633066639992
 0.15248718390411423
 0.26195044007731577

Optimal flux modes (OFMs)

We now wish to find the optimal flux modes (OFMs) of this optimal solution.

In our model, we have a fixed ATP maintenance flux

pruned_model.reactions["ATPM"].lower_bound

8.39

As a result of this fixed flux, standard EFM theory does not hold. Instead, we investigate the optimal flux modes (OFMs). These are the flux modes that are able to produce the optimal rate of the objective, whilst still being constrained to the fixed ATPM.

To find the OFMs, we must augment the stoichiometric matrix, and can then use the same find_efms function, since this augmented matrix behaves as required for the algorithm to work.

N = A.stoichiometry(pruned_model)

55×52 SparseArrays.SparseMatrixCSC{Float64, Int64} with 205 stored entries:
⎡⠀⠀⠀⠠⠀⠐⠀⠀⠀⠀⠀⠀⠁⠀⡀⠀⠀⠀⠨⡰⠀⠀⠀⠀⠀⠀⎤
⎢⡐⠀⢆⠀⠀⠀⠠⠀⠀⠀⠀⠈⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⠀⠀⠀⠀⎥
⎢⢡⠒⠀⣈⠈⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠀⢁⢀⠀⠀⢅⠀⠀⠀⠀⎥
⎢⠐⠄⠀⠚⡠⠀⠀⠀⠀⠀⠀⠀⠀⠊⡀⡁⠀⡐⠐⠀⢀⠐⠀⠀⠀⠀⎥
⎢⠂⠀⠀⢐⠑⠀⠀⠁⠀⠀⠀⠄⠀⠀⠀⠀⠀⠂⠀⠀⠀⠂⠀⠀⡀⡠⎥
⎢⠀⠀⠈⠠⠀⠀⠃⠐⠀⠀⠀⡀⠀⠀⠀⠀⠀⢠⠄⠀⠀⠀⠀⠀⠄⠄⎥
⎢⠀⠀⠀⢘⠀⠀⠀⠀⠄⠀⠀⠑⠱⡀⠀⠀⠀⠀⠂⠀⠀⠀⠀⠀⠉⠉⎥
⎢⡀⠒⣀⣚⢐⡒⡀⠀⠠⠀⠀⢀⡀⣙⠰⢀⠀⢀⠀⡂⣐⢀⠀⣀⠀⠀⎥
⎢⡄⠐⢡⢠⠀⠁⠁⠀⠀⠁⠀⠀⡄⠀⠀⢪⠀⡄⠀⠀⠁⠀⠀⠉⠀⠀⎥
⎢⠀⠀⠀⠘⠀⠀⠀⠀⠀⢀⠀⠘⠀⠼⠃⠃⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⢠⠠⢁⠀⠀⠀⠀⠂⠀⢀⠀⠀⠀⠘⠀⠀⠀⢠⢀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠩⠀⡀⠀⠀⠀⠀⠐⠀⠡⠁⠀⢀⠀⠄⠀⠀⠋⠡⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠐⠀⠁⠀⠀⠀⠀⠀⠀⠀⠀⠄⠈⠀⠀⠀⠀⠀⠀⠴⠀⣐⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠄⠛⠠⠄⎦

Take the fixed fluxes (here only ATPM) and the objective flux

atpm_idx = findfirst(x -> x == "ATPM", A.reactions(pruned_model))
biomass_idx = findfirst(x -> x == "BIOMASS_Ecoli_core_w_GAM", A.reactions(pruned_model))
fixed_fluxes = [atpm_idx, biomass_idx]
flux_values = [
    pruned_solution.tree.fluxes["ATPM"],
    pruned_solution.tree.fluxes["BIOMASS_Ecoli_core_w_GAM"],
]

2-element Vector{Float64}:
 8.390000014263562
 0.7407243746482269

Calculate the OFMs, using the get_ofms function

OFMs = get_ofms(Matrix(N), fixed_fluxes, flux_values)
OFM_dicts = [
    Dict(A.reactions(pruned_model) .=> OFMs[:, 1]),
    Dict(A.reactions(pruned_model) .=> OFMs[:, 2]),
]

2-element Vector{Dict{String, Float64}}:
 Dict("NH4t" => 4.039021870081851, "ACALD" => 48.44081880661824, "PTAr" => 0.0, "PGL" => 6.350267100077981, "ALCD2x" => 48.44081880661824, "PGK" => 57.72965060958195, "PDH" => 52.01607314573287, "CO2t" => 57.04288800562686, "PYK" => 23.206660196868704, "EX_nh4_e" => 4.039021870081851…)
 Dict("NH4t" => 4.039021870081861, "ACALD" => 0.0, "PTAr" => 24.22040940330914, "PGL" => 6.3502671000779864, "ALCD2x" => 0.0, "PGK" => 33.50924120627286, "PDH" => 27.795663742423752, "CO2t" => 32.82247860231776, "PYK" => 11.096455495214123, "EX_nh4_e" => 4.039021870081861…)

Scale the OFMs to have one unit flux through biomass

OFM_dicts = [
    Dict(x => y / OFM_dicts[1]["BIOMASS_Ecoli_core_w_GAM"] for (x, y) in OFM_dicts[1]),
    Dict(x => y / OFM_dicts[2]["BIOMASS_Ecoli_core_w_GAM"] for (x, y) in OFM_dicts[2]),
]

2-element Vector{Dict{String, Float64}}:
 Dict("ACALD" => 65.39655027502367, "PTAr" => 0.0, "ALCD2x" => 65.39655027502367, "PDH" => 70.22325027502372, "CO2t" => 77.0096002750237, "PYK" => 31.3296834708452, "EX_nh4_e" => 5.452799999999999, "CS" => 1.0789000000000004, "PGM" => 76.44075027502372, "TKT1" => 2.678783333333333…)
 Dict("ACALD" => 0.0, "PTAr" => 32.698275137511864, "ALCD2x" => 0.0, "PDH" => 37.52497513751188, "CO2t" => 44.3113251375119, "PYK" => 14.980545902089256, "EX_nh4_e" => 5.452800000000012, "CS" => 1.0788999999999962, "PGM" => 43.74247513751189, "TKT1" => 2.678783333333335…)

We see that one EFM is releasing ethanol, and the other acetate

OFM_dicts[1]["EX_etoh_e"]
OFM_dicts[1]["EX_ac_e"]


OFM_dicts[2]["EX_etoh_e"]
OFM_dicts[2]["EX_ac_e"]

32.698275137511864

Differentiate OFM usage

The weighted sum of these two OFMs is equal to the whole optimal flux solution, and can be easily calculated using two non-shared reactions

M = [
    OFM_dicts[1]["EX_etoh_e"] OFM_dicts[2]["EX_etoh_e"]
    OFM_dicts[1]["EX_ac_e"] OFM_dicts[2]["EX_ac_e"]
]

v = [
    pruned_solution.tree.fluxes["EX_etoh_e"]
    pruned_solution.tree.fluxes["EX_ac_e"]
]

λ = M \ v

2-element Vector{Float64}:
 0.6375419750388028
 0.10318239960943508

λ tells us the weighting of the two OFMs, so here we have 0.6 units of flux through OFM₁ and 0.1 units of flux through OFM₂ to produce our optimal biomass.

It is of interest to see how changes in the kinetic parameters affect this optimal weighting.

parameters = Ex.(collect(keys(parameter_values)))
rid_gcounts = Dict(
    rid => [v.gene_product_stoichiometry for (k, v) in d][1] for
    (rid, d) in pruned_reaction_isozymes
)
rid_pid =
    Dict(rid => [iso.kcat_forward for (k, iso) in v][1] for (rid, v) in parameter_isozymes)
sens = differentiate_ofm(
    OFM_dicts,
    parameters,
    rid_pid,
    parameter_values,
    rid_gcounts,
    capacity,
    gene_product_molar_masses,
    T.Optimizer,
)

2×33 Matrix{Float64}:
  2.34957e-6  -3.24672e-6   3.25167e-5   8.97662e-7   2.81279e-6   1.09438e-6  …   5.32262e-6   1.64894e-6   8.33172e-6  -4.56087e-5   0.000106976
 -9.72454e-7   6.36629e-6  -1.34582e-5  -3.71529e-7  -1.16417e-6  -4.52948e-7     -2.20296e-6  -6.82471e-7  -3.44838e-6   8.94313e-5  -4.42758e-5

To get control coefficients instead of sensititivities, we scale these derivatives.

scaled_sens = Matrix(undef, size(sens, 1), size(sens, 2))
for (i, col) in enumerate(eachcol(sens))
    scaled_sens[:, i] = collect(values(parameter_values))[i] .* col ./ λ
end

We plot the control coefficients to get a visual overview of the system.

sens_perm = sortperm(scaled_sens[1, :])
scaled_sens[:, sens_perm]
parameters[sens_perm]


f, a, hm = heatmap(
    scaled_sens[:, sens_perm]';
    colormap = CairoMakie.Reverse(:RdBu),
    axis = (
        yticks = (1:2, ["Ethanol producing", "Acetate producing"]),
        yticklabelrotation = pi / 2,
        xticklabelrotation = pi / 2,
        xlabel = "Parameters",
        ylabel = "Control coefficient, param/λ * ∂λ/∂param",
        xticks = (1:length(parameters), string.(parameters[sens_perm])),
        title = "Control coefficients of OFM weightings in optimal solution",
        ylabelsize = 23,
        xlabelsize = 23,
        titlesize = 25,
    ),
);
Colorbar(f[:, end+1], hm)
f

We see that most kinetic parameters have little effect on the optimal OFM weightings, see the reactions from RPE to PGI. Those parameters that do affect optimal weightings always increase the use of one OFM and decrease the use of the other.

Increasing the turnover number of PGK is predicted to decrease the optimal flux through the acetate producing OFM, and increase the optimal flux through the ethanol producing OFM.

This page was generated using Literate.jl.