Gap-filling#

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

Functions#

Function	Summary
`canExchange`	Check which metabolites a model can produce or consume.
`checkProduction`	Check which metabolites can be produced from a model.
`fillGaps`	Use template model(s) to fill gaps in a model.
`findLeakMetabolite`	Find a metabolite that can be freely made or consumed.
`fitTasks`	Fill gaps in a model so it can perform a list of tasks.
`ftINITFillGaps`	Variant of fillGaps specially adapted to speed up generation of ftINIT models.
`gapFillFastCore`	L1-norm LP subroutine for fastGapFill.
`gapFillFastLP`	LP-based gap-filling (fastGapFill / swiftGapFill).
`gapFillMILP`	Objective-based gap-filling via a global growth-floor MILP.
`gapFillSwiftCore`	SWIFTCORE LP subroutine for swiftGapFill.
`gapFillTopological`	BFS metabolite-producibility pre-screening.
`gapReport`	Perform a gap analysis and summarize the results.

Reference#

canExchange#

function canExchange=canExchange(model,direction,varargin)

Check which metabolites a model can produce or consume.

Adds an exchange reaction for each of the specified metabolites and checks which of them can carry a flux. This is the unified replacement for canProduce ('produce') and canConsume ('consume').

This is a less advanced but faster version of checkProduction, which additionally reports which metabolites must be connected in order to make the remaining ones producible.

Parameters

model : struct: a model structure.
direction : char: 'produce' to check for net production of each metabolite, or 'consume' to check for net consumption.

Name-Value Arguments

mets : 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 check for (default model.mets).

Returns

canExchange : logical: vector with true if the corresponding metabolite could be produced (direction 'produce') or consumed (direction 'consume').

Examples

produced = canExchange(model, 'produce');
consumed = canExchange(model, 'consume', mets);

See also

checkProduction, findLeakMetabolite

checkProduction#

function [notProduced, notProducedNames, neededForProductionMat,minToConnect,model]=checkProduction(model,varargin)

Check which metabolites can be produced from a model.

Checks which metabolites that can be produced from a model using the specified constraints.

The function is intended to be used to identify which metabolites must be connected in order to have a fully connected network. It does so by first identifying which metabolites could have a net production in the network. Then it calculates which other metabolites must be able to have net production in order to have production of all metabolites in the network. So, if a network contains the equations A[external]->B, C->D, and D->E it will identify that production of C will connect the metabolites D and E.

Parameters

model : struct: a model structure.

Name-Value Arguments

checkNeededForProduction : logical: for each of the metabolites that could not be produced, include an artificial production reaction and calculate which new metabolites that could be produced as an effect of this (default false).
excretionFromCompartments : cell: cell array with compartment ids from which metabolites can be excreted (default model.comps).
printDetails : logical: print details to the screen (default true).

Returns

notProduced : double: indices of metabolites that could not be produced.
notProducedNames : cell: cell array with names and compartments for metabolites that could not be produced.
neededForProductionMat : logical: matrix where n x m is true if metabolite n allows for production of metabolite m.
minToConnect : cell: the minimal number of metabolites that need to be connected in order to be able to produce all other metabolites, and which metabolites each of them connects.
model : struct: updated model structure with excretion reactions added.

Examples

[notProduced, notProducedNames, neededForProductionMat, minToConnect, model] = ...
    checkProduction(model, checkNeededForProduction, excretionFromCompartments, printDetails);

fillGaps#

function [newConnected, cannotConnect, addedRxns, newModel, exitFlag]=fillGaps(model,models,varargin)

Use template model(s) to fill gaps in a model.

This method works by merging the model with the reference model(s) and checking which reactions can carry flux. All reactions that cannot carry flux are removed (cannotConnect). If useModelConstraints is false it then solves the MILP problem of minimizing the number of active reactions from the reference models that are required to have flux in all the reactions in model. This requires that the input model has exchange reactions present for the nutrients that are needed for its metabolism. If useModelConstraints is true then the problem is to include as few reactions as possible from the reference models in order to satisfy the model constraints.

