PyCoMo discretized growth

This tutorial show-cases the use of discretized growth for optimization across the full solution space. With this, linearizing the model with fixed abundance or growth rate is no longer necessary. However, this approach ustilizes binary variable, using a MILP to solve. This is faster than non-linear solvers, but slower than LPs.

The implementation is based on an idea by Salvy and Hatzimanikatis:

Salvy, P., Hatzimanikatis, V. The ETFL formulation allows multi-omics integration in thermodynamics-compliant metabolism and expression models. Nat Commun 11, 30 (2020). https://doi.org/10.1038/s41467-019-13818-7

The expected runtime for this notebook is less than 10 minutes.

[1]:

import cobra
import pycomo

2026-03-15 15:39:37,627 - PyCoMo - INFO - Logger initialized.

[2]:

cobra.Configuration().solver="glpk"

[3]:

pycomo.configure_logger(level="info")

Base Model

[4]:

def construct_single_org_choice(is_one=True):
    model_id = 1 if is_one else 2
    other_id = 2 if is_one else 1
    model = cobra.Model()

    model.name = f"Org{model_id}"

    model.add_metabolites([cobra.Metabolite(i) for i in [
        "S",
        f"X{model_id}",
        "Y1",
        "Y2",
        "A",
        "P",
        "BM"]])
    model.add_reactions([cobra.Reaction(i) for i in [
        f"EX_X{model_id}",
        "EX_Y1",
        "EX_Y2",
        f"TP_X{model_id}",
        "TP_Y1",
        "TP_Y2",
        "v4",
        "v5",
        "bio"]])

    for internal_met in [
        "S",
        "A",
        "P",
        "BM"
    ]:
        model.metabolites.get_by_id(internal_met).compartment = "i"

    for external_met in [
        f"X{model_id}",
        "Y1",
        "Y2",
    ]:
        model.metabolites.get_by_id(external_met).compartment = "e"

    for met in model.metabolites:
        met.name = met.id

    for rxn in model.reactions:
        rxn.name = rxn.id

    model.reactions.get_by_id(f"EX_X{model_id}").add_metabolites({f"X{model_id}": -1})
    model.reactions.get_by_id(f"TP_X{model_id}").add_metabolites({f"X{model_id}": -1, "S": 1})

    model.reactions.get_by_id(f"EX_Y{model_id}").add_metabolites({f"Y{model_id}": -1})
    model.reactions.get_by_id(f"EX_Y{other_id}").add_metabolites({f"Y{other_id}": -1})
    model.reactions.get_by_id(f"TP_Y{model_id}").add_metabolites({f"Y{model_id}": -1, "A": 1})
    model.reactions.get_by_id(f"TP_Y{other_id}").add_metabolites({f"Y{other_id}": -1, "P": 1})

    model.reactions.get_by_id("v4").add_metabolites({"S": -1, "BM": 1, "P": 1})
    model.reactions.get_by_id("v5").add_metabolites({"A": -1, "BM": 1})

    model.reactions.get_by_id("bio").add_metabolites({"BM": -1})

    model.reactions.get_by_id(f"EX_X{model_id}").bounds = (-1000.,0.)
    model.reactions.get_by_id(f"TP_X{model_id}").bounds = (0.,1.)

    model.reactions.get_by_id(f"EX_Y{model_id}").bounds = (0.,1000.)
    model.reactions.get_by_id(f"EX_Y{other_id}").bounds = (0.,1000.)
    model.reactions.get_by_id(f"TP_Y{model_id}").bounds = (0.,1.)
    model.reactions.get_by_id(f"TP_Y{other_id}").bounds = (-1.,0.)

    model.reactions.get_by_id("v4").bounds = (0.,1000.)
    model.reactions.get_by_id("v5").bounds = (0.,1000.)

    model.reactions.get_by_id("bio").bounds = (0.,1000.)

    model.objective = "bio"

    return model

[5]:

def create_toy_pycomo_model():
    model1 = construct_single_org_choice(is_one=True)
    model2 = construct_single_org_choice(is_one=False)

    single_org1 = pycomo.SingleOrganismModel(model1, name="org1")
    single_org2 = pycomo.SingleOrganismModel(model2, name="org2")

    com_model = pycomo.CommunityModel([single_org1, single_org2], name="SymmetricToy")

    com_model.medium = {'EX_X1_medium': 1000.0,
     'EX_Y1_medium': 0.,
     'EX_Y2_medium': 0.,
     'EX_X2_medium': 1000.0}
    com_model.apply_medium()
    return com_model

