Using CoCoMiCo

CoCoMiCo reasons about cooperation and competition emerging from metabolic exchange in microbial communities. A community is a collection of genome-scale metabolic models that define, for each organism in the community, a set of biochemical reactions that produce and consume biomolecules.

Communities are conveniently built from SBML files.

from cocomico.community import Community
from pathlib import Path

# From a directory of SBML files
c0 = Community(sbml_dir=Path("tests/data_test/sbml"))

# From a single file in a directory of SBML files
c1o1 = Community(models=[Path("Com1Org1.sbml")], sbml_dir=Path("tests/data_test/sbml"))

# From a subset of files in a directory of SBML files
sbml_files = [Path(f"Com2Org{ i+1 }.sbml") for i in range(4)]
c2 = Community(models=sbml_files, sbml_dir=Path("tests/data_test/sbml"))

CoCoMiCo makes no taxonomic assumption about species, strains, or type of organism. CoCoMiCo presumes that every metabolic model it is provided defines a separate taxonomic unit, identified by a taxon.

c2.taxa

{'Com2Org1', 'Com2Org2', 'Com2Org3', 'Com2Org4'}

Strains of the same organism need only have different taxon identifiers but may otherwise contain the same reactions on the same biomolecules. In the examples, Com1Org2 and Com2Org2 have the same network.

# Community with two taxa for Org2
c21 = Community(models=[Path("Com1Org2.sbml"), *sbml_files], sbml_dir=Path("tests/data_test/sbml"))
c21.taxa

{'Com1Org2', 'Com2Org1', 'Com2Org2', 'Com2Org3', 'Com2Org4'}

Metabolites are biomolecules with provenance

CoCoMiCo models the exchange of biomolecules across taxa and a key element of that reasoning is the provenance of a biomolecule, the means by which it is produced. In CoCoMiCo an exchanged biomolecule is a metabolite, formally defined as a biomolecule with provenance.

By inspecting the metabolite-reaction relations in a community, we can retrieve the set of metabolites that could be produced by reactions in the community, and the set of biomolecules that could serve as reactants. For produced metabolites, we can retrieve their taxa and biomolecules.

c1o1.products

{Metabolite(provenance="Com1Org1", biomolecule="M_C_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_D_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_F2_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_G_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_H_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_I_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_J_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_X_c")}

c1o1.products.biomolecules

{'M_C_c', 'M_D_c', 'M_F2_c', 'M_G_c', 'M_H_c', 'M_I_c', 'M_J_c', 'M_X_c'}

c1o1.reactants

{Metabolite(provenance="Com1Org1", biomolecule="M_B_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_C_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_D_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_F2_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_G_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_H_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_I_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_K_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_X1_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_X_c")}

For cooperation and competition it is only necessary to explicitly represent the last taxon in which a metabolite was produced; in what follows, provenance means this immediate taxonomic provenance. The complete provenance, comprising the record of the taxon, the reaction that produced the biomolecule, and recursively the complete provenance of the reactants, is not materialized but used implicitly in the reasoning.

We can retrieve the taxa of a set of metabolites

c2.products.taxa

{'Com2Org1', 'Com2Org2', 'Com2Org3', 'Com2Org4'}

We can project sets of metabolites on selected taxa

from cocomico.base import Taxon
c0.products.select([Taxon("Com1Org1"), Taxon("Com1Org2")])

{Metabolite(provenance="Com1Org2", biomolecule="M_A_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_C_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_D_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_F2_c"),
 Metabolite(provenance="Com1Org2", biomolecule="M_F_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_G_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_H_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_I_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_J_c"),
 Metabolite(provenance="Com1Org2", biomolecule="M_M_c"),
 Metabolite(provenance="Com1Org2", biomolecule="M_O_c"),
 Metabolite(provenance="Com1Org2", biomolecule="M_P_c"),
 Metabolite(provenance="Com1Org2", biomolecule="M_Q_c"),
 Metabolite(provenance="Com1Org2", biomolecule="M_S_c"),
 Metabolite(provenance="Com1Org1", biomolecule="M_X_c"),
 Metabolite(provenance="Com1Org2", biomolecule="M_Y_c")}

