Manipulation

MATLAB functions in RAVEN/src/manipulation of the RAVEN toolbox.

Functions

Function	Summary
addExchangeRxns	Add exchange reactions for some metabolites.
addGenesRaven	Add genes to a model.
addMets	Add metabolites to a model.
addRxns	Add reactions to a model.
addRxnsGenesMets	Copy reactions from a source model into another model.
addTransport	Add transport reactions between compartments.
changeGrRules	Change the grRules of multiple reactions.
changeRxns	Modify the equations of reactions in a model.
closeModel	Add boundary metabolites to exchange reactions.
contractModel	Contract a model by merging identical reactions.
convertToIrrev	Convert a model to irreversible form.
copyToComps	Copy reactions to new compartment(s).
deleteUnusedGenes	Delete all genes not associated to any reaction.
expandModel	Expand reactions that use several gene associations.
findDuplicateRxns	Find reactions that share identical stoichiometry.
findPotentialErrors	Identify reactions with non-DNF GPR rules.
generateNewIds	Generate a list of new metabolite or reaction ids.
mergeCompartments	Merge all compartments in a model.
mergeModels	Merge models into one model structure.
permuteModel	Change the order of the reactions or metabolites in a model.
removeBadRxns	Remove reactions that enable production or consumption from nothing.
removeGenes	Delete a set of genes from a model.
removeLowScoreGenes	Remove low-score genes from a model.
removeMets	Delete a set of metabolites from a model.
removeReactions	Delete a set of reactions from a model.
replaceMets	Replace a metabolite with another already in the model.
reverseRxns	Reverse the direction of reactions in a model.
setExchangeBounds	Define exchange flux bounds for a set of metabolites.
setParam	Set parameters for reactions.
simplifyModel	Simplify a model by deleting reactions and metabolites.
sortIdentifiers	Sort model identifiers alphabetically.
sortModel	Sort the metabolites and reactions of a model.
standardizeGrRules	Standardise gene-reaction rules in a model.

Reference

addExchangeRxns

function [model, addedRxns]=addExchangeRxns(model,reactionType,varargin)

Add exchange reactions for some metabolites.

The function builds the stoichiometric columns directly from the metabolite list, without the equation parsing of addRxns, which makes it faster for adding exchange reactions. New reactions are named "metName exchange (TYPE)" and have the ids "EXC_TYPE_METID", where TYPE is IN, OUT or BOTH according to reactionType.

Parameters

model : struct

a model structure.

reactionType : char

the type of reactions to add:

• 'in' : input reactions

• 'out' : output reactions

• 'both' : reversible input and output reactions. The positive direction corresponds to output

Name-Value Arguments

mets : cell or logical or double

a cell array of metabolite ids, a logical vector with one element per metabolite of the model, or a vector of metabolite indexes to add reactions for (default model.mets).

Returns

model : struct

updated model structure.

addedRxns : cell

ids of the added reactions.

Examples

[model, addedRxns] = addExchangeRxns(model, reactionType, mets);

addGenesRaven

function newModel=addGenesRaven(model,genesToAdd)

Add genes to a model.

The function checks neither MIRIAM formats nor forbidden characters in gene ids.

Parameters

model : struct

a model structure.

genesToAdd : struct

the genes to add, which can have the following fields:

• genes : cell array with unique strings that identify each gene. The characters allowed in SBML ids (mainly a-z, 0-9 and "_") are the valid ones, but the function does not check them; they matter only for export to SBML

• geneShortNames : cell array of gene abbreviations (optional, default "")

• geneMiriams : cell array with MIRIAM structures (optional, default [])

• proteins : cell array of the protein names associated with each gene (optional, default "")

Returns

newModel : struct

an updated model structure.

Examples

newModel = addGenesRaven(model, genesToAdd);

addMets

function newModel=addMets(model,metsToAdd,varargin)

Add metabolites to a model.

The function checks neither MIRIAM formats nor forbidden characters in metabolite ids.

Parameters

model : struct

a model structure.

metsToAdd : struct

the metabolite structure, which can have the following fields:

• mets : cell array with unique strings that identify each metabolite (optional, default: ids numbered with the prefix defined below)

• metNames : cell array with the names of each metabolite (optional if mets is given, default: the ids)

• compartments : cell array with the compartment of each metabolite. Each entry must match model.comps. A single string, or a cell array with one entry, assigns all metabolites to that compartment

• b : Nx1 or Nx2 matrix with equality constraints for each metabolite (optional, default 0)

• unconstrained : vector describing if each metabolite is an exchange metabolite (1) or not (0) (optional, default 0)

• inchis : cell array with InChI strings (optional, default "")

• metSmiles : cell array with SMILES strings (optional, default "")

• metFormulas : cell array with the formulas (optional, default "")

• metMiriams : cell array with MIRIAM structures (optional, default [])

• metCharges : metabolite charge (optional, default NaN)

• metDeltaG : Gibbs free energy of formation at biochemical standard condition in kJ/mol (optional, default NaN)

• metNotes : cell array with metabolite notes as strings (optional, default "")

Name-Value Arguments

copyInfo : logical

true to copy the available annotation from the metabolite with the same name in another compartment of the model (default true).

prefix : char

prefix of the generated metabolite ids when metsToAdd.mets is not specified. If ids with the prefix exist in the model, the numbering starts at the highest existing number plus 1 (default "m_").

Returns

newModel : struct

an updated model structure.

Examples

newModel = addMets(model, metsToAdd, copyInfo, prefix);

Notes

When multiple metabolites are added at once, metMiriams is a cell array with one MIRIAM structure per metabolite (example with ChEBI and KEGG):

metsToAdd.metMiriams{1} = struct("name",{{"chebi";"kegg.compound"}},...
    "value",{{"CHEBI:18072";"C11821"}});
