Using positive lower bound produces unexpected behaviour

[behzadk]

Hello, I came across this issue while trying to use lower bound constraints to force cross feeding. I am trying to ensure a strain is outputting a non-essential amino acid that it can synthesise by setting the lower bound constraint to positive. When I try this using COMETSPy, the simulation seemingly ignores the positive lower bound.

I have provided some code using the Cobra e_coli 'text_book' model. I run three simulations with different lower bound constraints

lower bound negative: e_coli.change_bounds("EX_glc__D_e", -1000, 1000) lower bound zero: e_coli.change_bounds("EX_glc__D_e", 0.0, 1000) lower bound positive: e_coli.change_bounds("EX_glc__D_e", 1.0, 1000)

I would expect to see growth only in the first instance, no growth in the second, and no growth or an error in the third.

The image below shows the biomass of the three simulations (left) and the glucose concentrations (right). We have this unexpected growth when the lower bound is set to positive.

Any suggestions on how I could tackle this would be great, I'll continue to have a look into where the problem might be!

import cobra
import cobra.test
import cometspy as c
import os
import gurobipy
import matplotlib.pyplot as plt

def set_layout_metabolites(layout):
    layout.set_specific_metabolite("glc__D_e", 0.011)
    layout.set_specific_metabolite("o2_e", 10.0)

    # Add the rest of nutrients unlimited (ammonia, phosphate, water and protons)
    layout.set_specific_metabolite("nh4_e", 1000)
    layout.set_specific_metabolite("pi_e", 1000)
    layout.set_specific_metabolite("h2o_e", 1000)
    layout.set_specific_metabolite("h_e", 1000)

def generate_sim_params():
    # Set the parameters that are different from the default
    sim_params = c.params()
    sim_params.set_param("defaultVmax", 18.5)
    sim_params.set_param("defaultKm", 0.000015)
    sim_params.set_param("maxCycles", 1000)
    sim_params.set_param("timeStep", 0.1)
    sim_params.set_param("spaceWidth", 1)
    sim_params.set_param("maxSpaceBiomass", 10)
    sim_params.set_param("minSpaceBiomass", 1e-11)
    sim_params.set_param("writeMediaLog", True)

    return sim_params

def build_model():
    # create the model using CobraPy functionality
    e_coli_cobra = cobra.test.create_test_model("textbook")

    # use the loaded model to build a comets model
    e_coli = c.model(e_coli_cobra)

    # set the model's initial biomass
    e_coli.initial_pop = [0, 0, 5e-6]

    return e_coli

def plot_glucose(experiment, ax):
    media = experiment.media.copy()
    media = media.loc[media["metabolite"] == "glc__D_e"]
    media.groupby("metabolite").plot(x="cycle", ax=ax, y="conc_mmol", alpha=1.0)

def main():
    # Create empty 1x1 layout
    test_tube = c.layout()

    e_coli = build_model()

    # add it to the test_tube
    test_tube.add_model(e_coli)
    set_layout_metabolites(test_tube)

    sim_params = generate_sim_params()

    ## Run three experiments with different lower
    ## bound constraints on glucose
    # Run with open exchange
    e_coli.change_bounds("EX_glc__D_e", -1000, 1000)
    experiment_open = c.comets(test_tube, sim_params)
    experiment_open.run(delete_files=True)

    # Run with uptake blocked
    experiment_zero = c.comets(test_tube, sim_params)
    e_coli.change_bounds("EX_glc__D_e", 0.0, 1000)
    experiment_zero.run(delete_files=True)

    # Run with forced production
    experiment_invl = c.comets(test_tube, sim_params)
    e_coli.change_bounds("EX_glc__D_e", 1.0, 1000)
    experiment_invl.run(delete_files=True)

    # plotting
    fig, axes = plt.subplots(figsize=(15, 7), ncols=2)

    ax = axes[0]
    experiment_open.total_biomass.plot(
        x="cycle", ax=ax, label="EX_glc__D_e lower bound = -1000.0"
    )
    experiment_zero.total_biomass.plot(
        x="cycle", ax=ax, label="EX_glc__D_e lower bound = 0.0"
    )
    experiment_invl.total_biomass.plot(
        x="cycle", ax=ax, label="EX_glc__D_e lower bound = 1.0"
    )
    ax.set_ylabel("Biomass (gr.)")
    ax.legend(
        [
            "EX_glc__D_e lower bound = -1000.0",
            "EX_glc__D_e lower bound = 0.0",
            "EX_glc__D_e lower bound = 1.0",
        ]
    )
    ax = axes[1]
    plot_glucose(experiment_open, ax)
    plot_glucose(experiment_zero, ax)
    plot_glucose(experiment_invl, ax)

    ax.set_ylabel("Glucose concentration (mmol)")
    ax.legend(
        [
            "EX_glc__D_e lower bound = -1000.0",
            "EX_glc__D_e lower bound = 0.0",
            "EX_glc__D_e lower bound = 1.0",
        ]
    )
    plt.show()