Find maximum growth rate

The maximum growth rate is 2.0 at 50%/50% abundance of the two organisms. It can be found quickly by the discretized growth rate:

[6]:

com_model = create_toy_pycomo_model()

2026-03-15 15:39:37,924 - PyCoMo - INFO - No community model generated yet. Generating now:
2026-03-15 15:39:37,930 - PyCoMo - INFO - Identified biomass reaction from objective: bio
2026-03-15 15:39:37,930 - PyCoMo - INFO - Note: no products in the objective function, adding biomass to it.
2026-03-15 15:39:37,931 - PyCoMo - INFO - Biomass reaction bio is a boundary reaction.
2026-03-15 15:39:37,931 - PyCoMo - INFO - Moving reactant BM out of the boundary compartment into compartment 'bio', to avoid further boundary reactions being added to it.
2026-03-15 15:39:37,956 - PyCoMo - INFO - Identified biomass reaction from objective: bio
2026-03-15 15:39:37,956 - PyCoMo - INFO - Note: no products in the objective function, adding biomass to it.
2026-03-15 15:39:37,957 - PyCoMo - INFO - Biomass reaction bio is a boundary reaction.
2026-03-15 15:39:37,958 - PyCoMo - INFO - Moving reactant BM out of the boundary compartment into compartment 'bio', to avoid further boundary reactions being added to it.
2026-03-15 15:39:37,994 - PyCoMo - WARNING - No annotation overlap found for matching several metabolites (2). Please make sure that the matched metabolites are indeed representing the same substance in all models! The list of metaboliteswithout annotation overlap can be accessed via 'model.no_annotation_overlap'
2026-03-15 15:39:38,007 - PyCoMo - WARNING - There are 18 metabolites without elements in the model. Mass balance checks may be unreliable.
2026-03-15 15:39:38,008 - PyCoMo - INFO - Generated community model.

[7]:

pycomo.helper.discretized_growth.relax_model_linearisation(com_model)
f_vars = pycomo.helper.discretized_growth.init_f_variables(com_model)
pycomo.helper.discretized_growth.discretize_growth_of_model(com_model, f_vars, mu_lb=0., mu_ub=100., n_mu_bins=1000)

[0 <= 1 <= 1, 0 <= 2 <= 1, 0 <= 4 <= 1, 0 <= 8 <= 1, 0 <= 16 <= 1, 0 <= 32 <= 1, 0 <= 64 <= 1, 0 <= 128 <= 1, 0 <= 256 <= 1, 0 <= 512 <= 1]

[8]:

com_model.model.objective = "community_biomass"
sol = com_model.model.optimize("maximize")
sol

[8]:

Optimal solution with objective value 2.000

	fluxes	reduced_costs
org1_TF_X1_org1_e	-0.50	NaN
org1_TF_Y1_org1_e	-0.50	NaN
org1_TF_Y2_org1_e	0.50	NaN
org1_TP_X1	0.50	NaN
org1_TP_Y1	0.50	NaN
...	...	...
SK_org2_bio_ub	4.99	NaN
SK_org2_to_community_biomass_ub	4.99	NaN
f_final	1.00	NaN
abundance_reaction	0.00	NaN
community_biomass	2.00	NaN

55 rows × 2 columns

Find maximum of flux

flux v4 of organism 1 can reach its maximum flux of 1.0 at a growth rate 1.0 - not at maximum growth. See if the optimum is found:

[9]:

com_model.model.objective = "org1_v4"
sol = com_model.model.optimize("maximize")
sol

[9]:

Optimal solution with objective value 1.000

	fluxes	reduced_costs
org1_TF_X1_org1_e	-1.0	NaN
org1_TF_Y1_org1_e	0.0	NaN
org1_TF_Y2_org1_e	1.0	NaN
org1_TP_X1	1.0	NaN
org1_TP_Y1	0.0	NaN
...	...	...
SK_org2_bio_ub	-0.0	NaN
SK_org2_to_community_biomass_ub	-0.0	NaN
f_final	1.0	NaN
abundance_reaction	0.0	NaN
community_biomass	1.0	NaN

55 rows × 2 columns

[ ]: