Search code examples
pythoncbmpy

How to set a flux ratio as a constraint?

In some datasets, I sometimes observe fixed flux ratios which I would like to incorporate into my simulations. How could I do this in CBMPy?

For example, I have the model from here and would now like to constrain the ratio of succinate efflux and pyruvate efflux to 2.0. I know how to set constraints on individual reactions:

import cbmpy

# downloaded from http://bigg.ucsd.edu/models/e_coli_core
ecoli = cbmpy.CBRead.readSBML3FBC('e_coli_core.xml')

ecoli.setReactionBounds('R_EX_pyr_e', 1.0, 1000.0)
ecoli.setReactionBounds('R_EX_succ_e', 2.0, 1000.0)

# solve the model
cbmpy.doFBA(ecoli)

# get all reaction values
solution = ecoli.getReactionValues()
print(solution['R_EX_pyr_e'])
print(solution['R_EX_succ_e'])

For this case the ratio is correct, but how can I add it as a constraint that it will be fulfilled for all conditions?

Solution

  • That is indeed a common approach in Flux Balance Analysis (FBA) and you can use the function addUserConstraint to accomplish this.

    The entire code sample could look like this (explanation below):

    import cbmpy as cbm

# downloaded from http://bigg.ucsd.edu/models/e_coli_core
ecoli = cbm.CBRead.readSBML3FBC('e_coli_core.xml')

# make a clone of the original model
ecoli_ratio = ecoli.clone()

# add the desired user constraint; explanation follows below
ecoli_ratio.addUserConstraint("pyr_succ_ratio", fluxes=[(1.0, 'R_EX_pyr_e' ),(-0.5, 'R_EX_succ_e')], operator='=', rhs=0.0)

# now we have to set only one flux bound; if you think it is naturally excreted, this step is not needed
ecoli_ratio.setReactionBounds('R_EX_succ_e', 4.0, cbm.INF)

cbm.doFBA(ecoli_ratio)
solution = ecoli_ratio.getReactionValues()
print("{}: {}".format("succinate excretion rate", solution['R_EX_succ_e']))
print("{}: {}".format("pyruvate excretion rate", solution['R_EX_pyr_e']))

    This will print

    succinate excretion rate: 4.0
pyruvate excretion rate: 2.0

    As you can see the ratio is 2.0 as desired.

    Bit more explanation:

    The constraint is 

    J_succ / J_pyr = 2.0

    which can be rewritten to 

    J_succ = 2.0 J_pyr

    and finally

    J_pyr - 0.5 J_succ = 0

    That's exactly what we pass to fluxes in addUserConstraint:

    fluxes=[(1.0, 'R_EX_pyr_e' ),(-0.5, 'R_EX_succ_e')], operator='=', rhs=0.0)

    You can check the user defined constraints by printing:

    print(ecoli_ratio.user_constraints)
{'pyr_succ_ratio': {'operator': 'E', 'rhs': 0.0, 'fluxes': [(1.0, 'R_EX_pyr_e'), (-0.5, 'R_EX_succ_e')]}}

    As this is a dictionary you can delete the constraint by simply doing:

    del ecoli_ratio.user_constraints['pyr_succ_ratio']
print(ecoli_ratio.user_constraints)
{}

    but I highly recommend to create a clone everytime you introduce major changes to a model.