Screening through many model variants

screen is a function that you can use to run many different simulations on many different variants of the model efficiently in parallel.

using COBREXA

download_model(
    "http://bigg.ucsd.edu/static/models/e_coli_core.json",
    "e_coli_core.json",
    "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
)

import JSONFBCModels, HiGHS
import ConstraintTrees as C

model = load_model("e_coli_core.json")

JSONFBCModels.JSONFBCModel(#= 95 reactions, 72 metabolites =#)

Screening is done as follows: We specify a range of values to examine (in the example below to an actual range of -10 to 0), but any vector-like list of things can be used), and write a short "analysis" function that takes one of the values from the range as an argument and runs an analysis for that given value.

Here, we solve 11 different optimization problems for different bounds of the oxygen exchange:

screen(-10:0) do o2_bound
    c = flux_balance_constraints(model)
    c.fluxes.EX_o2_e.bound = C.EqualTo(o2_bound)
    optimized_values(
        c,
        objective = c.objective.value,
        optimizer = HiGHS.Optimizer,
        output = c.objective,
    )
end

11-element Vector{Float64}:
 0.5590505997120256
 0.5265539556675367
 0.4940573116230476
 0.4615606675785586
 0.4276455515060357
 0.3916484512109149
 0.35565135091579403
 0.31965425062067315
 0.28365715032555233
 0.24766005003043148
 0.2116629497353106

Screening in multiple dimensions

The simplest way to screen through multi-dimensional arrays is to use an iterator product. In the result, we receive a whole matrix of results instead of a vector.

screen(Iterators.product(-10:2:0, 0:2:10)) do (o2_bound, co2_bound)
    c = flux_balance_constraints(model)
    c.fluxes.EX_o2_e.bound = C.EqualTo(o2_bound)
    c.fluxes.EX_co2_e.bound = C.EqualTo(co2_bound)
    optimized_values(
        c,
        objective = c.objective.value,
        optimizer = HiGHS.Optimizer,
        output = c.objective,
    )
end

6×6 Matrix{Float64}:
 0.411525  0.505191  0.549994  0.559051  0.559051  0.556413
 0.451394  0.494057  0.492882  0.49012   0.487358  0.484596
 0.42659   0.423828  0.421066  0.418304  0.415542  0.412779
 0.354774  0.352012  0.349249  0.346487  0.343725  0.340963
 0.282957  0.280195  0.277433  0.274671  0.271909  0.269146
 0.211141  0.208378  0.205616  0.202854  0.200092  0.19733

Screening through non-numeric properties

If the problem at hand can not be expressed with mere ranges, we can specify vectors with any values. The following code examines the inter-dependency of blocking selected transport reactions together with selected exchanges, with 2 different bounds that implement the block. Because of 3-dimensional input, the result is expectably a 3-dimensional array:

screen(
    Iterators.product(
        [:H2Ot, :CO2t, :O2t, :NH4t], # a set of transport reactions
        [:EX_h2o_e, :EX_co2_e, :EX_o2_e, :EX_nh4_e], # a set of exchanges
        [C.Between(-0.1, 0.1), C.Between(-1, 1)], # bounds
    ),
) do (transport, exchange, bound)
    c = flux_balance_constraints(model)
    c.fluxes[transport].bound = 5 * bound
    c.fluxes[exchange].bound = 3 * bound
    optimized_values(
        c,
        objective = c.objective.value,
        optimizer = HiGHS.Optimizer,
        output = c.objective,
    )
end

4×4×2 Array{Float64, 3}:
[:, :, 1] =
 0.390761  0.394089  0.204711  0.0550176
 0.390029  0.46905   0.222462  0.0550176
 0.212193  0.229509  0.222462  0.0550176
 0.091696  0.091696  0.091696  0.0550176

[:, :, 2] =
 0.454097  0.494176  0.319654  0.494176
 0.454097  0.530039  0.319654  0.550176
 0.391648  0.391648  0.319654  0.391648
 0.454097  0.530039  0.319654  0.550176

Annotating the results

For reasons of tidiness, it is advisable to annotate all values with their actual meaning directly in the arrays.

With screen, the simplest (and usually sufficiently effective) way to do that is to return Pairs with annotation keys instead of plain values:

screen(
    Iterators.product(
        [:H2Ot, :CO2t, :O2t, :NH4t],
        [:EX_h2o_e, :EX_co2_e, :EX_o2_e, :EX_nh4_e],
        [C.Between(-0.1, 0.1), C.Between(-1, 1)],
    ),
) do (transport, exchange, bound)
    c = flux_balance_constraints(model)
    c.fluxes[transport].bound = 5 * bound
    c.fluxes[exchange].bound = 3 * bound
    (transport, exchange, bound.upper) => optimized_values(
        c,
        objective = c.objective.value,
        optimizer = HiGHS.Optimizer,
        output = c.objective,
    )
end

4×4×2 Array{Pair{Tuple{Symbol, Symbol, Float64}, Float64}, 3}:
[:, :, 1] =
 (:H2Ot, :EX_h2o_e, 0.1)=>0.390761  …  (:H2Ot, :EX_nh4_e, 0.1)=>0.0550176
 (:CO2t, :EX_h2o_e, 0.1)=>0.390029     (:CO2t, :EX_nh4_e, 0.1)=>0.0550176
  (:O2t, :EX_h2o_e, 0.1)=>0.212193      (:O2t, :EX_nh4_e, 0.1)=>0.0550176
 (:NH4t, :EX_h2o_e, 0.1)=>0.091696     (:NH4t, :EX_nh4_e, 0.1)=>0.0550176

[:, :, 2] =
 (:H2Ot, :EX_h2o_e, 1.0)=>0.454097  …  (:H2Ot, :EX_nh4_e, 1.0)=>0.494176
 (:CO2t, :EX_h2o_e, 1.0)=>0.454097     (:CO2t, :EX_nh4_e, 1.0)=>0.550176
  (:O2t, :EX_h2o_e, 1.0)=>0.391648      (:O2t, :EX_nh4_e, 1.0)=>0.391648
 (:NH4t, :EX_h2o_e, 1.0)=>0.454097     (:NH4t, :EX_nh4_e, 1.0)=>0.550176

Notably, this approach makes various indexing and ordering errors quite improbable.

This page was generated using Literate.jl.