# Analysis functions

## Common analysis functions

`COBREXA.envelope_lattice`

— Method```
envelope_lattice(model::MetabolicModel, ridxs::Vector{Int64}; samples, ranges, reaction_samples) -> Base.Generator{_A, COBREXA.var"#452#453"} where _A
```

Create a lattice (list of "tick" vectors) for reactions at indexes `ridxs`

in a model. Arguments `samples`

, `ranges`

, and `reaction_samples`

may be optionally specified to customize the lattice creation process.

`COBREXA.envelope_lattice`

— Method```
envelope_lattice(model::MetabolicModel, rids::Vector{String}; kwargs...) -> Base.Generator{_A, COBREXA.var"#452#453"} where _A
```

Version of `envelope_lattice`

that works on string reaction IDs instead of integer indexes.

`COBREXA.objective_envelope`

— Method```
objective_envelope(model::MetabolicModel, ridxs::Vector{Int64}, optimizer; modifications, lattice_args, lattice, analysis, kwargs...) -> NamedTuple{(:lattice, :values), _A} where _A<:Tuple{Any, Array}
```

Compute an array of objective values for the `model`

for rates of reactions specified `ridxs`

fixed to a regular range of values between their respective lower and upper bounds.

This can be used to compute a "production envelope" of a metabolite; but generalizes to any specifiable objective and to multiple dimensions of the examined space. For example, to retrieve a production envelope of any metabolite, set the objective coefficient vector of the `model`

to a vector that contains a single `1`

for the exchange reaction that "outputs" this metabolite, and run `objective_envelope`

with the exchange reaction of the "parameter" metabolite specified in `ridxs`

.

Returns a named tuple that contains `lattice`

with reference values of the metabolites, and an N-dimensional array `values`

with the computed objective values, where N is the number of specified reactions. Because of the increasing dimensionality, the computation gets very voluminous with increasing length of `ridxs`

. The `lattice`

for computing the optima can be supplied in the argument; by default it is created by `envelope_lattice`

called on the model and reaction indexes. Additional arguments for the call to `envelope_lattice`

can be optionally specified in `lattice_args`

.

`kwargs`

are internally forwarded to `screen_optmodel_modifications`

. `modifications`

are appended to the list of modifications after the lattice bounds are set. By default, this returns the objective values for all points in the lattice; alternate outputs can be implemented via the `analysis`

argument.

**Example**

```
julia> m = load_model("e_coli_core.xml");
julia> envelope = objective_envelope(m, ["R_EX_gln__L_e", "R_EX_fum_e"],
Tulip.Optimizer;
lattice_args=(samples=6,));
julia> envelope.lattice # the reaction rates for which the optima were computed
2-element Vector{Vector{Float64}}:
[0.0, 200.0, 400.0, 600.0, 800.0, 1000.0]
[0.0, 200.0, 400.0, 600.0, 800.0, 1000.0]
julia> envelope.values # the computed flux objective values for each reaction rate combination
6×6 Matrix{Float64}:
0.873922 9.25815 17.4538 19.56 20.4121 20.4121
13.0354 17.508 19.9369 21.894 22.6825 22.6825
16.6666 18.6097 20.2847 21.894 22.6825 22.6825
16.6666 18.6097 20.2847 21.894 22.6825 22.6825
16.6666 18.6097 20.2847 21.894 22.6825 22.6825
16.6666 18.6097 20.2847 21.894 22.6825 22.6825
```

`COBREXA.objective_envelope`

— Method```
objective_envelope(model::MetabolicModel, rids::Vector{String}, args...; kwargs...) -> NamedTuple{(:lattice, :values), _A} where _A<:Tuple{Any, Array}
```

Version of `objective_envelope`

that works on string reaction IDs instead of integer indexes.

`COBREXA.flux_balance_analysis`

— Method```
flux_balance_analysis(model::MetabolicModel, optimizer; modifications) -> JuMP.Model
```

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 `change_optimizer_attribute`

, `change_objective`

, and `change_sense`

.

Returns an optimized `JuMP`

model.

**Example**

```
model = load_model("e_coli_core.json")
solution = flux_balance_analysis(model, GLPK.optimizer)
value.(solution[:x]) # extract flux steady state from the optimizer
biomass_reaction_id = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM")
modified_solution = flux_balance_analysis(model, GLPK.optimizer;
modifications=[change_objective(biomass_reaction_id)])
```

`COBREXA.flux_balance_analysis_dict`

— Method```
flux_balance_analysis_dict(model::MetabolicModel, args...; kwargs...) -> Union{Nothing, Dict{String, Float64}}
```

A variant of FBA that returns a dictionary assigning fluxes to reactions, if the solution is found. Arguments are passed to `flux_balance_analysis`