metsToAdd.metMiriams{2} = struct("name",{{"chebi";"kegg.compound"}},...
    "value",{{"CHEBI:31132";"C12248"}});

addRxns

function newModel=addRxns(model,rxnsToAdd,varargin)

Add reactions to a model.

The function checks that grRules has one entry per reaction and that every gene of a rule is in model.genes (an error otherwise, unless allowNewGenes is true). standardizeGrRules then writes the rules in the standard format.

New metabolites receive the available annotation of the metabolite with the same name in another compartment, when one exists.

Parameters

model : struct

a model structure.

rxnsToAdd : struct

the reaction structure, which can have the following fields:

• rxns : cell array with unique strings that identify each reaction

• equations : cell array with equation strings. Decimal coefficients are expressed as "1.2". Reversibility is indicated by "<=>" or "=>"

• mets : (alternative to equations) cell array with the metabolites involved in each reaction as nested arrays, for example {{"met1","met2"},{"met1","met3","met4"}}. For a single added reaction, it can be a string array: {"met1","met2"}

• stoichCoeffs : (alternative to equations) cell array with the corresponding stoichiometries as nested vectors, for example {[-1,+2],[-1,-1,+1]}. For a single added reaction, it can be a vector: [-1,+2]

• rxnNames : cell array with the names of each reaction (optional, default "")

• lb : vector with the lower bounds (optional, default model.annotation.defaultLB or -inf for reversible reactions and 0 for irreversible reactions when "equations" is used. When "mets" and "stoichCoeffs" are used, the default for all reactions is model.annotation.defaultLB or -inf)

• ub : vector with the upper bounds (optional, default model.annotation.defaultUB or inf)

• c : vector with the objective function coefficients (optional, default 0)

• eccodes : cell array with the EC-numbers for each reaction. Delimit several EC-numbers with ";" (optional, default "")

• subSystems : cell array with the subsystems for each reaction (optional, default "")

• grRules : cell array with the gene-reaction relationship for each reaction. "(A and B) or (C)" means that the reaction is catalysed by a complex of A and B or by C on its own. All genes must be present in model.genes, unless allowNewGenes is true; addGenesRaven adds genes (optional, default "")

• rxnMiriams : cell array with Miriam structures (optional, default [])

• rxnComps : cell array with compartments (as in model.comps) (optional, default {})

• rxnNotes : cell array with reaction notes (optional, default "")

• rxnDeltaG : Gibbs free energy at biochemical standard condition in kJ/mol (optional, default NaN)

• rxnReferences : cell array with reaction references (optional, default "")

• rxnConfidenceScores : vector with reaction confidence scores (optional, default NaN)

Name-Value Arguments

eqnType : double or char

the interpretation of the equation strings (default 1). The values are the integers 1, 2 and 3, or the string aliases "id", "name" and "name[comp]":

• 1 or "id" : the metabolites are matched to model.mets. New metabolites (if allowed) are added to "compartment"

• 2 or "name" : the metabolites are matched to model.metNames and all metabolites are assigned to "compartment". New metabolites get ids generated by generateNewIds with the prefix "m_" (or the prefix given in allowNewMets), numbered from the highest existing number with that prefix plus 1

• 3 or "name[comp]" : the metabolites are written as "metNames[comps]". Only compartments in model.comps are allowed. New metabolites get ids generated in the same way as for eqnType 2

compartment : char

the compartment of the metabolites for eqnType=2 and of new metabolites for eqnType=1. Must match model.comps. Not used for eqnType=3.

allowNewMets : logical or char

true to allow the function to add new metabolites. A string also allows it and is used as the prefix of the new metabolite ids. The function calls addMets with copyInfo true, so a new metabolite takes the charge and deltaG from a metabolite with the same name in another compartment, and InChI, SMILES, formula and MIRIAM annotation from it when its own are empty. When the annotation of new metabolites is to be set explicitly, or when they are to be exchange (unconstrained) metabolites, adding them with addMets first gives that control (default false).

allowNewGenes : logical

true to allow the function to add new genes (default false).

Returns

newModel : struct

an updated model structure.

Examples

newModel = addRxns(model, rxnsToAdd, eqnType, compartment, ...
                   allowNewMets, allowNewGenes);

addRxnsGenesMets

function model=addRxnsGenesMets(model,sourceModel,rxns,varargin)

Copy reactions from a source model into another model.

Copies reactions from a source model to a draft model, together with their new metabolites and genes. Reactions that are already in the draft model are skipped.

The draft model and the source model must follow the same metabolite and compartment naming convention, because metabolites are matched only by metaboliteName[compartment]. The function can copy additional reactions from a source model to a draft model that getModelFromHomology created from the same models.

Parameters

model : struct

the draft model that receives the reactions.

sourceModel : struct

the model that the reactions and metabolites are copied from.

rxns : cell or char

reaction ids from the source model. A string adds one reaction.

Name-Value Arguments

addGene : logical or char or cell

three options (default false):

• false : no genes are annotated to the new reactions

• true : grRules are copied from the sourceModel and new genes are added when required

• string or cell array : the grRules of the new reactions, as a string or a cell array, and new genes are added when required

rxnNote : cell or char

strings that explain why the reactions were copied, stored in rxnNotes of the new reactions. A string applies the same note to every new reaction (default "Added via addRxnsGenesMets()").

confidence : double

integer confidence score for all new reactions, following doi:10.1038/nprot.2009.203 (default 0):

• 4 : biochemical data: direct evidence from enzyme assays

• 3 : genetic data: knockout/-in or overexpression analysis

• 2 : physiological data: indirect evidence, for example secretion products or defined medium requirement; sequence data: genome annotation

• 1 : modelling data: required for a functional model, hypothetical reaction

• 0 : no evidence