Metabolic scope determined from seeds

The nutritional environment of a community provides starting point of all metabolic exchanges and determines which reactions may take place.

from cocomico.base import Seeds, Biomolecule

seeds = Seeds([Biomolecule(f"M_{m_}_c") for m_ in ["E", "F", "X1"]])

The metabolic potential, or scope, of a community is evaluated with respect to sets of seeds.

c2.scope(seeds)

{Metabolite(provenance="Com2Org3", biomolecule="M_A2_c"),
 Metabolite(provenance="Com2Org4", biomolecule="M_A3_c"),
 Metabolite(provenance="Com2Org2", biomolecule="M_A_c"),
 Metabolite(provenance="Com2Org3", biomolecule="M_B2_c"),
 Metabolite(provenance="Com2Org4", biomolecule="M_B3_c"),
 Metabolite(provenance="Com2Org3", biomolecule="M_B_c"),
 Metabolite(provenance="Com2Org3", biomolecule="M_C2_c"),
 Metabolite(provenance="Com2Org1", biomolecule="M_C_c"),
 Metabolite(provenance="Com2Org3", biomolecule="M_C_c"),
 Metabolite(provenance="Com2Org4", biomolecule="M_D3_c"),
 Metabolite(provenance="Com2Org1", biomolecule="M_D_c"),
 Metabolite(provenance="Com2Org4", biomolecule="M_E3_c"),
 Metabolite(provenance="Com2Org1", biomolecule="M_F2_c"),
 Metabolite(provenance="Com2Org4", biomolecule="M_F3_c"),
 Metabolite(provenance="Com2Org2", biomolecule="M_F_c"),
 Metabolite(provenance="Com2Org4", biomolecule="M_G3_c"),
 Metabolite(provenance="Com2Org1", biomolecule="M_G_c"),
 Metabolite(provenance="Com2Org1", biomolecule="M_H_c"),
 Metabolite(provenance="Com2Org1", biomolecule="M_I_c"),
 Metabolite(provenance="Com2Org1", biomolecule="M_J_c"),
 Metabolite(provenance="Com2Org4", biomolecule="M_K_c"),
 Metabolite(provenance="Com2Org2", biomolecule="M_M_c"),
 Metabolite(provenance="Com2Org4", biomolecule="M_N_c"),
 Metabolite(provenance="Com2Org2", biomolecule="M_O_c"),
 Metabolite(provenance="Com2Org2", biomolecule="M_P_c"),
 Metabolite(provenance="Com2Org4", biomolecule="M_P_c"),
 Metabolite(provenance="Com2Org2", biomolecule="M_Q_c"),
 Metabolite(provenance="Com2Org2", biomolecule="M_S_c"),
 Metabolite(provenance="Com2Org3", biomolecule="M_V_c"),
 Metabolite(provenance="Com2Org3", biomolecule="M_W_c"),
 Metabolite(provenance="Com2Org1", biomolecule="M_X_c"),
 Metabolite(provenance="Com2Org2", biomolecule="M_Y_c"),
 Metabolite(provenance="Com2Org4", biomolecule="M_Z_c")}

The metabolite set defining the scope can be inspected and projected

c2.scope(seeds).taxa

{'Com2Org1', 'Com2Org2', 'Com2Org3', 'Com2Org4'}

c2.scope(seeds).select([Taxon("Com2Org1"), Taxon("Com2Org2")]).biomolecules

{'M_A_c',
 'M_C_c',
 'M_D_c',
 'M_F2_c',
 'M_F_c',
 'M_G_c',
 'M_H_c',
 'M_I_c',
 'M_J_c',
 'M_M_c',
 'M_O_c',
 'M_P_c',
 'M_Q_c',
 'M_S_c',
 'M_X_c',
 'M_Y_c'}

The metabolic capability of a taxon, without any contributions from the community, can be extracted by specifying an optional choice to scope:

