flux_variability_analysis.jl 6.39 KB
Newer Older
Sylvain Arreckx's avatar
Sylvain Arreckx committed
1
"""
2
    flux_variability_analysis(
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
3
4
5
6
7
        model::LM,
        reactions::Vector{Int},
        optimizer,
        workers = [myid()];
        gamma::AbstractFloat = 1.0,
8
    )::Matrix{Float64} where {LM<:MetabolicModel}
Sylvain Arreckx's avatar
Sylvain Arreckx committed
9

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
10
11
12
13
# Flux variability analysis (FVA)

FVA solves the pair of optimization problems in `model` for each flux xᵢ listed
in `reactions`
Sylvain Arreckx's avatar
Sylvain Arreckx committed
14
15
16
17
18
19
20
```
min/max xᵢ
s.t. S x = b
     xₗ ≤ x ≤ xᵤ
     cᵀx ≥ γ Z₀
```
where Z₀:= cᵀx₀ is the objective value of an optimal solution to the associated
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
21
FBA problem.
Sylvain Arreckx's avatar
Sylvain Arreckx committed
22

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
23
24
25
26
27
Internally uses the specified JuMP-compatible `optimizer`, and the work is
scheduled in parallel on `workers`.

Returns a matrix of minima and maxima of size (length(reactions),2).
"""
28
function flux_variability_analysis(
29
30
31
32
33
    model::LM,
    reactions::Vector{Int},
    optimizer,
    workers = [myid()];
    gamma::AbstractFloat = 1.0,
34
)::Matrix{Float64} where {LM<:MetabolicModel}
Sylvain Arreckx's avatar
Sylvain Arreckx committed
35

36
    if any(reactions .< 1) || any(reactions .> n_reactions(model))
37
        throw(DomainError(reactions, "Index exceeds number of reactions."))
38
39
    end

St. Elmo's avatar
St. Elmo committed
40
    optimization_model = flux_balance_analysis(model, optimizer)
St. Elmo's avatar
St. Elmo committed
41
    Z0 = COBREXA.JuMP.objective_value(optimization_model)
42
    optimization_model = nothing # we won't need this one anymore, so free the memory
43

44
45
46
    # store a JuMP optimization model at all workers
    save_model = :(
        begin
St. Elmo's avatar
St. Elmo committed
47
            optmodel = COBREXA.make_optimization_model($model, $optimizer)
St. Elmo's avatar
St. Elmo committed
48
49
50
51
52
53
54
            COBREXA._FVA_add_constraint(
                optmodel,
                $(objective(model)),
                optmodel[:x],
                $Z0,
                $gamma,
            )
55
56
            optmodel
        end
57
    )
58
59
    map(fetch, save_at.(workers, :cobrexa_parfva_model, Ref(save_model)))
    save_model = nothing # this has some volume, free it again
60

61
    # schedule FVA parts parallely using pmap
62
    fluxes = dpmap(
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
63
        rid -> :(COBREXA._FVA_optimize_reaction(cobrexa_parfva_model, $rid)),
64
        CachingPool(workers),
65
66
        [-reactions reactions],
    )
67

68
    # free the data on workers
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
69
    map(fetch, remove_from.(workers, :cobrexa_parfva_model))
70

71
72
    return fluxes
end
73

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
74
"""
75
    flux_variability_analysis(
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
76
77
78
        model::LM,
        optimizer;
        gamma::AbstractFloat = 1.0,
79
    ) where {LM<:MetabolicModel}
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
80
81
82

A simpler version of FVA that maximizes and minimizes all reactions in the model.
"""
83
function flux_variability_analysis(
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
84
85
86
    model::LM,
    optimizer;
    gamma::AbstractFloat = 1.0,
87
) where {LM<:MetabolicModel}
88
89
    n = n_reactions(model)
    return flux_variability_analysis(model, collect(1:n), optimizer; gamma = gamma)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
90
91
92
93
94
95
96
97
98
99
end