Returns

model : struct

an updated model structure.

Examples

newModel = addRxnsGenesMets(model, sourceModel, rxns, addGene, ...
                            rxnNote, confidence);

addTransport

function [model, addedRxns]=addTransport(model,fromComp,toComps,varargin)

Add transport reactions between compartments.

The function builds the stoichiometric columns directly from the metabolite names, without the equation parsing of addRxns, which makes it faster for adding transport reactions. New reaction names have the format "metaboliteName transport, fromCompName-toCompName", where the compartment names are the entries of model.compNames. New reaction ids are numbered sequentially with the prefix tr_ (tr_0001, tr_0002, ...).

Parameters

model : struct

a model structure.

fromComp : char

the id of the compartment to transport from. Must match one entry of model.comps.

toComps : cell

the ids of the compartments to transport to. Each must match an entry of model.comps.

Name-Value Arguments

metNames : cell

the metabolite names to add transport for (default all metabolites in fromComp).

isRev : logical

true if the transport reactions should be reversible (default true).

onlyToExisting : logical

true to add transport of a metabolite only if the metabolite exists in the target compartment. If false, the missing metabolites are first added with addMets (default true).

prefix : char

prefix of the reaction ids (default "tr_").

Returns

model : struct

updated model structure.

addedRxns : cell

ids of the added reactions.

Examples

[model, addedRxns] = addTransport(model, fromComp, toComps, ...
                     metNames, isRev, onlyToExisting, prefix);

changeGrRules

function model = changeGrRules(model,rxns,grRules,varargin)

Change the grRules of multiple reactions.

Genes in grRules that are not in the model are added to the model.

Parameters

model : struct

a model structure to change the gene association.

rxns : char or cell

reaction IDs.

grRules : char or cell

the replacement or additional gene association. ' and ' indicates subunits, ' or ' indicates isoenzymes, and brackets '()' separate different instances. A single rule applies to all reactions in rxns.

Name-Value Arguments

replace : logical

true to replace the old gene association with the new association. False to combine them as "(old) or (new)" (default true).

Returns

model : struct

an updated model structure.

Examples

model = changeGrRules(model, rxns, grRules, replace);

changeRxns

function model=changeRxns(model,rxns,equations,varargin)

Modify the equations of reactions in a model.

Parameters

model : struct

a model structure.

rxns : cell

cell array with reaction ids.

equations : cell or struct

cell array with equations, or a structure with the fields "mets" and "stoichCoeffs" as in addRxns, for example:

• equations.mets = {{"met1","met2"},{"met1","met3"}}

• equations.stoichCoeffs = {[-1,+2],[-1,+1]}

Name-Value Arguments

eqnType : double

the interpretation of the equation strings (default 1):

• 1 : the metabolites are matched to model.mets. New metabolites (if allowed) are added to "compartment".

• 2 : the metabolites are matched to model.metNames and all metabolites are assigned to "compartment". New metabolites get ids generated by generateNewIds with the prefix "m_", numbered from the highest existing number with that prefix plus 1.

• 3 : the metabolites are written as "metNames[comps]". Only compartments in model.comps are allowed. New metabolites get ids generated in the same way as for eqnType 2.

compartment : char

the compartment of the metabolites for eqnType=2 and of new metabolites for eqnType=1. Must match model.comps. Not used for eqnType=3.

allowNewMets : logical

true to allow the function to add new metabolites. The function calls addMets with copyInfo true, so a new metabolite takes the charge and deltaG from a metabolite with the same name in another compartment, and InChI, SMILES, formula and MIRIAM annotation from it when its own are empty. When the annotation of new metabolites is to be set explicitly, or when they are to be exchange (unconstrained) metabolites, adding them with addMets first gives that control (default false).

Returns

model : struct

an updated model structure.

Examples

model = changeRxns(model, rxns, equations, eqnType, compartment, allowNewMets);

Notes

The function keeps the bounds of the reactions. Changing an irreversible reaction to a reversible one, or the reverse, changes the model.rev field but not model.lb and model.ub, so the bounds still define the reversibility of the reaction. setParam changes the bounds.

New metabolites receive the annotation of the metabolite with the same name in another compartment, when one exists (see allowNewMets).

closeModel

function closedModel=closeModel(model)

Add boundary metabolites to exchange reactions.

Adds a compartment "b" (boundary) and a boundary metabolite in it for every reaction that has metabolites on one side only (exchange, sink and demand reactions). The boundary metabolite takes part in the reaction with the coefficient 1 or -1 and is marked as unconstrained in model.unconstrained; all other metabolites are marked as constrained. A reaction with one metabolite gets the boundary metabolite "_b" with the annotation of that metabolite; a reaction with several metabolites gets "_b" without annotation.

Parameters

model : struct

a model structure.

Returns

closedModel : struct

an updated model structure with boundary metabolites added.

Examples

closedModel = closeModel(model);

contractModel

function [reducedModel, removedRxns, indexedDuplicateRxns]=contractModel(model,varargin)

Contract a model by merging identical reactions.

simplifyModel runs contractModel for its deleteDuplicates option. If the duplicated reactions have "EXP*" suffixes (where * is a digit), the model has passed through expandModel, and the function removes these suffixes.

The gene associations and EC numbers of merged reactions are combined. The merged reaction keeps the widest bounds and the largest objective coefficient, and a warning is given when the duplicates differ in these. findDuplicateRxns lists duplicate reactions without changing the model.

Parameters

model : struct

a model structure.

Name-Value Arguments

distReverse : logical

true to treat reactions with the same metabolites but different reversibility as different reactions (default true).

mets : char or cell

metabolite ids as a string or a cell array of strings; only the reactions that involve these metabolites are checked for duplication (default: all reactions). replaceMets uses this option.

Returns