if __name__ == "__main__":
    main()

[dukovski]

Thank you for noticing this. Let me take a look and I'll get back to you.

Thanks, Ilija

[dukovski]

The way the exchange reactions are done, is that it checks the lowest of the absolute values of the lower bounds set by the user explicitly and the one calculated according to a model (e.g. Michaelis-Menten Monod curve), and sets it as a negative value exchange lower bound. This is probably something we should change in the code. For now, you can fix the bounds of the reaction upstream that produces glc_e, GLCpts. pep_c + glc__D_e → g6p_c + pyr_c. (Pay attention to the direction of the reaction)

Let me know if this is helpful, otherwise we can work on modifying the code and include this feature.

[behzadk]

Hi Ilija, thanks for the quick reply. Ok, I understood. Ideally I would be working directly with the exchange reactions but I should be able to go with the fix you suggested.

Can you point me to where in the code the exchange reactions are done? I can try and make a work around

[dukovski]

The code is the the class FBACell. Lines 762-769.

[dukovski]

// Start of modified code corrected lb 9/19/13 Ilija D. updated by DJORDJE if(media[j]/(cParams.getTimeStep()biomass[i])<calcMichaelisMentenRate(media[j]/(cParams.getSpaceVolume()), km, vMax, hill)) { rates[j] = Math.min(Math.abs(lb[i][j]),Math.abs(media[j]/(cParams.getTimeStep()biomass[i]))); } else rates[j] = Math.min(Math.abs(lb[i][j]), Math.abs(calcMichaelisMentenRate(media[j]/(cParams.getSpaceVolume()), km, vMax, hill)));

[dukovski]

This is for Monod exchange, the other ones are below that. You could actually contribute the code when done.

[behzadk]

Ok great, I'll do that. How do I compile from the src that's on GitHub? I can't find the instructions for this and would like to test the fix before I submit it

edit: Want I would like to implement was a linear rate in the case that production is forced (lb > 0). We can't have a michaelis-menten dynamic for production since we don't model intracellular concentrations.

So something like this in each of case MONOD, case PSEUDO_MONOD, default would work.

if (lb[i][j] > 0.0)
{
    rates[j] = -lb[i][j];
    continue;
}

Note we use -lb here because the rate is multiplied by -1 when writing the lower bounds later on

[dukovski]

I use Eclipse for coding and debugging. https://www.eclipse.org/downloads/packages/release/oxygen/3a/eclipse-ide-java-developers

I'm not sue how familiar you are with Eclipse and IDEs. I also use the egit plugin in Eclipse to have git integrated with it. From Eclispe you can get comets source code with egit. Let me know if you need more instructions, I'm not sure how familiar you are with these tools.

[dukovski]

Yeah we have to be careful. This would work for your case, but it would be nice to have it written in a way that would satisfy any general case.

[behzadk]

Ok thanks, trying to build with Eclipse now. I haven't done much Java so I am a bit out of my depth!

I've got the dependencies and I think I have them added to added to the build libraries. I'm getting one error stopping the build:

Have you seen this before? The file aidaref.jar is causing the error, when I open the file it is empty. Below is the first section of my Java build path, aidaref.jar is there. Any ideas?

[sulheim]

Hey, my workaround for this previously was to set secretion of a particular metabolote as the first objective (and then biomass as the second objective). You then set an upper bound for secretion rather than a lower bound.

However, I do think that it should be possible to enforce this constraint in tje way you have tried here. Mayne this can be a separate parameter, or just have that MM-calculated bounds are ignored if the lower bound is > 0

[dukovski]

Hey, my workaround for this previously was to set secretion of a particular metabolote as the first objective (and then biomass as the second objective). You then set an upper bound for secretion rather than a lower bound.

However, I do think that it should be possible to enforce this constraint in tje way you have tried here. Mayne this can be a separate parameter, or just have that MM-calculated bounds are ignored if the lower bound is > 0

Good point, but yes, we should change the code, and I'm happy to have behzadk getting involved into coding.

[dukovski]

I've got the dependencies and I think I have them added to added to the build libraries. I'm getting one error stopping the build:

You are missing a colt library. I think it is best if you install colt from scratch again. BTW, here's a screenshot of my build path in Eclipse.

[dukovski]

behzadk, BTW would you mind telling us who you are and your affiliation by emailing to comets@bu.edu. If you contribute with code we will credit you for that. Thanks.

[behzadk]

Thanks Ilija, I sent you an email yesterday.

Managed to get the build working just now and have been able to run it with cometspy.

Can you tell me the export settings you use to produce the executable that you have in the COMETS releases?

[dukovski]

So I use these settings when I debug:

[dukovski]

[dukovski]

When I export I just do the base .jar file, no libraries. I link them later when I install it. BTW, for some reason I did not get the email. Can you please send it again at dukovski@bu.edu. Thanks.

[behzadk]

For reference here is the commit that addresses this in the main COMETS repo:

https://github.com/segrelab/comets/commit/cd96c97d71bb22c77971ee83fb0cf46ffd480e37