diff --git a/docs/src/examples/06-thermodynamic-models.jl b/docs/src/examples/06-thermodynamic-models.jl index d7b67bb9f..1ac8afd0f 100644 --- a/docs/src/examples/06-thermodynamic-models.jl +++ b/docs/src/examples/06-thermodynamic-models.jl @@ -112,10 +112,11 @@ m = build_max_min_driving_force_model( R = 8.31446261815324e-3, # kJ/K/mol ) -m *= :metabolite_ratio_constraints^log_ratio_constraints( - Dict("atp" => ("atp_c", "adp_c", 10.0), "nadh" => ("nadh_c", "nad_c", 0.13)), - m.log_metabolite_concentrations, -) +m *= + :metabolite_ratio_constraints^log_ratio_constraints( + Dict("atp" => ("atp_c", "adp_c", 10.0), "nadh" => ("nadh_c", "nad_c", 0.13)), + m.log_metabolite_concentrations, + ) # solve the model mmdf_solution = optimized_constraints( diff --git a/docs/src/examples/07-loopless-models.jl b/docs/src/examples/07-loopless-models.jl index 6f515a400..1dd71143b 100644 --- a/docs/src/examples/07-loopless-models.jl +++ b/docs/src/examples/07-loopless-models.jl @@ -33,10 +33,14 @@ sol = loopless_flux_balance_analysis( optimizer = GLPK.Optimizer, max_flux_bound = 1000.0, # needs to be an order of magnitude bigger, big M method heuristic strict_inequality_tolerance = 1.0, # heuristic from paper - modifications = [set_optimizer_attribute("tol_int", 1e-9)] + modifications = [set_optimizer_attribute("tol_int", 1e-9)], ) sol.pseudo_gibbs_free_energy_reaction -@test isapprox(mmdf_solution.max_min_driving_force, 0.8739215069684292, atol = TEST_TOLERANCE) #src +@test isapprox( + mmdf_solution.max_min_driving_force, + 0.8739215069684292, + atol = TEST_TOLERANCE, +) #src diff --git a/docs/src/examples/08-community-models.jl b/docs/src/examples/08-community-models.jl index 0b434d2aa..1594d71e6 100644 --- a/docs/src/examples/08-community-models.jl +++ b/docs/src/examples/08-community-models.jl @@ -23,3 +23,51 @@ model = load_model("e_coli_core.json") # numbers, molar masses of enzymes, and capacity bounds. import AbstractFBCModels as A + +m1 = fbc_model_constraints(model) +m2 = fbc_model_constraints(model) + +lbs, ubs = A.bounds(model) +env_ex_rxns = Dict(rid => (lbs[i], ubs[i]) for (i, rid) in enumerate(A.reactions(model)) if startswith(rid, "EX_")) + + +m = build_community_environment(env_ex_rxns) +m += :bug1^m1 +m += :bug2^m2 + +member_abundances = [(:bug1, 0.2), (:bug2, 0.8)] +m *= + :environmental_exchange_balances^link_environmental_exchanges( + m, + [(:bug1, 0.2), (:bug2, 0.8)], + ) + +m *= + :equal_growth_rate_constraint^equal_growth_rate_constraints( + [(:bug1, m.bug1.fluxes.:BIOMASS_Ecoli_core_w_GAM.value), (:bug2, m.bug2.fluxes.:BIOMASS_Ecoli_core_w_GAM.value)] + ) + +m.bug1.fluxes.EX_glc__D_e.bound = (-1000.0, 1000.0) +m.bug2.fluxes.EX_glc__D_e.bound = (-1000.0, 1000.0) +m.bug1.fluxes.CYTBD.bound = (-10.0, 10.0) # respiration limited + +m *= :objective^C.Constraint(m.bug1.fluxes.:BIOMASS_Ecoli_core_w_GAM.value) + +using Gurobi + +sol = optimized_constraints(m; objective = m.objective.value, optimizer=Gurobi.Optimizer) + +Dict(k => v for (k, v) in sol.bug1.fluxes if startswith(string(k), "EX_")) +Dict(k => v for (k, v) in sol.bug2.fluxes if startswith(string(k), "EX_")) + +# exchange cytosolic metabolites +mets = [:EX_akg_e, :EX_succ_e, :EX_pyr_e, :EX_acald_e, :EX_fum_e, :EX_mal__L_e] +for met in mets + m.bug1.fluxes[met].bound = (-1000.0, 1000.0) + m.bug2.fluxes[met].bound = (-1000.0, 1000.0) +end + +sol = optimized_constraints(m; objective = m.objective.value, optimizer=Tulip.Optimizer, modifications=[set_optimizer_attribute("IPM_IterationsLimit", 100000)]) + +Dict(k => v for (k, v) in sol.bug1.fluxes if startswith(string(k), "EX_")) +Dict(k => v for (k, v) in sol.bug2.fluxes if startswith(string(k), "EX_")) diff --git a/src/COBREXA.jl b/src/COBREXA.jl index c9d73692d..ce9bd2493 100644 --- a/src/COBREXA.jl +++ b/src/COBREXA.jl @@ -37,6 +37,7 @@ include("builders/objectives.jl") include("builders/enzymes.jl") include("builders/thermodynamic.jl") include("builders/loopless.jl") +include("builders/communities.jl") include("analysis/modifications.jl") include("analysis/flux_balance.jl") diff --git a/src/analysis/loopless_flux_balance.jl b/src/analysis/loopless_flux_balance.jl index e69de29bb..8b1378917 100644 --- a/src/analysis/loopless_flux_balance.jl +++ b/src/analysis/loopless_flux_balance.jl @@ -0,0 +1 @@ + diff --git a/src/builders/communities.jl b/src/builders/communities.jl new file mode 100644 index 000000000..686fe2797 --- /dev/null +++ b/src/builders/communities.jl @@ -0,0 +1,51 @@ +""" +$(TYPEDSIGNATURES) + +Helper function to create environmental exchange rections. +""" +function environment_exchange_variables(env_ex_rxns = Dict{String,Tuple{Float64,Float64}}()) + rids = collect(keys(env_ex_rxns)) + lbs_ubs = collect(values(env_ex_rxns)) + C.variables(; keys = Symbol.(rids), bounds = lbs_ubs) +end + +export environment_exchange_variables + +function build_community_environment(env_ex_rxns = Dict{String,Tuple{Float64,Float64}}()) + C.ConstraintTree( + :environmental_exchange_reactions => environment_exchange_variables(env_ex_rxns), + ) +end + +export build_community_environment + +function link_environmental_exchanges( + m::C.ConstraintTree, + member_abundances::Vector{Tuple{Symbol,Float64}}; + on = m.:environmental_exchange_reactions, + member_fluxes_id = :fluxes, +) + C.ConstraintTree( + rid => C.Constraint( + value = -rxn.value + sum( + abundance * m[member][member_fluxes_id][rid].value for + (member, abundance) in member_abundances if + haskey(m[member][member_fluxes_id], rid); + init = zero(C.LinearValue), + ), + bound = 0.0, + ) for (rid, rxn) in on + ) +end + +export link_environmental_exchanges + +function equal_growth_rate_constraints(member_biomasses::Vector{Tuple{Symbol,C.LinearValue}}) + C.ConstraintTree( + Symbol(bid1, :_, bid2) => C.Constraint(value = bval1 - bval2, bound = 0.0) for + ((bid1, bval1), (bid2, bval2)) in + zip(member_biomasses[1:end-1], member_biomasses[2:end]) + ) +end + +export equal_growth_rate_constraints diff --git a/src/builders/loopless.jl b/src/builders/loopless.jl index a9d0a4112..aa0383053 100644 --- a/src/builders/loopless.jl +++ b/src/builders/loopless.jl @@ -35,65 +35,84 @@ function loopless_flux_balance_analysis( model; max_flux_bound = 1000.0, # needs to be an order of magnitude bigger, big M method heuristic strict_inequality_tolerance = 10.0, # heuristic from paper - modifications=[], + modifications = [], optimizer, ) - - m = fbc_model_constraints(model) - # find all internal reactions - internal_reactions = [(i, Symbol(rid)) for (i, rid) in enumerate(A.reactions(model)) if !is_boundary(model, rid)] - internal_reaction_ids = last.(internal_reactions) - internal_reaction_idxs = first.(internal_reactions) + m = fbc_model_constraints(model) - # add loopless variables: flux direction setters (binary) and pseudo gibbs free energy of reaction - m += :loopless_binary_variables^C.variables(keys=internal_reaction_ids, bounds=Ref(C.Binary)) - m += :pseudo_gibbs_free_energy_reaction^C.variables(keys=internal_reaction_ids) + # find all internal reactions + internal_reactions = [ + (i, Symbol(rid)) for + (i, rid) in enumerate(A.reactions(model)) if !is_boundary(model, rid) + ] + internal_reaction_ids = last.(internal_reactions) + internal_reaction_idxs = first.(internal_reactions) - # add -1000 * (1-a) ≤ v ≤ 1000 * a which need to be split into forward and backward components - m *= :loopless_reaction_directions^:backward^C.ConstraintTree( + # add loopless variables: flux direction setters (binary) and pseudo gibbs free energy of reaction + m += + :loopless_binary_variables^C.variables( + keys = internal_reaction_ids, + bounds = Ref(C.Binary), + ) + m += :pseudo_gibbs_free_energy_reaction^C.variables(keys = internal_reaction_ids) + + # add -1000 * (1-a) ≤ v ≤ 1000 * a which need to be split into forward and backward components + m *= + :loopless_reaction_directions^:backward^C.ConstraintTree( rid => C.Constraint( - value = m.fluxes[rid].value + max_flux_bound * (1 - m.loopless_binary_variables[rid].value), + value = m.fluxes[rid].value + + max_flux_bound * (1 - m.loopless_binary_variables[rid].value), bound = (0, Inf), ) for rid in internal_reaction_ids ) - m *= :loopless_reaction_directions^:forward^C.ConstraintTree( + m *= + :loopless_reaction_directions^:forward^C.ConstraintTree( rid => C.Constraint( - value = max_flux_bound * m.loopless_binary_variables[rid].value - m.fluxes[rid].value , + value = max_flux_bound * m.loopless_binary_variables[rid].value - + m.fluxes[rid].value, bound = (0, Inf), ) for rid in internal_reaction_ids ) - # add -1000*a + 1 * (1-a) ≤ Gibbs ≤ -1 * a + 1000 * (1 - a) which also need to be split - m *= :loopless_pseudo_gibbs_sign^:backward^C.ConstraintTree( + # add -1000*a + 1 * (1-a) ≤ Gibbs ≤ -1 * a + 1000 * (1 - a) which also need to be split + m *= + :loopless_pseudo_gibbs_sign^:backward^C.ConstraintTree( rid => C.Constraint( - value = m.pseudo_gibbs_free_energy_reaction[rid].value + max_flux_bound * m.loopless_binary_variables[rid].value - strict_inequality_tolerance * (1 - m.loopless_binary_variables[rid].value), + value = m.pseudo_gibbs_free_energy_reaction[rid].value + + max_flux_bound * m.loopless_binary_variables[rid].value - + strict_inequality_tolerance * + (1 - m.loopless_binary_variables[rid].value), bound = (0, Inf), ) for rid in internal_reaction_ids ) - m *= :loopless_pseudo_gibbs_sign^:forward^C.ConstraintTree( + m *= + :loopless_pseudo_gibbs_sign^:forward^C.ConstraintTree( rid => C.Constraint( - value = -strict_inequality_tolerance * m.loopless_binary_variables[rid].value + max_flux_bound * (1 - m.loopless_binary_variables[rid].value) - m.pseudo_gibbs_free_energy_reaction[rid].value, + value = -strict_inequality_tolerance * + m.loopless_binary_variables[rid].value + + max_flux_bound * (1 - m.loopless_binary_variables[rid].value) - + m.pseudo_gibbs_free_energy_reaction[rid].value, bound = (0, Inf), ) for rid in internal_reaction_ids ) - # add N_int' * Gibbs = 0 where N_int = nullspace - N_int = nullspace(Array(A.stoichiometry(model)[:, internal_reaction_idxs])) # no sparse nullspace function - m *= :loopless_condition^C.ConstraintTree( + # add N_int' * Gibbs = 0 where N_int = nullspace + N_int = nullspace(Array(A.stoichiometry(model)[:, internal_reaction_idxs])) # no sparse nullspace function + m *= + :loopless_condition^C.ConstraintTree( Symbol(:nullspace_vector, i) => C.Constraint( - value = sum(col[j] * m.pseudo_gibbs_free_energy_reaction[internal_reaction_ids[j]].value for j in eachindex(col)), - bound = 0, + value = sum( + col[j] * + m.pseudo_gibbs_free_energy_reaction[internal_reaction_ids[j]].value + for j in eachindex(col) + ), + bound = 0, ) for (i, col) in enumerate(eachcol(N_int)) ) - # solve - optimized_constraints( - m; - objective = m.objective.value, - optimizer, - modifications, - ) + # solve + optimized_constraints(m; objective = m.objective.value, optimizer, modifications) end export loopless_flux_balance_analysis diff --git a/src/misc/utils.jl b/src/misc/utils.jl index a1a7e0bf6..7c559b045 100644 --- a/src/misc/utils.jl +++ b/src/misc/utils.jl @@ -8,6 +8,7 @@ reaction, otherwise return false. Checks if on boundary by inspecting the number of metabolites in the reaction stoichiometry. Boundary reactions have only one metabolite, e.g. an exchange reaction, or a sink/demand reaction. """ -is_boundary(model::A.AbstractFBCModel, rxn_id::String) = length(keys(A.reaction_stoichiometry(model, rxn_id))) == 1 +is_boundary(model::A.AbstractFBCModel, rxn_id::String) = + length(keys(A.reaction_stoichiometry(model, rxn_id))) == 1 export is_boundary