reducedModel : struct

a model structure without duplicate reactions.

removedRxns : cell

cell array with the ids of the removed duplicate reactions.

indexedDuplicateRxns : cell

cell array with one entry per reaction of reducedModel, holding the ids of the duplicates removed into that reaction (multiple ids separated by semicolon).

Examples

[reducedModel, removedRxns, indexedDuplicateRxns] = contractModel(model, distReverse, mets);

Notes

The function may not work for grRules strings with nested expressions of "and" and "or".

convertToIrrev

function [irrevModel,matchRev,rev2irrev,irrev2rev]=convertToIrrev(model,varargin)

Convert a model to irreversible form.

Reversible reactions are split into one forward and one reverse reaction. The reverse reactions are named 'rxnID_REV'. A warning is shown if some reaction ids already end with '_REV'.

Parameters

model : struct

a model structure.

Name-Value Arguments

rxns : cell

cell array with the ids of the reactions to split, if reversible (default model.rxns).

Returns

irrevModel : struct

a model structure where reversible reactions have been split into one forward and one reverse reaction.

matchRev : double

matching forward reaction to its backward reaction.

rev2irrev : cell

forward and backward reactions for reversible reactions.

irrev2rev : double

matching all reactions back to original model.

Examples

[irrevModel,matchRev,rev2irrev,irrev2rev]=convertToIrrev(model,rxns);

copyToComps

function model=copyToComps(model,toComps,varargin)

Copy reactions to new compartment(s).

Parameters

model : struct

a model structure.

toComps : cell

cell array of compartment ids. An id that is not in model.comps is added as a new compartment (see compNames and compOutside).

Name-Value Arguments

rxns : cell or logical or double

a cell array of reaction ids, a logical vector with one element per reaction of the model, or a vector of reaction indexes to copy (default model.rxns).

deleteOriginal : logical

true to remove the original reactions, which moves the reactions (default false).

compNames : cell

cell array of compartment names. Used when new compartments are added (default toComps).

compOutside : cell

cell array with the id (as in comps) of the compartment that surrounds each new compartment. Used when new compartments are added (default all {''}).

Returns

model : struct

an updated model structure.

Examples

model=copyToComps(model,toComps,rxns,deleteOriginal,compNames,compOutside);

Notes

New reactions and metabolites have the ids "id_toComps(i)".

deleteUnusedGenes

function reducedModel=deleteUnusedGenes(model,varargin)

Delete all genes not associated to any reaction.

Parameters

model : struct

a model structure.

Name-Value Arguments

verbose : double

0 for silent; 1 for printing the number of deleted genes; 2 for printing the list of deleted genes (default 1).

Returns

reducedModel : struct

an updated model structure.

Examples

reducedModel=deleteUnusedGenes(model);

expandModel

function [newModel, rxnToCheck]=expandModel(model)

Expand reactions that use several gene associations.

Each reaction with several gene associations (isozymes) is split into one reaction per association.

Parameters

model : struct

a model structure.

Returns

newModel : struct

model structure with separate reactions for isozymes. The ids of split reactions are id_EXP_1, id_EXP_2 and so on, numbered per isozyme.

rxnToCheck : cell

cell array with the original ids of the reactions whose grRule is not in disjunctive normal form and was expanded by distributivity, and of the reactions whose grRule cannot be parsed and is left unchanged.

Examples

[newModel, rxnToCheck]=expandModel(model);

Notes

The isozymes are the AND-clauses of the grRule in disjunctive normal form, so nested expressions are expanded: "g1 and (g2 or g3)" yields "g1 and g2" and "g1 and g3". Reactions listed in rxnToCheck because of distributivity are expanded in the same way; they are listed because a rule that needs distributivity can indicate that the GPR differs from what its author intended. standardizeGrRules and findPotentialErrors inspect them.

findDuplicateRxns

function pairs = findDuplicateRxns(model, varargin)

Find reactions that share identical stoichiometry.

Only the stoichiometry is compared; bounds, GPRs and annotations are ignored. With the default ignoreDirection, A→B and B→A are duplicates. The pairs identify reactions that could be merged (contractModel merges duplicate reactions).

Parameters

model : struct

a model structure.

Name-Value Arguments

ignoreDirection : logical

true to treat A→B and B→A as duplicates (default true).

Returns

pairs : double

Nx2 numeric array of reaction-index pairs (i, j) where reactions i and j share the same (possibly negated) stoichiometry, with i < j. Empty if the model has no duplicates.

Examples

pairs = findDuplicateRxns(model);
pairs = findDuplicateRxns(model, false);

findPotentialErrors

function issues=findPotentialErrors(model)

Identify reactions with non-DNF GPR rules.

Detects grRules that are not in disjunctive normal form (DNF), that is, rules in which an OR operator is nested inside an AND context. The RAVEN format requires DNF: groups of genes that must all be present together (enzyme complexes) are written as conjunctions, and the reaction is catalysed by at least one complex (a disjunction of conjunctions).

Parameters

model : struct

a model structure with grRules and rxns fields.

Returns

issues : struct

struct array with one element per problematic reaction and the fields:

• index : double, 1-based index into model.rxns.

• rxn : char, reaction id.

• grRule : char, the problematic grRule string.

• reason : char, short explanation.

Notes

Detection walks the parse tree from parseGrRule, so gene ids that contain substrings such as "and" or "or" are not misread as operators.

A rule is non-DNF when an AND operator has an OR anywhere beneath it. Bracketing alone does not make a rule non-DNF: "((G1 and G2) or G3)" and "(G1 or G2) or (G3 and G4)" are both in DNF. standardizeGrRules attempts an automatic repair.

Rules that cannot be parsed at all are reported as issues too, with the parser's message as the reason.

Examples

