add basic pipe-able pfba

src/analysis/flux_balance_analysis.jl

@@ -82,3 +82,16 @@ function flux_balance_analysis(
     optimize!(opt_model)
     return opt_model
 end
+
+"""
+$(TYPEDSIGNATURES)
+
+Pipe friendly variant of [`flux_balance_analysis`](@ref). Forwards all arguments. 
+
+# Example
+```
+model |> flux_balance_analysis(Tulip.Optimizer; modifications = [...])
+```
+"""
+flux_balance_analysis(optimizer; kwargs...) =
+    model -> flux_balance_analysis(model, optimizer; kwargs...)

src/reconstruction/pipes/parsimonious.jl

@@ -0,0 +1,9 @@
+"""
+$(TYPEDSIGNATURES)
+
+Pipe-able version of the [`ParsimoniousModel`](@ref) wrapper.
+"""
+with_parsimonious_solution(semantics::Vector{Symbol}) =
+    model -> ParsimoniousModel(model, semantics)
+
+with_parsimonious_solution(semantic::Symbol) = with_parsimonious_solution([semantic])

src/types/wrappers/ParsimoniousModel.jl

@@ -0,0 +1,48 @@
+
+"""
+$(TYPEDEF)
+
+A wrapper that adds a quadratic objective that minimizes the sum of a set of
+squared variables. If the bounds of the growth rate of the model are fixed, this
+corresponds to finding a parsimonious solution (i.e. pFBA).
+
+This is used to implement [`parsimonious_flux_balance_analysis`](@ref).
+
+# Example
+```
+model |> with_changed_bound("biomass", lower_bound = 0.1) |> with_parsimonious_solution(:enzymes) |> flux_balance_analysis(Clarabel.Optimizer)
+```
+"""
+struct ParsimoniousModel <: AbstractModelWrapper
+    inner::AbstractMetabolicModel
+    var_ids::Vector{String}
+end
+
+function ParsimoniousModel(model::AbstractMetabolicModel, semantics::Vector{Symbol})
+    var_ids = vcat([@eval $(Symbol(sem, :s))($(model)) for sem in semantics]...)
+    ParsimoniousModel(model, var_ids)
+end
+
+ParsimoniousModel(model::AbstractMetabolicModel, semantic::Symbol) =
+    ParsimoniousModel(model, [semantic])
+
+Accessors.unwrap_model(m::ParsimoniousModel) = m.inner
+
+"""
+$(TYPEDSIGNATURES)
+
+Return a negative, uniformly weighted, quadratic-only objective representing the
+squared sum of `model.var_ids`. 
+"""
+function Accessors.objective(model::ParsimoniousModel)::SparseMat
+    obj = spzeros(n_variables(model), n_variables(model) + 1) # + 1 for QP solver formulation 
+
+    idxs = indexin(model.var_ids, variables(model))
+    j = findfirst(isnothing, idxs)
+    isnothing(j) || throw(DomainError(model.var_ids[j], "No variable with this ID."))
+
+    for i in idxs
+        obj[i, i] = -1.0 # negative because objective will be maximized
+    end
+    obj
+end

test/analysis/parsimonious_flux_balance_analysis.jl

@@ -14,4 +14,16 @@
     # The used optimizer doesn't really converge to the same answer everytime
     # here, we therefore tolerate a wide range of results.
     @test isapprox(d["biomass1"], 10.0, atol = QP_TEST_TOLERANCE)
+
+    sol =
+        model |>
+        with_changed_bound("biomass1", lower_bound = 10.0) |>
+        with_parsimonious_solution(:reaction) |>
+        flux_balance_analysis(Clarabel.Optimizer)
+
+    @test isapprox(
+        values_dict(:reaction, model, sol)["biomass1"],
+        10.0,
+        atol = QP_TEST_TOLERANCE,
+    )
 end

test/reconstruction/enzyme_constrained.jl

@@ -72,19 +72,34 @@
     @test isapprox(prot_mass, total_gene_product_mass, atol = TEST_TOLERANCE)
     @test isapprox(prot_mass, mass_groups["uncategorized"], atol = TEST_TOLERANCE)
 
-    # test enzyme objective
+    # test pFBA
     growth_lb = rxn_fluxes["BIOMASS_Ecoli_core_w_GAM"] * 0.9
-    opt_model = flux_balance_analysis(
-        gm,
-        Tulip.Optimizer;
-        modifications = [
-            change_objective(genes(gm); weights = [], sense = MIN_SENSE),
-            change_constraint("BIOMASS_Ecoli_core_w_GAM", lower_bound = growth_lb),
-            change_optimizer_attribute("IPM_IterationsLimit", 1000),
-        ],
+
+    gm_changed_bound =
+        model |>
+        with_changed_bound("BIOMASS_Ecoli_core_w_GAM", lower_bound = growth_lb) |>
+        with_enzyme_constraints(
+            gene_product_mass_group_bound = Dict(
+                "uncategorized" => total_gene_product_mass,
+            ),
+        )
+
+    pfba_sol =
+        gm_changed_bound |>
+        with_parsimonious_solution(:enzyme) |>
+        flux_balance_analysis(Clarabel.Optimizer)
+
+    @test isapprox(
+        values_dict(:reaction, gm_changed_bound, pfba_sol)["BIOMASS_Ecoli_core_w_GAM"],
+        0.7315450597991255,
+        atol = QP_TEST_TOLERANCE,
+    )
+
+    @test isapprox(
+        values_dict(:enzyme_group, gm_changed_bound, pfba_sol)["uncategorized"],
+        91.425,
+        atol = QP_TEST_TOLERANCE,
     )
-    mass_groups_min = values_dict(:enzyme_group, gm, opt_model)
-    @test mass_groups_min["uncategorized"] < mass_groups["uncategorized"]
 end
 
 @testset "GECKO small model" begin