Example: Metabolic modeling
In this example we demonstrate the use of ConstraintTree structure for solving the metabolic modeling tasks. At the same time, we show how to export the structure to JuMP, and use value trees to find useful information about the result.
First, let's import some packages:
import ConstraintTrees as CWe will need a constraint-based metabolic model; for this test we will use the usual "E. Coli core metabolism" model as available from BiGG:
import Downloads: download
download("http://bigg.ucsd.edu/static/models/e_coli_core.xml", "e_coli_core.xml")
import SBML
ecoli = SBML.readSBML("e_coli_core.xml")SBML.Model with 95 reactions, 72 species, and 5 parameters.Allocating and constraining variables
Let's first build the constrained representation of the problem. First, we will need a variable for each of the reactions in the model.
c = C.variables(keys = Symbol.(keys(ecoli.reactions)))ConstraintTree with 95 elements:
:R_ACALD => Constraint(LinearValue(#= ... =#))
:R_ACALDt => Constraint(LinearValue(#= ... =#))
:R_ACKr => Constraint(LinearValue(#= ... =#))
:R_ACONTa => Constraint(LinearValue(#= ... =#))
:R_ACONTb => Constraint(LinearValue(#= ... =#))
:R_ACt2r => Constraint(LinearValue(#= ... =#))
:R_ADK1 => Constraint(LinearValue(#= ... =#))
:R_AKGDH => Constraint(LinearValue(#= ... =#))
:R_AKGt2r => Constraint(LinearValue(#= ... =#))
:R_ALCD2x => Constraint(LinearValue(#= ... =#))
:R_ATPM => Constraint(LinearValue(#= ... =#))
:R_ATPS4r => Constraint(LinearValue(#= ... =#))
:R_BIOMASS_Ecoli_core_w_GAM => Constraint(LinearValue(#= ... =#))
:R_CO2t => Constraint(LinearValue(#= ... =#))
:R_CS => Constraint(LinearValue(#= ... =#))
:R_CYTBD => Constraint(LinearValue(#= ... =#))
:R_D_LACt2 => Constraint(LinearValue(#= ... =#))
:R_ENO => Constraint(LinearValue(#= ... =#))
:R_ETOHt2r => Constraint(LinearValue(#= ... =#))
⋮ => ⋮By default, Julia shows relatively long namespace prefixes before all identifiers, which substantially clutters the output. To improve the pretty-printing, all type names are marked as exported from the package, and you can import them via using ConstraintTrees.
The above operation returns a ConstraintTree. You can browse these as a dictionary:
c[:R_PFK]Constraint(LinearValue([49], [1.0]), nothing)...or much more conveniently using the record dot syntax as properties:
c.R_PFKConstraint(LinearValue([49], [1.0]), nothing)The individual LinearValues in constraints behave like sparse vectors that refer to variables: The first field represents the referenced variable indexes, and the second field represents the coefficients. Compared to the sparse vectors, information about the total number of variables is not stored explicitly.
Operator ^ is used to name individual constraints and directories in the hierarchy. Let us name our constraints as "fluxes" (which is a common name in metabolic modeling) and explore the result:
c = :fluxes^cConstraintTree with 1 element:
:fluxes => ConstraintTree(#= 95 elements =#)We can see that there is now only a single "top-level directory" in the constraint system, which can be explored with the dot access again:
c.fluxes.R_PFKConstraint(LinearValue([49], [1.0]), nothing)Indexing via values is again possible via the usual bracket notation, and can be freely combined with the dot notation:
c[:fluxes][:R_PFK]Constraint(LinearValue([49], [1.0]), nothing)Adding single-variable constraints
Each element in the constraint tree consists of a linear combination of the variables, which can be freely used to construct (and constraint) new linear combinations of variables. As the simplest use, we can constraint the variables (using Constraints) to their valid bounds as defined by the model:
rxn_constraints =
let rxn_bounds = Symbol.(keys(ecoli.reactions)) .=> zip(SBML.flux_bounds(ecoli)...)
C.ConstraintTree(
r => C.Constraint(value = c.fluxes[r].value, bound = (lb, ub)) for
(r, ((lb, _), (ub, _))) in rxn_bounds # SBML units are ignored for simplicity
)
endConstraintTree with 95 elements:
:R_ACALD => Constraint(LinearValue(#= ... =#), Between(-1000.0, 1000.0))
:R_ACALDt => Constraint(LinearValue(#= ... =#), Between(-1000.0, 1000.0))
:R_ACKr => Constraint(LinearValue(#= ... =#), Between(-1000.0, 1000.0))
:R_ACONTa => Constraint(LinearValue(#= ... =#), Between(-1000.0, 1000.0))
:R_ACONTb => Constraint(LinearValue(#= ... =#), Between(-1000.0, 1000.0))
:R_ACt2r => Constraint(LinearValue(#= ... =#), Between(-1000.0, 1000.0))
:R_ADK1 => Constraint(LinearValue(#= ... =#), Between(-1000.0, 1000.0))
:R_AKGDH => Constraint(LinearValue(#= ... =#), Between(0.0, 1000.0))
:R_AKGt2r => Constraint(LinearValue(#= ... =#), Between(-1000.0, 1000.0))
:R_ALCD2x => Constraint(LinearValue(#= ... =#), Between(-1000.0, 1000.0))
:R_ATPM => Constraint(LinearValue(#= ... =#), Between(8.39, 1000.0))
:R_ATPS4r => Constraint(LinearValue(#= ... =#), Between(-1000.0, 1000.0))
:R_BIOMASS_Ecoli_core_w_GAM => Constraint(LinearValue(#= ... =#), Between(0.0, 1000.0))
:R_CO2t => Constraint(LinearValue(#= ... =#), Between(-1000.0, 1000.0))
:R_CS => Constraint(LinearValue(#= ... =#), Between(0.0, 1000.0))
:R_CYTBD => Constraint(LinearValue(#= ... =#), Between(0.0, 1000.0))
:R_D_LACt2 => Constraint(LinearValue(#= ... =#), Between(-1000.0, 1000.0))
:R_ENO => Constraint(LinearValue(#= ... =#), Between(-1000.0, 1000.0))
:R_ETOHt2r => Constraint(LinearValue(#= ... =#), Between(-1000.0, 1000.0))
⋮ => ⋮Note that in the example we use a simplified Dict-like construction of the ConstraintTree. You might equivalently write the code as a product (using prod()) of constraints that are individually labeled using the ^ operator, but the direct dictionary construction is faster because it skips many intermediate steps, and looks much more like idiomatic Julia code.
To combine the constraint trees, we can make a nice directory for the constraints and add them to the tree using operator *. Making "products" of constraint trees combines the trees in a way that they share their variables. In particular, using the values from c.fluxes in the constraints within rxn_constraints here will constraint precisely the same variables (and thus values) as the ones in the original system.
c = c * :constraints^rxn_constraintsConstraintTree with 2 elements:
:constraints => ConstraintTree(#= 95 elements =#)
:fluxes => ConstraintTree(#= 95 elements =#)Our model representation now contains 2 "directories":
collect(keys(c))2-element Vector{Symbol}:
:constraints
:fluxesValue and constraint arithmetics
Values may be combined additively and multiplied by real constants; which allows us to easily create more complex linear combination of any values already occurring in the model:
3 * c.fluxes.R_PFK.value - c.fluxes.R_ACALD.value / 2LinearValue([49, 73], [3.0, -0.5])For simplicity, you can also scale whole constraints, but it is impossible to add them together because the meaning of the bounds would get broken:
(3 * c.fluxes.R_PFK, -c.fluxes.R_ACALD / 2)(Constraint(LinearValue([49], [3.0]), nothing), Constraint(LinearValue([73], [-0.5]), nothing))To process constraints in bulk, you may use C.value for easier access to values when making new constraints:
sum(C.value.(values(c.fluxes)))LinearValue([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 86, 87, 88, 89, 90, 91, 92, 93, 94, 95], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])Notably, ConstraintTrees provide their own implementation of sum which typically works faster when adding many Values together. The basic interface and results are otherwise the same as with the sum from Base:
C.sum(C.value.(values(c.fluxes)))LinearValue([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 86, 87, 88, 89, 90, 91, 92, 93, 94, 95], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])Since the sum from Base package is usually implemented as a left fold, it does not behave optimally when the temporary sub-results grow during the computation (and thus their addition becomes gradually slower). In turn, using the Base.sum for summing up LinearValues and QuadraticValues may take time quadratic in the number of added items. sum from ConstraintTrees uses a different addition order which reduces the amount of large items added together (implemented by "pairwise" preduce), and in works in almost-linear time in most cases.
Affine values
To simplify various modeling goals (mainly calculation of various kinds of "distances"), the values support inclusion of an affine element – the variable with index 0 is assumed to be the "affine unit", and its assigned value is fixed at 1.0.
To demonstrate, let's make a small system with 2 variables.
system = C.variables(keys = [:x, :y])ConstraintTree with 2 elements:
:x => Constraint(LinearValue(#= ... =#))
:y => Constraint(LinearValue(#= ... =#))To add an affine element to a LinearValue, simply add it as a Real number using +, as in the linear transformations below:
system =
:original_coords^system *
:transformed_coords^C.ConstraintTree(
:xt => C.Constraint(1 + system.x.value + 4 + system.y.value),
:yt => C.Constraint(0.1 * (3 - system.y.value)),
)
C.pretty(system)┬─original_coords
│ ╰─┬─x: 1.0*x[1]
│ ╰─y: 1.0*x[2]
╰─transformed_coords
╰─┬─xt: 5.0 + 1.0*x[1] + 1.0*x[2]
╰─yt: 0.30000000000000004 + -0.1*x[2]Adding combined constraints
Metabolic modeling relies on the fact that the total rates of any metabolite getting created and consumed by the reaction equals to zero (which corresponds to conservation of mass). We can now add corresponding "stoichiometric" network constraints by following the reactants and products in the SBML structure:
stoi_constraints = C.ConstraintTree(
Symbol(m) => C.Constraint(
value = -C.sum(
(
sr.stoichiometry * c.fluxes[Symbol(rid)].value for
(rid, r) in ecoli.reactions for sr in r.reactants if sr.species == m
),
init = zero(C.LinearValue), # sometimes the sums are empty
) + C.sum(
(
sr.stoichiometry * c.fluxes[Symbol(rid)].value for
(rid, r) in ecoli.reactions for sr in r.products if sr.species == m
),
init = zero(C.LinearValue),
),
bound = 0.0,
) for m in keys(ecoli.species)
);Let's have a closer look at one of the constraints:
stoi_constraints.M_acald_cConstraint(LinearValue([41, 73, 94], [1.0, -1.0, 1.0]), EqualTo(0.0))Again, we can label the stoichiometry properly and add it to the bigger model representation:
c = c * :stoichiometry^stoi_constraintsConstraintTree with 3 elements:
:constraints => ConstraintTree(#= 95 elements =#)
:fluxes => ConstraintTree(#= 95 elements =#)
:stoichiometry => ConstraintTree(#= 72 elements =#)Saving the objective
Constraint based models typically optimize a certain linear formula. Constraint trees do not support setting objectives (they are not constraints), but we can save the objective as a harmless unconstrained "constraint" that can be used later to refer to the objective more easily. We can save that information into the constraint system immediately:
c *=
:objective^C.Constraint(
C.sum(
c.fluxes[Symbol(rid)].value * coeff for
(rid, coeff) in (keys(ecoli.reactions) .=> SBML.flux_objective(ecoli)) if
coeff != 0.0;
init = 0.0,
),
)ConstraintTree with 4 elements:
:constraints => ConstraintTree(#= 95 elements =#)
:fluxes => ConstraintTree(#= 95 elements =#)
:objective => Constraint(LinearValue(#= ... =#))
:stoichiometry => ConstraintTree(#= 72 elements =#)Constrained system solutions and value trees
To aid exploration of variable assignments in the constraint trees, we can convert them to value trees. These have the very same structure as constraint trees, but carry only the "solved" constraint values instead of full constraints.
Let's demonstrate this quickly on the example of system with affine variables from above. First, let's assume that someone solved the system (in some way) and produced a solution of variables as follows:
solution = [1.0, 5.0] # corresponds to :x and :y in order given in `variables`2-element Vector{Float64}:
1.0
5.0A value tree for this solution is constructed in a straightforward manner:
st = C.substitute_values(system, solution)Tree{Float64} with 2 elements:
:original_coords => Tree{Float64}(#= 2 elements =#)
:transformed_coords => Tree{Float64}(#= 2 elements =#)We can now check the values of the "original" coordinates:
st.original_coordsTree{Float64} with 2 elements:
:x => 1.0
:y => 5.0The more complex constraints get their values computed automatically from the variable assignment:
st.transformed_coordsTree{Float64} with 2 elements:
:xt => 11.0
:yt => -0.2Solving the constraint system using JuMP
We can make a small function that throws our model into JuMP, optimizes it, and gives us back a variable assignment vector. This vector can then be used to determine and browse the values of constraints and variables using a Float64-valued tree.
import JuMP
function optimized_vars(cs::C.ConstraintTree, objective::C.LinearValue, optimizer)
model = JuMP.Model(optimizer)
JuMP.@variable(model, x[1:C.variable_count(cs)])
JuMP.@objective(model, JuMP.MAX_SENSE, C.substitute(objective, x))
C.traverse(cs) do c
b = c.bound
if b isa C.EqualTo
JuMP.@constraint(model, C.substitute(c.value, x) == b.equal_to)
elseif b isa C.Between
val = C.substitute(c.value, x)
isinf(b.lower) || JuMP.@constraint(model, val >= b.lower)
isinf(b.upper) || JuMP.@constraint(model, val <= b.upper)
end
end
JuMP.optimize!(model)
JuMP.value.(model[:x])
endoptimized_vars (generic function with 1 method)With this in hand, we can use an external linear problem solver to find the optimum of the constrained system:
import GLPK
optimal_variable_assignment = optimized_vars(c, c.objective.value, GLPK.Optimizer)95-element Vector{Float64}:
0.0
6.00724957535031
7.477381962160286
-5.064375661482085
0.22346172933182415
-3.2148950476847733
2.5043094703687165
21.799492655998776
4.959984944574683
1.4969837572615767
⋮
0.0
29.175827135565836
0.0
0.0
-6.400958585906926e-15
-21.799492655998776
-1.776356839400254e-15
0.0
3.2148950476847733To explore the solution more easily, we can make a tree with values that correspond to ones in our constraint tree:
result = C.substitute_values(c, optimal_variable_assignment)
result.fluxes.R_BIOMASS_Ecoli_core_w_GAM0.8739215069684383result.fluxes.R_PFK7.477381962160286Sometimes it is unnecessary to recover the values for all constraints, so we are better off selecting just the right subtree:
C.substitute_values(c.fluxes, optimal_variable_assignment)Tree{Float64} with 95 elements:
:R_ACALD => -1.98196e-15
:R_ACALDt => -1.98196e-15
:R_ACKr => -2.38143e-14
:R_ACONTa => 6.00725
:R_ACONTb => 6.00725
:R_ACt2r => -2.38143e-14
:R_ADK1 => 8.88178e-16
:R_AKGDH => 5.06438
:R_AKGt2r => 0.0
:R_ALCD2x => 0.0
:R_ATPM => 8.39
:R_ATPS4r => 45.514
:R_BIOMASS_Ecoli_core_w_GAM => 0.873922
:R_CO2t => -22.8098
:R_CS => 6.00725
:R_CYTBD => 43.599
:R_D_LACt2 => 1.95151e-14
:R_ENO => 14.7161
:R_ETOHt2r => 0.0
⋮ => ⋮C.substitute_values(c.objective, optimal_variable_assignment)0.8739215069684383We'll save the result for future use at the end of this example:
result_single_organism = resultTree{Float64} with 4 elements:
:constraints => Tree{Float64}(#= 95 elements =#)
:fluxes => Tree{Float64}(#= 95 elements =#)
:objective => 0.873922
:stoichiometry => Tree{Float64}(#= 72 elements =#)Combining and extending constraint systems
Constraint trees can be extended with new variables from another constraint trees using the + operator. Contrary to the * operator, adding the constraint trees does not share the variables between operands, and the resulting constraint tree will basically contain two disconnected trees that solve independently. The user is expected to create additional constraints to connect the independent parts.
Here, we demonstrate this by creating a community of two slightly different E. Coli species: First, we disable functionality of a different reaction in each of the models to create a diverse group of differently handicapped organisms:
c =
:community^(
:species1^(c * :handicap^C.Constraint(c.fluxes.R_PFK.value, 0)) +
:species2^(c * :handicap^C.Constraint(c.fluxes.R_ACALD.value, 0))
)ConstraintTree with 1 element:
:community => ConstraintTree(#= 2 elements =#)We can create additional variables that represent total community intake of oxygen, and total community production of biomass:
c += :exchanges^C.variables(keys = [:oxygen, :biomass], bounds = [(-10.0, 10.0), nothing])ConstraintTree with 2 elements:
:community => ConstraintTree(#= 2 elements =#)
:exchanges => ConstraintTree(#= 2 elements =#)These can be constrained so that the total influx (or outflux) of each of the registered metabolites is in fact equal to total consumption or production by each of the species:
c *=
:exchange_constraints^C.ConstraintTree(
:oxygen => C.Constraint(
value = c.exchanges.oxygen.value - c.community.species1.fluxes.R_EX_o2_e.value -
c.community.species2.fluxes.R_EX_o2_e.value,
bound = 0.0,
),
:biomass => C.Constraint(
value = c.exchanges.biomass.value -
c.community.species1.fluxes.R_BIOMASS_Ecoli_core_w_GAM.value -
c.community.species2.fluxes.R_BIOMASS_Ecoli_core_w_GAM.value,
bound = 0.0,
),
)ConstraintTree with 3 elements:
:community => ConstraintTree(#= 2 elements =#)
:exchange_constraints => ConstraintTree(#= 2 elements =#)
:exchanges => ConstraintTree(#= 2 elements =#)Let's see how much biomass are the two species capable of producing together:
result =
C.substitute_values(c, optimized_vars(c, c.exchanges.biomass.value, GLPK.Optimizer))
result.exchangesTree{Float64} with 2 elements:
:biomass => 0.469658
:oxygen => -10.0Finally, we can iterate over all species in the small community and see how much biomass was actually contributed by each:
Dict(k => v.fluxes.R_BIOMASS_Ecoli_core_w_GAM for (k, v) in result.community)Dict{Symbol, Float64} with 2 entries:
:species2 => 0.442232
:species1 => 0.027425Modifying constraint systems in-place
Constraint trees can be modified in-place in a way that allows you to easily change small values in the trees without reconstructing them from the ground up.
Although in-place modification is extremely convenient and looks much easier than rebuilding the tree, it may be very detrimental to the robustness and efficiency of the programs, for several reasons:
- changing any data breaks assumptions on anything that was already derived from the data
- for efficiency, the tree structures are not copied by default if there's no need to do it, and only shared by references; which means that a naive change at a single place of the tree may easily change values also in other parts of any trees, including completely different trees
- the "convenient way" of making sure that the above problem never happens is to copy-on-write the whole tree structure, which is typically quite detrimental to memory use and program efficiency
Only use the in-place modifications if:
- there is code that explicitly makes sure there is no false sharing via references, e.g. using a deep copy
- the in-place modifications are the last thing happening to the constraint tree before it is used by the solver
- the in-place modification code is not a part of a re-usable library
- you are using a suitable wrapper interface such as Accessors.jl
Now, if you are completely sure that ignoring the robustness guidelines will help your code, you can do the in-place tree modifications quite easily using both dot-access and array-index syntax.
You can thus, e.g., set a single bound:
c.exchanges.oxygen.bound = C.Between(-20.0, 20.0)Between(-20.0, 20.0)...or rebuild a whole constraint (using a tuple shortcut for Between):
c.exchanges.biomass = C.Constraint(c.exchanges.biomass.value, (-20, 20))Constraint(LinearValue([192], [1.0]), Between(-20.0, 20.0))...or even add new constraints, here using the index syntax for demonstration:
c[:exchanges][:production_is_zero] = C.Constraint(c.exchanges.biomass.value, 0)Constraint(LinearValue([192], [1.0]), EqualTo(0.0))...or remove some constraints (this erases the constraint that was added just above):
delete!(c.exchanges, :production_is_zero)DataStructures.SortedDict{Symbol, Union{Constraint, ConstraintTree}, Base.Order.ForwardOrdering} with 2 entries:
:biomass => Constraint(LinearValue(#= ... =#), Between(-20.0, 20.0))
:oxygen => Constraint(LinearValue(#= ... =#), Between(-20.0, 20.0))In the end, the flux optimization yields an expectably different result:
result_with_more_oxygen =
C.substitute_values(c, optimized_vars(c, c.exchanges.biomass.value, GLPK.Optimizer))
result.exchangesTree{Float64} with 2 elements:
:biomass => 0.469658
:oxygen => -10.0Alternative: Using Accessors.jl
Accessors.jl implement a "lensy" way to update immutable data structures. That comes with a nice outcome of doing the right amount of shallow copies for you automatically, thus avoiding much of the technical danger of in-place modifications. (You still lose the equational reasoning on your code, but that may not be an issue at all in usual codebases.)
Accessors interface is used simply through macros @set (which sets a deeply nested field in a structure, returning a modified copy), or with @reset which automatically "assigns" the result back to the original variable:
using Accessors
c = @set c.exchanges.biomass.bound = C.Between(-50, 50)ConstraintTree with 3 elements:
:community => ConstraintTree(#= 2 elements =#)
:exchange_constraints => ConstraintTree(#= 2 elements =#)
:exchanges => ConstraintTree(#= 2 elements =#)The above code is equivalent to:
@reset c.exchanges.biomass.bound = C.Between(-50, 50)ConstraintTree with 3 elements:
:community => ConstraintTree(#= 2 elements =#)
:exchange_constraints => ConstraintTree(#= 2 elements =#)
:exchanges => ConstraintTree(#= 2 elements =#)...and it is also possible to use string and symbol indexes to pick the individual tree items:
@reset c[:exchanges]["biomass"].bound = C.Between(-50, 50)ConstraintTree with 3 elements:
:community => ConstraintTree(#= 2 elements =#)
:exchange_constraints => ConstraintTree(#= 2 elements =#)
:exchanges => ConstraintTree(#= 2 elements =#)All of these operations give us:
c.exchanges.biomassConstraint(LinearValue([192], [1.0]), Between(-50.0, 50.0))Seeing the differences between the results
ConstraintTrees.jl defines its own version of zip function that can apply a function to the contents of several trees, "zipping" them over the same keys in the structure. This is vaguely similar but otherwise not related to the zip from Julia base (similarly, ConstraintTrees.jl have their own specific map).
In practice, this allows you to create combined trees with various nice properties very quickly. For example, you can find how much the values have changed between our two communities:
C.zip((x, y) -> y - x, result, result_with_more_oxygen, Float64)Tree{Float64} with 3 elements:
:community => Tree{Float64}(#= 2 elements =#)
:exchange_constraints => Tree{Float64}(#= 2 elements =#)
:exchanges => Tree{Float64}(#= 2 elements =#)The result is again a Tree, with the contained type specified by the last argument (Float64 in this case). We can explore it right away as the other result trees. Also, it is possible to call this kind of function using the Julia do notation, making the syntax a bit neater:
difference = C.zip(result, result_with_more_oxygen, Float64) do x, y
y - x
endTree{Float64} with 3 elements:
:community => Tree{Float64}(#= 2 elements =#)
:exchange_constraints => Tree{Float64}(#= 2 elements =#)
:exchanges => Tree{Float64}(#= 2 elements =#)Exploring the difference works as expected:
difference.community.species1.fluxesTree{Float64} with 95 elements:
:R_ACALD => 0.0
:R_ACALDt => 0.0
:R_ACKr => 0.614082
:R_ACONTa => 0.272252
:R_ACONTb => 0.272252
:R_ACt2r => 0.614082
:R_ADK1 => 1.65141
:R_AKGDH => 1.70003e-16
:R_AKGt2r => 0.0
:R_ALCD2x => 0.0
:R_ATPM => 0.0
:R_ATPS4r => 16.5237
:R_BIOMASS_Ecoli_core_w_GAM => 0.252343
:R_CO2t => -6.96009
:R_CS => 0.272252
:R_CYTBD => 14.6871
:R_D_LACt2 => 0.0
:R_ENO => 2.17284
:R_ETOHt2r => 0.0
⋮ => ⋮For convenience in special cases, zip is also overloaded for 3 arguments. We can, for a completely artificial example, check if the absolute flux change was bigger in the first or in the second organism in the community when compared to the original single-organism flux (which we luckily saved above):
changes = C.zip(
result.community.species1,
result.community.species2,
result_single_organism,
Bool,
) do s1, s2, orig
abs(s1 - orig) > abs(s2 - orig)
end
changes.fluxesTree{Bool} with 95 elements:
:R_ACALD => false
:R_ACALDt => false
:R_ACKr => false
:R_ACONTa => true
:R_ACONTb => true
:R_ACt2r => false
:R_ADK1 => true
:R_AKGDH => false
:R_AKGt2r => false
:R_ALCD2x => false
:R_ATPM => false
:R_ATPS4r => false
:R_BIOMASS_Ecoli_core_w_GAM => true
:R_CO2t => false
:R_CS => true
:R_CYTBD => true
:R_D_LACt2 => false
:R_ENO => true
:R_ETOHt2r => false
⋮ => ⋮More high-level functions like zip are described in an example on functional tree processing.
Simplifying the modified constraint systems
We may notice that some constraints that we placed on the community make solving of some variables quite trivial. For example, we knocked out PFK in species1:
c.community.species1.handicapConstraint(LinearValue([49], [1.0]), EqualTo(0.0))Obviously, the variable at index 49 must be equal to zero in the whole model, and we could completely remove it from the model. Unfortunately, just removing the constraint will not help: The variable is also referenced from elsewhere (e.g., from stoichiometry) and removing the constraints actually creates additional feasible solutions for the model, leaving us with possibly invalid solutions!
Instead of that, we can substitute a literal zero value into the variable in the whole model, and then use a pruning function to get rid of the (now unused) variable index.
First, let's create a list of all variables in the model that we can use for substitution:
vars = [C.variable(; idx).value for idx = 1:C.variable_count(c)]192-element Vector{LinearValue}:
LinearValue([1], [1.0])
LinearValue([2], [1.0])
LinearValue([3], [1.0])
LinearValue([4], [1.0])
LinearValue([5], [1.0])
LinearValue([6], [1.0])
LinearValue([7], [1.0])
LinearValue([8], [1.0])
LinearValue([9], [1.0])
LinearValue([10], [1.0])
⋮
LinearValue([184], [1.0])
LinearValue([185], [1.0])
LinearValue([186], [1.0])
LinearValue([187], [1.0])
LinearValue([188], [1.0])
LinearValue([189], [1.0])
LinearValue([190], [1.0])
LinearValue([191], [1.0])
LinearValue([192], [1.0])Now we use a bit of the knowledge about the model structure – the handicaps constraint single variables, so we can substitute for them directly. (If the handicaps constrained larger linear combinations of variables, we would have to resort to algebra.)
vars[c.community.species1.handicap.value.idxs[1]] = zero(C.LinearValue)
vars[c.community.species2.handicap.value.idxs[1]] = zero(C.LinearValue)LinearValue(Int64[], Float64[])substitute can be used to feed these new variables into all variables in the existing constraint system, which we immediately follow by pruning of the variables (pruning is described closer in the functional tree processing example).
c_simplified = C.prune_variables(C.substitute(c, vars))ConstraintTree with 3 elements:
:community => ConstraintTree(#= 2 elements =#)
:exchange_constraints => ConstraintTree(#= 2 elements =#)
:exchanges => ConstraintTree(#= 2 elements =#)The result contains exactly 2 variables less than the original community:
(C.variable_count(c), C.variable_count(c_simplified))(192, 190)The constraints that were substituted for are, at this point, roughly equivalent to saying 0 == 0, and most solvers will simply drop them as tautologies:
c_simplified.community.species1.handicapConstraint(LinearValue(Int64[], Float64[]), EqualTo(0.0))Finally, the variables are properly substituted for in the other equations, including the stoichiometry: We can compare e.g. the original balance of metabolite f6p (which is consumed by PFK):
c.community.species1.stoichiometry.M_f6p_cConstraint(LinearValue([10, 42, 49, 72, 77, 82, 93], [1.0, 1.0, -1.0, 1.0, 1.0, -0.0709, 1.0]), EqualTo(0.0))...to the simplified stoichiometry, which no longer refers to the PFK reaction:
c_simplified.community.species1.stoichiometry.M_f6p_cConstraint(LinearValue([10, 42, 71, 76, 81, 92], [1.0, 1.0, 1.0, 1.0, -0.0709, 1.0]), EqualTo(0.0))This page was generated using Literate.jl.