diff --git a/.gitignore b/.gitignore
index 46e71d415..caac77fde 100644
--- a/.gitignore
+++ b/.gitignore
@@ -20,3 +20,8 @@ Manifest.toml
 
 # ignore container files
 *.sif
+
+# ignore models
+*.xml
+*.json
+*.mat
\ No newline at end of file
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 68f70fccf..4f7aa3de7 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -71,9 +71,9 @@ variables:
   before_script:
     - docker login -u $CI_USER_NAME -p $GITLAB_ACCESS_TOKEN $CI_REGISTRY
 
-.global_julia18: &global_julia18
+.global_julia19: &global_julia19
   variables:
-    JULIA_VER: "v1.8.3"
+    JULIA_VER: "v1.9.4"
 
 .global_julia16: &global_julia16
   variables:
@@ -120,13 +120,11 @@ variables:
 # any available docker and current julia
 #
 
-docker:julia1.8:
+docker:julia1.9:
   stage: test
   image: $CI_REGISTRY/r3/docker/julia-custom
   script:
     - julia --check-bounds=yes --inline=yes --project=@. -e "import Pkg; Pkg.test(; coverage = true)"
-  after_script:
-    - julia --project=test/coverage test/coverage/coverage-summary.jl
   <<: *global_trigger_pull_request
 
 #
@@ -134,12 +132,12 @@ docker:julia1.8:
 # built & deployed
 #
 
-linux:julia1.8:
+linux:julia1.9:
   stage: test
   tags:
     - slave01
   <<: *global_trigger_full_tests
-  <<: *global_julia18
+  <<: *global_julia19
   <<: *global_env_linux
 
 linux:julia1.6:
@@ -154,22 +152,22 @@ linux:julia1.6:
 # Additional platform&environment compatibility tests
 #
 
-windows8:julia1.8:
+windows8:julia1.9:
   stage: test-compat
   <<: *global_trigger_compat_tests
-  <<: *global_julia18
+  <<: *global_julia19
   <<: *global_env_win8
 
-windows10:julia1.8:
+windows10:julia1.9:
   stage: test-compat
   <<: *global_trigger_compat_tests
-  <<: *global_julia18
+  <<: *global_julia19
   <<: *global_env_win10
 
-mac:julia1.8:
+mac:julia1.9:
   stage: test-compat
   <<: *global_trigger_compat_tests
-  <<: *global_julia18
+  <<: *global_julia19
   <<: *global_env_mac
 
 windows8:julia1.6:
diff --git a/Project.toml b/Project.toml
index dc4027e93..14daa42be 100644
--- a/Project.toml
+++ b/Project.toml
@@ -18,20 +18,32 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 
 [compat]
 AbstractFBCModels = "0.2"
+Aqua = "0.7"
 Clarabel = "0.3"
-ConstraintTrees = "0.4"
+ConstraintTrees = "0.5"
 DistributedData = "0.2"
 DocStringExtensions = "0.8, 0.9"
+Downloads = "1"
+GLPK = "1"
+JSONFBCModels = "0.1"
 JuMP = "1"
+SBMLFBCModels = "0.1"
+SHA = "0.7, 1"
 StableRNGs = "1.0"
+Test = "1"
+Tulip = "0.9"
 julia = "1.5"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
+Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
+JSONFBCModels = "475c1105-d6ed-49c1-9b32-c11adca6d3e8"
+SBMLFBCModels = "3e8f9d1a-ffc1-486d-82d6-6c7276635980"
+SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa"
 
 [targets]