.

This function is kept for backwards compatibility, use `flux_dict`

instead.

`COBREXA.flux_balance_analysis_vec`

— Method```
flux_balance_analysis_vec(model::MetabolicModel, args...; kwargs...) -> Union{Nothing, Vector{Float64}}
```

A variant of FBA that returns a vector of fluxes in the same order as reactions of the model, if the solution is found.

Arguments are passed to `flux_balance_analysis`

.

This function is kept for backwards compatibility, use `flux_vector`

instead.

`COBREXA._max_variability_flux`

— Method```
_max_variability_flux(opt_model, flux, sense, ret) -> Any
```

Internal helper for maximizing reactions in optimization model.

`COBREXA.flux_variability_analysis`

— Method```
flux_variability_analysis(model::MetabolicModel, optimizer; kwargs...) -> Array
```

A simpler version of `flux_variability_analysis`

that maximizes and minimizes all declared fluxes in the model. Arguments are forwarded.

`COBREXA.flux_variability_analysis`

— Method```
flux_variability_analysis(model::MetabolicModel, fluxes::SparseArrays.SparseMatrixCSC{Float64, Int64}, optimizer; modifications, workers, optimal_objective_value, bounds, ret) -> Array
```

Flux variability analysis solves a pair of optimization problems in `model`

for each flux `f`

described in `fluxes`

:

```
min,max fᵀxᵢ
s.t. S x = b
xₗ ≤ x ≤ xᵤ
cᵀx ≥ bounds(Z₀)[1]
cᵀx ≤ bounds(Z₀)[2]
```

where Z₀:= cᵀx₀ is the objective value of an optimal solution of the associated FBA problem (see `flux_balance_analysis`

for a related analysis, also for explanation of the `modifications`

argument).

The `bounds`

is a user-supplied function that specifies the objective bounds for the variability optimizations, by default it restricts the flux objective value to the precise optimum reached in FBA. It can return `-Inf`

and `Inf`

in first and second pair to remove the limit. Use `gamma_bounds`

and `objective_bounds`

for simple bounds.

`optimizer`

must be set to a `JuMP`

-compatible optimizer. The computation of the individual optimization problems is transparently distributed to `workers`

(see `Distributed.workers()`

). The value of Z₀ can be optionally supplied in argument `optimal_objective_value`

, which prevents this function from calling the non-parallelizable FBA. Separating the single-threaded FBA and multithreaded variability computation can be used to improve resource allocation efficiency in many common use-cases.

`ret`

is a function used to extract results from optimized JuMP models of the individual fluxes. By default, it calls and returns the value of `JuMP.objective_value`

. More information can be extracted e.g. by setting it to a function that returns a more elaborate data structure; such as `m -> (JuMP.objective_value(m), JuMP.value.(m[:x]))`

.

Returns a matrix of extracted `ret`

values for minima and maxima, of total size (`size(fluxes,2)`

,2). The optimizer result status is checked with `is_solved`

; `nothing`

is returned if the optimization failed for any reason.

**Example**

```
model = load_model("e_coli_core.json")
flux_variability_analysis(model, GLPK.optimizer)
```

`COBREXA.flux_variability_analysis`

— Method```
flux_variability_analysis(model::MetabolicModel, flux_indexes::Vector{Int64}, optimizer; kwargs...) -> Array
```

An overload of `flux_variability_analysis`

that explores the fluxes specified by integer indexes

`COBREXA.flux_variability_analysis_dict`

— Method```
flux_variability_analysis_dict(model::MetabolicModel, optimizer; kwargs...) -> Tuple{Any, Any}
```

A variant of `flux_variability_analysis`

that returns the individual maximized and minimized fluxes as two dictionaries (of dictionaries). All keyword arguments except `ret`

are passed through.

**Example**

```
mins, maxs = flux_variability_analysis_dict(
model,
Tulip.Optimizer;
bounds = objective_bounds(0.99),
modifications = [
change_optimizer_attribute("IPM_IterationsLimit", 500),
change_constraint("EX_glc__D_e"; lb = -10, ub = -10),
change_constraint("EX_o2_e"; lb = 0, ub = 0),
],
)
```

`COBREXA.reaction_variability_analysis`

— Method```
reaction_variability_analysis(model::MetabolicModel, optimizer; kwargs...) -> Array
```

Shortcut for `reaction_variability_analysis`

that examines all reactions.

`COBREXA.reaction_variability_analysis`

— Method```
reaction_variability_analysis(model::MetabolicModel, reaction_indexes::Vector{Int64}, optimizer; kwargs...) -> Array
```

A variant for `flux_variability_analysis`

that examines actual reactions (selected by their indexes in `reactions`

argument) instead of whole fluxes. This may be useful for models where the sets of reactions and fluxes differ.

