parsimonious_flux_balance_analysis.jl 3.64 KB
Newer Older
St. Elmo's avatar
St. Elmo committed
1
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
2
    parsimonious_flux_balance_analysis(
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
3
4
5
6
7
8
        model::MetabolicModel,
        optimizer;
        modifications = [],
        qp_modifications = [],
        relax_bounds=[1.0, 0.999999, 0.99999, 0.9999, 0.999, 0.99],
    )
St. Elmo's avatar
St. Elmo committed
9
10

Run parsimonious flux balance analysis (pFBA) on the `model`.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
11
12
13
14
15
16
17
18
19
20
21
22
23
24

pFBA gets the model optimum by standard FBA (using
[`flux_balance_analysis`](@ref) 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.

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 `nothing` if the
optimization failed.
St. Elmo's avatar
St. Elmo committed
25
26
27
28
29

# Example
```
optimizer = Gurobi.Optimizer
atts = Dict("OutputFlag" => 0)
30
model = load_model(StandardModel, "iJO1366.json")
St. Elmo's avatar
St. Elmo committed
31
biomass = findfirst(model.reactions, "BIOMASS_Ec_iJO1366_WT_53p95M")
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
32
sol = pfba(model, biomass, Gurobi.optimizer)
St. Elmo's avatar
St. Elmo committed
33
34
35
36
37
```
"""
function parsimonious_flux_balance_analysis(
    model::MetabolicModel,
    optimizer;
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
38
39
40
    modifications = [],
    qp_modifications = [],
    relax_bounds = [1.0, 0.999999, 0.99999, 0.9999, 0.999, 0.99],
St. Elmo's avatar
St. Elmo committed
41
42
43
)
    # Run FBA
    opt_model = flux_balance_analysis(model, optimizer; modifications = modifications)
St. Elmo's avatar
St. Elmo committed
44
    is_solved(opt_model) || return nothing # FBA failed
St. Elmo's avatar
St. Elmo committed
45

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
46
47
48
49
50
51
52
    # get the objective
    Z = objective_value(opt_model)
    original_objective = COBREXA.JuMP.objective_function(opt_model)

    # prepare the model for pFBA
    for mod in qp_modifications
        mod(model, opt_model)
St. Elmo's avatar
St. Elmo committed
53
54
    end

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
55
    # add the minimization constraint for total flux
St. Elmo's avatar
St. Elmo committed
56
    v = opt_model[:x] # fluxes
57
58
59
60
61

    if solver_name(opt_model) in ["Gurobi", ]
        set_start_value.(all_variables(opt_model), value.(all_variables(opt_model))) # warm start
    end
    
St. Elmo's avatar
St. Elmo committed
62
63
    @objective(opt_model, Min, sum(dot(v, v)))

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
64
65
66
67
68
    for rb in relax_bounds
        lb, ub = objective_bounds(rb)(Z)
        @_models_log @info "pFBA step relaxed to [$lb,$ub]"
        @constraint(opt_model, pfba_constraint, lb <= original_objective <= ub)

St. Elmo's avatar
St. Elmo committed
69
        optimize!(opt_model)
St. Elmo's avatar
St. Elmo committed
70
        is_solved(opt_model) && break
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
71
72

        COBREXA.JuMP.delete(opt_model, pfba_constraint)
St. Elmo's avatar
St. Elmo committed
73
        COBREXA.JuMP.unregister(opt_model, :pfba_constraint)
St. Elmo's avatar
St. Elmo committed
74
75
    end

St. Elmo's avatar
St. Elmo committed
76
    is_solved(opt_model) || return nothing # pFBA failed
St. Elmo's avatar
St. Elmo committed
77
78
79
80
81

    return opt_model
end

"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
82
    parsimonious_flux_balance_analysis_vec(args...; kwargs...)
St. Elmo's avatar
St. Elmo committed
83
84
85

Perform parsimonious flux balance analysis on `model` using `optimizer`. 
Returns a vector of fluxes in the same order as the reactions in `model`. 
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
86
Arguments are forwarded to [`parsimonious_flux_balance_analysis`](@ref) internally.
St. Elmo's avatar
St. Elmo committed
87
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
88
89
90
91
function parsimonious_flux_balance_analysis_vec(args...; kwargs...)
    opt_model = parsimonious_flux_balance_analysis(args...; kwargs...)

    isnothing(opt_model) && return nothing
St. Elmo's avatar
St. Elmo committed
92
93
94
95
96

    return value.(opt_model[:x])
end

"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
97
    parsimonious_flux_balance_analysis_dict(model::MetabolicModel, args...; kwargs...)
St. Elmo's avatar
St. Elmo committed
98
99

Perform parsimonious flux balance analysis on `model` using `optimizer`. 
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
100
101
Returns a dictionary mapping the reaction IDs to fluxes. 
Arguments are forwarded to [`parsimonious_flux_balance_analysis`](@ref) internally.
St. Elmo's avatar
St. Elmo committed
102
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
103
104
105
106
function parsimonious_flux_balance_analysis_dict(model::MetabolicModel, args...; kwargs...)
    opt_fluxes = parsimonious_flux_balance_analysis_vec(model, args...; kwargs...)

    isnothing(opt_fluxes) && return nothing
St. Elmo's avatar
St. Elmo committed
107

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
108
    return Dict(zip(reactions(model), opt_fluxes))
St. Elmo's avatar
St. Elmo committed
109
end