flux_balance_analysis.jl 3.18 KB
Newer Older
1
@testset "Flux balance analysis with CoreModel" begin
Sylvain Arreckx's avatar
Sylvain Arreckx committed
2
    cp = test_simpleLP()
3
    lp = flux_balance_analysis(cp, Tulip.Optimizer)
4
    @test termination_status(lp) == MOI.OPTIMAL
St. Elmo's avatar
St. Elmo committed
5
    sol = JuMP.value.(lp[:x])
6
    @test sol  [1.0, 2.0]
Sylvain Arreckx's avatar
Sylvain Arreckx committed
7
8
9

    # test the maximization of the objective
    cp = test_simpleLP2()
10
    lp = flux_balance_analysis(cp, Tulip.Optimizer)
11
    @test termination_status(lp) == MOI.OPTIMAL
St. Elmo's avatar
St. Elmo committed
12
    sol = JuMP.value.(lp[:x])
13
    @test sol  [-1.0, 2.0]
Sylvain Arreckx's avatar
Sylvain Arreckx committed
14
15

    # test with a more biologically meaningfull model
16
    cp = load_model(CoreModel, model_paths["iJR904.mat"])
17
    expected_optimum = 0.9219480950504393
Sylvain Arreckx's avatar
Sylvain Arreckx committed
18

19
    lp = flux_balance_analysis(cp, Tulip.Optimizer)
20
    @test termination_status(lp) == MOI.OPTIMAL
St. Elmo's avatar
St. Elmo committed
21
    sol = JuMP.value.(lp[:x])
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
22
23
    @test isapprox(objective_value(lp), expected_optimum, atol = TEST_TOLERANCE)
    @test isapprox(cp.c' * sol, expected_optimum, atol = TEST_TOLERANCE)
24
25

    # test the "nicer output" variants
26
    fluxes_vec = flux_balance_analysis_vec(cp, Tulip.Optimizer)
27
    @test all(fluxes_vec .== sol)
28
    fluxes_dict = flux_balance_analysis_dict(cp, Tulip.Optimizer)
29
30
    rxns = reactions(cp)
    @test all([fluxes_dict[rxns[i]] == sol[i] for i in eachindex(rxns)])
31
end
St. Elmo's avatar
St. Elmo committed
32

St. Elmo's avatar
St. Elmo committed
33
@testset "Flux balance analysis with StandardModel" begin
St. Elmo's avatar
St. Elmo committed
34

35
    model = load_model(StandardModel, model_paths["e_coli_core.json"])
St. Elmo's avatar
St. Elmo committed
36

St. Elmo's avatar
St. Elmo committed
37
38
39
40
    sol = flux_balance_analysis_dict(
        model,
        Tulip.Optimizer;
        modifications = [
St. Elmo's avatar
St. Elmo committed
41
42
            change_objective("BIOMASS_Ecoli_core_w_GAM"),
            change_constraint("EX_glc__D_e", -12, -12),
St. Elmo's avatar
St. Elmo committed
43
            change_sense(MAX_SENSE),
44
            change_optimizer_attribute("IPM_IterationsLimit", 110),
St. Elmo's avatar
St. Elmo committed
45
46
        ],
    )
St. Elmo's avatar
format    
St. Elmo committed
47
48
49
    @test isapprox(
        sol["BIOMASS_Ecoli_core_w_GAM"],
        1.0572509997013568,
50
        atol = TEST_TOLERANCE,
St. Elmo's avatar
format    
St. Elmo committed
51
    )
St. Elmo's avatar
St. Elmo committed
52

St. Elmo's avatar
St. Elmo committed
53
54
    pfl_frac = 0.8
    biomass_frac = 0.2
St. Elmo's avatar
St. Elmo committed
55
56
57
    sol_multi = flux_balance_analysis_dict(
        model,
        Tulip.Optimizer;
cylon-x's avatar
cylon-x committed
58
59
60
61
62
63
        modifications = [
            change_objective(
                ["BIOMASS_Ecoli_core_w_GAM", "PFL"];
                weights = [biomass_frac, pfl_frac],
            ),
        ],
St. Elmo's avatar
St. Elmo committed
64
65
66
67
    )
    @test isapprox(
        biomass_frac * sol_multi["BIOMASS_Ecoli_core_w_GAM"] + pfl_frac * sol_multi["PFL"],
        31.999999998962604,
68
        atol = TEST_TOLERANCE,
St. Elmo's avatar
St. Elmo committed
69
    )
St. Elmo's avatar
St. Elmo committed
70
end
St. Elmo's avatar
St. Elmo committed
71
72
73

@testset "Flux balance analysis with CoreModelCoupled" begin

74
    model = load_model(CoreModel, model_paths["e_coli_core.json"])
St. Elmo's avatar
St. Elmo committed
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91

    # assume coupling constraints of the form:
    # -γ ≤ vᵢ/μ  ≤ γ
    # I.e., enforces that the ratio between any reaction flux
    # and the growth rate is bounded by γ.
    γ = 40

    # construct coupling bounds
    nr = n_reactions(model)
    biomass_index = first(indexin(["BIOMASS_Ecoli_core_w_GAM"], reactions(model)))

    Cf = sparse(1.0I, nr, nr)
    Cf[:, biomass_index] .= -γ

    Cb = sparse(1.0I, nr, nr)
    Cb[:, biomass_index] .= γ

cylon-x's avatar
cylon-x committed
92
    C = [Cf; Cb]
St. Elmo's avatar
St. Elmo committed
93

cylon-x's avatar
cylon-x committed
94
    clb = spzeros(2 * nr)
St. Elmo's avatar
St. Elmo committed
95
    clb[1:nr] .= -1000.0
cylon-x's avatar
cylon-x committed
96
    cub = spzeros(2 * nr)
St. Elmo's avatar
St. Elmo committed
97
98
99
100
101
102
103
    cub[nr+1:end] .= 1000

    cmodel = CoreModelCoupled(model, C, clb, cub) # construct

    dc = flux_balance_analysis_dict(cmodel, Tulip.Optimizer)
    @test isapprox(dc["BIOMASS_Ecoli_core_w_GAM"], 0.665585699298256, atol = TEST_TOLERANCE)
end