The intended use is that the user can attempt a general gap-filling using useModelConstraints=false, or a more targeted gap-filling by setting constraints in the model structure and then using useModelConstraints=true. For example, to include reactions so that all biomass components can be synthesized, the user could define a biomass equation and set its lower bound to >0. Running this function with useModelConstraints=true would then give the smallest set of reactions that have to be included in order for the model to produce biomass.

Parameters

model : struct: a model structure that may contain gaps to be filled.
models : cell or struct: a cell array of reference models or a model structure. The gaps will be filled using reactions from these models.

Name-Value Arguments

allowNetProduction : logical

true if net production of all metabolites is allowed. A reaction can be unable to carry flux because one of the reactants is unavailable or because one of the products cannot be further processed. If true, only the first type of unconnectivity is considered (default false).

useModelConstraints : logical

true if the constraints specified in the model structure should be used. If false then reactions are included from the template model(s) so that as many reactions as possible in model can carry flux (default false).

supressWarnings : logical

false if warnings should be displayed (default false).

rxnScores : double or cell

array with scores for each of the reactions in the reference model(s). If more than one model is supplied, then rxnScores should be a cell array of vectors. The solver will try to maximize the sum of the scores for the included reactions (default is -1 for all reactions).

algorithm : char

gap-filling algorithm to use (default 'reference'):

'reference' : MILP-based connectivity maximisation
'fastLP' : L1-norm LP relaxation via FASTCORE (gapFillFastLP)
'swiftLP' : SWIFTCORE single-LP variant (gapFillFastLP)
'gapfillMILP' : growth-floor MILP adding database reactions (gapFillMILP). Uses only models{1} as the database.
'topological' : BFS metabolite producibility pre-screening (gapFillTopological); no model modification.

epsilon : double

minimum flux threshold used by fastLP and swiftLP (default 1e-4).

minGrowth : double

minimum growth rate for gapfillMILP (default [] = 10% of max FBA).

verbose : logical

print progress messages for new algorithms (default true).

Returns

newConnected : cell

cell array with the reactions that could be connected. This is not calculated if useModelConstraints is true.

cannotConnect : cell

cell array with reactions that could not be connected. This is not calculated if useModelConstraints is true.

addedRxns : cell

cell array with the reactions that were added from "models".

newModel : struct

the model with reactions added to fill gaps.

exitFlag : double

exit status:

1 : optimal solution found
-1 : no solution found, either because the problem is infeasible or because the solver reached its time limit before finding one
-2 : a solution was found but is not proven optimal, because the solver stopped on its time limit

Examples

[newConnected, cannotConnect, addedRxns, newModel, exitFlag]=...
    fillGaps(model,models,allowNetProduction,useModelConstraints,...
    supressWarnings,rxnScores,params);

Notes

The 'reference' and 'gapfillMILP' algorithms solve a MILP with the solver set by setRavenSolver. RAVEN's bundled GLPK solves it for small models (with a RAVEN:glpkMILP warning); for larger models Gurobi is much faster and more reliable (setRavenSolver('gurobi')).

findLeakMetabolite#

function [solution, metabolite]=findLeakMetabolite(model, direction, varargin)

Find a metabolite that can be freely made or consumed.

Tests whether the model can produce (excrete) or consume (take up) any metabolite without balanced stoichiometric constraints — a sign that the model contains a stoichiometric leak, futile cycle, or unconstrained exchange. This is the unified replacement for makeSomething ('produce') and consumeSomething ('consume').

Parameters

model : struct: a model structure.
direction : char: 'produce' to look for freely excreted metabolites, or 'consume' to look for freely consumed metabolites.

Name-Value Arguments