"""
    _FVA_add_constraint(model, c, x, Z0, gamma)

Internal helper function for adding constraints to a model. Exists mainly
because for avoiding namespace problems on remote workers.
"""
function _FVA_add_constraint(model, c, x, Z0, gamma)
St. Elmo's avatar
St. Elmo committed
100
    COBREXA.JuMP.@constraint(model, c' * x  gamma * Z0)
101
102
end

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
103
104
105
106
107
108
109
"""
    _FVA_get_opt(model, rid)

Helper for creating the optimized model on a remote worker, for avoiding
namespace problems.
"""
function _FVA_optimize_reaction(model, rid)
110
    sense = rid > 0 ? MOI.MAX_SENSE : MOI.MIN_SENSE
St. Elmo's avatar
St. Elmo committed
111
    var = COBREXA.JuMP.all_variables(model)[abs(rid)]
112

St. Elmo's avatar
St. Elmo committed
113
114
115
    COBREXA.JuMP.@objective(model, sense, var)
    COBREXA.JuMP.optimize!(model)
    return COBREXA.JuMP.objective_value(model)
116
end
117
118

"""
119
    fva(model::StandardModel, optimizer; optimum_bound=1.0-_constants.tolerance, modifications)
120

St. Elmo's avatar
St. Elmo committed
121
122
123
Run flux variability analysis (FVA) on the `model` (of type `StandardModel`). 
Optionally specifying problem modifications like in [`flux_balance_analysis`](@ref).
This algorithm runs FBA on the model to determine the optimum of the objective.
St. Elmo's avatar
St. Elmo committed
124
125
This optimum then constrains subsequent problems, where `optimum_bound` can be used to 
relax this constraint as a fraction of the FBA optimum, e.g. 
126
Note, the `optimizer` must be set to perform the analysis, any JuMP solver will work.
St. Elmo's avatar
St. Elmo committed
127
128
129
Note, this function only runs serially. 
Consider using a different model type for parallel implementations. 
Returns two dictionaries (`fva_max` and `fva_min`) that map each reaction `id` to dictionaries of the resultant flux distributions when that `id` is optimized.
130

131
See also: [`CoreModel`](@ref)
St. Elmo's avatar
St. Elmo committed
132

133
134
135
136
137
138
139
140
# Example
```
optimizer = Gurobi.Optimizer
model = Cobraread_model("iJO1366.json")
biomass = findfirst(model.reactions, "BIOMASS_Ec_iJO1366_WT_53p95M")
fva_max, fva_min = fva(model, biomass, optimizer; solver_attributes=atts)
```
"""
St. Elmo's avatar
St. Elmo committed
141
142
function flux_variability_analysis(
    model::StandardModel,
143
    optimizer;
144
    optimum_bound = 1.0 - _constants.tolerance,
St. Elmo's avatar
St. Elmo committed
145
    modifications = [(model, opt_model) -> nothing],
146
)
St. Elmo's avatar
St. Elmo committed
147
    # Run FBA
St. Elmo's avatar
St. Elmo committed
148
    opt_model = flux_balance_analysis(model, optimizer; modifications = modifications)
149

St. Elmo's avatar
St. Elmo committed
150
151
    COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ||
        (return nothing, nothing)
152

St. Elmo's avatar
St. Elmo committed
153
154
    fva_min = Dict{String,Union{Nothing,Dict{String,Float64}}}()
    fva_max = Dict{String,Union{Nothing,Dict{String,Float64}}}()
155

156
    # Now do FVA
St. Elmo's avatar
St. Elmo committed
157
158
159
    v = opt_model[:x]

    λ = COBREXA.JuMP.objective_value(opt_model) # objective value
St. Elmo's avatar
St. Elmo committed
160
161
    λmin = min(optimum_bound * λ, λ * 1.0 / optimum_bound)
    λmax = max(optimum_bound * λ, λ * 1.0 / optimum_bound)
St. Elmo's avatar
St. Elmo committed
162

St. Elmo's avatar
St. Elmo committed
163
164
    COBREXA.JuMP.@constraint(
        opt_model,
St. Elmo's avatar
St. Elmo committed
165
        λmin <= COBREXA.JuMP.objective_function(opt_model) <= λmax # in case there is a negative bound
St. Elmo's avatar
St. Elmo committed
166
    )
167
168

    for i = 1:length(v)
St. Elmo's avatar
St. Elmo committed
169
170
        COBREXA.JuMP.@objective(opt_model, Max, v[i])
        COBREXA.JuMP.optimize!(opt_model)
171
        status = (
St. Elmo's avatar
St. Elmo committed
172
173
            COBREXA.JuMP.termination_status(opt_model) == MOI.OPTIMAL ||
            COBREXA.JuMP.termination_status(opt_model) == MOI.LOCALLY_SOLVED
174
175
        )
        if status
St. Elmo's avatar
St. Elmo committed
176
            fva_max[model.reactions[i].id] =
St. Elmo's avatar
St. Elmo committed
177
                Dict(zip(reactions(model), value.(opt_model[:x])))
178
        else
St. Elmo's avatar
St. Elmo committed
179
            @warn "Error maximizing index: $i with error $(termination_status(opt_model))"
St. Elmo's avatar
St. Elmo committed
180
            fva_max[model.reactions[i].id] = nothing
181
182
        end

St. Elmo's avatar
St. Elmo committed
183
184
        @objective(opt_model, Min, v[i])
        optimize!(opt_model)
185
        status = (
St. Elmo's avatar
St. Elmo committed
186
187
            termination_status(opt_model) == MOI.OPTIMAL ||
            termination_status(opt_model) == MOI.LOCALLY_SOLVED
188
189
        )
        if status
St. Elmo's avatar
St. Elmo committed
190
            fva_min[model.reactions[i].id] =
St. Elmo's avatar
St. Elmo committed
191
                Dict(zip(reactions(model), value.(opt_model[:x])))
192
        else
St. Elmo's avatar
St. Elmo committed
193
            @warn "Error minimizing index: $i with error $(termination_status(opt_model))"
St. Elmo's avatar
St. Elmo committed
194
            fva_min[model.reactions[i].id] = nothing
195
196
197
198
199
        end
    end

    return fva_max, fva_min
end