7. Building a model from scratch#

Add metabolites, a reaction, a gene association and the exchanges that let material in and out, and end with a model that carries flux. This is how a hand-built model starts, and how you test an idea before an automated reconstruction is involved.

Functions on this page#

MATLAB

Python

addMets

Metabolite cobrapy

add metabolites

addRxns

add_reactions_from_equations

add reactions from equation strings

addGenesRaven

auto-created from the GPR cobrapy

add genes

addExchangeRxns

Model.add_boundary cobrapy

add exchange reactions

constructEquations

Reaction.reaction cobrapy

read the equations back

The important part is the same in both: write the reaction as an equation string and let the toolbox derive the stoichiometry, rather than filling in a column of the stoichiometric matrix by hand.

Setup#

empty.xml is RAVEN's starter model: four metabolites (sucrose, glucose, fructose and water) and one reaction, invertase, which splits sucrose into glucose and fructose. Everything is in one compartment, e.

model = importModel('empty.xml');
fprintf('%d reactions, %d metabolites\n', numel(model.rxns), numel(model.mets));
eqn = constructEquations(model);
fprintf('%s\n', eqn{1});

The model contains 0 errors and 1 warnings.

[Warning: The following fields have prefixes removed from all entries. If this is undesired, run importModel with removePrefix as false. Example: importModel('filename.xml',[],false); model.rxns (R_ prefix) model.mets (M_ prefix)]
1 reactions, 4 metabolites
sucrose[e] + H2O[e] => glucose[e] + fructose[e]

from cobra.io import read_sbml_model

model = read_sbml_model("empty.xml")
print(len(model.reactions), "reactions,", len(model.metabolites), "metabolites")
for reaction in model.reactions:
    print(f"  {reaction.id}: {reaction.reaction}")

1 reactions, 4 metabolites
  r1: m1 + m4 --> m2 + m3

7.1 Add metabolites#

We are going to add hexokinase, so the model needs ATP, ADP and glucose-6-phosphate first.

metsToAdd.mets = {'m5', 'm6', 'm7'};
metsToAdd.metNames = {'ATP', 'ADP', 'glucose-6-phosphate'};
metsToAdd.compartments = 'e';
metsToAdd.metFormulas = {'C10H12N5O13P3', 'C10H12N5O10P2', 'C6H11O9P'};
model = addMets(model, metsToAdd);
fprintf('%d metabolites\n', numel(model.mets));

7 metabolites

from cobra import Metabolite

model.add_metabolites([
    Metabolite("m5", name="ATP", formula="C10H12N5O13P3", compartment="e"),
    Metabolite("m6", name="ADP", formula="C10H12N5O10P2", compartment="e"),
    Metabolite("m7", name="glucose-6-phosphate", formula="C6H11O9P", compartment="e"),
])
print(len(model.metabolites), "metabolites")

7 metabolites

Give every metabolite a formula and a compartment as you add it. Both are what 9. Quality control checks against, and adding them later means revisiting every reaction written in between.

7.2 Add a reaction#

The equation string is the same in both toolboxes: <=> for a reversible reaction, => for an irreversible one.

rxnsToAdd.rxns = {'HEX1'};
rxnsToAdd.rxnNames = {'hexokinase'};
rxnsToAdd.equations = {'m2 + m5 => m7 + m6'};
rxnsToAdd.grRules = {'YFR053C'};
model = addRxns(model, rxnsToAdd, 'eqnType', 1, 'allowNewGenes', true);

idx = getIndexes(model, 'HEX1', 'rxns');
eqn = constructEquations(model, {'HEX1'});
fprintf('%s\n', eqn{1});
fprintf('bounds: [%g %g], genes: %s\n', model.lb(idx), model.ub(idx), ...
    model.grRules{idx});

New genes added to the model:
YFR053C
glucose[e] + ATP[e] => ADP[e] + glucose-6-phosphate[e]
bounds: [0 Inf], genes: YFR053C