ignoreMets : 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 for metabolites to exclude from this analysis (default []).
isNames : logical: true if the supplied mets represent metabolite names (as opposed to IDs). This is a way to delete metabolites in several compartments at once without knowing the exact IDs. This only works if ignoreMets is a cell array (default false).
minNrFluxes : logical: solves the MILP problem of minimizing the number of fluxes instead of the sum. Slower, but can be used if the sum gives too many fluxes (default false).
allowExcretion : logical: allow for excretion of all metabolites. Only used when direction is 'produce' (default true).
params : struct: obsolete option.
ignoreIntBounds : logical: true if internal bounds (including reversibility) should be ignored. Exchange reactions are not affected. This can be used to find unbalanced solutions which are not possible using the default constraints (default false).

Returns

solution : double: flux vector for the solution.
metabolite : double: the index of the metabolite(s) which was produced or consumed. If possible only one metabolite is reported, but there are situations where metabolites can only be exchanged in pairs (or more).

Examples

[solution, metabolite] = findLeakMetabolite(model, 'produce');
[solution, metabolite] = findLeakMetabolite(model, 'consume', ignoreMets);

Notes

Works by forcing at least 1 unit of "any metabolite" to be produced or consumed, then minimising for the sum of fluxes. If more than one metabolite is produced/consumed, it tries to find a single metabolite that satisfies the constraint on its own.

fitTasks#

function [outModel, addedRxns, failedTasks]=fitTasks(model,refModel,inputFile,varargin)

Fill gaps in a model so it can perform a list of tasks.

Fills gaps in a model by including reactions from a reference model, so that the resulting model can perform all the tasks in a task list. The gap-filling is done in a task-by-task manner, rather than solving for all tasks at once. This means that the order of the tasks could influence the result.

Parameters

model : struct: a model structure.
refModel : struct: reference model from which to include reactions.
inputFile : char, optional: a task list in Excel format. See the function parseTaskList for details (optional if taskStructure is supplied).

Name-Value Arguments

printOutput : logical

true if the results of the test should be displayed (default true).

rxnScores : double

scores for each of the reactions in the reference model. Only negative scores are allowed. The solver will try to maximize the sum of the scores for the included reactions (default is -1 for all reactions).

taskStructure : struct

structure with the tasks, as from parseTaskList. If supplied then inputFile is ignored.

params : struct

solver parameters, forwarded to the gap-filling back-end and from there to optimizeProb for each per-task MILP (default []).

gapFillMode : char

which gap-filling back-end to use (default 'merge'):

'merge' : refModel holds only the candidate reactions, and fillGaps merges it with the model for every task.
'preMerged' : refModel already contains the model's own reactions, so no merge is needed per task. The task constraints are written into refModel as well, and ftINITFillGaps does the gap-filling.

verbose : logical

if true, the MILP progression will be shown. Only read when gapFillMode is 'preMerged' (default false).

resolveTies : logical

if true, pin the gap-fill MILP's degenerate optimum to a canonical answer instead of relying on the solver Seed alone. Only read when gapFillMode is 'preMerged' -- see ftINITFillGaps' resolveTies (default false).

Returns

outModel : struct: model structure with reactions added to perform the tasks.
addedRxns : logical: MxN matrix with the added reactions (M) from refModel for each task (N). An element is true if the corresponding reaction is added in the corresponding task. Failed tasks and SHOULD FAIL tasks are ignored.
failedTasks : logical: Nx1, true for each task that could not be gap-filled, either because no feasible solution exists using the reference model or because the attempt threw. Such tasks add no reactions, so they are otherwise indistinguishable from tasks that already worked.

Examples

[outModel, addedRxns]=fitTasks(model,refModel,inputFile,printOutput,...
    rxnScores,taskStructure,params);

ftINITFillGaps#

function [addedRxns, newModel, exitFlag]=ftINITFillGaps(tModel, origModel, tRefModel,allowNetProduction,supressWarnings,rxnScores,params,verbose,resolveTies)

