fix a few errors to make enzymes work

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),
     )

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
 # </details>
 # ```
 
-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,20 +108,20 @@ 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
 
 # With all the parameters specified, we can directly use the enzyme constrained
 # 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,
 )
 

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

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),
+)

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(

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);

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);

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