parsimonious_flux_balance_analysis.jl 4.32 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

St. Elmo's avatar
St. Elmo committed
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
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,
27
Schramm, Gunnar, Purvine, Samuel O, Lopez-Ferrer, Daniel, Weitz, Karl K, Eils,
St. Elmo's avatar
St. Elmo committed
28
29
Roland, König, Rainer, Smith, Richard D, Palsson, Bernhard Ø, (2010) Omic data
from evolved E. coli are consistent with computed optimal growth from
30
genome-scale models. Molecular Systems Biology, 6. 390. doi:
St. Elmo's avatar
St. Elmo committed
31
accession:10.1038/msb.2010.47" for more details.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
32
33
34
35
36
37

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
38
using the callback in `qp_modifications`, which are applied after the FBA. See
St. Elmo's avatar
St. Elmo committed
39
the documentation of [`flux_balance_analysis`](@ref) for usage examples of modifications.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
40
41
42
43
44
45
46

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
47
48
49

# Example
```
50
51
52
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
St. Elmo's avatar
St. Elmo committed
53
54
55
56
57
```
"""
function parsimonious_flux_balance_analysis(
    model::MetabolicModel,
    optimizer;
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
58
59
60
    modifications = [],
    qp_modifications = [],
    relax_bounds = [1.0, 0.999999, 0.99999, 0.9999, 0.999, 0.99],
St. Elmo's avatar
St. Elmo committed
61
62
63
)
    # Run FBA
    opt_model = flux_balance_analysis(model, optimizer; modifications = modifications)
St. Elmo's avatar
St. Elmo committed
64
    is_solved(opt_model) || return nothing # FBA failed
St. Elmo's avatar
St. Elmo committed
65

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
66
67
    # get the objective
    Z = objective_value(opt_model)
68
    original_objective = objective_function(opt_model)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
69
70
71
72

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

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
75
    # add the minimization constraint for total flux
St. Elmo's avatar
St. Elmo committed
76
77
78
    v = opt_model[:x] # fluxes
    @objective(opt_model, Min, sum(dot(v, v)))

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
79
80
81
82
83
    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
84
        optimize!(opt_model)
St. Elmo's avatar
St. Elmo committed
85
        is_solved(opt_model) && break
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
86

St. Elmo's avatar
St. Elmo committed
87
88
        delete(opt_model, pfba_constraint)
        unregister(opt_model, :pfba_constraint)
St. Elmo's avatar
St. Elmo committed
89
90
    end

St. Elmo's avatar
St. Elmo committed
91
    is_solved(opt_model) || return nothing # pFBA failed
St. Elmo's avatar
St. Elmo committed
92
93
94
95
96

    return opt_model
end

"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
97
    parsimonious_flux_balance_analysis_vec(args...; kwargs...)
St. Elmo's avatar
St. Elmo committed
98

99
100
101
102
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`](@ref)
internally.
St. Elmo's avatar
St. Elmo committed
103
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
104
105
106
107
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
108
109
110
111
112

    return value.(opt_model[:x])
end

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

115
116
117
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`](@ref) internally.
St. Elmo's avatar
St. Elmo committed
118
"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
119
120
121
122
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
123

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