Variant of fillGaps specially adapted to speed up generation of ftINIT models.

  tModel              model that contains the task-specific rxns.
  origModel           model without task-specific rxns.
  tRefModel           reference tModel - the full tModel, containing all rxns
                      used for gap-filling of tModel + the task-specific rxns
  allowNetProduction  true if net production of all metabolites is
                      allowed. A reaction can be unable to carry flux because one of
                      the reactants is unavailable or because one of the
                      products can't be further processed. If this
                      parameter is true, only the first type of
                      unconnectivity is considered (optional, default false)
  supressWarnings     false if warnings should be displayed (optional, default
                      false)
  rxnScores           scores for each of the reactions in the
                      reference tModel.
                      The solver will try to maximize the sum of the
                      scores for the included reactions
  params              *obsolete option*
  verbose             if true, the MILP progression will be shown.
  resolveTies         if true, pin the gap-fill MILP's degenerate optimum to a
                      canonical answer instead of relying on the solver Seed
                      alone -- see getMinNrFluxes' resolveTies (optional,
                      default false)

  addedRxns           the rxns added
  newModel            the tModel with reactions added to fill gaps
  exitFlag            1: optimal solution found
                     -1: no solution found, either because the problem is
                         infeasible or because the solver reached its time
                         limit before finding one
                     -2: a solution was found but is not proven optimal,
                         because the solver stopped before reaching
                         optimality

  This method works by merging the tModel to the reference model and
  checking which reactions can carry flux. All reactions that can't
  carry flux are removed. It then solves the MILP problem of including
  as few additional reactions as possible from the reference model so
  that the constraints already set on tModel (e.g. a required biomass
  flux, or exchange bounds) can be satisfied.

  See also: fillGaps, the tINIT-side counterpart, which merges a
  reference model with the model being gap-filled on every call instead
  of expecting them pre-merged.

Usage: [addedRxns, newModel, exitFlag]=...
          ftINITFillGaps(tModel,origModel,tRefModel,allowNetProduction,...
          supressWarnings,rxnScores,params,verbose,resolveTies)

gapFillFastCore#

function activeRxns = gapFillFastCore(model, coreRxns, epsilon)

L1-norm LP subroutine for fastGapFill.

Finds the minimum-L1-flux extension of the model that activates all reactions in the core set. Any reaction with non-zero flux in the solution is "consistent" with the core. Used internally by gapFillFastLP.

Based on FASTCORE (Vlassis et al. 2014, PLoS Comput Biol 10:e1003424).

Parameters

model : struct: a RAVEN model structure.
coreRxns : cell or logical or double: core reaction set: cell array of reaction IDs, a logical vector, or a vector of reaction indices into model.rxns. All core reactions are forced to carry flux >= epsilon in the LP solution.
epsilon : double: minimum flux threshold for a reaction to be considered active.

Returns

activeRxns : logical: logical vector (length numel(model.rxns)) where true indicates the reaction carries flux in the L1-minimal solution. Empty if the LP is infeasible (core reactions cannot all carry flux simultaneously).

See also

gapFillSwiftCore, gapFillFastLP

gapFillFastLP#

function [addedRxns, newModel, cannotConnect] = gapFillFastLP(model, universalModel, varargin)

LP-based gap-filling (fastGapFill / swiftGapFill).

Merges the draft model with a universal reaction database, then for each blocked draft reaction calls an LP subroutine (FASTCORE or SWIFTCORE) to find the minimum set of database reactions that makes the blocked reaction carry flux. The union of all activated database reactions across all blocked reactions is returned as addedRxns.

Based on fastGapFill (Thiele et al. 2014, PLoS Comput Biol 10:e1003515) and swiftGapFill (Tefagh & Boyd 2020, BMC Bioinformatics 21:23).

Parameters

model : struct: draft RAVEN model to be gap-filled.
universalModel : struct: universal reaction database (RAVEN model). The caller is responsible for building this, e.g. via getKEGGModelForOrganism with broad scope.

Optional name-value parameters

'epsilon' (default 1e-4) minimum flux threshold for a reaction to be considered active. 'variant' (default 'fast') 'fast' — use FASTCORE L1-norm LP (gapFillFastCore) 'swift' — use SWIFTCORE single-LP (gapFillSwiftCore); faster but stochastic (results may vary between runs/solvers). 'verbose' (default true) print progress messages.