`COBREXA.make_gecko_model`

— Method```
make_gecko_model(model::MetabolicModel; reaction_isozymes, gene_product_bounds, gene_product_molar_mass, gene_product_mass_group, gene_product_mass_group_bound)
```

Wrap a model into a `GeckoModel`

, following the structure given by GECKO algorithm (see `GeckoModel`

documentation for details).

**Arguments**

`reaction_isozymes`

is a function that returns a vector of`Isozyme`

s for each reaction, or empty vector if the reaction is not enzymatic.`gene_product_bounds`

is a function that returns lower and upper bound for concentration for a given gene product (specified by the same string gene ID as in`reaction_isozymes`

), as`Tuple{Float64,Float64}`

.`gene_product_molar_mass`

is a function that returns a numeric molar mass of a given gene product specified by string gene ID.`gene_product_mass_group`

is a function that returns a string group identifier for a given gene product, again specified by string gene ID. By default, all gene products belong to group`"uncategorized"`

which is the behavior of original GECKO.`gene_product_mass_group_bound`

is a function that returns the maximum mass for a given mass group.

Alternatively, all function arguments may be also supplied as dictionaries that provide the same data lookup.

`COBREXA._mmdf_add_df_bound`

— Method```
_mmdf_add_df_bound(lb, ub) -> COBREXA.var"#538#539"
```

Helper function to add a new constraint on the driving force.

`COBREXA._mmdf_concen_objective`

— Method```
_mmdf_concen_objective(midx, sense) -> COBREXA.var"#536#537"
```

Helper function to change the objective to optimizing some concentration.

`COBREXA._mmdf_dgr_objective`

— Method```
_mmdf_dgr_objective(ridx, sense) -> COBREXA.var"#534#535"
```

Helper function to change the objective to optimizing some dG.

`COBREXA.max_min_driving_force`

— Method```
max_min_driving_force(model::MetabolicModel, reaction_standard_gibbs_free_energies::Dict{String, Float64}, optimizer; flux_solution, proton_ids, water_ids, constant_concentrations, concentration_ratios, concentration_lb, concentration_ub, T, R, small_flux_tol, modifications, ignore_reaction_ids) -> Union{Nothing, NamedTuple{(:mmdf, :dg_reactions, :concentrations), _A} where _A<:Tuple{Any, Dict, Dict}}
```

Perform a max-min driving force analysis on the `model`

, as defined by Noor, et al., "Pathway thermodynamics highlights kinetic obstacles in central metabolism.", PLoS computational biology, 2014.

The function uses the supplied `optimizer`

and `reaction_standard_gibbs_free_energies`

. Optionally, `flux_solution`

can be used to set the directions of each reaction in `model`

(all reactions are assumed to proceed forward and are active by default). The supplied `flux_solution`

should be free of internal cycles i.e. thermodynamically consistent. This optional input is important if a reaction in `model`

normally runs in reverse (negative flux). Note, reactions in `flux_solution`

that are smaller than `small_flux_tol`

are also ignored in the analysis function (for numerical stability).

The max-min driving force algorithm returns the Gibbs free energy of the reactions, the concentrations of metabolites and the actual maximum minimum driving force. The optimization problem solved is:

```
max min -ΔᵣG
s.t. ΔᵣG = ΔrG⁰ + R T S' ln(C)
ΔᵣG ≤ 0
ln(Cₗ) ≤ ln(C) ≤ ln(Cᵤ)
```

where `ΔrG`

are the Gibbs energies dissipated by the reactions, R is the gas constant, T is the temperature, S is the stoichiometry of the model, and C is the vector of metabolite concentrations (and their respective lower and upper bounds).

In case no feasible solution exists, `nothing`

is returned.

Reactions specified in `ignore_reaction_ids`

are internally ignored when calculating the max-min driving force. This should include water and proton importers.

Since biochemical thermodynamics are assumed, the `proton_ids`

and `water_ids`

need to be specified so that they can be ignored in the calculations. Effectively this assumes an aqueous environment at constant pH is used.

`constant_concentrations`

is used to fix the concentrations of certain metabolites (such as CO₂). `concentration_ratios`

