Analysis functions

envelope_lattice(model::MetabolicModel, ridxs::Vector{Int64}; samples, ranges, reaction_samples) -> Base.Generator{_A, COBREXA.var"#458#459"} 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.

envelope_lattice(model::MetabolicModel, rids::Vector{String}; kwargs...) -> Base.Generator{_A, COBREXA.var"#458#459"} where _A

Version of envelope_lattice that works on string reaction IDs instead of integer indexes.

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

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.

make_expression_limited_model(model::MetabolicModel; relative_expression, bounding_function)

Build an ExpressionLimitedModel.

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

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.

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.

_max_variability_flux(opt_model, flux, sense, ret) -> Any

Internal helper for maximizing reactions in optimization model.

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.

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)

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

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

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

Shortcut for reaction_variability_analysis that examines all reactions.

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.

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 Isozymes 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.

_mmdf_add_df_bound(lb, ub) -> COBREXA.var"#545#546"

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

_mmdf_concen_objective(midx, sense) -> COBREXA.var"#543#544"

Helper function to change the objective to optimizing some concentration.

_mmdf_dgr_objective(ridx, sense) -> COBREXA.var"#541#542"

Helper function to change the objective to optimizing some dG.

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.

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.

minimize_metabolic_adjustment(flux_ref_dict::Dict{String, Float64}) -> COBREXA.var"#550#552"{Dict{String, Float64}}

Overload of minimize_metabolic_adjustment that works with a dictionary of fluxes.

minimize_metabolic_adjustment(flux_ref::Vector{Float64}) -> COBREXA.var"#548#549"{Vector{Float64}}

An optimization model modification that implements the MOMA in minimize_metabolic_adjustment_analysis.

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

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.

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.

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

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.

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.

_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.

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

The actual implementation of screen.

_screen_optmodel_item() -> Any

Internal helper for screen_optmodel_modifications that computes one item of the screening task.

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

The actual implementation of screen_optmodel_modifications.

_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.

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

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

A variant of optimize_objective directly usable in screen_optmodel_modifications.

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.

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.

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

A shortcut for screen that only works with model variants.

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.

_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.

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

_maximize_warmup_reaction(opt_model, rid, ret) -> Any

A helper function for finding warmup points from reaction variability.

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.

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.

add_crowding_constraints(weights::Dict{Int64, Float64}) -> COBREXA.var"#577#579"{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.

add_crowding_constraints(weights::Dict{String, Float64}) -> COBREXA.var"#581#582"{Dict{String, Float64}}

Variant of add_crowding_constraints that takes a dictinary of reactions ids instead of reaction indices mapped to weights.

change_constraint(id::String; lb, ub) -> COBREXA.var"#586#587"{Nothing, Nothing, String}

Change the lower and upper bounds (lb and ub respectively) of reaction id if supplied.

change_objective(new_objective::Union{String, Vector{String}}; weights, sense) -> Union{COBREXA.var"#589#591"{MathOptInterface.OptimizationSense, Vector{String}}, COBREXA.var"#589#591"{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.

constrain_objective_value(tolerance) -> COBREXA.var"#583#584"

Limit the objective value to tolerance-times the current objective value, as with objective_bounds.

_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.

knockout(gene_id::String) -> COBREXA.var"#593#594"{Vector{String}}

A helper variant of knockout for a single gene.

knockout(gene_ids::Vector{String}) -> COBREXA.var"#593#594"{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).

add_loopless_constraints(; max_flux_bound, strict_inequality_tolerance) -> COBREXA.var"#598#603"{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.

add_moment_constraints(kcats::Dict{String, Float64}, protein_mass_fraction::Float64) -> COBREXA.var"#608#614"{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) ],
)

change_optimizer(optimizer) -> COBREXA.var"#622#623"

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.

change_optimizer_attribute(attribute_key, value) -> COBREXA.var"#624#625"

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.

change_sense(objective_sense) -> COBREXA.var"#620#621"

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.

silence

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