-test = ["Aqua", "Clarabel", "GLPK", "Test", "Tulip"]
+test = ["Aqua", "Clarabel", "Downloads", "GLPK", "SHA", "Test", "Tulip", "JSONFBCModels"]
diff --git a/README.md b/README.md
index 47237b005..b9fd793f1 100644
--- a/README.md
+++ b/README.md
@@ -115,7 +115,7 @@ download("http://bigg.ucsd.edu/static/models/e_coli_core.xml", "e_coli_core.xml"
 model = load_model("e_coli_core.xml")
 
 # run a FBA
-fluxes = flux_balance_analysis_dict(model, Tulip.Optimizer)
+fluxes = flux_balance_dict(model, Tulip.Optimizer)
 ```
 
 The variable `fluxes` will now contain a dictionary of the computed optimal
diff --git a/docs/Project.toml b/docs/Project.toml
index 5b3dcde68..c8ea85a2a 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,15 +1,10 @@
 [deps]
+AbstractFBCModels = "5a4f3dfa-1789-40f8-8221-69268c29937c"
 COBREXA = "babc4406-5200-4a30-9033-bf5ae714c842"
-CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
-Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
-ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-Escher = "8cc96de1-1b23-48cb-9272-618d67962629"
-GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+JSONFBCModels = "475c1105-d6ed-49c1-9b32-c11adca6d3e8"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa"
 
diff --git a/docs/make.jl b/docs/make.jl
index 0b98cc0c3..1e716e90c 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -25,6 +25,7 @@ find_mds(path) =
         filter(x -> endswith(x, ".md"), readdir(joinpath(@__DIR__, "src", path))),
     )
 
+#TODO migrate this to Documenter-1, and make all checks strict
 # build the docs
 makedocs(
     modules = [COBREXA],
@@ -51,10 +52,11 @@ makedocs(
             "Contents" => "concepts.md"
             find_mds("concepts")
         ],
-        "Reference" => [
-            "Contents" => "reference.md"
-            find_mds("reference")
-        ],
+        "Reference" => "reference.md",
+        #[ # TODO re-add this when the reference gets bigger
+        #"Contents" => "reference.md"
+        #find_mds("reference")
+        #],
     ],
 )
 
diff --git a/docs/src/concepts/1_screen.md b/docs/src/concepts/1_screen.md
index a00204092..cb482a4b6 100644
--- a/docs/src/concepts/1_screen.md
+++ b/docs/src/concepts/1_screen.md
@@ -24,7 +24,7 @@ screen_variants(
         [with_changed_bound("O2t", lb = 0, ub = 0)],  # disable O2 transport
         [with_changed_bound("CO2t", lb = 0, ub = 0), with_changed_bound("O2t", lb = 0, ub = 0)],  # disable both transports
     ],
-    m -> flux_balance_analysis_dict(m, Tulip.Optimizer)["BIOMASS_Ecoli_core_w_GAM"],
+    m -> flux_balance_dict(m, Tulip.Optimizer)["BIOMASS_Ecoli_core_w_GAM"],
 )
 ```
 The call specifies a model (the `m` that we have loaded) that is being tested,
@@ -89,7 +89,7 @@ res = screen_variants(m,
             ["EX_h2o_e", "EX_co2_e", "EX_o2_e", "EX_nh4_e"], # and this set of exchanges
         )
     ],
-    m -> flux_balance_analysis_dict(m, Tulip.Optimizer)["BIOMASS_Ecoli_core_w_GAM"],
+    m -> flux_balance_dict(m, Tulip.Optimizer)["BIOMASS_Ecoli_core_w_GAM"],
 )
 ```
 
@@ -199,7 +199,7 @@ screen_variants(
         [with_disabled_oxygen_transport],
         [with_disabled_reaction("NH4t")],
     ],
-    m -> flux_balance_analysis_dict(m, Tulip.Optimizer)["BIOMASS_Ecoli_core_w_GAM"],
+    m -> flux_balance_dict(m, Tulip.Optimizer)["BIOMASS_Ecoli_core_w_GAM"],
 )
 ```
 
@@ -222,7 +222,7 @@ That should get you the results for all new variants of the model:
 
 Some analysis functions may take additional arguments, which you might want to
 vary for the analysis. `modifications` argument of
-[`flux_balance_analysis_dict`](@ref) is one example of such argument, allowing
+[`flux_balance_dict`](@ref) is one example of such argument, allowing
 you to specify details of the optimization procedure.
 
 [`screen`](@ref) function allows you to do precisely that -- apart from