is used to specify additional constraints on metabolite pair concentrations (typically, this is done with various cofactors such as the ATP/ADP ratio. For example, you can fix the concentration of ATP to be always 5× higher than of ADP by specifying `Dict(("ATP","ADP") => 5.0)`

`concentration_lb`

and `concentration_ub`

set the `Cₗ`

and `Cᵤ`

in the optimization problems.

`T`

and `R`

can be specified in the corresponding units; defaults are K and kJ/K/mol.

`COBREXA.max_min_driving_force_variability`

— Method```
max_min_driving_force_variability(model::MetabolicModel, reaction_standard_gibbs_free_energies::Dict{String, Float64}, optimizer; workers, optimal_objective_value, bounds, modifications, kwargs...) -> Array
```

Perform a variant of flux variability analysis on a max min driving force type problem. Arguments are forwarded to `max_min_driving_force`

. Calls `screen`

internally and possibly distributes computation across `workers`

. If `optimal_objective_value = nothing`

, the function first performs regular max min driving force analysis to find the max min driving force of the model and sets this to `optimal_objective_value`

. Then iteratively maximizes and minimizes the driving force across each reaction, and then the concentrations while staying close to the original max min driving force as specified in `bounds`

.

The `bounds`

is a user-supplied function that specifies the max min driving force bounds for the variability optimizations, by default it restricts the flux objective value to the precise optimum reached in the normal max min driving force analysis. It can return `-Inf`

and `Inf`

in first and second pair to remove the limit. Use `gamma_bounds`

and `objective_bounds`

for simple bounds.

Returns a matrix of solutions to `max_min_driving_force`

additionally constrained as described above, where the rows are in the order of the reactions and then the metabolites of the `model`

. For the reaction rows the first column is the maximum dG of that reaction, and the second column is the minimum dG of that reaction subject to the above constraints. For the metabolite rows, the first column is the maximum concentration, and the second column is the minimum concentration subject to the constraints above.

`COBREXA.minimize_metabolic_adjustment`

— Method```
minimize_metabolic_adjustment(flux_ref_dict::Dict{String, Float64}) -> COBREXA.var"#543#545"{Dict{String, Float64}}
```

Overload of `minimize_metabolic_adjustment`

that works with a dictionary of fluxes.

`COBREXA.minimize_metabolic_adjustment`

— Method```
minimize_metabolic_adjustment(flux_ref::Vector{Float64}) -> COBREXA.var"#541#542"{Vector{Float64}}
```

An optimization model modification that implements the MOMA in `minimize_metabolic_adjustment_analysis`

.

`COBREXA.minimize_metabolic_adjustment_analysis`

— Method```
minimize_metabolic_adjustment_analysis(model::MetabolicModel, flux_ref::Union{Dict{String, Float64}, Vector{Float64}}, optimizer; modifications, kwargs...) -> JuMP.Model
```

Run minimization of metabolic adjustment (MOMA) on `model`

with respect to `flux_ref`

, which is a vector of fluxes in the order of `reactions(model)`

. MOMA finds the shortest Euclidian distance between `flux_ref`

and `model`

with `modifications`

:

```
min Σᵢ (xᵢ - flux_refᵢ)²
s.t. S x = b
xₗ ≤ x ≤ xᵤ
```

Because the problem has a quadratic objective, a QP solver is required. See "Daniel, Vitkup & Church, Analysis of Optimality in Natural and Perturbed Metabolic Networks, Proceedings of the National Academy of Sciences, 2002" for more details.

Additional arguments are passed to `flux_balance_analysis`

.

Returns an optimized model that contains the resultant nearest flux.

**Example**

```
model = load_model("e_coli_core.json")
flux_ref = flux_balance_analysis_vec(model, Gurobi.Optimizer)
optmodel = minimize_metabolic_adjustment(
model,
flux_ref,
Gurobi.Optimizer;
modifications = [change_constraint("PFL"; lb=0, ub=0)], # find flux of mutant that is closest to the wild type (reference) model
)
value.(solution[:x]) # extract the flux from the optimizer
```

`COBREXA.minimize_metabolic_adjustment_analysis_dict`

— Method```
minimize_metabolic_adjustment_analysis_dict(model::MetabolicModel, args...; kwargs...) -> Union{Nothing, Dict{String, Float64}}
```

Perform minimization of metabolic adjustment (MOMA) and return a dictionary mapping the reaction IDs to fluxes. Arguments are forwarded to `minimize_metabolic_adjustment`

internally.

This function is kept for backwards compatibility, use `flux_vector`

instead.

`COBREXA.minimize_metabolic_adjustment_analysis_vec`

— Method```
minimize_metabolic_adjustment_analysis_vec(model::MetabolicModel, args...; kwargs...) -> Union{Nothing, Vector{Float64}}
```

Perform minimization of metabolic adjustment (MOMA) and return a vector of fluxes in the same order as the reactions in `model`

. Arguments are forwarded to `minimize_metabolic_adjustment`

internally.

This function is kept for backwards compatibility, use `flux_vector`

instead.

`COBREXA.parsimonious_flux_balance_analysis`

— Method```
parsimonious_flux_balance_analysis(model::MetabolicModel, optimizer; modifications, qp_modifications, relax_bounds) -> JuMP.Model
```

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_analysis`

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_analysis`

for usage examples of modifications.

Thhe 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).