eqnType says how the equation is written: 1 or "id" matches metabolites against model.mets, 2 or "name" against model.metNames with compartment deciding where they go, and 3 or "name[comp]" reads the compartment from the token itself.

allowNewGenes is needed because YFR053C is not in the model yet. Both allowNewGenes and allowNewMets default to false, so addRxns refuses anything it does not recognise until told otherwise, and a mistyped identifier is an error rather than a new entity. allowNewMets also accepts a string, used as the prefix for the ids it generates. Adding metabolites with addMets first carries more annotation than addRxns can infer.

from raven_toolbox.manipulation import add_reactions_from_equations

add_reactions_from_equations(model, [
    {
        "id": "HEX1",
        "name": "hexokinase",
        "equation": "m2 + m5 => m7 + m6",
        "gene_reaction_rule": "YFR053C",
    },
])

hex1 = model.reactions.get_by_id("HEX1")
print(hex1.reaction)
print("bounds:", hex1.bounds, "genes:", hex1.gene_reaction_rule)

m2 + m5 --> m6 + m7
bounds: (0.0, 1000.0) genes: YFR053C

The arrow sets the bounds, so an irreversible reaction needs no bounds key; passing bounds overrides it. <=> is reversible, and -->, -> and => are all accepted for an irreversible reaction.

allow_new_genes and allow_new_mets both default to True here, the opposite of the MATLAB defaults, so a gene named in the rule is created without a separate step and without a message. That removes the addGenesRaven call MATLAB needs, and it means a typo becomes a new gene rather than an error. Set them to False once the model's metabolites and genes are all defined.

7.3 Add exchanges#

Without boundary reactions nothing can enter or leave, so there is no net uptake, secretion or growth. Internal cycles can still carry flux (see 18.3), and under the steady-state assumption a metabolite that appears in only one reaction blocks that reaction. Exchange reactions are the boundary.

model = addExchangeRxns(model, 'both', 'mets', model.mets);
fprintf('%d reactions\n', numel(model.rxns));

NOTE: The exchange reactions are assigned to the first compartment
9 reactions

for metabolite in list(model.metabolites):
    model.add_boundary(metabolite, type="exchange")

print(len(model.reactions), "reactions")
print(sorted(rxn.id for rxn in model.exchanges))

9 reactions
['EX_m1', 'EX_m2', 'EX_m3', 'EX_m4', 'EX_m5', 'EX_m6', 'EX_m7']

add_boundary works out that e is the external compartment. It refuses when nothing looks external, which is the usual reason it fails on a small hand-built model: give a compartment a recognisable name, or build the exchange as an ordinary reaction with a single metabolite. type="demand" and type="sink" add the other two kinds of boundary reaction, which cobrapy tracks in separate collections.

7.4 Does it carry flux?#

Solving the model shows whether the reactions assembled so far can carry flux at all.

model = setParam(model, 'obj', 'HEX1', 1);
sol = solveLP(model);
fprintf('HEX1 flux: %.2f\n', sol.x(getIndexes(model, 'HEX1', 'rxns')));

HEX1 flux: 1000.00

model.objective = "HEX1"
solution = model.optimize()
print(f"HEX1 flux: {solution.fluxes['HEX1']:.2f}")

HEX1 flux: 1000.00

Warning

What can go wrong

  • A typo creates a metabolite. The defaults differ: add_reactions_from_equations creates unrecognised metabolites and genes, while addRxns refuses them unless allowNewMets and allowNewGenes are set. So m7 and M7 become two different things in Python and an error in MATLAB. Set allow_new_mets=False and allow_new_genes=False once everything is defined.

  • No exchange reactions. Nothing can be taken up or secreted, so growth is zero, with no error to explain why.

  • A model with no external compartment. add_boundary cannot guess one, and RAVEN's addExchangeRxns adds exchanges for internal metabolites without warning, opening the model to uptake and secretion of intermediates.

See also#