diff --git a/Project.toml b/Project.toml index 8d802fe5d90e99641433993a40684a486eed5269..c362e0b11d8d4f1995dd2130efbe3e907d79d92a 100644 --- a/Project.toml +++ b/Project.toml @@ -45,6 +45,7 @@ OSQP = "ab2f91bb-94b4-55e3-9ba0-7f65df51de79" SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa" +GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6" [targets] -test = ["Aqua", "Downloads", "OSQP", "SHA", "Test", "Tulip"] +test = ["Aqua", "Downloads", "OSQP", "SHA", "Test", "Tulip", "GLPK"] diff --git a/src/analysis/modifications/crowding.jl b/src/analysis/modifications/crowding.jl index 6e0ff7ef38542049bd9b04cf895a77b221918be0..76213583911faa24aaae3280055dfa406e65d29b 100644 --- a/src/analysis/modifications/crowding.jl +++ b/src/analysis/modifications/crowding.jl @@ -22,22 +22,6 @@ add_crowding_constraint(weights::Dict{Int64,Float64}) = @constraint(opt_model, sum(w * crowding_slack[i] for (i, w) in enumerate(ws)) <= 1) end -""" - add_crowding_constraint(weight::Float64; kwargs) - -Variant of [`add_crowding_constraint`](@ref) that takes a single weight and assigns it to -each internal reaction flux, where internal reactions are identified with -[`find_internal_reactions`](@ref) and `kwargs` are passed to this function. -""" -add_crowding_constraint(weight::Float64; kwargs...) = - (model, opt_model) -> begin - idxs = find_internal_reactions(model; kwargs...) - add_crowding_constraint(Dict(zip(idxs, fill(weight, length(idxs)))))( - model, - opt_model, - ) - end - """ add_crowding_constraint(weights::Dict{String, Float64}) diff --git a/src/analysis/modifications/loopless.jl b/src/analysis/modifications/loopless.jl new file mode 100644 index 0000000000000000000000000000000000000000..aaf035453df029418855659bb3ab7b24554c5c5e --- /dev/null +++ b/src/analysis/modifications/loopless.jl @@ -0,0 +1,60 @@ +""" + loopless(; + max_flux_bound = _constants.default_reaction_bound, + strict_inequality_tolerance = _constants.loopless_strict_inequality_tolerance, + ) + +Add quasi-thermodynamic constraints to the model to ensure that no thermodynamically +infeasible internal cycles can occur. Adds the following constraints to the problem: +``` +-max_flux_bound × (1 - yᵢ) ≤ xᵢ ≤ max_flux_bound × yᵢ +-max_flux_bound × yᵢ + strict_inequality_tolerance × (1 - yᵢ) ≤ Gᵢ +Gᵢ ≤ -strict_inequality_tolerance × yᵢ + max_flux_bound × (1 - yᵢ) +Nᵢₙₜ' × G = 0 +yᵢ ∈ {0, 1} +Gᵢ ∈ ℠+i ∈ internal reactions +Nᵢₙₜ is the nullspace of the internal stoichiometric matrix +``` +Note, this modification introduces binary variables, so an optimization solver capable of +handing mixed integer problems needs to be used. The arguments `max_flux_bound` and +`strict_inequality_tolerance` implement the "big-M" method of indicator constraints. + +For more details about the algorithm, see `Schellenberger, Lewis, and, Palsson. "Elimination +of thermodynamically infeasible loops in steady-state metabolic models.", Biophysical +journal, 2011`. +""" +loopless(; + max_flux_bound = _constants.default_reaction_bound, # needs to be an order of magnitude bigger, big M method heuristic + strict_inequality_tolerance = _constants.loopless_strict_inequality_tolerance, +) = + (model, opt_model) -> begin + + internal_rxn_idxs = [ + ridx for (ridx, rid) in enumerate(reactions(model)) if + !is_boundary(reaction_stoichiometry(model, rid)) + ] + + N_int = nullspace(Array(stoichiometry(model)[:, internal_rxn_idxs])) # no sparse nullspace function + + y = @variable(opt_model, y[1:length(internal_rxn_idxs)], Bin) + G = @variable(opt_model, G[1:length(internal_rxn_idxs)]) # approx ΔG for internal reactions + + x = opt_model[:x] + for (cidx, ridx) in enumerate(internal_rxn_idxs) + @constraint(opt_model, -max_flux_bound * (1 - y[cidx]) <= x[ridx]) + @constraint(opt_model, x[ridx] <= max_flux_bound * y[cidx]) + + @constraint( + opt_model, + -max_flux_bound * y[cidx] + strict_inequality_tolerance * (1 - y[cidx]) <= G[cidx] + ) + @constraint( + opt_model, + G[cidx] <= + -strict_inequality_tolerance * y[cidx] + max_flux_bound * (1 - y[cidx]) + ) + end + + @constraint(opt_model, N_int' * G .== 0) + end diff --git a/src/base/constants.jl b/src/base/constants.jl index 5379f4c9c76ec501baab59fb7301ff4164895882..864a688479e36e2eee56eea6b8a2e4cf6a6569b6 100644 --- a/src/base/constants.jl +++ b/src/base/constants.jl @@ -9,7 +9,8 @@ const _constants = ( tolerance = 1e-6, sampling_keep_iters = 100, sampling_size = 1000, - exchange_suffixes = ["_e", "[e]", "(e)"], + loopless_strict_inequality_tolerance = 1, + extracellular_suffixes = ["_e", "[e]", "(e)"], exchange_prefixes = ["EX_", "Exch_", "Ex_", "R_EX_", "R_Ex", "R_Exch_"], biomass_strings = ["BIOMASS", "biomass", "Biomass"], keynames = ( diff --git a/src/base/utils/Reaction.jl b/src/base/utils/Reaction.jl index 8009acbc7b06f8a159a6b6ab86c07ca9c52df598..809e2638756372d6f122ada0fd9d9e68db0c3e70 100644 --- a/src/base/utils/Reaction.jl +++ b/src/base/utils/Reaction.jl @@ -39,17 +39,20 @@ function check_duplicate_reaction( end """ - is_boundary(reaction::Reaction) + is_boundary(rxn_dict::Dict{String, Float64}) -Return true if reaction is a boundary reaction, otherwise return false. -Checks if on boundary by inspecting the number of metabolites in reaction equation. +Return true if the reaction denoted by `rxn_dict` is a boundary reaction, otherwise return false. +Checks if on boundary by inspecting the number of metabolites in `rxn_dict`. Boundary reactions have only one metabolite, e.g. an exchange reaction, or a sink/demand reaction. """ -function is_boundary(rxn::Reaction)::Bool - length(keys(rxn.metabolites)) == 1 ? true : false +function is_boundary(rxn_dict::Dict{String,Float64})::Bool + length(keys(rxn_dict)) == 1 end -is_boundary(model::StandardModel, rxn_id::String) = is_boundary(model.reactions[rxn_id]) +is_boundary(model::MetabolicModel, rxn_id::String) = + is_boundary(reaction_stoichiometry(model, rxn_id)) + +is_boundary(rxn::Reaction) = is_boundary(rxn.metabolites) is_boundary(model::StandardModel, rxn::Reaction) = is_boundary(rxn) # for consistency with functions below diff --git a/src/base/utils/looks_like.jl b/src/base/utils/looks_like.jl index dadccfa323a2a2a95f784a085962e10e5375396f..b3f05577b1223876d3b8e1a22344fb1e95e06eba 100644 --- a/src/base/utils/looks_like.jl +++ b/src/base/utils/looks_like.jl @@ -99,101 +99,40 @@ find_biomass_reaction_ids(m::MetabolicModel; kwargs...) = filter(id -> looks_like_biomass_reaction(id, kwargs...), reactions(m)) """ - looks_like_internal_reaction( - rxn_id::String; - exchange_prefixes = _constants.exchange_prefixes, - biomass_strings = _constants.biomass_strings, - )::Bool - -A predicate that matches reaction identifiers that look like -internal reactions, i.e. reactions that are neither exchange nor biomass reactions. -Exchange reactions are identified based on matching prefixes in -the set `exchange_prefixes` and biomass reactions are identified by looking for -occurences of `biomass_strings` in the reaction id. - -Also see [`find_internal_reactions`](@ref). - -# Example -``` -findall(looks_like_internal_reaction, reactions(model)) # returns indices -filter(looks_like_internal_reaction, reactions(model)) # returns Strings - -# to use the optional arguments you need to expand the function's arguments -# using an anonymous function -findall(x -> looks_like_internal_reaction(x; exchange_prefixes=["EX_", "R_EX_"]), reactions(model)) # returns indices -filter(x -> looks_like_internal_reaction(x; exchange_prefixes=["EX_", "R_EX_"]), reactions(model)) # returns Strings -``` -""" -function looks_like_internal_reaction( - rxn_id::String; - exchange_prefixes = _constants.exchange_prefixes, - biomass_strings = _constants.biomass_strings, -)::Bool - !looks_like_biomass_reaction( - rxn_id; - exchange_prefixes = exchange_prefixes, - biomass_strings = biomass_strings, - ) && - !looks_like_exchange_reaction( - rxn_id; - exchange_prefixes = exchange_prefixes, - biomass_strings = biomass_strings, - ) -end - -""" - find_internal_reactions(m::MetabolicModel; kwargs...) - -Shortcut for finding internal reaction indices in a model; arguments are -forwarded to [`looks_like_internal_reaction`](@ref). -""" -find_internal_reactions(m::MetabolicModel; kwargs...) = - findall(id -> looks_like_internal_reaction(id; kwargs...), reactions(m)) - -""" - find_internal_reaction_ids(m::MetabolicModel; kwargs...) - -Shortcut for finding internal reaction identifiers in a model; arguments are -forwarded to [`looks_like_internal_reaction`](@ref). -""" -find_internal_reaction_ids(m::MetabolicModel; kwargs...) = - filter(id -> looks_like_internal_reaction(id, kwargs...), reactions(m)) - -""" - looks_like_exchange_metabolite(rxn_id::String; - exchange_suffixes = _constants.exchange_suffixes, + looks_like_extracellular_metabolite(rxn_id::String; + extracellular_suffixes = _constants.extracellular_suffixes, )::Bool -A predicate that matches metabolite identifiers that look like involved in -exchange reactions. Exchange metabolites are identified by `exchange_suffixes` -at the end of the metabolite id. +A predicate that matches metabolite identifiers that look like they are extracellular +metabolites. Extracellular metabolites are identified by `extracellular_suffixes` at the end of the +metabolite id. # Example ``` -filter(looks_like_exchange_metabolite, metabolites(model)) # returns strings -findall(looks_like_exchange_metabolite, metabolites(model)) # returns indices +filter(looks_like_extracellular_metabolite, metabolites(model)) # returns strings +findall(looks_like_extracellular_metabolite, metabolites(model)) # returns indices ``` """ -function looks_like_exchange_metabolite( +function looks_like_extracellular_metabolite( met_id::String; - exchange_suffixes = _constants.exchange_suffixes, + extracellular_suffixes = _constants.extracellular_suffixes, )::Bool - any(endswith(met_id, x) for x in exchange_suffixes) + any(endswith(met_id, x) for x in extracellular_suffixes) end """ - find_exchange_metabolites(m::MetabolicModel; kwargs...) + find_extracellular_metabolites(m::MetabolicModel; kwargs...) -Shortcut for finding exchange metabolite indexes in a model; arguments are -forwarded to [`looks_like_exchange_metabolite`](@ref). +Shortcut for finding extracellular metabolite indexes in a model; arguments are +forwarded to [`looks_like_extracellular_metabolite`](@ref). """ -find_exchange_metabolites(m::MetabolicModel; kwargs...) = - findall(id -> looks_like_exchange_metabolite(id; kwargs...), metabolites(m)) +find_extracellular_metabolites(m::MetabolicModel; kwargs...) = + findall(id -> looks_like_extracellular_metabolite(id; kwargs...), metabolites(m)) """ - find_exchange_metabolite_ids(m::MetabolicModel; kwargs...) + find_extracellular_metabolite_ids(m::MetabolicModel; kwargs...) -Shortcut for finding exchange metabolite identifiers in a model; arguments are -forwarded to [`looks_like_exchange_metabolite`](@ref). +Shortcut for finding extracellular metabolite identifiers in a model; arguments are +forwarded to [`looks_like_extracellular_metabolite`](@ref). """ -find_exchange_metabolite_ids(m::MetabolicModel; kwargs...) = - findall(id -> looks_like_exchange_metabolite(id; kwargs...), metabolites(m)) +find_extracellular_metabolite_ids(m::MetabolicModel; kwargs...) = + findall(id -> looks_like_extracellular_metabolite(id; kwargs...), metabolites(m)) diff --git a/test/analysis/fba_variants.jl b/test/analysis/fba_with_crowding.jl similarity index 72% rename from test/analysis/fba_variants.jl rename to test/analysis/fba_with_crowding.jl index 22c90ee10c4f0149aba1ceb3ac8a7206ff1855db..0eca14f571da61e7b55f6fa19abc4a7c17db0d6c 100644 --- a/test/analysis/fba_variants.jl +++ b/test/analysis/fba_with_crowding.jl @@ -1,13 +1,16 @@ @testset "FBA with crowding constraints" begin model = load_model(model_paths["e_coli_core.json"]) - idxs = find_internal_reactions(model) + rid_weight = Dict( + rid => 0.004 for rid in reactions(model) if + !looks_like_biomass_reaction(rid) && !looks_like_exchange_reaction(rid) + ) sol = flux_balance_analysis_dict( model, Tulip.Optimizer; modifications = [ change_optimizer_attribute("IPM_IterationsLimit", 1000), - add_crowding_constraint(0.004), + add_crowding_constraint(rid_weight), change_constraint("EX_glc__D_e"; lb = -6), ], ) diff --git a/test/analysis/loopless.jl b/test/analysis/loopless.jl new file mode 100644 index 0000000000000000000000000000000000000000..5e1532fd9fa01029f66a4cd6d680b0f05727a1c1 --- /dev/null +++ b/test/analysis/loopless.jl @@ -0,0 +1,12 @@ +@testset "Loopless FBA" begin + + model = load_model(model_paths["e_coli_core.json"]) + + sol = flux_balance_analysis_dict(model, GLPK.Optimizer; modifications = [loopless()]) + + @test isapprox( + sol["BIOMASS_Ecoli_core_w_GAM"], + 0.8739215069684292, + atol = TEST_TOLERANCE, + ) +end diff --git a/test/base/utils/looks_like.jl b/test/base/utils/looks_like.jl index 7e56bf5ec45678de6a64e15577877f307d53745c..c09966d330760d43f2f3fd087e74a15c934e5fda 100644 --- a/test/base/utils/looks_like.jl +++ b/test/base/utils/looks_like.jl @@ -15,8 +15,6 @@ ["m1"; "m2"; "m3"], ) @test filter(looks_like_exchange_reaction, reactions(cp)) == ["EX_m1"] - @test find_internal_reactions(cp) == [2, 3] - @test find_internal_reaction_ids(cp) == ["r2", "r3"] cp = CoreModel( [-1.0 0 0; 0 0 -1; 0 -1 0], @@ -43,7 +41,7 @@ x -> looks_like_exchange_reaction(x; exclude_biomass = true), reactions(cp), ) == ["EX_m1(e)", "EX_m3(e)"] - @test filter(looks_like_exchange_metabolite, metabolites(cp)) == ["m1[e]", "m3[e]"] + @test filter(looks_like_extracellular_metabolite, metabolites(cp)) == ["m1[e]", "m3[e]"] @test filter(looks_like_biomass_reaction, reactions(cp)) == ["EX_biomass(e)", "biomass1"] @test filter( @@ -55,32 +53,26 @@ end @testset "Looks like functions, basic" begin model = load_model(model_paths["e_coli_core.json"]) @test length(filter(looks_like_exchange_reaction, reactions(model))) == 20 - @test length(filter(looks_like_exchange_metabolite, metabolites(model))) == 20 @test length(filter(looks_like_biomass_reaction, reactions(model))) == 1 - @test length(filter(looks_like_internal_reaction, reactions(model))) == 74 + @test length(filter(looks_like_extracellular_metabolite, metabolites(model))) == 20 model = load_model(model_paths["e_coli_core.xml"]) @test length(filter(looks_like_exchange_reaction, reactions(model))) == 20 - @test length(filter(looks_like_exchange_metabolite, metabolites(model))) == 20 @test length(filter(looks_like_biomass_reaction, reactions(model))) == 1 - @test length(filter(looks_like_internal_reaction, reactions(model))) == 74 + @test length(filter(looks_like_extracellular_metabolite, metabolites(model))) == 20 model = load_model(model_paths["e_coli_core.mat"]) @test length(filter(looks_like_exchange_reaction, reactions(model))) == 20 - @test length(filter(looks_like_exchange_metabolite, metabolites(model))) == 20 @test length(filter(looks_like_biomass_reaction, reactions(model))) == 1 - @test length(filter(looks_like_internal_reaction, reactions(model))) == 74 + @test length(filter(looks_like_extracellular_metabolite, metabolites(model))) == 20 model = convert(StandardModel, model) @test length(filter(looks_like_exchange_reaction, reactions(model))) == 20 - @test length(filter(looks_like_exchange_metabolite, metabolites(model))) == 20 @test length(filter(looks_like_biomass_reaction, reactions(model))) == 1 - @test length(filter(looks_like_internal_reaction, reactions(model))) == 74 + @test length(filter(looks_like_extracellular_metabolite, metabolites(model))) == 20 model = convert(CoreModelCoupled, model) @test length(filter(looks_like_exchange_reaction, reactions(model))) == 20 - @test length(filter(looks_like_exchange_metabolite, metabolites(model))) == 20 @test length(filter(looks_like_biomass_reaction, reactions(model))) == 1 - @test length(filter(looks_like_internal_reaction, reactions(model))) == 74 - + @test length(filter(looks_like_extracellular_metabolite, metabolites(model))) == 20 end diff --git a/test/runtests.jl b/test/runtests.jl index e9a0f20691a7ff871070c93f56fd026cbd64e74f..847c7f35e63f1f90cdac9ff15c25c957dd479a20 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,6 +13,7 @@ using SHA using SparseArrays using Statistics using Tulip +using GLPK # for MILPs # tolerance for comparing analysis results (should be a bit bigger than the # error tolerance in computations)