Model

Aim

Our model refers to a standard existing Saccharomyces cerevisiae model and aims to predict optimal flux distributions within metabolic reactions by means of flux balance analysis (FBA). Additionally, this model contains enzyme constraints that define a GECKO model. Mainly, our model tries to reconstruct all the existing metabolic pathways with their fluxes and enzyme levels to optimize the phosphate uptake into the cell. These mentioned specific variations need to be adapted to the model. In addition to that, phosphate levels need to be introduced to the model by collecting experimental data on phosphate levels in cells under different conditions. Those phosphate levels can easily be added as mass balance equations that have been constructed for different phosphate species (e.g. ATP, Pi) based on known biochemical reactions and transport reactions.

Model structure

We used the S. cerevisiae model and added aforementioned requirements to it. It was planned to use the experimentally collected data to receive exact values that relate more directly to our experiments.

Contribution to future iGEM teams

Future iGEM teams could benefit from our model as this is an existing S. cerevisiae simulation containing all the important metabolic reactions and gene products. For every need there is, the model can be adjusted by adding the desired variable, parameter and/or constraint, and if that is not enough, experimental data extracted directly from the laboratory phase can be inserted too.

Creating the code

To analyze the model data, we used the programming language ‘Julia’ and its ‘COBREXA’.jl library which is a toolkit for working with large constraint-based metabolic models. In addition, we chose ‘Gurobi’ as our optimization solver as it can handle linear programming (LP), quadratic programming (QP), quadratically constrained programming (QCP), mixed-integer linear programming (MILP), mixed-integer quadratic programming (MIQP), and mixed-integer quadratically constrained programming (MIQCP).

Understanding the code (step by step)

We need to use COBREXA as our library, Gurobi as our mathematical optimization solver and the package ‘Distributions’ to calculate probability distributions and statistics. We need to load the model “iND750.json”, which is a Saccharomyces cerevisiae model from literature ¹. Basically, this is just an enzyme constrained model, also called a GECKO model, bounding gene product concentrations and supplying molar mass of gene products. After establishing such an enzyme constrained model, we can perform an FBA to receive optimal fluxes from reactions to enzymes. The enzymes have been mapped by creating the variable “reaction_isozymes” that holds the dictionary storing all the isozymes involved in all the reactions.
Finally, after creating all the required variables, we are able to solve the whole model to receive all the fluxes listed then in one dictionary to have an overview of all the solutions. Furthermore, other variables like substrate concentrations can then be added to this model. Lastly, this model uses by default glucose as the carbon source which can also be modified with other carbon sources and the model can then be resolved using another carbon source. All in all, this GECKO model explores the metabolic behavior of the organism under different conditions, focusing on the carbon source used for metabolism.

The code

The complete model code:

using COBREXA, Gurobi, Distributions
model = load_model(StandardModel, "iND750.json") # this is a saccharomyces cerevisiae model 

# remove all non-directional bounds on the model, except for ATPM!
"ATPM" in reactions(model)
for rid in reactions(model)
    rid == "ATPM" && continue
    lb = model.reactions[rid].lb
    ub = model.reactions[rid].ub
    model.reactions[rid].lb = lb < 0 ? -1000.0 : 0.0
    model.reactions[rid].ub = ub > 0 ? 1000.0 : 0.0
end

#=
Enzyme constrained models are called GECKO models in COBREXA 1.x. The most
important jargon to keep in mind is that an enzyme can be composed of multiple
gene products (these are what the genes in a model encode for). 

To create a GECKO model, these dictionaries are required:

1. reaction_isozymes = Dict{String,Vector{Isozyme}}. This maps reactions to
enzymes (sometimes complexes). Multiple isozymes may catalyze the same reaction.
2. gene_product_bounds = Dict{String,Tuple{Float64,Float64}}. This bounds the
actual gene product concentrations.
3. gene_product_molar_mass = Dict{String,Float64}. This supplies the molar mass of a gene product.
4. gene_product_mass_group = Dict{String,String}. This maps gene products to specific capacity bound IDs.
5. gene_product_mass_group_bound = Dict{String,Float64}. This is the numeric bound of each capacity ID.

Isozymes looks like this: 
mutable struct Isozyme
    gene_product_count::Dict{String,Int} 
    kcat_forward::Float64 
    kcat_reverse::Float64
end
=#

kcat() = 10^rand(Normal(1,1)) # realistic for kcats with units of 1/s
mw() = rand(Uniform(10, 100)) #  realistic for proteins with units of kDa

# first create all the isozymes
reaction_isozymes = Dict{String,Vector{Isozyme}}()
for rid in reactions(model)
    isos = Isozyme[]
    grrs = reaction_gene_association(model, rid)
    isnothing(grrs) && continue # only assign isozymes to reactions involving genes    
    for grr in grrs
        kfor = kcat() * 3600 # change unit to 1/h
        kback = kcat() * 3600
        push!(isos, Isozyme(Dict(k => 1 for k in grr), kfor, kback))
    end
    reaction_isozymes[rid] = isos
end

gene_product_mass_group = Dict{String,String}()
gene_product_molar_mass = Dict{String,Float64}()
gene_product_bounds = Dict{String,Tuple{Float64,Float64}}()
for gid in genes(model)
    gene_product_mass_group[gid] = "capacity_bound_"*(rand() < 0.5 ? "A" : "B")
    gene_product_molar_mass[gid] = mw()
    gene_product_bounds[gid] = (0.0, 10000.0) # unit in mmol/gDW
end
gene_product_mass_group_bound = Dict("capacity_bound_A" => 0.5, "capacity_bound_B" => 0.3) # units g/gDW

gmodel = make_gecko_model(
    model;
    reaction_isozymes,
    gene_product_bounds,
    gene_product_mass_group,
    gene_product_molar_mass,
    gene_product_mass_group_bound,
)

# solve it
opt_model = flux_balance_analysis(gmodel, Gurobi.Optimizer) # uses glucose by default
sol = flux_dict(gmodel, opt_model)
flux_summary(sol)

# now change the carbon source
opt_model = flux_balance_analysis(gmodel, Gurobi.Optimizer; modifications=[
    change_constraint("EX_glc__D_e"; lb = 0, ub = 0)
    change_constraint("EX_xyl__D_e"; lb = -12, ub = -12)
])
sol = flux_dict(gmodel, opt_model)
flux_summary(sol)

# etc.

---

1. Natalie C. Duarte, Markus J. Herrgård, and Bernhard Ø. Palsson. Reconstruction and Validation of Saccharomyces cerevisiae iND750, a Fully Compartmentalized Genome-Scale Metabolic Model. Genome Research, 2004.↩