flux_variability_analysis.jl 5.5 KB
Newer Older
Sylvain Arreckx's avatar
Sylvain Arreckx committed
1
"""
2
    flux_variability_analysis(
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
3
        model::MetabolicModel,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
4
        reactions::Vector{Int},
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
5
6
7
        optimizer;
        modifications = [],
        workers = [myid()],
8
9
        optimal_objective_value = nothing,
        bounds = z -> (z, Inf),
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
10
11
12
13
14
        ret = objective_value,
    )::Matrix{Float64}

Flux variability analysis solves a pair of optimization problems in `model` for
each flux listed in `reactions`:
Sylvain Arreckx's avatar
Sylvain Arreckx committed
15
```
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
16
 min,max xᵢ
Sylvain Arreckx's avatar
Sylvain Arreckx committed
17
s.t. S x = b
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
18
19
20
    xₗ ≤ x ≤ xᵤ
     cᵀx ≥ bounds(Z₀)[1]
     cᵀx ≤ bounds(Z₀)[2]
Sylvain Arreckx's avatar
Sylvain Arreckx committed
21
```
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
22
where Z₀:= cᵀx₀ is the objective value of an optimal solution of the associated
23
24
FBA problem (see [`flux_balance_analysis`](@ref) for a related analysis, also
for explanation of the `modifications` argument).
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
25
26
27
28

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
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
29
30
first and second pair to remove the limit. Use [`gamma_bounds`](@ref) and
[`objective_bounds`](@ref) for simple bounds.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
31
32
33

`optimizer` must be set to a `JuMP`-compatible optimizer. The computation of
the individual optimization problems is transparently distributed to `workers`
34
35
36
37
38
(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.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
39
40

`ret` is a function used to extract results from optimized JuMP models of the
41
42
43
44
individual reactions. 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]))`.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
45
46

Returns a matrix of extracted `ret` values for minima and maxima, of total size
47
48
49
(`length(reactions)`,2). The optimizer result status is checked with
[`is_solved`](@ref); `nothing` is returned if the optimization failed for any
reason.
50
51
52
53
54
55

# Example
```
model = load_model("e_coli_core.json")
flux_variability_analysis(model, [1, 2, 3, 42], GLPK.optimizer)
```
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
56
"""
57
function flux_variability_analysis(
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
58
    model::MetabolicModel,
59
    reactions::Vector{Int},
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
60
61
62
    optimizer;
    modifications = [],
    workers = [myid()],
63
64
    optimal_objective_value = nothing,
    bounds = z -> (z, Inf),
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
65
66
    ret = objective_value,
)
67
    if any(reactions .< 1) || any(reactions .> n_reactions(model))
68
        throw(DomainError(reactions, "Index exceeds number of reactions."))
69
70
    end

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
71
    Z = bounds(
72
        isnothing(optimal_objective_value) ?
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
73
74
        objective_value(
            flux_balance_analysis(model, optimizer; modifications = modifications),
75
        ) : optimal_objective_value,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
76
    )
77

78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
    return screen_optmodel_modifications(
        model,
        optimizer;
        common_modifications = vcat(
            modifications,
            [
                (model, opt_model) -> begin
                    Z[1] > -Inf && @constraint(
                        opt_model,
                        objective(model)' * opt_model[:x] >= Z[1]
                    )
                    Z[2] < Inf && @constraint(
                        opt_model,
                        objective(model)' * opt_model[:x] <= Z[2]
                    )
                end,
            ],
        ),
        args = [-reactions reactions],
        analysis = (_, opt_model, ridx) -> _max_variability_flux(opt_model, ridx, ret),
98
    )
99
end
100

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
101
"""
102
    flux_variability_analysis(
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
103
        model::MetabolicModel,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
104
        optimizer;
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
105
106
        kwargs...
    )
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
107

108
109
A simpler version of [`flux_variability_analysis`](@ref) that maximizes and
minimizes all reactions in the model. Arguments are forwarded.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
110
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
111
function flux_variability_analysis(model::MetabolicModel, optimizer; kwargs...)
112
    n = n_reactions(model)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
113
    return flux_variability_analysis(model, collect(1:n), optimizer; kwargs...)
114
115
end

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
116
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
117
118
119
120
121
    flux_variability_analysis_dict(
        model::MetabolicModel,
        optimizer;
        kwargs...
    )
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
122

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
123
124
125
A variant of [`flux_variability_analysis`](@ref) that returns the individual
maximized and minimized fluxes of all reactions as two dictionaries (of
dictionaries). All keyword arguments except `ret` are passed through.
St. Elmo's avatar
St. Elmo committed
126
127
128
129
130
131
132
133
134

# Example
```
mins, maxs = flux_variability_analysis_dict(
    model,
    Tulip.Optimizer;
    bounds = objective_bounds(0.99),
    modifications = [
        change_optimizer_attribute("IPM_IterationsLimit", 500),
135
136
        change_constraint("EX_glc__D_e"; lb = -10, ub = -10),
        change_constraint("EX_o2_e"; lb = 0, ub = 0),
St. Elmo's avatar
St. Elmo committed
137
138
139
    ],
)
```
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
140
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
141
142
143
144
145
function flux_variability_analysis_dict(model::MetabolicModel, optimizer; kwargs...)
    vs = flux_variability_analysis(
        model,
        optimizer;
        kwargs...,
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
146
        ret = m -> JuMP.value.(m[:x]),
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
147
148
    )
    rxns = reactions(model)
149

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
150
    return (
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
151
152
        Dict(zip(rxns, [Dict(zip(rxns, fluxes)) for fluxes in vs[:, 1]])),
        Dict(zip(rxns, [Dict(zip(rxns, fluxes)) for fluxes in vs[:, 2]])),
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
153
    )
154
end
155

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
156
"""
157
    _max_variability_flux(opt_model, rid, ret)
158

159
Internal helper for maximizing reactions in optimization model.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
160
"""
161
function _max_variability_flux(opt_model, rid, ret)
St. Elmo's avatar
St. Elmo committed
162
    sense = rid > 0 ? MAX_SENSE : MIN_SENSE
163
    var = all_variables(opt_model)[abs(rid)]
164

165
166
    @objective(opt_model, sense, var)
    optimize!(opt_model)
167

168
169
    if is_solved(opt_model)
        return ret(opt_model)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
170
171
172
    else
        return nothing
    end
173
end