c2.scope(seeds, choice=Taxon("Com2Org2")).biomolecules

{'M_A_c', 'M_F_c', 'M_M_c', 'M_O_c', 'M_P_c', 'M_Q_c', 'M_S_c'}

Note that the scope will be empty if the community members do not have the nutrients that they need.

no_growth = Seeds({Biomolecule("M_unneeded_c"), Biomolecule("M_another_c")})
c21.scope(no_growth)

MetaboliteSet()

Metabolic exchanges

A metabolite is exchanged between a producer taxon $P$ and a consumer taxon $C$, with respect to seeds $S$, if the producer could synthesize the biomolecule but the consumer could not. Formally, the metabolite is in the $S$ scope of the community with provenance $P$, but not producable from $S$ by $C$ alone.

The exchanges are a map from biomolecules to pairs of producer and consumer taxa.

ex = c21.exchange(seeds)
for b,x in sorted(ex.items()):
    print(f"{ b }: { ', '.join(str(i) for i in x) }")

M_A_c: Com1Org2--Com2Org3, Com2Org2--Com2Org3
M_B_c: Com2Org3--Com2Org1
M_C_c: Com2Org3--Com1Org2, Com2Org1--Com1Org2, Com2Org1--Com2Org3, Com2Org3--Com2Org2, Com2Org1--Com2Org2
M_K_c: Com2Org4--Com2Org1
M_N_c: Com2Org4--Com1Org2, Com2Org4--Com2Org2
M_P_c: Com2Org2--Com2Org4, Com1Org2--Com2Org4
M_Y_c: Com1Org2--Com2Org3, Com2Org2--Com2Org3
M_Z_c: Com2Org4--Com2Org3

Monopsonist, polyopsonist metabolites

A metabolite is produced in polyopsony if there are at least two consumer taxons. Polyopsonist metabolites are those for which taxa may be in competition. The cardinality of each metabolite’s competition is reported in the resulting Dict.

c21.polyopsonist(seeds)

{'M_C_c': 3, 'M_F_c': 2, 'M_X1_c': 2, 'M_N_c': 2}

The polyopsonist biomolecules can be queried for membership:

Biomolecule('M_X1_c') in c21.polyopsonist(seeds)

True

The degree of competition can be queried as well:

c21.polyopsonist(seeds)[Biomolecule('M_C_c')]

3

A metabolite is produced in monopsony if there is exactly one consumer taxon, who dictates how the biomolecule is consumed. A monopsonist always has exactly one exchange.

c21.monopsonist(seeds)

{'M_B_c': 1, 'M_A_c': 1, 'M_K_c': 1, 'M_P_c': 1, 'M_Z_c': 1, 'M_Y_c': 1}

Evaluating cooperation and competition potential

CoCoMiCo computes numerical measures of competition and cooperation potential to aid in ranking communities. The competition potential of a community is calculated as the ratio of polyopsonist metabolites to taxa.

import cocomico.score as score
score.competition(c21, seeds)

(1.8, {'competition': 1.8, 'number of polyopsonist metabolites': 4})

The cooperation potential of a community is calculated from the metabolic exchanges. It is the sum, over all exchanged metabolites, of a bonus determined for each metabolite from the number of its producer and consumer taxa.

score.cooperation(c21, seeds)

(19.25,
 {'cooperation': 19.25,
  'number of exchanged metabolites': 8,
  'coop bonus producers': 10.0,
  'coop bonus consumers': 9.25})

Evaluating added value

The added value of a community are the metabolic capabilities that arise from the interaction of its members, beyond the aggregated capabilities of the individual members. CoCoMiCo evaluates this added value using two numerical measures: delta is the difference between the number of metabolites in the community scope and the sum of the individual scopes, rho is the difference between the number of activated×taxon tuples in the community scope and the sum of the individual activated×taxon tuples.

import cocomico.score as score
score.delta(c21, seeds)

(20,
 {'added value community': 20,
  'sum community scope': 41,
  'sum individual scope': 21})