**Example**

```
model = load_model("e_coli_core.json")
optmodel = parsimonious_flux_balance_analysis(model, biomass, Gurobi.Optimizer)
value.(solution[:x]) # extract the flux from the optimizer
```

`COBREXA.parsimonious_flux_balance_analysis_dict`

— Method```
parsimonious_flux_balance_analysis_dict(model::MetabolicModel, args...; kwargs...) -> Union{Nothing, Dict{String, Float64}}
```

Perform parsimonious flux balance analysis on `model`

using `optimizer`

. Returns a dictionary mapping the reaction IDs to fluxes. Arguments are forwarded to `parsimonious_flux_balance_analysis`

internally.

This function is kept for backwards compatibility, use `flux_dict`

instead.

`COBREXA.parsimonious_flux_balance_analysis_vec`

— Method```
parsimonious_flux_balance_analysis_vec(model::MetabolicModel, args...; kwargs...) -> Union{Nothing, Vector{Float64}}
```

Perform parsimonious flux balance analysis on `model`

using `optimizer`

. Returns a vector of fluxes in the same order as the reactions in `model`

. Arguments are forwarded to `parsimonious_flux_balance_analysis`

internally.

This function is kept for backwards compatibility, use `flux_vector`

instead.

`COBREXA._screen_args`

— Method```
_screen_args(argtuple, kwargtuple, modsname) -> Union{Tuple{}, NamedTuple}
```

Internal helper to check the presence and shape of modification and argument arrays in `screen`

and pals.

`COBREXA._screen_impl`

— Method```
_screen_impl(model::MetabolicModel; variants, analysis, args, workers)
```

The actual implementation of `screen`

.

`COBREXA._screen_optmodel_item`

— Method```
_screen_optmodel_item() -> Any
```

Internal helper for `screen_optmodel_modifications`

that computes one item of the screening task.

`COBREXA._screen_optmodel_modifications_impl`

— Method```
_screen_optmodel_modifications_impl(model::MetabolicModel, optimizer; common_modifications, modifications, args, analysis, workers)
```

The actual implementation of `screen_optmodel_modifications`

.

`COBREXA._screen_optmodel_prepare`

— Method```
_screen_optmodel_prepare(model, optimizer, common_modifications) -> Tuple{Any, JuMP.Model}
```

Internal helper for `screen_optmodel_modifications`

that creates the model and applies the modifications.

`COBREXA.screen`

— Method```
screen(args...; kwargs...)
```

Take an array of model-modifying function vectors in `variants`

, and execute the function `analysis`

on all variants of the `model`

specified by `variants`

. The computation is distributed over worker IDs in `workers`

. If `args`

are supplied (as an array of the same size as the `variants`

), they are forwarded as arguments to the corresponding analysis function calls.

The array of variants must contain vectors of single-parameter functions, these are applied to model in order. The functions must *not* modify the model, but rather return a modified copy. The copy should be made as shallow as possible, to increase memory efficiency of the process. Variant generators that modify the argument model in-place will cause unpredictable results. Refer to the definition of `screen_variant`

for details.

The function `analysis`

will receive a single argument (the modified model), together with arguments from `args`

expanded by `...`

. Supply an array of tuples or vectors to pass in multiple arguments to each function. If the argument values should be left intact (not expanded to multiple arguments), they must be wrapped in single-item tuple or vector.

The modification and analysis functions are transferred to `workers`

as-is; all packages required to run them (e.g. the optimization solvers) must be loaded there. Typically, you want to use the macro `@everywhere using MyFavoriteSolver`

from `Distributed`

package for loading the solvers.

**Return value**

The results of running `analysis`

are collected in to the resulting array, in a way that preserves the shape of the `variants`

, similarly as with `pmap`

.

The results of `analysis`

function must be serializable, preferably made only from pure Julia structures, because they may be transferred over the network between the computation nodes. For that reason, functions that return whole JuMP models that contain pointers to allocated C structures (such as `flux_balance_analysis`

used with `GLPK`

or `Gurobi`

otimizers) will generally not work in this context.

Note: this function is a thin argument-handling wrapper around `_screen_impl`

.

**Example**

```
function reverse_reaction(i::Int)
(model::CoreModel) -> begin
mod = copy(model)
mod.S[:,i] .*= -1 # this is unrealistic but sufficient for demonstration
mod
end
end
m = load_model(CoreModel, "e_coli_core.xml")
screen(m,
variants = [
[reverse_reaction(5)],
[reverse_reaction(3), reverse_reaction(6)]
],
analysis = mod -> mod.S[:,3]) # observe the changes in S
screen(m,
variants = [
[reverse_reaction(5)],
[reverse_reaction(3), reverse_reaction(6)]
],
analysis = mod -> flux_balance_analysis_vec(mod, GLPK.Optimizer))
```

