Bio Systems Biology Context Specific Models

Version Compatibility

Reference examples tested with: COBRApy 0.29+, corda 0.5+, numpy 1.26+, pandas 2.2+, Python 3.10+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Note: COBRApy core does NOT implement GIMME/iMAT/INIT. Real Python implementations live in troppo (multi-method) and corda (CORDA); the reference multi-method implementations are the MATLAB COBRA Toolbox createTissueSpecificModel and RAVEN (INIT/tINIT). Do not expect a cobra.flux_analysis.gimme(); it does not exist.

Context-Specific Models

"Build a liver-specific metabolic model from my expression data" -> Prune/constrain a generic genome-scale model to the reactions an extraction algorithm judges active in that context, given omics data mapped through GPR rules and a threshold.

• Python: corda.CORDA (CORDA), troppo (GIMME/iMAT/tINIT/FASTCORE/CORDA); MATLAB: COBRA Toolbox createTissueSpecificModel, RAVEN (COBRApy for downstream FBA)

The governing principle: the threshold and method are the experiment, not the data

Two facts govern every context-specific model:

• Expression is NOT flux. A highly transcribed gene need not carry high flux (post-transcriptional regulation, allostery, kinetics decouple mRNA from enzyme activity), and an absent transcript does not prove its reaction is off. Absence is a moderately strong constraint; presence is a weak one. mRNA -> enzyme -> flux is a lossy chain, and every extraction method encodes a DIFFERENT guess about that chain.
• Method choice and thresholding dominate the result more than the biology does. Systematic evaluations found that no extraction method reliably beats the others, and that expression-integration methods often fail to beat parsimonious FBA, which uses NO expression at all; the choice of method x threshold x objective changes the model content MORE than the input data does (Machado & Herrgard 2014; Opdam 2017; Richelle 2019). The single on/off threshold is the highest-leverage hidden decision. Consequence: report the method, the threshold strategy, and the objective as first-class methods, and treat a reaction whose inclusion flips with a plausible threshold change as a hypothesis, not a finding.

A corollary trap: GIMME-family methods require a protected objective (usually biomass). For a differentiated, non-proliferating tissue (hepatocyte, neuron), forcing a growth objective is a category error - those cells are not making copies of themselves. Use an objective-free method (iMAT) or a task-based one (tINIT) for non-growing tissue, or define a genuine maintenance/functional task instead of biomass.

Decision: which extraction method

Honest tooling reality: the most complete, best-validated implementations are MATLAB (COBRA Toolbox / RAVEN). In Python, troppo is the multi-method option and corda is the most turnkey native implementation. Steering a user to "just use COBRApy" for iMAT/GIMME sends them into reimplementing an algorithm.

Map Expression Through GPR Rules

Goal: Convert per-gene expression into a per-reaction activity score that respects enzyme logic.

Approach: Evaluate the GPR with min for AND (a complex is limited by its scarcest subunit) and max for OR (any isozyme suffices). This min/max convention is standard but lossy - it discards the quantitative contribution of all but the limiting/dominant gene.

def reaction_activity(rxn, gene_expr, default=0.0):
    '''Aggregate gene expression to a reaction score: min over AND (complex), max over OR (isozyme).'''
    if not rxn.genes:
        return default
    values = [gene_expr.get(g.id, default) for g in rxn.genes]
    return max(values)   # simplified OR; a full parser applies min within each AND-clause first

CORDA in Python (a real, turnkey extraction method)

Goal: Reconstruct a context-specific model that keeps as many high-confidence reactions as possible while excluding absent ones, rescuing reactions that high-confidence ones depend on.

Approach: Translate expression into CORDA's five confidence classes (-1 absent, 0 unknown, 1 low, 2 medium, 3 high) via the GPR, then let CORDA build a "concise but not minimal" model.

from corda import CORDA, reaction_confidence

# gene_conf maps gene id -> confidence in {-1, 0, 1, 2, 3}; derive it from expression quantiles.
gene_conf = {g.id: 2 for g in model.genes}
rxn_conf = {r.id: reaction_confidence(r, gene_conf) for r in model.reactions}   # pass the Reaction, not its GPR string

opt = CORDA(model, rxn_conf)
opt.build()
context_model = opt.cobra_model('liver')   # verify the exact accessor for the installed corda version