From e8f7ce1e923a7b22f3ce59d777535e84b682e0f9 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Fri, 5 Jan 2024 17:05:26 +0100 Subject: [PATCH] fix a few errors to make enzymes work --- .../examples/03-parsimonious-flux-balance.jl | 4 +-- .../examples/05-enzyme-constrained-models.jl | 18 +++++------ src/COBREXA.jl | 2 +- src/builders/objectives.jl | 14 ++++---- src/frontend/enzyme_constrained.jl | 32 ++++++++++++------- src/frontend/moma.jl | 16 ++++------ src/frontend/parsimonious.jl | 10 +++--- src/solver.jl | 1 + 8 files changed, 52 insertions(+), 45 deletions(-) diff --git a/docs/src/examples/03-parsimonious-flux-balance.jl b/docs/src/examples/03-parsimonious-flux-balance.jl index 5fa08f335..a39a5cfc6 100644 --- a/docs/src/examples/03-parsimonious-flux-balance.jl +++ b/docs/src/examples/03-parsimonious-flux-balance.jl @@ -54,7 +54,7 @@ model |> parsimonious_flux_balance_analysis(Clarabel.Optimizer; settings = [sile # the quadratic objective (this approach is much more flexible). ctmodel = fbc_model_constraints(model) -ctmodel *= :l2objective^squared_sum_objective(ctmodel.fluxes) +ctmodel *= :l2objective^squared_sum_value(ctmodel.fluxes) ctmodel.objective.bound = 0.3 # set growth rate # TODO currently breaks opt_model = optimization_model( @@ -90,7 +90,7 @@ minimize_metabolic_adjustment(ref_sol, Clarabel.Optimizer; settings = [silence]) ctmodel = fbc_model_constraints(model) ctmodel *= - :minoxphospho^squared_sum_error_objective( + :minoxphospho^squared_sum_error_value( ctmodel.fluxes, Dict(:ATPS4r => 33.0, :CYTBD => 22.0), ) diff --git a/docs/src/examples/05-enzyme-constrained-models.jl b/docs/src/examples/05-enzyme-constrained-models.jl index 72e321cda..ab6d75a78 100644 --- a/docs/src/examples/05-enzyme-constrained-models.jl +++ b/docs/src/examples/05-enzyme-constrained-models.jl @@ -64,10 +64,10 @@ for rid in A.reactions(model) for (i, grr) in enumerate(grrs) d = get!(reaction_isozymes, rid, Dict{String,Isozyme}()) # each isozyme gets a unique name - d["isozyme_"*string(i)] = Isozyme( # SimpleIsozyme struct is defined by COBREXA + d["isozyme_"*string(i)] = Isozyme( gene_product_stoichiometry = Dict(grr .=> fill(1.0, size(grr))), # assume subunit stoichiometry of 1 for all isozymes - kcat_forward = ecoli_core_reaction_kcats[rid] * 3600.0, # forward reaction turnover number units = 1/h - kcat_backward = ecoli_core_reaction_kcats[rid] * 3600.0, # reverse reaction turnover number units = 1/h + kcat_forward = ecoli_core_reaction_kcats[rid] * 3.6, # forward reaction turnover number units = 1/h + kcat_reverse = ecoli_core_reaction_kcats[rid] * 3.6, # reverse reaction turnover number units = 1/h ) end end @@ -90,7 +90,7 @@ end # # ``` -gene_molar_masses = ecoli_core_gene_product_masses +gene_product_molar_masses = ecoli_core_gene_product_masses #!!! warning "Molar mass units" # Take care with the units of the molar masses. In literature they are @@ -108,7 +108,7 @@ gene_molar_masses = ecoli_core_gene_product_masses # The capacity limitation usually denotes an upper bound of protein available to # the cell. -total_enzyme_capacity = 0.1 # g enzyme/gDW +total_enzyme_capacity = 100.0 # mg enzyme/gDW ### Running a basic enzyme constrained model @@ -116,12 +116,12 @@ total_enzyme_capacity = 0.1 # g enzyme/gDW # convenience function to run enzyme constrained FBA in one shot: ec_solution = enzyme_constrained_flux_balance_analysis( - model, + model; reaction_isozymes, - gene_molar_masses, - [("total_proteome_bound", A.genes(model), total_enzyme_capacity)]; + gene_product_molar_masses, + capacity = total_enzyme_capacity, settings = [set_optimizer_attribute("IPM_IterationsLimit", 10_000)], - unconstrain_reactions = ["EX_glc__D_e"], + #unconstrain_reactions = ["EX_glc__D_e"], optimizer = Tulip.Optimizer, ) diff --git a/src/COBREXA.jl b/src/COBREXA.jl index 7e5d2adc3..c83cfa28c 100644 --- a/src/COBREXA.jl +++ b/src/COBREXA.jl @@ -55,8 +55,8 @@ include("builders/enzymes.jl") include("builders/fbc.jl") include("builders/knockouts.jl") include("builders/loopless.jl") +include("builders/mmdf.jl") include("builders/objectives.jl") -include("builders/thermodynamic.jl") include("builders/unsigned.jl") # simplified front-ends for the above diff --git a/src/builders/objectives.jl b/src/builders/objectives.jl index 9b1d487e7..6861e56d1 100644 --- a/src/builders/objectives.jl +++ b/src/builders/objectives.jl @@ -19,15 +19,14 @@ $(TYPEDSIGNATURES) TODO """ -squared_sum_objective(x::C.ConstraintTree) = - squared_sum_error_objective(x, Dict(keys(x) .=> 0.0)) +squared_sum_value(x::C.ConstraintTree) = squared_sum_error_value(x, Dict(keys(x) .=> 0.0)) """ $(TYPEDSIGNATURES) TODO useful for L1 parsimonious stuff """ -function sum_objective(x...) +function sum_value(x...) res = zero(C.LinearValue) for ct in x C.map(ct) do c @@ -42,8 +41,7 @@ $(TYPEDSIGNATURES) TODO """ -squared_sum_error_objective(constraints::C.ConstraintTree, target::Dict{Symbol,Float64}) = - sum( - (C.squared(C.value(c) - target[k]) for (k, c) in constraints if haskey(target, k)), - init = zero(C.LinearValue), - ) +squared_sum_error_value(constraints::C.ConstraintTree, target::Dict{Symbol,Float64}) = sum( + (C.squared(C.value(c) - target[k]) for (k, c) in constraints if haskey(target, k)), + init = zero(C.LinearValue), +) diff --git a/src/frontend/enzyme_constrained.jl b/src/frontend/enzyme_constrained.jl index 7223ca286..d7b29dec8 100644 --- a/src/frontend/enzyme_constrained.jl +++ b/src/frontend/enzyme_constrained.jl @@ -73,16 +73,16 @@ function enzyme_constrained_flux_balance_analysis( :isozyme_reverse_amounts^isozyme_amounts + :gene_product_amounts^C.variables( keys = Symbol.(A.genes(model)), - bounds = Between(0, Inf), + bounds = C.Between(0, Inf), ) # connect all parts with constraints constraints = constraints * :directional_flux_balance^sign_split_constraints( - constraints.fluxes_forward, - constraints.fluxes_reverse, - constraints.fluxes, + positive = constraints.fluxes_forward, + negative = constraints.fluxes_reverse, + signed = constraints.fluxes, ) * :isozyme_flux_forward_balance^isozyme_flux_constraints( constraints.isozyme_forward_amounts, @@ -104,18 +104,28 @@ function enzyme_constrained_flux_balance_analysis( constraints.gene_product_amounts, (constraints.isozyme_forward_amounts, constraints.isozyme_reverse_amounts), (rid, isozyme) -> maybemap( - x -> x.gene_product_stoichiometry, + x -> [(Symbol(k), v) for (k, v) in x.gene_product_stoichiometry], maybeget(reaction_isozymes, string(rid), string(isozyme)), ), ) * - :gene_product_capacity_limits^C.ConstraintTree( - Symbol(id) => C.Constraint( + :gene_product_capacity^( + capacity isa Float64 ? + C.Constraint( value = sum( - constraints.gene_product_amounts[gp].value * - gene_product_molar_masses[gp] for gp in gps + gpa.value * gene_product_molar_masses[String(gp)] for + (gp, gpa) in constraints.gene_product_amounts ), - bound = C.Between(0, limit), - ) for (id, gps, limit) in capacity_limits + bound = C.Between(0, capacity), + ) : + C.ConstraintTree( + Symbol(id) => C.Constraint( + value = sum( + constraints.gene_product_amounts[Symbol(gp)].value * + gene_product_molar_masses[gp] for gp in gps + ), + bound = C.Between(0, limit), + ) for (id, gps, limit) in capacity_limits + ) ) optimized_constraints( diff --git a/src/frontend/moma.jl b/src/frontend/moma.jl index 4e2b621bd..f97b4a4f3 100644 --- a/src/frontend/moma.jl +++ b/src/frontend/moma.jl @@ -29,7 +29,7 @@ similar perturbation). MOMA solution then gives an expectable "easiest" adjustment of the organism towards a somewhat working state. Reference fluxes that do not exist in the model are ignored (internally, the -objective is constructed via [`squared_sum_error_objective`](@ref)). +objective is constructed via [`squared_sum_error_value`](@ref)). Additional parameters are forwarded to [`optimized_constraints`](@ref). """ @@ -40,7 +40,7 @@ function minimization_of_metabolic_adjustment_analysis( kwargs..., ) constraints = fbc_model_constraints(model) - objective = squared_sum_error_objective(constraints.fluxes, reference_fluxes) + objective = squared_sum_error_value(constraints.fluxes, reference_fluxes) optimized_constraints( constraints * :minimal_adjustment_objective^C.Constraint(objective); optimizer, @@ -122,15 +122,13 @@ function linear_minimization_of_metabolic_adjustment_analysis( constraints *= :reference_diff^difference constraints *= :reference_directional_diff_balance^sign_split_constraints( - constraints.reference_positive_diff, - constraints.reference_negative_diff, - difference, + positive = constraints.reference_positive_diff, + negative = constraints.reference_negative_diff, + signed = difference, ) - objective = sum_objective( - constraints.reference_positive_diff, - constraints.reference_negative_diff, - ) + objective = + sum_value(constraints.reference_positive_diff, constraints.reference_negative_diff) optimized_constraints( constraints * :linear_minimal_adjustment_objective^C.Constraint(objective); diff --git a/src/frontend/parsimonious.jl b/src/frontend/parsimonious.jl index 0075e505d..a995e2148 100644 --- a/src/frontend/parsimonious.jl +++ b/src/frontend/parsimonious.jl @@ -105,7 +105,7 @@ function parsimonious_flux_balance_analysis( kwargs..., ) constraints = fbc_model_constraints(model) - parsimonious_objective = squared_sum_objective(constraints.fluxes) + parsimonious_objective = squared_sum_value(constraints.fluxes) parsimonious_optimized_constraints( constraints * :parsimonious_objective^C.Constraint(parsimonious_objective); optimizer, @@ -151,12 +151,12 @@ function linear_parsimonious_flux_balance_analysis( :fluxes_reverse^unsigned_negative_contribution_variables(ct.fluxes) constraints *= :directional_flux_balance^sign_split_constraints( - ct.fluxes_forward, - ct.fluxes_reverse, - ct.fluxes, + positive = ct.fluxes_forward, + negative = ct.fluxes_reverse, + signed = ct.fluxes, ) - parsimonious_objective = sum_objective(ct.fluxes_forward, ct.fluxes_reverse) + parsimonious_objective = sum_value(ct.fluxes_forward, ct.fluxes_reverse) parsimonious_optimized_constraints( constraints * :parsimonious_objective^C.Constraint(parsimonious_objective); diff --git a/src/solver.jl b/src/solver.jl index 4a987a361..ad8032eef 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -72,6 +72,7 @@ function optimization_model( boolean = J.@variable(model, binary = true) J.@constraint(model, C.substitute(v, x) == b.a + boolean * (b.b - b.a)) end + add_constraint(::C.Value, _::Nothing) = nothing function add_constraint(c::C.Constraint) add_constraint(c.value, c.bound) end