`COBREXA.screen_optimize_objective`

— Method```
screen_optimize_objective(_, optmodel) -> Union{Nothing, Float64}
```

A variant of `optimize_objective`

directly usable in `screen_optmodel_modifications`

.

`COBREXA.screen_optmodel_modifications`

— Method```
screen_optmodel_modifications(args...; kwargs...)
```

Screen multiple modifications of the same optimization model.

This function is potentially more efficient than `screen`

because it avoids making variants of the model structure and remaking of the optimization model. On the other hand, modification functions need to keep the optimization model in a recoverable state (one that leaves the model usable for the next modification), which limits the possible spectrum of modifications applied.

Internally, `model`

is distributed to `workers`

and transformed into the optimization model using `make_optimization_model`

. `common_modifications`

are applied to the models at that point. Next, vectors of functions in `modifications`

are consecutively applied, and the result of `analysis`

function called on model are collected to an array of the same extent as `modifications`

. Calls of `analysis`

are optionally supplied with extra arguments from `args`

expanded with `...`

, just like in `screen`

.

Both the modification functions (in vectors) and the analysis function here have 2 base parameters (as opposed to 1 with `screen`

): first is the `model`

(carried through as-is), second is the prepared JuMP optimization model that may be modified and acted upon. As an example, you can use modification `change_constraint`

and analysis `screen_optimize_objective`

.

Note: This function is a thin argument-handling wrapper around `_screen_optmodel_modifications_impl`

.

`COBREXA.screen_variant`

— Function```
screen_variant(model::MetabolicModel, variant::Vector, analysis) -> Any
screen_variant(model::MetabolicModel, variant::Vector, analysis, args) -> Any
```

Helper function for `screen`

that applies all single-argument functions in `variant`

to the `model`

(in order from "first" to "last"), and executes `analysis`

on the result.

Can be used to test model variants locally.

`COBREXA.screen_variants`

— Method```
screen_variants(model, variants, analysis; workers) -> Array
```

A shortcut for `screen`

that only works with model variants.

`COBREXA.make_smoment_model`

— Method```
make_smoment_model(model::MetabolicModel; reaction_isozyme, gene_product_molar_mass, total_enzyme_capacity)
```

Construct a model with a structure given by sMOMENT algorithm; returns a `SMomentModel`

(see the documentation for details).

**Arguments**

`reaction_isozyme`

parameter is a function that returns a single`Isozyme`

for each reaction, or`nothing`

if the reaction is not enzymatic. If the reaction has multiple isozymes, use`smoment_isozyme_speed`

to select the fastest one, as recommended by the sMOMENT paper.`gene_product_molar_mass`

parameter is a function that returns a molar mass of each gene product as specified by sMOMENT.`total_enzyme_capacity`

is the maximum "enzyme capacity" in the model.

Alternatively, all function arguments also accept dictionaries that are used to provide the same data lookup.

## Sampling

`COBREXA._affine_hit_and_run_chain`

— Method```
_affine_hit_and_run_chain(warmup, lbs, ubs, C, cl, cu, iters, seed) -> Matrix{Float64}
```

Internal helper function for computing a single affine hit-and-run chain.

`COBREXA.affine_hit_and_run`

— Method```
affine_hit_and_run(m::MetabolicModel, warmup_points::Matrix{Float64}; sample_iters, workers, chains, seed) -> Any
```

Run a hit-and-run style sampling that starts from `warmup_points`

and uses their affine combinations for generating the run directions to sample the space delimited by `lbs`

and `ubs`

. The reaction rate vectors in `warmup_points`

should be organized in columns, i.e. `warmup_points[:,1]`

is the first set of reaction rates.

There are total `chains`

of hit-and-run runs, each on a batch of `size(warmup_points, 2)`

points. The runs are scheduled on `workers`

, for good load balancing `chains`

should be ideally much greater than `length(workers)`

.

Each run continues for `maximum(sample_iters)`

iterations; the numbers in `sample_iters`

represent the iterations at which the whole "current" batch of points is collected for output. For example, `sample_iters=[1,4,5]`

causes the process run for 5 iterations, returning the sample batch that was produced by 1st, 4th and last (5th) iteration.

Returns a matrix of sampled reaction rates (in columns), with all collected samples horizontally concatenated. The total number of samples (columns) will be `size(warmup_points,2) * chains * length(sample_iters)`

.

**Example**

```
warmup_points = warmup_from_variability(model, GLPK.Optimizer)
samples = affine_hit_and_run(model, warmup_points, sample_iters = 101:105)
# convert the result to flux (for models where the distinction matters):
fluxes = reaction_flux(model)' * samples
```