score.rho(c21, seeds)

(37,
 {'added activated community': 37,
  'sum community activated': 68,
  'sum individual activated': 31})

Identifying activated reactions

The metabolic scope of a community defines a set of metabolites, all of which may be available for consumption of community members. When all of the reactants of a reaction are available, we say that the reaction is activated. CoCoMiCo predicts the set of activated reactions in the community:

c2.activated(seeds)

{Reaction(name='R_A2_B2', taxon='Com2Org3'),
 Reaction(name='R_A2_B2rev', taxon='Com2Org3'),
 Reaction(name='R_A3_P', taxon='Com2Org4'),
 Reaction(name='R_A3_Prev', taxon='Com2Org4'),
 Reaction(name='R_AZ_W', taxon='Com2Org3'),
 Reaction(name='R_A_S', taxon='Com2Org2'),
 Reaction(name='R_A_Srev', taxon='Com2Org2'),
 Reaction(name='R_B3_K', taxon='Com2Org4'),
 Reaction(name='R_B3_Krev', taxon='Com2Org4'),
 Reaction(name='R_B_JD', taxon='Com2Org1'),
 Reaction(name='R_C2_B2', taxon='Com2Org3'),
 Reaction(name='R_C2_B2rev', taxon='Com2Org3'),
 Reaction(name='R_C_C2', taxon='Com2Org3'),
 Reaction(name='R_C_C2rev', taxon='Com2Org3'),
 Reaction(name='R_C_PY', taxon='Com2Org2'),
 Reaction(name='R_D3_E3', taxon='Com2Org4'),
 Reaction(name='R_D3_E3rev', taxon='Com2Org4'),
 Reaction(name='R_D_X', taxon='Com2Org1'),
 Reaction(name='R_D_Xrev', taxon='Com2Org1'),
 Reaction(name='R_E3_F3', taxon='Com2Org4'),
 Reaction(name='R_E3_F3rev', taxon='Com2Org4'),
 Reaction(name='R_E_V', taxon='Com2Org3'),
 Reaction(name='R_F2_C', taxon='Com2Org1'),
 Reaction(name='R_F2_Crev', taxon='Com2Org1'),
 Reaction(name='R_F_M', taxon='Com2Org2'),
 Reaction(name='R_F_Mrev', taxon='Com2Org2'),
 Reaction(name='R_G3_A3', taxon='Com2Org4'),
 Reaction(name='R_G3_A3rev', taxon='Com2Org4'),
 Reaction(name='R_G3_F3', taxon='Com2Org4'),
 Reaction(name='R_G_H', taxon='Com2Org1'),
 Reaction(name='R_G_Hrev', taxon='Com2Org1'),
 Reaction(name='R_H_C', taxon='Com2Org1'),
 Reaction(name='R_I_H', taxon='Com2Org1'),
 Reaction(name='R_I_Hrev', taxon='Com2Org1'),
 Reaction(name='R_K_H', taxon='Com2Org1'),
 Reaction(name='R_M_O', taxon='Com2Org2'),
 Reaction(name='R_M_Q', taxon='Com2Org2'),
 Reaction(name='R_M_Qrev', taxon='Com2Org2'),
 Reaction(name='R_N_D3', taxon='Com2Org4'),
 Reaction(name='R_N_D3rev', taxon='Com2Org4'),
 Reaction(name='R_N_M', taxon='Com2Org2'),
 Reaction(name='R_O_P', taxon='Com2Org2'),
 Reaction(name='R_O_Prev', taxon='Com2Org2'),
 Reaction(name='R_Q_A', taxon='Com2Org2'),
 Reaction(name='R_Q_Arev', taxon='Com2Org2'),
 Reaction(name='R_V_B', taxon='Com2Org3'),
 Reaction(name='R_V_Brev', taxon='Com2Org3'),
 Reaction(name='R_W_A2', taxon='Com2Org3'),
 Reaction(name='R_W_A2rev', taxon='Com2Org3'),
 Reaction(name='R_W_B', taxon='Com2Org3'),
 Reaction(name='R_X1_F2', taxon='Com2Org1'),
 Reaction(name='R_X1_ZB3', taxon='Com2Org4'),
 Reaction(name='R_X_G', taxon='Com2Org1'),
 Reaction(name='R_X_Grev', taxon='Com2Org1'),
 Reaction(name='R_Y_W', taxon='Com2Org3')}

