fba.jl 3.96 KB
Newer Older
1
"""
2
    flux_balance_analysis(model::M, optimizer) where {M<:MetabolicModel}
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
3
4

Flux balance analysis solves the following problem for the input `model`:
Sylvain Arreckx's avatar
Sylvain Arreckx committed
5
6
7
8
9
```
max cᵀx
s.t. S x = b
     xₗ ≤ x ≤ xᵤ
```
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
10

11
Returns a solved model from [`optimize_model`](@ref).
Sylvain Arreckx's avatar
Sylvain Arreckx committed
12
"""
13
14
flux_balance_analysis(model::M, optimizer) where {M<:MetabolicModel} =
    optimize_model(model, optimizer; sense = MOI.MAX_SENSE)
15
16

"""
St. Elmo's avatar
St. Elmo committed
17
    fba(model::CobraModel, optimizer; objective_func::Union{Reaction, Array{Reaction, 1}}=Reaction[], weights=Float64[], solver_attributes=Dict{Any, Any}(), constraints=Dict{String, Tuple{Float64,Float64}}())
18

St. Elmo's avatar
St. Elmo committed
19
Run flux balance analysis (FBA) on the `model` optionally specifying `objective_rxn(s)` and their `weights` (empty `weights` mean equal weighting per reaction).
20
Optionally also specify any additional flux constraints with `constraints`, a dictionary mapping reaction `id`s to tuples of (lb, ub) flux constraints.
21
Note, the `optimizer` must be set to perform the analysis, any JuMP solver will work.
22
23
24
The `solver_attributes` can also be specified in the form of a dictionary where each (key, value) pair will be passed to `set_optimizer_attribute(cbmodel, key, value)`.
This function builds the optimization problem from the model, and hence uses the constraints implied by the model object.
Returns a dictionary of reaction `id`s mapped to fluxes if solved successfully, otherwise an empty dictionary.
25

26
27
28
29
30
31
32
33
# Example
```
optimizer = Gurobi.Optimizer
atts = Dict("OutputFlag" => 0)
model = CobraTools.read_model("iJO1366.json")
biomass = findfirst(model.reactions, "BIOMASS_Ec_iJO1366_WT_53p95M")
sol = fba(model, biomass, optimizer; solver_attributes=atts)
```
34
"""
35
36
37
function fba(
    model::CobraModel,
    optimizer;
St. Elmo's avatar
St. Elmo committed
38
    objective_func::Union{Reaction,Array{Reaction,1}} = Reaction[],
39
40
41
    weights = Float64[],
    solver_attributes = Dict{Any,Any}(),
    constraints = Dict{String,Tuple{Float64,Float64}}(),
St. Elmo's avatar
St. Elmo committed
42
    sense = MOI.MAX_SENSE,
43
)
St. Elmo's avatar
St. Elmo committed
44
    # get core optimization problem
45
    cbm, v, mb, lbcons, ubcons = make_optimization_model(model, optimizer, sense = sense)
46

St. Elmo's avatar
St. Elmo committed
47
    # modify core optimization problem according to user specifications
48
49
50
51
52
    if !isempty(solver_attributes) # set other attributes
        for (k, val) in solver_attributes
            set_optimizer_attribute(cbm, k, val)
        end
    end
53

54
55
56
    # set additional constraints
    for (rxnid, con) in constraints
        ind = model.reactions[findfirst(model.reactions, rxnid)]
57
        set_bound(ind, lbcons, ubcons; lb = con[1], ub = con[2])
58
    end
59

St. Elmo's avatar
St. Elmo committed
60
    # if an objective function is supplied, modify the default objective
St. Elmo's avatar
St. Elmo committed
61
    if typeof(objective_func) == Reaction || !isempty(objective_func)
St. Elmo's avatar
St. Elmo committed
62
        # ensure that an array of objective indices are fed in
63
64
        if typeof(objective_func) == Reaction
            objective_indices = [model[objective_func]]
St. Elmo's avatar
St. Elmo committed
65
        else
66
            objective_indices = [model[rxn] for rxn in objective_func]
St. Elmo's avatar
St. Elmo committed
67
        end
68

St. Elmo's avatar
St. Elmo committed
69
70
71
72
        if isempty(weights)
            weights = ones(length(objective_indices))
        end
        opt_weights = zeros(length(model.reactions))
73

St. Elmo's avatar
St. Elmo committed
74
        # update the objective function tracker
75
        # don't update model objective function - silly thing to do
St. Elmo's avatar
St. Elmo committed
76
77
78
        wcounter = 1
        for i in eachindex(model.reactions)
            if i in objective_indices
79
                # model.reactions[i].objective_coefficient = weights[wcounter]
St. Elmo's avatar
St. Elmo committed
80
81
                opt_weights[i] = weights[wcounter]
                wcounter += 1
St. Elmo's avatar
St. Elmo committed
82
                # else
83
                # model.reactions[i].objective_coefficient = 0.0
St. Elmo's avatar
St. Elmo committed
84
            end
85
86
        end

87
        @objective(cbm, sense, sum(opt_weights[i] * v[i] for i in objective_indices))
88
    else # use default objective
89
        # automatically assigned by make_optimization_model
St. Elmo's avatar
St. Elmo committed
90
    end
St. Elmo's avatar
St. Elmo committed
91

92
93
94
95
96
97
98
99
100
101
102
103
104
    optimize!(cbm)

    status = (
        termination_status(cbm) == MOI.OPTIMAL ||
        termination_status(cbm) == MOI.LOCALLY_SOLVED
    )

    if status
        return map_fluxes(v, model)
    else
        @warn "Optimization issues occurred."
        return Dict{String,Float64}()
    end
105
end