`COBREXA._maximize_warmup_reaction`

— Method```
_maximize_warmup_reaction(opt_model, rid, ret) -> Any
```

A helper function for finding warmup points from reaction variability.

`COBREXA.warmup_from_variability`

— Function```
warmup_from_variability(model::MetabolicModel, optimizer, n_points::Int64) -> Matrix{Float64}
warmup_from_variability(model::MetabolicModel, optimizer, n_points::Int64, seed; kwargs...) -> Matrix{Float64}
```

Generates FVA-like warmup points for samplers, by selecting random points by minimizing and maximizing reactions. Can not return more than 2 times the number of reactions in the model.

`COBREXA.warmup_from_variability`

— Function```
warmup_from_variability(model::MetabolicModel, optimizer) -> Matrix{Float64}
warmup_from_variability(model::MetabolicModel, optimizer, min_reactions::AbstractVector{Int64}) -> Matrix{Float64}
warmup_from_variability(model::MetabolicModel, optimizer, min_reactions::AbstractVector{Int64}, max_reactions::AbstractVector{Int64}; modifications, workers) -> Matrix{Float64}
```

Generate FVA-like warmup points for samplers, by minimizing and maximizing the specified reactions. The result is returned as a matrix, each point occupies as single column in the result.

!!! warning Limited effect of modifications in `warmup_from_variability`

Modifications of the optimization model applied in `modifications`

parameter that change the semantics of the model have an effect on the warmup points, but do not automatically carry to the subsequent sampling. Users are expected to manually transplant any semantic changes to the actual sampling functions, such as `affine_hit_and_run`

.

## Analysis modifiers

`COBREXA.add_crowding_constraints`

— Method```
add_crowding_constraints(weights::Dict{Int64, Float64}) -> COBREXA.var"#570#572"{Dict{Int64, Float64}}
```

Adds a molecular crowding constraint to the optimization problem: `∑ wᵢ × vᵢ ≤ 1`

where `wᵢ`

is a weight and `vᵢ`

is a flux index in the model's reactions specified in `weights`

as `vᵢ => wᵢ`

pairs.

See Beg, Qasim K., et al. "Intracellular crowding defines the mode and sequence of substrate uptake by Escherichia coli and constrains its metabolic activity." Proceedings of the National Academy of Sciences 104.31 (2007) for more details.

`COBREXA.add_crowding_constraints`

— Method```
add_crowding_constraints(weights::Dict{String, Float64}) -> COBREXA.var"#574#575"{Dict{String, Float64}}
```

Variant of `add_crowding_constraints`

that takes a dictinary of reactions `ids`

instead of reaction indices mapped to weights.

`COBREXA.change_constraint`

— Method```
change_constraint(id::String; lb, ub) -> COBREXA.var"#579#580"{Nothing, Nothing, String}
```

Change the lower and upper bounds (`lb`

and `ub`

respectively) of reaction `id`

if supplied.

`COBREXA.change_objective`

— Method```
change_objective(new_objective::Union{String, Vector{String}}; weights, sense) -> Union{COBREXA.var"#582#584"{MathOptInterface.OptimizationSense, Vector{String}}, COBREXA.var"#582#584"{MathOptInterface.OptimizationSense, String}}
```

Modification that changes the objective function used in a constraint based analysis function. `new_objective`

can be a single reaction identifier, or an array of reactions identifiers.

Optionally, the objective can be weighted by a vector of `weights`

, and a optimization `sense`

can be set to either `MAX_SENSE`

or `MIN_SENSE`

.

`COBREXA.constrain_objective_value`

— Method```
constrain_objective_value(tolerance) -> COBREXA.var"#576#577"
```

Limit the objective value to `tolerance`

-times the current objective value, as with `objective_bounds`

.

`COBREXA._do_knockout`

— Method```
_do_knockout(model::MetabolicModel, opt_model, gene_ids::Vector{String})
```

Internal helper for knockouts on generic MetabolicModels. This can be overloaded so that the knockouts may work differently (more efficiently) with other models.

`COBREXA.knockout`

— Method```
knockout(gene_id::String) -> COBREXA.var"#586#587"{Vector{String}}
```

A helper variant of `knockout`

for a single gene.

`COBREXA.knockout`

— Method```
knockout(gene_ids::Vector{String}) -> COBREXA.var"#586#587"{Vector{String}}
```

A modification that zeroes the bounds of all reactions that would be knocked out by the combination of specified genes (effectively disabling the reactions).

A slightly counter-intuitive behavior may occur if knocking out multiple genes: Because this only changes the reaction bounds, multiple gene knockouts *must* be specified in a single call to `knockout`