As with scope, we can identify the reactions that would activated in the metabolic network of a taxon, without any contributions from the community:

c21.activated(seeds, choice=Taxon("Com2Org2"))

{Reaction(name='R_A_S', taxon='Com2Org2'),
 Reaction(name='R_A_Srev', taxon='Com2Org2'),
 Reaction(name='R_F_M', taxon='Com2Org2'),
 Reaction(name='R_F_Mrev', taxon='Com2Org2'),
 Reaction(name='R_M_O', taxon='Com2Org2'),
 Reaction(name='R_M_Q', taxon='Com2Org2'),
 Reaction(name='R_M_Qrev', taxon='Com2Org2'),
 Reaction(name='R_O_P', taxon='Com2Org2'),
 Reaction(name='R_O_Prev', taxon='Com2Org2'),
 Reaction(name='R_Q_A', taxon='Com2Org2'),
 Reaction(name='R_Q_Arev', taxon='Com2Org2')}

Seed consumption and production

Seed biomolecules are provided in the environment, but not all seeds will be consumed by community members. CoCoMiCo can identify which seeds are in fact consumed by activated reactions.

c21.consumed_seeds(seeds)

{'M_E_c', 'M_F_c', 'M_X1_c'}

Seed biomolecules are not not necessarily exhausted from the environment, because community members may also produce them. CoCoMiCo can identify seeds that are produced by activated reactions.

c21.produced_seeds(seeds)

{'M_F_c'}

As with scope and activated, identifying which seeds are produced and consumed can be restricted to a single taxon. Note that the set of produced or consumed keys could be empty.

[c21.consumed_seeds(seeds, choice=Taxon("Com2Org1")), c21.produced_seeds(seeds, choice=Taxon("Com2Org1"))]

[{'M_X1_c'}, set()]

unneeded = Seeds({Biomolecule("M_unneeded_c"), Biomolecule("M_another_c")})
c21.consumed_seeds(no_growth)

set()

c21.consumed_seeds(seeds.union(unneeded))

{'M_E_c', 'M_F_c', 'M_X1_c'}

CoCoMiCo internals

Constructing models

Genome-scale metabolic models are collections of reactions, defined as bipartite relations between metabolites and reactions. Models in CoCoMiCo are usually created from SBML files.

from cocomico.base import Taxon, Biomolecule, Metabolite, Reaction

T = [ Taxon(t) for t in ["Org1", "Org2", "Org3", "Org4"] ]
B = [ Biomolecule(b) for b in ["A", "B", "C", "D"] ]
M = [ Metabolite(biomolecule=b, provenance=t) for t in T for b in B]
R = { str(x := Reaction(r, t)): x for r in ["R1", "R2", "R2rev", "EX"] for t in T }

from cocomico.model import Model

g1 = Model(
    biomolecule=M,
    tuples={
        (M[0], R["R1.Org1"]), (M[1], R["R1.Org1"]), (R["R1.Org1"], M[2]),
        (M[2], R["R2.Org1"]), (R["R2.Org1"], M[3]),
        (M[3], R["R2rev.Org1"]), (R["R2rev.Org1"], M[2]),
        (R["EX.Org1"], M[0])
    },
)

g1.reactions

{Reaction(name='EX', taxon='Org1'),
 Reaction(name='R1', taxon='Org1'),
 Reaction(name='R2', taxon='Org1'),
 Reaction(name='R2rev', taxon='Org1')}

Internally, models are encoded as reactant-reaction and product-reaction relations

g1.relations