issues = findPotentialErrors(model);
if ~isempty(issues)
    fprintf("%d reactions with non-DNF grRules\n", numel(issues));
end
indexes2check = vertcat(issues.index);

generateNewIds

function newIds=generateNewIds(model,type,prefix,varargin)

Generate a list of new metabolite or reaction ids.

The ids are numbered sequentially with a defined prefix. The highest existing number among the ids of that type with the prefix is read from the model. Suffixes are read with str2double. A suffix counts as a number when the result is a finite integer value below 2^52; text, Inf, NaN, fractions (r_2.5) and numbers of 2^52 or more (r_1e300, r_1e400) do not count.

Parameters

model : struct

model structure.

type : char

type of id, 'rxns' or 'mets'.

prefix : char

prefix to be used in all ids, for example 's_' or 'r_'.

Name-Value Arguments

quantity : double

number of new ids to generate (default 1).

numLength : double

length of the numerical part of the id; 4 gives ids like r_0001 and 6 gives ids like r_000001. If the prefix is used in the model, the length of the existing ids applies instead (default 4).

Returns

newIds : cell

cell array with the generated ids.

Examples

newIds = generateNewIds(model, type, prefix, quantity, numLength);

mergeCompartments

function [model, deletedRxns, duplicateRxns]=mergeCompartments(model,varargin)

Merge all compartments in a model.

Parameters

model : struct

a model structure.

Name-Value Arguments

keepUnconstrained : logical

true to keep the unconstrained metabolites in a separate compartment "b" (name "Unconstrained"). If they are merged, an exchange reaction between an unconstrained metabolite and a metabolite of the same name has no net stoichiometry after merging and is deleted (default false).

deleteRxnsWithOneMet : logical

true to delete reactions that have one metabolite after merging and had more before. Such reactions arise from reactions such as A[c] + B[c] => A[m]. Some models balance hydrogen around each membrane with reactions like this (default false).

distReverse : logical

true to treat reactions with the same metabolites but different reversibility as different reactions (default true).

Returns

model : struct

a model with all reactions in one compartment.

deletedRxns : cell

the ids of the reactions deleted because they have one metabolite after merging.

duplicateRxns : cell

one entry per reaction of the returned model: '' when no duplicate was merged into it, otherwise the ids of the identical reactions from different compartments that were removed into it, joined by ";".

Examples

[model, deletedRxns, duplicateRxns] = mergeCompartments(model);

Notes

Merges all compartments into one compartment "s" (name "System"). The merged model can be used to test whether the network can synthesise all metabolites.

Metabolites with the same name are merged into one. The merged metabolite keeps the id and the other fields of one of the copies: the copy that is already in compartment "s" if there is one, otherwise the copy whose id comes first in character code order (for example ACCOA_c before ACCOA_m). With keepUnconstrained only the copies that are not unconstrained are candidates. If the metabolite ids contain the compartment, the ids of merged metabolites can name a compartment that no longer exists.

mergeModels

function model=mergeModels(models,varargin)

Merge models into one model structure.

Merges models into one model structure. Reactions are added without a check for duplicate reactions, so duplicates can appear. A reaction id that already exists in the merged model is renamed by appending "_" and the id of its source model. Metabolites are matched by name and compartment (metaboliteName[comp]) and genes by name.

Parameters

models : cell

a cell array with model structures.

Name-Value Arguments

metParam : char

the metabolite field used for matching: "metNames" (name and compartment) or "mets" (id) (default "metNames").

supressWarnings : logical

true to suppress warnings (default false).

copyToComps : logical

internal flag set by copyToComps when it calls mergeModels to duplicate a model into a new compartment. Left at the default when calling mergeModels directly (default false).

Returns

model : struct

the merged model structure. In addition to the usual fields, it has the fields rxnFrom, metFrom and geneFrom, which hold the id of the model that each reaction, metabolite and gene was taken from. If a model already has these fields, its entries are not modified.

Examples

model = mergeModels(models);

permuteModel

function newModel=permuteModel(model, indexes, type)

Change the order of the reactions or metabolites in a model.

Parameters

model : struct

a model structure.

indexes : double

a vector with the same length as the number of items in the model, which gives the new order of items.

type : char

'rxns' for reactions, 'mets' for metabolites, 'genes' for genes, 'comps' for compartments.

Returns

newModel : struct

an updated model structure.

Examples

newModel = permuteModel(model, indexes, type);

removeBadRxns

function [newModel, removedRxns]=removeBadRxns(model,varargin)

Remove reactions that enable production or consumption from nothing.

Iteratively removes reactions that enable production or consumption of a metabolite without any uptake or excretion.

Parameters

model : struct

a model structure. The model must not allow any uptake or excretion; the function warns about open exchange reactions when the model has no unconstrained field. importModel("filename", false) imports a model without uptake and excretion.

Name-Value Arguments

rxnRules : double

which reactions may be removed (default 1):

• 1 : remove only unbalanced reactions

• 2 : also remove reactions that could not be checked for mass balance

• 3 : any reaction can be removed

ignoreMets : cell or logical or double

a cell array of metabolite ids, a logical vector with one element per metabolite of the model, or a vector of metabolite indexes to exclude from the analysis (default []).

isNames : logical

true if ignoreMets holds metabolite names (false: ids), which excludes the metabolites of that name in all compartments. Applies only if ignoreMets is a cell array (default false).

balanceElements : cell

a cell array with the elements for which to balance the reactions, chosen from the elements defined in parseFormulas (default {"C";"P";"S";"N";"O"}).

refModel : struct

a reference model, intended as a copy of model with uptake and excretion allowed and some objectives (such as production of biomass) constrained to a non-zero flux. A supplied reference model triggers a warning that the feature is not supported, followed by an error.

ignoreIntBounds : logical