Returns

addedRxns : cell: reaction IDs from universalModel that should be added to rescue blocked reactions.
newModel : struct: copy of model with addedRxns incorporated.
cannotConnect : cell: blocked draft reaction IDs that remain blocked even after augmenting with the full universalModel (cannot be rescued by the database).

See also

gapFillFastCore, gapFillSwiftCore, gapFillMILP, fillGaps

gapFillMILP#

function [addedRxns, newModel, exitFlag] = gapFillMILP(model, universalModel, varargin)

Objective-based gap-filling via a global growth-floor MILP.

Finds the smallest set of reactions from a universal database that makes the model's objective (biomass) >= minGrowth, following the database repair of Kumar et al. 2007 (BMC Bioinformatics 8:212, doi:10.1186/1471-2105-8-212). Reactions already in the draft model keep their bounds: a reaction with the wrong direction is a curation issue and is not changed here.

The universal model should contain all candidate reactions, including exchange and transport reactions if those repair classes are desired.

This is a single-level MILP. It is computationally more expensive than gapFillFastLP but directly optimises for objective feasibility.

Parameters

model : struct: draft RAVEN model to be gap-filled.
universalModel : struct: universal reaction database (RAVEN model).

Optional name-value parameters

'minGrowth' (default []) minimum objective value to achieve. If empty, runs FBA on the merged model and uses 10% of the maximum as the floor. 'bigM' (default 1000) Big-M constant for coupling constraints. Must be larger than any expected flux magnitude. 'params' (default []) solver parameter struct passed to optimizeProb. Fields such as params.TimeLimit and params.intTol can be used to tune the MILP. 'verbose' (default true) print progress messages.

Returns

addedRxns : cell: reaction IDs from universalModel that are selected by the MILP.
newModel : struct: model with addedRxns incorporated.
exitFlag : double: 1 = optimal, -1 = infeasible, -2 = timeout or other failure.

Notes

The MILP is solved with the solver set by setRavenSolver. RAVEN's bundled GLPK solves it for small models (with a RAVEN:glpkMILP warning); for larger models Gurobi is much faster and more reliable (setRavenSolver('gurobi')).

See also

gapFillFastLP, gapFillTopological, fillGaps

gapFillSwiftCore#

function activeRxns = gapFillSwiftCore(model, coreRxns, epsilon)

SWIFTCORE LP subroutine for swiftGapFill.

Single-LP alternative to gapFillFastCore. Maximises the sum of non-core fluxes while enforcing that every core reaction carries at least epsilon flux. Any reaction with non-zero flux in the solution is "consistent" with the core.

Because the maximisation objective makes the LP degenerate (many optima exist that differ only in which non-core reactions carry flux), the returned active set can vary between LP solvers and between runs with slightly different numerical conditions. This is the expected behaviour of SWIFTCORE (Tefagh & Boyd 2020).

Based on SWIFTCORE (Tefagh & Boyd 2020, BMC Bioinformatics 21:23).

Parameters

model : struct: a RAVEN model structure.
coreRxns : cell or logical or double: core reaction set: cell array of reaction IDs, a logical vector, or a vector of reaction indices into model.rxns. All core reactions are forced to carry flux >= epsilon in the LP solution.
epsilon : double: minimum flux threshold for a reaction to be considered active.

Returns

activeRxns : logical: logical vector (length numel(model.rxns)) where true indicates the reaction carries flux in the solution. Empty on infeasibility.

See also

gapFillFastCore, gapFillFastLP

gapFillTopological#

function result = gapFillTopological(model, universalModel, varargin)

BFS metabolite-producibility pre-screening.

Identifies which metabolites in the draft model cannot be produced from a given set of seed metabolites (medium), using a graph-based fixed-point "scope" computation (no solver required). For each unreachable metabolite, finds candidate reactions in the universal database that produce it directly.

