Making adjustments to the model

Typically, we do not need to solve the models as they come from the authors (someone else already did that!), but we want to perform various perturbations in the model structure and conditions, and explore how the model behaves in the changed conditions.

With COBREXA, there are 2 different approaches that one can take:

  1. We can change the model structure, and use the changed metabolic model. This is better for doing simple and small, but systematic modifications, such as removing metabolites, adding reactions, etc.
  2. We can intercept the pipeline that converts the metabolic model to constraints and to the optimizer representation, and make modifications along that way. This is better suited to making global model adjustments, such as using combined objectives, adding reaction-coupling constraints, and combining multiple models into a bigger one.

Here we demonstrate the first, "modeling" approach. The main advantage of this approach is that the modified model is still a FBC model, and we can export, save and share it via the AbstractFBCModels interace. The main disadvantage is that the "common" FBC model interface does not easily express various complicated constructions (communities, reaction coupling, enzyme constraints, etc.) – see the example about modifying the constraints for more details.

Getting the base model

using COBREXA

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

import JSONFBCModels
[ Info: using cached `e_coli_core.json'

For applying the modifications, we will use the canonical model as exported from package AbstractFBCModels. There are other possibilities, but the canonical one is easiest to use for common purposes.

import AbstractFBCModels.CanonicalModel as CM

We can now load the model:

model = convert(CM.Model, load_model("e_coli_core.json"))
AbstractFBCModels.CanonicalModel.Model(
  reactions = Dict{String, AbstractFBCModels.CanonicalModel.Reaction}("ACALD" =…
  metabolites = Dict{String, AbstractFBCModels.CanonicalModel.Metabolite}("glu_…
  genes = Dict{String, AbstractFBCModels.CanonicalModel.Gene}("b4301" => Abstra…
  couplings = Dict{String, AbstractFBCModels.CanonicalModel.Coupling}(),
)

The canonical model is quite easy to work with, made basically of the most accessible Julia structures possible. For example, we can look at a reaction as such:

model.reactions["PFK"]
AbstractFBCModels.CanonicalModel.Reaction(
  name = "Phosphofructokinase",
  lower_bound = 0.0,
  upper_bound = 1000.0,
  stoichiometry = Dict("adp_c" => 1.0, "atp_c" => -1.0, "f6p_c" => -1.0, "fdp_c…
  objective_coefficient = 0.0,
  gene_association_dnf = [["b1723"], ["b3916"]],
  annotations = Dict("bigg.reaction" => ["PFK"], "metanetx.reaction" => ["MNXR1…
  notes = Dict("original_bigg_ids" => ["PFK"]),
)

model.reactions["CS"].stoichiometry
Dict{String, Float64} with 6 entries:
  "coa_c"   => 1.0
  "h2o_c"   => -1.0
  "accoa_c" => -1.0
  "cit_c"   => 1.0
  "oaa_c"   => -1.0
  "h_c"     => 1.0
Create custom model types!

For some applications, CanonicalModel might be too restrictive. Creating a custom model type that perfectly fits a use-case can be done simply by overloading several functions. The documentation of AbstractFBCModels describes the process closer. Further, all model types that adhere to the AbstractFBCModels' interface will "just work" with all analysis functions in COBREXA!

Running FBA on modified models

Since the canonical model is completely mutable, we can change it in any way we like and feed the result directly into flux_balance_analysis. Let's first find a "original" solution, so that we have a base solution for comparing:

import HiGHS

base_solution = flux_balance_analysis(model, optimizer = HiGHS.Optimizer)
base_solution.objective
0.8739215069684304

Now, for example, we can limit the intake of glucose by the model:

model.reactions["EX_glc__D_e"]
AbstractFBCModels.CanonicalModel.Reaction(
  name = "D-Glucose exchange",
  lower_bound = -10.0,
  upper_bound = 1000.0,
  stoichiometry = Dict("glc__D_e" => -1.0),
  objective_coefficient = 0.0,
  gene_association_dnf = nothing,
  annotations = Dict("sabiork" => ["7002", "601"], "metanetx.reaction" => ["MNX…
  notes = Dict("original_bigg_ids" => ["EX_glc_e"]),
)

Since the original intake limit is 10 units, let's try limiting that to 5:

model.reactions["EX_glc__D_e"].lower_bound = -5.0
-5.0

...and solve the modified model:

low_glucose_solution = flux_balance_analysis(model, optimizer = HiGHS.Optimizer)
low_glucose_solution.objective
0.415597775092906

Preventing reference-based sharing problems with deepcopy

People often want to try different perturbations with a single base model. It would therefore look feasible to save the "unmodified" model in a single variable, and make copies of that with the modifications applied. Let's observe what happens:

base_model = convert(CM.Model, load_model("e_coli_core.json")) # load the base

modified_model = base_model # seemingly make a "copy" for modification

modified_model.reactions["EX_glc__D_e"].lower_bound = -123.0 # modify the glucose intake limit
-123.0

Surprisingly, the base model got modified too!

base_model.reactions["EX_glc__D_e"]
AbstractFBCModels.CanonicalModel.Reaction(
  name = "D-Glucose exchange",
  lower_bound = -123.0,
  upper_bound = 1000.0,
  stoichiometry = Dict("glc__D_e" => -1.0),
  objective_coefficient = 0.0,
  gene_association_dnf = nothing,
  annotations = Dict("sabiork" => ["7002", "601"], "metanetx.reaction" => ["MNX…
  notes = Dict("original_bigg_ids" => ["EX_glc_e"]),
)

This is because Julia uses reference-based sharing whenever anything mutable is copied using the = operator. While this is extremely useful in many scenarios for data processing efficiency and computational speed, it unfortunately breaks this simple use-case.

To fix this situation, we must always remember to make an actual copy of the model data, by either carefully copying the changed parts (e.g., using a similar approach as with the "shallow" copy()), or simply by copying the whole model structure as is with deepcopy(). Let's try again:

base_model = convert(CM.Model, load_model("e_coli_core.json"))
modified_model = deepcopy(base_model) # this forces an actual copy of the data
modified_model.reactions["EX_glc__D_e"].lower_bound = -123.0
-123.0

With deepcopy, the result works as intended:

(
    modified_model.reactions["EX_glc__D_e"].lower_bound,
    base_model.reactions["EX_glc__D_e"].lower_bound,
)
(-123.0, -10.0)
Avoid overwriting base models when using in-place modifications

Whenever changing a copy of the model, check that the base model is not inadvertently changed via a reference. Always use some copy mechanism such as copy or deepcopy to prevent the default reference-based sharing.

Observing the differences

We already have a base_solution and low_glucose_solution from above. What is the easiest way to see what has changed? We can quite easily compute squared distance between all dictionary entries using Julia function for merging dictionaries (called mergewith).

With that, we can extract the plain difference in fluxes:

flux_differences = mergewith(-, base_solution.fluxes, low_glucose_solution.fluxes)
ConstraintTrees.Tree{Float64} with 95 elements:
  :ACALD                    => 0.0
  :ACALDt                   => 0.0
  :ACKr                     => 0.0
  :ACONTa                   => 2.56379
  :ACONTb                   => 2.56379
  :ACt2r                    => 0.0
  :ADK1                     => 0.0
  :AKGDH                    => 2.06931
  :AKGt2r                   => 0.0
  :ALCD2x                   => 0.0
  :ATPM                     => 0.0
  :ATPS4r                   => 20.6429
  :BIOMASS_Ecoli_core_w_GAM => 0.458324
  :CO2t                     => -10.4958
  :CS                       => 2.56379
  :CYTBD                    => 19.9319
  :D_LACt2                  => 0.0
  :ENO                      => 7.13113
  :ETOHt2r                  => 0.0
  ⋮                         => ⋮

...and see what were the biggest directional differences:

sort(collect(flux_differences), by = last)
95-element Vector{Pair{Symbol, Union{Float64, ConstraintTrees.Tree{Float64}}}}:
        :H2Ot => -13.834417820645202
        :CO2t => -10.495804428171123
     :EX_o2_e => -9.965936361749836
         :PGK => -7.816778467951078
         :PGM => -7.131126165065294
 :EX_glc__D_e => -5.0
    :EX_nh4_e => -2.499147645170859
       :GLUDy => -2.381954266930288
      :SUCOAS => -2.0693072054856274
     :EX_pi_e => -1.6860355124504807
              ⋮
         :ENO => 7.131126165065294
        :GAPD => 7.816778467951078
      :EX_h_e => 9.193974061423127
         :O2t => 9.965936361749836
    :EX_co2_e => 10.495804428171123
    :EX_h2o_e => 13.834417820645202
      :NADH16 => 17.86256551801404
       :CYTBD => 19.931872723499673
      :ATPS4r => 20.642858106791966

...or compute the squared distance, to see the "absolute" changes:

flux_changes =
    mergewith((x, y) -> (x - y)^2, base_solution.fluxes, low_glucose_solution.fluxes)
ConstraintTrees.Tree{Float64} with 95 elements:
  :ACALD                    => 0.0
  :ACALDt                   => 0.0
  :ACKr                     => 0.0
  :ACONTa                   => 6.57303
  :ACONTb                   => 6.57303
  :ACt2r                    => 0.0
  :ADK1                     => 0.0
  :AKGDH                    => 4.28203
  :AKGt2r                   => 0.0
  :ALCD2x                   => 0.0
  :ATPM                     => 0.0
  :ATPS4r                   => 426.128
  :BIOMASS_Ecoli_core_w_GAM => 0.210061
  :CO2t                     => 110.162
  :CS                       => 6.57303
  :CYTBD                    => 397.28
  :D_LACt2                  => 0.0
  :ENO                      => 50.853
  :ETOHt2r                  => 0.0
  ⋮                         => ⋮

...and again see what changed most:

sort(collect(flux_changes), by = last)
95-element Vector{Pair{Symbol, Union{Float64, ConstraintTrees.Tree{Float64}}}}:
    :ACALD => 0.0
   :ACALDt => 0.0
     :ACKr => 0.0
    :ACt2r => 0.0
     :ADK1 => 0.0
   :AKGt2r => 0.0
   :ALCD2x => 0.0
     :ATPM => 0.0
  :D_LACt2 => 0.0
  :ETOHt2r => 0.0
           ⋮
  :EX_o2_e => 99.31988756644756
      :O2t => 99.31988756644756
     :CO2t => 110.16191059441655
 :EX_co2_e => 110.16191059441655
 :EX_h2o_e => 191.39111643618554
     :H2Ot => 191.39111643618554
   :NADH16 => 319.07124688534424
    :CYTBD => 397.27955026579025
   :ATPS4r => 426.12759081714677
Always use a uniquely defined flux solutions for flux comparisons

Since the usual flux balance allows a lot of freedom in the "solved" flux and the only value that is "reproducible" by the analysis is the objective, one should never compare the flux distributions directly. Typically, that may result in false-positive (and sometimes false-negative) differences. Use e.g. parsimonious FBA to obtain uniquely determined and safely comparable flux solutions.

Coupling constraints

Some model types support additional constraints over the reaction fluxes, which are historically called "coupling". These allow to e.g. place a bound on a total flux through several reactions.

Canonical model supports these as "couplings":

model.couplings["total_energy_intake"] = CM.Coupling(
    lower_bound = 0,
    upper_bound = 5,
    reaction_weights = Dict("EX_glc__D_e" => -1.0, "EX_fru_e" => -1.0, "EX_pyr_e" => -1.0),
)
AbstractFBCModels.CanonicalModel.Coupling(
  name = nothing,
  reaction_weights = Dict("EX_pyr_e" => -1.0, "EX_glc__D_e" => -1.0, "EX_fru_e"…
  lower_bound = 0.0,
  upper_bound = 5.0,
  annotations = Dict{String, Vector{String}}(),
  notes = Dict{String, Vector{String}}(),
)

The values of any coupling constraints can be inspected directly in the solved model:

solution_with_coupling = flux_balance_analysis(model, optimizer = HiGHS.Optimizer)

solution_with_coupling.coupling.total_energy_intake
5.0

This page was generated using Literate.jl.