@@ -242,7 +242,7 @@ iterations needed for Tulip solver to find a feasible solution:
 screen(m,
     args = [(i,) for i in 5:15],  # the iteration counts, packed in 1-tuples
     analysis = (m,a) -> # `args` elements get passed as the extra parameter here
-        flux_balance_analysis_vec(m,
+        flux_balance_vec(m,
             Tulip.Optimizer;
             modifications=[modify_optimizer_attribute("IPM_IterationsLimit", a)],
         ),
diff --git a/docs/src/examples/01-loading-and-saving.jl b/docs/src/examples/01-loading-and-saving.jl
index 584cee8cb..05e8adc68 100644
--- a/docs/src/examples/01-loading-and-saving.jl
+++ b/docs/src/examples/01-loading-and-saving.jl
@@ -4,3 +4,5 @@
 using COBREXA
 
 # TODO: download the models into a single directory that can get cached. Probably best have a fake mktempdir().
+#
+# TODO: demonstrate download_model here and explain how to get hashes (simply not fill them in for the first time)
diff --git a/docs/src/examples/02-flux-balance-analysis.jl b/docs/src/examples/02-flux-balance-analysis.jl
index 678930970..9bbef41cf 100644
--- a/docs/src/examples/02-flux-balance-analysis.jl
+++ b/docs/src/examples/02-flux-balance-analysis.jl
@@ -1,6 +1,51 @@
 
-# # Flux balance analysis
+# # Flux balance analysis (FBA)
+
+# Here we use [`flux_balance`](@ref) and several related functions to
+# find an optimal flux in the *E. coli* "core" model. We will need the model,
+# which we can download using [`download_model`](@ref):
 
 using COBREXA
 
-# TODO: run FBA on a FBC model
+download_model(
+    "http://bigg.ucsd.edu/static/models/e_coli_core.json",
+    "e_coli_core.json",
+    "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
+)
+
+# Additionally to COBREXA and the model format package, we will need a solver
+# -- let's use Tulip here:
+
+import JSONFBCModels
+import Tulip
+
+model = load_model("e_coli_core.json")
+
+# ## Running a FBA
+#
+# There are many possibilities on how to arrange the metabolic model into the
+# optimization framework and how to actually solve it. The "usual" assumed one
+# is captured in the default behavior of function
+# [`flux_balance`](@ref):
+
+solution = flux_balance(model, Tulip.Optimizer)
+
+@test isapprox(solution.objective, 0.8739, atol = TEST_TOLERANCE) #src
+
+# The result contains a tree of all optimized values in the model, including
+# fluxes, the objective value, and possibly others (given by what the model
+# contains).
+#
+# You can explore the dot notation to explore the solution, extracting e.g. the
+# value of the objective:
+
+solution.objective
+
+# ...or the value of the flux through the given reaction (note the solution is
+# not unique in FBA):
+
+solution.fluxes.PFK
+
+# ...or make a "table" of all fluxes through all reactions:
+
+collect(solution.fluxes)
diff --git a/docs/src/examples/02a-optimizer-parameters.jl b/docs/src/examples/02a-optimizer-parameters.jl
new file mode 100644
index 000000000..766d6c46a
--- /dev/null
+++ b/docs/src/examples/02a-optimizer-parameters.jl
@@ -0,0 +1,57 @@
+
+# # Changing optimizer parameters
+#
+# Many optimizers require fine-tuning to produce best results. You can pass in
+# additional optimizer settings via the `modifications` parameter of
+# [`flux_balance`](@ref). These include e.g.
+#
+# - [`set_optimizer_attribute`](@ref) (typically allowing you to tune e.g.
+#   iteration limits, tolerances, or floating-point precision)
+# - [`set_objective_sense`](@ref) (allowing you to change and reverse the
+#   optimization direction, if required)
+# - [`silence`](@ref) to disable the debug output of the optimizer
+# - and even [`set_optimizer`](@ref), which changes the optimizer
+#   implementation used (this is not quite useful in this case, but becomes
+#   beneficial with more complex, multi-stage optimization problems)
+#
+# To demonstrate this, we'll use the usual toy model:
+
+using COBREXA
+import JSONFBCModels, Tulip
+
+download_model(
+    "http://bigg.ucsd.edu/static/models/e_coli_core.json",
+    "e_coli_core.json",
+    "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
+)
+
+model = load_model("e_coli_core.json")
+
+# Running a FBA with a silent optimizer that has slightly increased iteration
+# limit for IPM algorithm may now look as follows:
+solution = flux_balance(
+    model,
+    Tulip.Optimizer;
+    modifications = [silence, set_optimizer_attribute("IPM_IterationsLimit", 1000)],
+)
+
+@test !isnothing(solution) #src
+
+# To see some of the effects of the configuration changes, you may e.g.
+# deliberately cripple the optimizer's possibilities to a few iterations, which
+# will cause it to fail and return no solution:
+
+solution = flux_balance(
+    model,
+    Tulip.Optimizer;
+    modifications = [silence, set_optimizer_attribute("IPM_IterationsLimit", 2)],
+)
+
+println(solution)
+
+@test isnothing(solution) #src
+
+# Applicable optimizer attributes are documented in the documentations of the
+# respective optimizers. To browse the possibilities, you may want to see the
+# [JuMP documentation page that summarizes the references to the available
+# optimizers](https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers).
diff --git a/docs/src/examples/02b-model-modifications.jl b/docs/src/examples/02b-model-modifications.jl
new file mode 100644
index 000000000..22307fe7c
--- /dev/null
+++ b/docs/src/examples/02b-model-modifications.jl
@@ -0,0 +1,28 @@
+
+# # Making adjustments to the model
+#
+# Typically, we do not need to solve the models as they come from the authors
+# (someone else already did that!), but we want to perform various
+# perturbations in the model structure and conditions, and explore how the
+# model behaves in the changed conditions.
+#
+# With COBREXA, there are 2 different approaches that one can take:
+# 1. We can change the model structure and use the changed metabolic model.
+# This is better for doing simple and small but systematic modifications, such
+# as removing metabolites, adding reactions, etc.
+# 2. We can intercept the pipeline that converts the metabolic model to
+# constraints and then to the optimizer representation, and make small
+# modifications along that way. This is better for various technical model
+# adjustments, such as using combined objectives or adding reaction-coupling
+# constraints.
+#
+# Here we demonstrate the first, "modelling" approach. The main advantage of
+# that approach is that the modified model is still a FBC model, and you can
+# export, save and share it via the AbstractFBCModels interace. The main
+# disadvantage is that the "common" FBC model interface does not easily express
+# various complicated constructions (communities, reaction coupling, enzyme
+# constraints, etc.) -- see the [example about modifying the
+# constraints](02c-constraint-modifications.md) for a closer look on how to
+# modify even such complex constructions.
+#
+# TODO here. :)
diff --git a/docs/src/examples/02c-constraint-modifications.jl b/docs/src/examples/02c-constraint-modifications.jl
new file mode 100644
index 000000000..584ba25ac
--- /dev/null
+++ b/docs/src/examples/02c-constraint-modifications.jl
@@ -0,0 +1,64 @@
+
+# # Making adjustments to the constraint system
+#
+# In the [previous example about model
+# adjustments](02b-model-modifications.md), we noted that some constraint
+# systems may be to complex to be changed within the limits of the usual FBC
+# model view, and we may require a sharper tool to do the changes we need. This
+# example shows how to do that by modifying the constraint systems that are
+# generated within COBREXA to represent the metabolic model contents.
+#
+# ## Background: Model-to-optimizer pipeline
+#
+# ## Background: Constraint trees
+#
+# ## Changing the model-to-optimizer pipeline
+#
+# TODO clean up the stuff below:
+
+using COBREXA
+
+download_model(
+    "http://bigg.ucsd.edu/static/models/e_coli_core.json",
+    "e_coli_core.json",
+    "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
+)
+
+import JSONFBCModels
+import Tulip
+
+model = load_model("e_coli_core.json")
+
+# ## Customizing the model
+
+# We can also modify the model. The most explicit way to do this is
+# to make a new constraint tree representation of the model.
+
+import ConstraintTrees as C
+
+ctmodel = fbc_model_constraints(model)
+
+fermentation = ctmodel.fluxes.EX_ac_e.value + ctmodel.fluxes.EX_etoh_e.value
+
+forced_mixed_fermentation =
+    ctmodel * :fermentation^C.Constraint(fermentation, (10.0, 1000.0)) # new modified model is created
+
+vt = flux_balance(forced_mixed_fermentation, Tulip.Optimizer; modifications = [silence])
+
+@test isapprox(vt.objective, 0.6337, atol = TEST_TOLERANCE) #src
+
+# Models that cannot be solved return `nothing`. In the example below, the
+# underlying model is modified.
+
+ctmodel.fluxes.ATPM.bound = (1000.0, 10000.0)
+
+vt = flux_balance(ctmodel, Tulip.Optimizer; modifications = [silence])
+
+@test isnothing(vt) #src
+
+# Models can also be piped into the analysis functions
+
+ctmodel.fluxes.ATPM.bound = (8.39, 10000.0) # revert
+vt = ctmodel |> flux_balance(Tulip.Optimizer; modifications = [silence])
+
+@test isapprox(vt.objective, 0.8739, atol = TEST_TOLERANCE) #src
diff --git a/docs/src/reference.md b/docs/src/reference.md
index f0b01a550..51abbb38c 100644
--- a/docs/src/reference.md
+++ b/docs/src/reference.md
@@ -1,6 +1,59 @@
 # API reference
 
-```@contents
-Pages = ["reference/types.md", "reference/io.md", "reference/builders.md", "reference/solver.md"]
-Depth = 2
+## Helper types
+
+```@autodocs
+Modules = [COBREXA]
+Pages = ["src/types.jl"]
+```
+
+## Model loading and saving
+
+```@autodocs
+Modules = [COBREXA]
+Pages = ["src/io.jl"]
+```
+
+## Solver interface
+
+```@autodocs
+Modules = [COBREXA]
+Pages = ["src/solver.jl"]
+```
+
+## Constraint system building
+
+```@autodocs
+Modules = [COBREXA]
+Pages = ["src/builders/core.jl"]
 ```
+
+### Genetic constraints
+
+```@autodocs
+Modules = [COBREXA]
+Pages = ["src/builders/genes.jl"]
+```
+
+### Objective function helpers
+
+```@autodocs
+Modules = [COBREXA]
+Pages = ["src/builders/objectives.jl"]
+```
+
+## Analysis functions
+
+```@autodocs
+Modules = [COBREXA]
+Pages = ["src/analysis/flux_balance.jl", "src/analysis/parsimonious_flux_balance.jl"]
+```
+
+### Analysis modifications
+
+```@autodocs
+Modules = [COBREXA]
+Pages = ["src/analysis/modifications.jl"]
+```
+
+## Distributed analysis
diff --git a/docs/src/reference/io.md b/docs/src/reference/io.md
deleted file mode 100644
index 5c494ea0a..000000000
--- a/docs/src/reference/io.md
+++ /dev/null
@@ -1,7 +0,0 @@
-
-# Model loading and saving
-
-```@autodocs
-Modules = [COBREXA]
-Pages = ["src/io.jl"]
-```
diff --git a/docs/src/reference/types.md b/docs/src/reference/types.md
deleted file mode 100644
index bb70365c6..000000000
--- a/docs/src/reference/types.md
+++ /dev/null
@@ -1,7 +0,0 @@
-
-# Helper types
-
-```@autodocs
-Modules = [COBREXA]
-Pages = ["src/types.jl"]
-```
diff --git a/src/COBREXA.jl b/src/COBREXA.jl
index 5baaa18cc..14e8a2257 100644
--- a/src/COBREXA.jl
+++ b/src/COBREXA.jl
@@ -10,10 +10,15 @@ import ConstraintTrees as C
 import JuMP as J
 
 include("types.jl")
+include("io.jl")
 include("solver.jl")
 
 include("builders/core.jl")
 include("builders/genes.jl")
 include("builders/objectives.jl")
 
+include("analysis/modifications.jl")
+include("analysis/flux_balance.jl")
+include("analysis/parsimonious_flux_balance.jl")
+
 end # module COBREXA
diff --git a/src/analysis/flux_balance.jl b/src/analysis/flux_balance.jl
new file mode 100644
index 000000000..f755b45fc
--- /dev/null
+++ b/src/analysis/flux_balance.jl
@@ -0,0 +1,65 @@
+"""
+$(TYPEDSIGNATURES)
+
+Run flux balance analysis (FBA) on the `model`, optionally specifying
+`modifications` to the problem.  Basically, FBA solves this optimization
+problem:
+```
+max cᵀx
+s.t. S x = b
+     xₗ ≤ x ≤ xᵤ
+```
+See "Orth, J., Thiele, I. & Palsson, B. What is flux balance analysis?. Nat
+Biotechnol 28, 245-248 (2010). https://doi.org/10.1038/nbt.1614" for more
+information.
+
+The `optimizer` must be set to a `JuMP`-compatible optimizer, such as
+`GLPK.Optimizer` or `Tulip.Optimizer`.
+
+Optionally, you may specify one or more modifications to be applied to the model
+before the analysis, such as [`set_objective_sense`](@ref),
+[`set_optimizer`](@ref), [`set_optimizer_attribute`](@ref), and
+[`silence`](@ref).
+
+Returns a [`C.ValueTree`](@ref).
+
+# Example
+```
+model = load_model("e_coli_core.json")
+solution = flux_balance(model, GLPK.optimizer)
+```
+"""
+function flux_balance(model::A.AbstractFBCModel, optimizer; modifications = [])
+    ctmodel = fbc_model_constraints(model)
+    flux_balance(ctmodel, optimizer; modifications)
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+A variant of [`flux_balance`](@ref) that takes in a
+[`C.ConstraintTree`](@ref) as the model to optimize. The objective is inferred
+from the field `objective` in `ctmodel`. All other arguments are forwarded.
+"""
+function flux_balance(ctmodel::C.ConstraintTree, optimizer; modifications = [])
+    opt_model = optimization_model(ctmodel; objective = ctmodel.objective.value, optimizer)
+
+    for mod in modifications
+        mod(opt_model)
+    end
+
+    J.optimize!(opt_model)
+
+    is_solved(opt_model) || return nothing
+
+    C.ValueTree(ctmodel, J.value.(opt_model[:x]))
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+Pipe-able variant of [`flux_balance`](@ref).
+"""
+flux_balance(optimizer; modifications = []) = m -> flux_balance(m, optimizer; modifications)
+
+export flux_balance
diff --git a/src/analysis/modifications.jl b/src/analysis/modifications.jl
new file mode 100644
index 000000000..585077480
--- /dev/null
+++ b/src/analysis/modifications.jl
@@ -0,0 +1,44 @@
+
+#TODO: at this point, consider renaming the whole thing to "settings"
+
+"""
+$(TYPEDSIGNATURES)
+
+Change the objective sense of optimization. Possible arguments are
+`JuMP.MAX_SENSE` and `JuMP.MIN_SENSE`.
+"""
+set_objective_sense(objective_sense) =
+    opt_model -> J.set_objective_sense(opt_model, objective_sense)
+
+export set_objective_sense
+
+"""
+$(TYPEDSIGNATURES)
+
+Change the JuMP optimizer used to run the optimization.
+"""
+set_optimizer(optimizer) = opt_model -> J.set_optimizer(opt_model, optimizer)
+
+export set_optimizer
+
+"""
+$(TYPEDSIGNATURES)
+
+Change a JuMP optimizer attribute. The attributes are optimizer-specific, refer
+to the JuMP documentation and the documentation of the specific optimizer for
+usable keys and values.
+"""
+set_optimizer_attribute(attribute_key, value) =
+    opt_model -> J.set_optimizer_attribute(opt_model, attribute_key, value)
+
+export set_optimizer_attribute
+
+"""
+    silence
+
+Modification that disable all output from the JuMP optimizer (shortcut for
+`set_silent` from JuMP).
+"""
+silence(opt_model) = J.set_silent(opt_model)
+
+export silence
diff --git a/src/analysis/parsimonious_flux_balance.jl b/src/analysis/parsimonious_flux_balance.jl
new file mode 100644
index 000000000..26d842a6e
--- /dev/null
+++ b/src/analysis/parsimonious_flux_balance.jl
@@ -0,0 +1,92 @@
+"""
+$(TYPEDSIGNATURES)
+
+Run parsimonious flux balance analysis (pFBA) on the `model`. In short, pFBA
+runs two consecutive optimization problems. The first is traditional FBA:
+```
+max cᵀx = μ
+s.t. S x = b
+     xₗ ≤ x ≤ xᵤ
+```
+And the second is a quadratic optimization problem:
+```
+min Σᵢ xᵢ²
+s.t. S x = b
+     xₗ ≤ x ≤ xᵤ
+     μ = μ⁰
+```
+Where the optimal solution of the FBA problem, μ⁰, has been added as an
+additional constraint. See "Lewis, Nathan E, Hixson, Kim K, Conrad, Tom M,
+Lerman, Joshua A, Charusanti, Pep, Polpitiya, Ashoka D, Adkins, Joshua N,
+Schramm, Gunnar, Purvine, Samuel O, Lopez-Ferrer, Daniel, Weitz, Karl K, Eils,
+Roland, König, Rainer, Smith, Richard D, Palsson, Bernhard Ø, (2010) Omic data
+from evolved E. coli are consistent with computed optimal growth from
+genome-scale models. Molecular Systems Biology, 6. 390. doi:
+accession:10.1038/msb.2010.47" for more details.
+
+pFBA gets the model optimum by standard FBA (using
+[`flux_balance`](@ref) with `optimizer` and `modifications`), then
+finds a minimal total flux through the model that still satisfies the (slightly
+relaxed) optimum. This is done using a quadratic problem optimizer. If the
+original optimizer does not support quadratic optimization, it can be changed
+using the callback in `qp_modifications`, which are applied after the FBA. See
+the documentation of [`flux_balance`](@ref) for usage examples of
+modifications.
+
+The optimum relaxation sequence can be specified in `relax` parameter, it
+defaults to multiplicative range of `[1.0, 0.999999, ..., 0.99]` of the original
+bound.
+
+Returns an optimized model that contains the pFBA solution (or an unsolved model
+if something went wrong).
+
+# Performance
+
+This implementation attempts to save time by executing all pFBA steps on a
+single instance of the optimization model problem, trading off possible
+flexibility. For slightly less performant but much more flexible use, one can
+construct parsimonious models directly using
+[`with_parsimonious_objective`](@ref).
+
+# Example
+```
+model = load_model("e_coli_core.json")
+parsimonious_flux_balance(model, biomass, Gurobi.Optimizer) |> values_vec
+```
+"""
+function parsimonious_flux_balance(
+    model::C.ConstraintTree,
+    optimizer;
+    modifications = [],
+    qp_modifications = [],
+    relax_bounds = [1.0, 0.999999, 0.99999, 0.9999, 0.999, 0.99],
+)
+    # Run FBA
+    opt_model = flux_balance(model, optimizer; modifications)
+    J.is_solved(opt_model) || return nothing # FBA failed
+
+    # get the objective
+    Z = J.objective_value(opt_model)
+    original_objective = J.objective_function(opt_model)
+
+    # prepare the model for pFBA
+    for mod in qp_modifications
+        mod(model, opt_model)
+    end
+
+    # add the minimization constraint for total flux
+    v = opt_model[:x] # fluxes
+    J.@objective(opt_model, Min, sum(dot(v, v)))
+
+    for rb in relax_bounds
+        # lb, ub = objective_bounds(rb)(Z)
+        J.@constraint(opt_model, pfba_constraint, lb <= original_objective <= ub)
+
+        J.optimize!(opt_model)
+        J.is_solved(opt_model) && break
+
+        J.delete(opt_model, pfba_constraint)
+        J.unregister(opt_model, :pfba_constraint)
+    end
+
+end
diff --git a/src/builders/core.jl b/src/builders/core.jl
index d85e1d71c..97cf6d62a 100644
--- a/src/builders/core.jl
+++ b/src/builders/core.jl
@@ -1,5 +1,4 @@
 
-import AbstractFBCModels as F
 import SparseArrays: sparse
 
 """
@@ -8,22 +7,22 @@ $(TYPEDSIGNATURES)
 A constraint tree that models the content of the given instance of
 `AbstractFBCModel`.
 """
-function fbc_model_constraints(model::F.AbstractFBCModel)
-    rxns = Symbol.(F.reactions(model))
-    mets = Symbol.(F.metabolites(model))
-    lbs, ubs = F.bounds(model)
-    stoi = F.stoichiometry(model)
-    bal = F.balance(model)
-    obj = F.objective(model)
+function fbc_model_constraints(model::A.AbstractFBCModel)
+    rxns = Symbol.(A.reactions(model))
+    mets = Symbol.(A.metabolites(model))
+    lbs, ubs = A.bounds(model)
+    stoi = A.stoichiometry(model)
+    bal = A.balance(model)
+    obj = A.objective(model)
 
     #TODO: is sparse() required below?
     return C.ConstraintTree(
-        :fluxes => C.variables(keys = rxns, bounds = zip(lbs, ubs)),
-        :balances => C.ConstraintTree(
-            m => C.Constraint(value = C.LinearValue(sparse(row)), bound = b) for
-            (m, row, b) in zip(mets, eachrow(stoi), bals)
-        ),
-        :objective => C.Constraint(value = C.Value(sparse(obj))),
+        :fluxes^C.variables(keys = rxns, bounds = zip(lbs, ubs)) *
+        :flux_stoichiometry^C.ConstraintTree(
+            met => C.Constraint(value = C.LinearValue(sparse(row)), bound = b) for
+            (met, row, b) in zip(mets, eachrow(stoi), bal)
+        ) *
+        :objective^C.Constraint(C.LinearValue(sparse(obj))),
     )
 end
 
@@ -70,11 +69,26 @@ sign_split_constraints(;
     signed::C.ConstraintTree,
 ) = C.ConstraintTree(
     k => C.Constraint(
-        value = s + (haskey(negative, k) ? C.value(negative[k]) : zero(C.Value)) -
-                (haskey(positive, k) ? C.value(positive[k]) : zero(C.Value)),
+        value = s.value +
+                (haskey(negative, k) ? negative[k].value : zero(typeof(s.value))) -
+                (haskey(positive, k) ? positive[k].value : zero(typeof(s.value))),
         bound = 0.0,
-    ) for (k, s) in C.elems(signed)
+    ) for (k, s) in signed
 )
 #TODO the example above might as well go to docs
 
 export sign_split_constraints
+
+function fluxes_in_direction(fluxes::C.ConstraintTree, direction = :forward)
+    keys = Symbol[]
+    for (id, flux) in fluxes
+        if direction == :forward
+            last(flux.bound) > 0 && push!(keys, id)
+        else
+            first(flux.bound) < 0 && push!(keys, id)
+        end
+    end
+    C.variables(; keys, bounds = Ref((0.0, Inf)))
+end
+
+export fluxes_in_direction
diff --git a/src/builders/genes.jl b/src/builders/genes.jl
index 2996a08f9..8cd1587dd 100644
--- a/src/builders/genes.jl
+++ b/src/builders/genes.jl
@@ -1,19 +1,42 @@
 
+"""
+$(TYPEDSIGNATURES)
+"""
 knockout_constraints(; fluxes::C.ConstraintTree, knockout_test) = C.ConstraintTree(
     rxn => C.Constraint(C.value(fluxes[rxn]), 0.0) for
     rxn in keys(fluxes) if knockout_test(rxn)
 )
 
-export knockout_constraints
-
-fbc_gene_knockout_constraints(;
+"""
+$(TYPEDSIGNATURES)
+"""
+gene_knockouts(;
     fluxes::C.ConstraintTree,
-    genes,
-    fbc_model::A.AbstractFBCModel,
+    ko_genes::Vector{String},
+    model::A.AbstractFBCModel,
 ) = knockout_constraints(;
     fluxes,
-    knockout_test = rxn ->
-        !A.reaction_gene_products_available(rxn, g -> not(g in genes)),
+    knockout_test = rxn -> begin
+        maybe_avail = A.reaction_gene_products_available(
+            model,
+            string(rxn),
+            g -> !(g in ko_genes), # not available if knocked out
+        )
+        isnothing(maybe_avail) ? false : !maybe_avail # negate here because of knockout_constraints
+    end,
 )
 
-export fbc_gene_knockout_constraints
+#TODO remove the bang from here, there's no side effect
+"""
+$(TYPEDSIGNATURES)
+"""
+knockout!(ctmodel::C.ConstraintTree, ko_genes::Vector{String}, model::A.AbstractFBCModel) =
+    ctmodel * :gene_knockouts^gene_knockouts(; fluxes = ctmodel.fluxes, ko_genes, model)
+
+"""
+$(TYPEDSIGNATURES)
+
+Pipe-able variant.
+"""
+knockout!(ko_genes::Vector{String}, model::A.AbstractFBCModel) =
+    m -> knockout!(m, ko_genes, model)
diff --git a/src/io.jl b/src/io.jl
index f517bdb58..b3b59f640 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -1,6 +1,6 @@
 
 """
-    $(TYPEDSIGNATURES)
+$(TYPEDSIGNATURES)
 
 Load a FBC model representation while guessing the correct model type. Uses
 `AbstractFBCModels.load`.
@@ -8,28 +8,39 @@ Load a FBC model representation while guessing the correct model type. Uses
 This overload almost always involves a search over types; do not use it in
 environments where performance is critical.
 """
-function load_fbc_model(path::String) where {A<:AbstractFBCModel}
+function load_model(path::String)
     A.load(path)
 end
 
 """
-    $(TYPEDSIGNATURES)
+$(TYPEDSIGNATURES)
 
 Load a FBC model representation. Uses `AbstractFBCModels.load`.
 """
-function load_fbc_model(model_type::Type{A}, path::String) where {A<:AbstractFBCModel}
+function load_model(model_type::Type{T}, path::String) where {T<:A.AbstractFBCModel}
     A.load(model_type, path)
 end
 
-export load_fbc_model
+export load_model
 
 """
-    $(TYPEDSIGNATURES)
+$(TYPEDSIGNATURES)
 
 Save a FBC model representation. Uses `AbstractFBCModels.save`.
 """
-function save_fbc_model(model::A, path::String) where {A<:AbstractFBCModel}
+function save_model(model::T, path::String) where {T<:A.AbstractFBCModel}
     A.save(model, path)
 end
 
-export save_fbc_model
+export save_model
+
+"""
+$(TYPEDSIGNATURES)
+
+Safely download a model with a known hash. All arguments are forwarded to
+`AbstractFBCModels.download_data_file` -- see the documentation in the
+AbstractFBCModels package for details.
+"""
+download_model(args...; kwargs...) = A.download_data_file(args...; kwargs...)
+
+export download_model
diff --git a/src/solver.jl b/src/solver.jl
index 7cd7d1027..c1a3b2b7b 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -47,36 +47,3 @@ is_solved(opt_model::J.Model) =
     J.termination_status(opt_model) in [J.MOI.OPTIMAL, J.MOI.LOCALLY_SOLVED]
 
 export is_solved
-
-"""
-$(TYPEDSIGNATURES)
-
-The optimized objective value of a JuMP model, if solved.
-"""
-optimized_objective_value(opt_model::J.Model)::Maybe{Float64} =
-    is_solved(opt_model) ? J.objective_value(opt_model) : nothing
-
-export optimized_objective_value
-
-"""
-$(TYPEDSIGNATURES)
-
-The optimized variable assignment of a JuMP model, if solved.
-"""
-optimized_variable_assignment(opt_model::J.Model)::Maybe{Vector{Float64}} =
-    is_solved(opt_model) ? J.value.(opt_model[:x]) : nothing
-
-export optimized_variable_assignment
-
-"""
-$(TYPEDSIGNATURES)
-
-Annotate a `ConstraintTree` with the values given by the optimization model,
-producing a `ValueTree` (if solved).
-"""
-solution(c::C.ConstraintTree, opt_model::J.Model)::Maybe{C.ValueTree} =
-    let vars = optimized_variable_assignment(opt_model)
-        isnothing(vars) ? nothing : C.ValueTree(c, vars)
-    end
-
-export solution
diff --git a/test/runtests.jl b/test/runtests.jl
index 0a4324e65..1bdc67ade 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,6 +6,10 @@ using Distributed
 import AbstractFBCModels as A
 using GLPK # for MILPs
 
+# testing constants
+const TEST_TOLERANCE = 1e-3
+const QP_TEST_TOLERANCE = 1e-2 # for Clarabel
+
 # helper functions for running tests en masse
 print_timing(fn, t) = @info "$(fn) done in $(round(t; digits = 2))s"
 
@@ -15,23 +19,30 @@ function run_test_file(path...)
     print_timing(fn, t)
 end
 
-function run_doc_ex(path...)
-    run_test_file("..", "docs", "src", "examples", path...)
+function run_doc_examples()
+    for dir in filter(endswith(".jl"), readdir("../docs/src/examples", join = true))
+        run_test_file(dir)
+    end
 end
 
 # set up the workers for Distributed, so that the tests that require more
 # workers do not unnecessarily load the stuff multiple times
 W = addprocs(2)
-t = @elapsed @everywhere using COBREXA, Tulip, JuMP
+t = @elapsed @everywhere begin
+    using COBREXA
+    import Tulip, JuMP
+end
 
 # load the test models
 run_test_file("data_static.jl")
 run_test_file("data_downloaded.jl")
 
-# import base files
+# TODO data_static and data_downloaded need to be interned into the demos.
+# Instead let's make a single "doc running directory" that runs all the
+# documentation, which doesn't get erased to improve the test caching.
+
 @testset "COBREXA test suite" begin
-    run_doc_ex("01-loading-and-saving.jl")
-    run_doc_ex("02-flux-balance-analysis.jl")
+    run_doc_examples()
     run_test_file("aqua.jl")
 end