true to ignore internal bounds (including reversibility). Exchange reactions are not affected. Ignoring the bounds finds unbalanced solutions that the default constraints do not allow (default false).

printReport : logical

true to print a report (default false).

seed : double

seed for the random number generator. A seed makes the choice among several candidate reactions reproducible (default: the current rng state is not changed).

Returns

newModel : struct

the model structure without the problematic reactions.

removedRxns : cell

a cell array with the ids of the removed reactions.

Notes

The function removes reactions that enable production or consumption of metabolites when no exchange reaction is used. Models inferred automatically from databases can contain unbalanced reactions that allow net production or consumption of metabolites without any consumption or excretion. One cause is a general compound that has a different meaning in different reactions (as DNA has in these two reactions):

dATP + dGTP + dCTP + dTTP <=> DNA + 4 PPi
0.25 dATP + 0.25 dGTP + 0.25 dCTP + 0.25 dTTP <=> DNA + PPi

Such reactions are elementally unbalanced, but the model can contain further unbalanced reactions that are not problematic. The function removes as few reactions as possible so that the model cannot produce or consume anything from nothing. It calls findLeakMetabolite repeatedly, selects at random one of the elementally unbalanced reactions that carry flux in the solution, removes it, and iterates until no metabolite can be produced or consumed. The 'produce' direction runs before the 'consume' direction.

Examples

[newModel, removedRxns] = removeBadRxns(model, rxnRules, ignoreMets, ...
    isNames, balanceElements, refModel, ignoreIntBounds, printReport);

removeGenes

function reducedModel = removeGenes(model,genesToRemove,varargin)

Delete a set of genes from a model.

Parameters

model : struct

a model structure.

genesToRemove : cell or logical or double

either a cell array of gene IDs, a logical vector with the same number of elements as genes in the model, or a vector of indexes to remove.

Name-Value Arguments

removeUnusedMets : logical

remove metabolites that no remaining reaction uses (default false).

removeBlockedRxns : logical

remove reactions that cannot carry flux after the genes are deleted. When false, the bounds of these reactions are set to 0 (default false).

standardizeRules : logical

write the gene rules in the standard format (default true).

Returns

reducedModel : struct

an updated model structure. Genes in genesToRemove that are not in the model are ignored. When none of them is in the model, no gene, reaction or bound changes (with standardizeRules true the gene rules are still written in the standard format). Otherwise, besides genesToRemove, every gene that no remaining reaction used before the removal is deleted as well, also when it was not requested (for example a gene that no reaction uses). A gene that only occurs in a rule part removed with another gene stays, because the reaction used it before the removal (for example G3 in (G1 and G3) or G2 after removing G1). A gene also stays when its reaction remains, even if that reaction is blocked (removeBlockedRxns false), such as the other subunit of a complex that lost one.

Examples

reducedModel = removeGenes(model, genesToRemove, removeUnusedMets, ...
                   removeBlockedRxns, standardizeRules);

removeLowScoreGenes

function [newModel,remGenes] = removeLowScoreGenes(model,geneScores,varargin)

Remove low-score genes from a model.

Removes genes from a model based on their scores, a step used by the tINIT package. The function distinguishes isozymes from subunits of an enzyme complex. Genes are removed from each grRule, subject to the following conditions:

1) At least one gene must remain associated with the reaction.
2) Genes involved in a complex (joined by ANDs) are not removed.

Parameters

model : struct

model structure from which genes are to be removed.

geneScores : double

a vector of scores associated with the model genes, with the same dimensions as model.genes. Genes with a score of zero or above and genes without a score (NaN) remain in the rules; a gene that no rule uses afterwards is removed from model.genes. Genes with a negative score are removed from the rules when the conditions above allow it.

If all genes associated with a reaction have a negative score, the least-negative gene remains; in a tie, the gene listed first in the rule is kept.

A negative-scoring gene that is a subunit in a complex is not removed, but the entire complex may be removed. Example cases:

• Original: G1 or (G2 and G3 and G4); Negative: G1, G2; New: G2 and G3 and G4.

• Original: G1 or (G2 and G3) or (G4 and G5); Negative: G1, G2; New: G4 and G5 [using default complexScoring].

• Original: (G1 and (G2 or G3) and G4); Negative: G2; New: (G1 and G3 and G4).

Name-Value Arguments

isozymeScoring : char

method used to combine the scores of multiple isozymes: "min", "max", "median" or "average" (default "max").

complexScoring : char

method used to combine the scores of enzyme complex subunits: "min", "max", "median" or "average" (default "min").

Returns

newModel : struct

model with updated genes, grRules and rxnGeneMat fields, and with the gene-associated fields (geneShortNames, proteins, geneMiriams, geneFrom, geneComps) trimmed to the remaining genes.

remGenes : cell

the genes removed from model.genes: negative-score genes removed from every grRule and genes that no rule uses, whatever their score. Negative-score genes that are removed from some grRules but not all are not included in this list.

removeMets

function reducedModel=removeMets(model,metsToRemove,varargin)

Delete a set of metabolites from a model.

Parameters

model : struct

a model structure.

metsToRemove : cell or logical or double

either a cell array of metabolite IDs, a logical vector with the same number of elements as metabolites in the model, or a vector of indexes to remove.

Name-Value Arguments

isNames : logical

true if the supplied mets are metabolite names, false if they are IDs. All metabolites with a listed name are deleted, which removes metabolites in several compartments at once. Requires metsToRemove to be a cell array (default false).

removeUnusedRxns : logical

remove reactions that involve no metabolite (default false).

removeUnusedGenes : logical

remove genes that no remaining reaction uses. Takes effect only when removeUnusedRxns is true (default false).

removeUnusedComps : logical

remove compartments that contain no metabolite (default false).

Returns

reducedModel : struct

an updated model structure.

Examples