, because the modifications have no way to remember which genes are already knocked out and which not.

In turn, having a reaction that can be catalyzed either by Gene1 or by Gene2, specifying `modifications = [knockout(["Gene1", "Gene2"])]`

does indeed disable the reaction, but `modifications = [knockout("Gene1"), knockout("Gene2")]`

does *not* disable the reaction (although reactions that depend either only on Gene1 or only on Gene2 are disabled).

`COBREXA.add_loopless_constraints`

— Method```
add_loopless_constraints(; max_flux_bound, strict_inequality_tolerance) -> COBREXA.var"#591#596"{Float64, Int64}
```

Add quasi-thermodynamic constraints to the model to ensure that no thermodynamically infeasible internal cycles can occur. Adds the following constraints to the problem:

```
-max_flux_bound × (1 - yᵢ) ≤ xᵢ ≤ max_flux_bound × yᵢ
-max_flux_bound × yᵢ + strict_inequality_tolerance × (1 - yᵢ) ≤ Gᵢ
Gᵢ ≤ -strict_inequality_tolerance × yᵢ + max_flux_bound × (1 - yᵢ)
Nᵢₙₜ' × G = 0
yᵢ ∈ {0, 1}
Gᵢ ∈ ℝ
i ∈ internal reactions
Nᵢₙₜ is the nullspace of the internal stoichiometric matrix
```

Note, this modification introduces binary variables, so an optimization solver capable of handing mixed integer problems needs to be used. The arguments `max_flux_bound`

and `strict_inequality_tolerance`

implement the "big-M" method of indicator constraints.

For more details about the algorithm, see `Schellenberger, Lewis, and, Palsson. "Elimination of thermodynamically infeasible loops in steady-state metabolic models.", Biophysical journal, 2011`

.

`COBREXA.add_moment_constraints`

— Method```
add_moment_constraints(kcats::Dict{String, Float64}, protein_mass_fraction::Float64) -> COBREXA.var"#601#607"{Dict{String, Float64}, Float64}
```

A modification that adds enzyme capacity constraints to the problem using a *modified* version of the MOMENT algorithm. Requires specific activities, `ksas`

[mmol product/g enzyme/h], for each reaction. Proteins are identified by their associated gene IDs. Adds a variable vector `y`

to the problem corresponding to the protein concentration [g enzyme/gDW cell] of each gene product in the order of `genes(model)`

. The total protein concentration [g protein/gDW cell] is constrained to be less than or equal to the `protein_mass_fraction`

. Reaction flux constraints are changed to the MOMENT constraints (see below) for all reactions that have a gene reaction rule, otherwise the flux bounds are left unaltered.

See Adadi, Roi, et al. "Prediction of microbial growth rate versus biomass yield by a metabolic network with kinetic parameters." PLoS computational biology (2012) for more details of the original algorithm.

Here, a streamlined version of the algorithm is implemented to ensure that the correct units are used. Specifically, this implementation uses specific activities instead of `kcats`

. Thus, for a reaction that can only proceed forward and is catalyzed by protein `a`

, the flux `x[i]`

is bounded by `x[i] <= ksas[i] * y[a]`

. If isozymes `a`

or `b`

catalyse the reaction, then `x[i] <= ksas[i] * (y[a] + y[b])`

. If a reaction is catalyzed by subunits `a`

and `b`

then `x[i] <= ksas[i] * min(y[a], y[b])`

. These rules are applied recursively in the model like in the original algorithm. The enzyme capacity constraint is then implemented by `sum(y) ≤ protein_mass_fraction`

. The major benefit of using `ksas`

instead of `kcats`

is that active site number and unit issues are prevented.

**Example**

```
flux_balance_analysis(
...,
modifications = [ add_moment_constraints(my_kcats, 0.6) ],
)
```

`COBREXA.change_optimizer`

— Method```
change_optimizer(optimizer) -> COBREXA.var"#615#616"
```

Change the JuMP optimizer used to run the optimization.

This may be used to try different approaches for reaching the optimum, and in problems that may require different optimizers for different parts, such as the `parsimonious_flux_balance_analysis`

.

`COBREXA.change_optimizer_attribute`

— Method```
change_optimizer_attribute(attribute_key, value) -> COBREXA.var"#617#618"
```

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.

`COBREXA.change_sense`

— Method```
change_sense(objective_sense) -> COBREXA.var"#613#614"
```

Change the objective sense of optimization. Possible arguments are `MAX_SENSE`

and `MIN_SENSE`

.

If you want to change the objective and sense at the same time, use `change_objective`

instead to do both at once.

`COBREXA.silence`

— Function`silence`

Modification that disable all output from the JuMP optimizer (shortcut for `set_silent`

from JuMP).