This is a computationally cheap pre-screening step that prunes the universal database before a more expensive gap-filling solve. It is inspired by the topological analysis used in Meneco (Prigent et al. 2017, PLoS Comput Biol).

Algorithm

Starting from seed metabolites (available in the medium), iteratively mark reactions as "fireable" when all their substrates are reachable, then mark their products as newly reachable. Continue until no new metabolite is reached (fixed point). Reversible reactions (lb < 0) are also considered in their reverse direction.

Parameters

model : struct: draft RAVEN model.
universalModel : struct: universal reaction database (RAVEN model).

Optional name-value parameters

'seeds' (default []) Metabolite IDs (cell array) available from the medium. Default: metabolites of every exchange reaction whose bounds let it supply the metabolite, i.e. getExchangeRxns 'uptake' plus 'reverse'. 'targets' (default []) Metabolite IDs (cell array) that should be produced. Default: substrates of the objective (biomass) reaction. 'verbose' (default true) print a summary of the results.

Returns

result : struct: .reachableMets — logical vector (length nMets) marking producible mets .blockedMets — cell array of metabolite IDs that cannot be reached .candidateRxns — n-by-1 cell array of universal rxn ID lists, one cell per blocked metabolite (same order as blockedMets) .pruningFraction — fraction of universal reactions pruned (not needed for any blocked metabolite)

See also

gapFillFastLP, gapFillMILP, fillGaps

gapReport#

function [noFluxRxns, noFluxRxnsRelaxed, subGraphs, notProducedMets, minToConnect,...
    neededForProductionMat, canProduceWithoutInput, canConsumeWithoutOutput, ...
    connectedFromTemplates, addedFromTemplates]=gapReport(model, varargin)

Perform a gap analysis and summarize the results.

Runs a series of connectivity checks on the model - which reactions are blocked, how many isolated sub-networks exist, which metabolites cannot be produced (and the minimal set of metabolites that would need to be connected to fix that), and which metabolites can be produced/consumed without any input/output, with every boundary reaction closed. If templateModels is supplied, it also attempts to gap-fill using those models. A formatted report is printed to the command window; the same results are returned as outputs for further use.

Parameters

model : struct: a model structure.

Name-Value Arguments

templateModels : cell: a cell array of template models to use for gap filling. If omitted, the automated gap-filling section is skipped.

Returns

noFluxRxns : cell: reactions that cannot carry flux.
noFluxRxnsRelaxed : cell: reactions that cannot carry flux even if the mass balance constraint is relaxed so that net production of all metabolites is allowed.
subGraphs : struct: the metabolites in each of the isolated sub-networks.
notProducedMets : cell: the metabolites that could not have net production.
minToConnect : struct: the minimal number of metabolites that need to be connected in order to be able to produce all other metabolites, and which metabolites each of them connects.
neededForProductionMat : double: matrix where n x m is true if metabolite n allows for production of metabolite m.
canProduceWithoutInput : cell: metabolites that could be produced with every boundary reaction closed.
canConsumeWithoutOutput : cell: metabolites that could be consumed with every boundary reaction closed.
connectedFromTemplates : cell: the reactions that could be connected using the template models.
addedFromTemplates : struct: the reactions that were added from the template models and which model they were added from.

Notes

The boundary reactions are those getExchangeRxns returns: exchanges, sinks and demands. For a model that keeps its boundary metabolites (an unconstrained field), those metabolites are deleted for the connectivity checks and kept, balanced, for the closed model.

The model must be feasible under its own bounds and with its boundary closed. Otherwise gapReport stops with a RAVEN:infeasibleModel error that names the reactions with a forced flux; set the bounds of ATP maintenance and any other forced flux to zero first.

Examples

[noFluxRxns, noFluxRxnsRelaxed, subGraphs, notProducedMets, ...
    minToConnect, neededForProductionMat, connectedFromTemplates, ...
    addedFromTemplates] = gapReport(model, templateModels);