reducedModel = removeMets(model, metsToRemove, isNames, ...
    removeUnusedRxns, removeUnusedGenes, removeUnusedComps);

removeReactions

function reducedModel=removeReactions(model,rxnsToRemove,varargin)

Delete a set of reactions from a model.

Parameters

model : struct

a model structure.

rxnsToRemove : cell or logical or double

either a cell array of reaction IDs, a logical vector with the same number of elements as reactions in the model, or a vector of indexes to remove.

Name-Value Arguments

removeUnusedMets : logical

remove metabolites that no remaining reaction uses (default false).

removeUnusedGenes : logical

remove genes that no remaining reaction uses (default false).

removeUnusedComps : logical

remove compartments that contain no metabolite. Takes effect only when removeUnusedMets is true (default false).

Returns

reducedModel : struct

an updated model structure.

Examples

reducedModel = removeReactions(model, rxnsToRemove, removeUnusedMets, ...
    removeUnusedGenes, removeUnusedComps);

replaceMets

function [model, removedRxns, idxDuplRxns]=replaceMets(model,metabolite,replacement,varargin)

Replace a metabolite with another already in the model.

Replaces the name and annotation (formula, miriams, charge, deltaG, InChI and SMILES) of a metabolite with those of a replacement metabolite that is already in the model. If this results in duplicate metabolites, the S matrix is updated to use one of them. With identifiers true the replacement metabolite is kept; by name, the metabolite with the lowest index among those of the same name and compartment is kept. At the end, contractModel is run to remove duplicate reactions that this creates.

Parameters

model : struct

a model structure.

metabolite : char

name of the metabolite to be replaced (its identifier if identifiers is true).

replacement : char

name of the replacement metabolite (its identifier if identifiers is true).

Name-Value Arguments

verbose : logical

whether to print the ids of reactions that involve the replaced metabolite (default false).

identifiers : logical

true if 'metabolite' and 'replacement' are metabolite identifiers, false if they are metabolite names (default false).

Returns

model : struct

model structure with the selected metabolites replaced.

removedRxns : cell

not assigned by the function; requesting this output raises an error.

idxDuplRxns : double

not assigned by the function; requesting this output raises an error.

Examples

[model, removedRxns, idxDuplRxns] = replaceMets(model, metabolite, ...
    replacement, verbose);

Notes

The function applies to a model that contains both 'oxygen' and 'o2' as metabolite names. If 'oxygen' and 'o2' are identifiers, 'identifiers' is set to true.

reverseRxns

function model=reverseRxns(model, rxns)

Reverse the direction of reactions in a model.

Negates the stoichiometric coefficients and the objective coefficient of the reactions, and sets the lower bound to the negated upper bound and the upper bound to the negated lower bound. A flux in the new positive direction equals a flux in the original negative direction. Reaction IDs in rxns that are not in the model are ignored.

Parameters

model : struct

a model structure.

rxns : cell

IDs of the reactions to reverse.

Returns

model : struct

an updated model structure.

Examples

model = reverseRxns(model, rxns);

setExchangeBounds

function [exchModel,unusedMets] = setExchangeBounds(model,varargin)

Define exchange flux bounds for a set of metabolites.

Parameters

model : struct

a model structure.

Name-Value Arguments

mets : cell or double

a cell array of metabolite names (case insensitive) or metabolite IDs, or a vector of indices in model.mets (default all exchanged metabolites).

lb : double

lower bound of exchange flux. Either a vector of bounds corresponding to each of the provided metabolites, or a single value that is applied to all (default model.annotation.defaultLB if it exists, otherwise -1000). Without mets, only a single value is accepted.

ub : double

upper bound of exchange flux. Either a vector of bounds corresponding to each of the provided metabolites, or a single value that is applied to all (default model.annotation.defaultUB if it exists, otherwise 1000). Without mets, only a single value is accepted.

closeOthers : logical

close the exchange reactions of all exchanged metabolites that are not in the provided list. This prevents import of these metabolites and leaves their export unchanged. Set to false, with a warning, when the exchange reactions are not all written in the same direction (default true).

mediaOnly : logical

only consider exchange reactions involving exchange to or from the extracellular (media) compartment. Reactions such as "sink" reactions that exchange metabolites directly with an intracellular compartment are therefore ignored, although "getExchangeRxns" identifies such reactions as exchange reactions. The extracellular compartment is the one named "extracellular" in the "compNames" field (case insensitive). The "compNames" and "metComps" fields are required; an error is thrown if either is missing or if no compartment has this name (default false).

Returns

exchModel : struct

a model structure with updated exchange flux bounds for the provided set of metabolites.

unusedMets : cell

metabolites in mets that were not used because they are not involved in any of the considered exchange reactions in the model.

Examples

exchModel = setExchangeBounds(model, mets, lb, ub, closeOthers, ...
    mediaOnly);

Notes

Exchange reactions involving more than one metabolite are ignored.

setParam

function model=setParam(model, paramType, rxnList, params, var)

Set parameters for reactions.

Parameters

model : struct

a model structure.

paramType : char

the type of parameter to set. The value is matched case-insensitively, except for the reset of the other objective coefficients under "obj" (below):

• "lb" : lower bound.

• "ub" : upper bound.

• "eq" : both upper and lower bound (equality constraint).

• "obj" : objective coefficient. The objective coefficients of all other reactions are set to 0 only when every element of paramType contains "obj" in lower case; with another capitalisation, or when paramType also holds other types, the other coefficients keep their values.

• "rev" : reversibility (only changes the model.rev field, does not affect model.lb and model.ub). A non-zero value sets the reaction reversible.

• "var" : variance around measured bound.

• "unc" : unconstrained, sets the lower and upper bound to model.annotation.defaultLB and model.annotation.defaultUB if these are defined, otherwise to -1000 and 1000.