{'reactant': {(Metabolite(provenance="Org1", biomolecule="A"),
   Reaction(name='R1', taxon='Org1'),
   'Org1'),
  (Metabolite(provenance="Org1", biomolecule="B"),
   Reaction(name='R1', taxon='Org1'),
   'Org1'),
  (Metabolite(provenance="Org1", biomolecule="C"),
   Reaction(name='R2', taxon='Org1'),
   'Org1'),
  (Metabolite(provenance="Org1", biomolecule="D"),
   Reaction(name='R2rev', taxon='Org1'),
   'Org1')},
 'product': {(Metabolite(provenance="Org1", biomolecule="A"),
   Reaction(name='EX', taxon='Org1'),
   'Org1'),
  (Metabolite(provenance="Org1", biomolecule="C"),
   Reaction(name='R1', taxon='Org1'),
   'Org1'),
  (Metabolite(provenance="Org1", biomolecule="C"),
   Reaction(name='R2rev', taxon='Org1'),
   'Org1'),
  (Metabolite(provenance="Org1", biomolecule="D"),
   Reaction(name='R2', taxon='Org1'),
   'Org1')}}

Some modeling tools generate reactions with no products or reactants, for example to represent exchanges with the environment. CoCoMiCo will interpret these reactions implicitly as seeds, in addition to those that are explicitly provided. To remove reactions of this kind, use the remove_seed_reactions method on a Model or a Community.

g1.remove_seed_reactions()
g1.reactions

{Reaction(name='R1', taxon='Org1'),
 Reaction(name='R2', taxon='Org1'),
 Reaction(name='R2rev', taxon='Org1')}

Constructing communities

A community is a collection of merged models. To create a community explicitly as a merge of models, provide a dictionary mapping Taxon to Model objects:

models = {
    taxon: Model(relations=model.relations, biomolecule=model.biomolecule)
    for taxon, model in c2.models.items()
}

merged = Community(merge=models)

c2.taxa == merged.taxa
c2.biomolecule == merged.biomolecule
c2.products == merged.products
c2.reactants == merged.reactants
c2.reactions == merged.reactions

True

Knowledge base

Biomolecules, seeds, models, and communities are encoded in a knowledge base that forms the factual basis for CoCoMiCo’s reasoning. The knowledge base for a community is a set of Atoms, that can be inspected using Community.lp_facts.

{atom.predicate for atom in c1o1.knowledge_base}

{'biomolecule', 'product', 'reactant', 'reaction', 'taxon'}

{atom.arg(0) for atom in c1o1.knowledge_base if atom.predicate == 'reaction'}

{'"R_B_JD"',
 '"R_D_X"',
 '"R_D_Xrev"',
 '"R_F2_C"',
 '"R_F2_Crev"',
 '"R_G_H"',
 '"R_G_Hrev"',
 '"R_H_C"',
 '"R_I_H"',
 '"R_I_Hrev"',
 '"R_K_H"',
 '"R_X1_F2"',
 '"R_X_G"',
 '"R_X_Grev"'}

{atom for atom in c1o1.knowledge_base if atom.predicate == 'product' and atom.arg(1) in _}

{product("M_C_c","R_F2_C","Com1Org1"),
 product("M_C_c","R_H_C","Com1Org1"),
 product("M_D_c","R_B_JD","Com1Org1"),
 product("M_D_c","R_D_Xrev","Com1Org1"),
 product("M_F2_c","R_F2_Crev","Com1Org1"),
 product("M_F2_c","R_X1_F2","Com1Org1"),
 product("M_G_c","R_G_Hrev","Com1Org1"),
 product("M_G_c","R_X_G","Com1Org1"),
 product("M_H_c","R_G_H","Com1Org1"),
 product("M_H_c","R_I_H","Com1Org1"),
 product("M_H_c","R_K_H","Com1Org1"),
 product("M_I_c","R_I_Hrev","Com1Org1"),
 product("M_J_c","R_B_JD","Com1Org1"),
 product("M_X_c","R_D_X","Com1Org1"),
 product("M_X_c","R_X_Grev","Com1Org1")}