rxnList : cell or double

a cell array of reaction IDs or a vector with their corresponding indexes. IDs that are not in the model are ignored with a warning.

params : double

a vector of the corresponding values, or a single value that is applied to all reactions in rxnList.

var : double, optional

percentage of variance around measured value, if "var" is set as paramType. Defining "var" as 5 results in lb and ub at 97.5% and 102.5% of the provided params value (if params value is negative, then lb and ub are 102.5% and 97.5%).

Returns

model : struct

an updated model structure.

Examples

model = setParam(model, paramType, rxnList, params, var);

simplifyModel

function [reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,varargin)

Simplify a model by deleting reactions and metabolites.

Deletes reactions and metabolites that cannot carry flux. The deleted reactions and metabolites identify different types of gaps.

Parameters

model : struct

a model structure.

Name-Value Arguments

deleteUnconstrained : logical

delete metabolites marked as unconstrained in model.unconstrained, and the compartments that are then left without metabolites (default true).

deleteDuplicates : logical

delete all but the first occurrence of duplicate reactions, as found by contractModel. Duplicate reactions have the same stoichiometry and reversibility. The kept reaction takes the minimum lower bound, the maximum upper bound and the maximum objective coefficient of the duplicates, with a warning when these differ (default false).

deleteZeroInterval : logical

delete reactions that are constrained to zero flux (default false).

deleteInaccessible : logical

delete dead end reactions (default false).

deleteMinMax : logical

delete reactions that cannot carry a flux, as determined by haveFlux under the bounds of the model. The run time grows with the number of reactions, since haveFlux solves linear programs (default false).

groupLinear : logical

group linearly dependent pathways. The grouped model is irreversible, and all gene information (genes, rxnGeneMat, grRules and the gene-associated fields) is deleted (default false).

constrainReversible : logical

check if there are reversible reactions which can only carry flux in one direction, and if so constrain them to be irreversible. A reversible reaction that is constrained to one direction is not split into two reactions when groupLinear converts the model to irreversible, so the metabolites that it connects can meet the condition for grouping (one producing and one consuming reaction) (default false).

reservedRxns : cell

cell array with reaction IDs that are not deleted by deleteZeroInterval, deleteInaccessible and deleteMinMax and not merged by groupLinear (default none).

suppressWarnings : logical

true if warnings should be suppressed (default false).

Returns

reducedModel : struct

an updated model structure.

deletedReactions : cell

a cell array with the IDs of all deleted reactions.

deletedMetabolites : cell

a cell array with the IDs of all deleted metabolites.

Examples

[reducedModel, deletedReactions, deletedMetabolites] = ...
    simplifyModel(model, deleteUnconstrained, deleteDuplicates, ...
    deleteZeroInterval, deleteInaccessible, deleteMinMax, ...
    groupLinear, constrainReversible, reservedRxns, suppressWarnings);

sortIdentifiers

function newModel = sortIdentifiers(model)

Sort model identifiers alphabetically.

Sorts reactions, metabolites, genes (if the genes field exists) and compartments alphabetically by their identifier.

Parameters

model : struct

a model structure.

Returns

newModel : struct

an updated model structure with alphabetically sorted identifiers.

Examples

newModel = sortIdentifiers(model);

sortModel

function model=sortModel(model,varargin)

Sort the metabolites and reactions of a model.

Depending on the options, the function sorts the metabolites by name and compartment, orients reversible reactions consistently, and reorders the reactions within each subsystem.

Parameters

model : struct

a model structure.

Name-Value Arguments

sortReversible : logical

orient each reversible reaction by changing the sign of its stoichiometric coefficients. The metabolite that is first in lexicographical order (metabolite name followed by compartment ID) becomes a reactant. Not applied when sortReactionOrder is true (default false).

sortMetName : logical

sort the metabolites of the model by metabolite name followed by compartment ID (default false).

sortReactionOrder : logical

reorder the reactions within each subsystem, such that reactions consuming a metabolite come after reactions producing it. The order is found by a random local search that also changes the direction of reversible reactions to match production and consumption direction. Requires the subSystems field. Subsystems with fewer than 2 or more than 250 reactions are not reordered (default false).

seed : double

seed for the random number generator used by the sortReactionOrder local search, which makes its result reproducible across runs (default: the current rng state is not changed).

Returns

model : struct

an updated model structure.

Examples

model = sortModel(model, sortReversible, sortMetName, sortReactionOrder);

standardizeGrRules

function [grRules,rxnGeneMat,indexes2check] = standardizeGrRules(model,varargin)

Standardise gene-reaction rules in a model.

The grRules are standardised according to the following:

• No overall containing brackets

• Only enzyme complexes are enclosed in brackets

• " and " and " or " strings are set to lowercase

A rxnGeneMat matrix consistent with the standardised grRules is created. grRules that are not in disjunctive normal form (an "or" nested inside an "and") or that cannot be parsed are not standardised, apart from the removal of double spaces, spaces inside brackets and uppercase operators. The function throws an error if the model has no grRules field.

Parameters

model : struct

a model structure.

Name-Value Arguments

embedded : logical

if true, suppresses the console warning about grRules that are not in disjunctive normal form or cannot be parsed, for use when this function is called inside another RAVEN function that handles reporting itself (default false).

Returns

grRules : cell

[nRxns x 1] cell array with the standardised grRules.

rxnGeneMat : double

[nRxns x nGenes] sparse matrix consistent with the standardised grRules.

indexes2check : double

indices of the grRules that are not in disjunctive normal form or cannot be parsed. These grRules need a manual check.

Examples

[grRules, rxnGeneMat, indexes2check] = standardizeGrRules(model, embedded);

Notes

If used in a model reconstruction or modification pipeline, run this function first.