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 adapted to speed up the generation of ftINIT models.
`gapFillFastCore`	L1-norm LP subroutine for fastGapFill.
`gapFillFastLP`	LP-based gap-filling (fastGapFill or 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 summarise 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 an absolute flux above 1e-5. The direction 'produce' corresponds to the deprecated function canProduce and 'consume' to canConsume.

canExchange only tests whether each added exchange reaction can carry flux. checkProduction also reports which metabolites must be connected in order to make the remaining ones producible, which requires further problems to be solved.

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 can be produced from a model under the current constraints.

The function identifies which metabolites must be connected in order to have a fully connected network. It first identifies which metabolites can have a net production in the network. It then calculates which other metabolites must be able to have net production in order to have production of all metabolites in the network. If a network contains the equations A[external]->B, C->D, and D->E, production of C connects the metabolites D and E.

Parameters

model : struct: a model structure.

Name-Value Arguments

checkNeededForProduction : logical: true to calculate, for each metabolite that could not be produced, which of the other such metabolites can be produced when an uptake reaction for it is added (default false).
excretionFromCompartments : cell: cell array with compartment ids from which metabolites can be excreted (default model.comps).
printDetails : logical: true to print the metabolites that could not be produced to the command window (default true).

Returns

notProduced : double: indices of metabolites that could not be produced, counted among the metabolites in the compartments of excretionFromCompartments, in model order.
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. The rows and columns follow the order of notProduced. It is empty if checkNeededForProduction is false.
minToConnect : cell: metabolites that need to be connected in order to be able to produce all other metabolites, each with the number of metabolites it connects. The metabolites are selected iteratively, each time the one that connects the most remaining metabolites, so the number of metabolites need not be the global minimum. It is empty if checkNeededForProduction is false.
model : struct: model structure with an excretion reaction added for each metabolite in the compartments of excretionFromCompartments. All added reactions are excretion reactions, also when checkNeededForProduction is true.

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.

The model is merged with the reference model(s). If useModelConstraints is false, the reactions of the model that cannot carry flux in the model are tested for flux in the merged model. The non-reversible reactions that carry flux in the merged model are the newConnected reactions. The other reactions are returned as cannotConnect: the reactions that cannot carry flux in the merged model and the reversible reactions that can, which are left out of the MILP target set. A MILP then minimises the number of active reactions from the reference models that are required for flux in the reactions that can be connected (newConnected). This requires that the input model has exchange reactions for the nutrients that its metabolism needs. If useModelConstraints is true, the flux test and the split into newConnected and cannotConnect are not run, and the MILP includes as few reactions as possible from the reference models to satisfy the model constraints.

A general gap-filling uses useModelConstraints=false. A targeted gap-filling sets constraints in the model structure and uses useModelConstraints=true. For example, a model with a biomass equation whose lower bound is set above 0 gets the smallest set of reactions from the reference models that allows it 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. Reactions from these models are used to fill the gaps.

Name-Value Arguments

allowNetProduction : logical

true to allow net production of all metabolites. 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 cause is considered (default false).

useModelConstraints : logical

true to use the constraints specified in the model structure. If false, reactions are included from the template model(s) so that as many reactions as possible in model can carry flux (default false).

supressWarnings : logical

true to suppress the warning and the note about missing exchange reactions (default false).

rxnScores : double or cell

scores for each of the reactions in the reference model(s). If more than one model is supplied, rxnScores is a cell array of vectors. The solver maximises the sum of the scores of the included reactions. Used by the 'reference' algorithm (default -1 for all reactions).

params : struct

solver parameters passed to getMinNrFluxes (algorithm 'reference') or to gapFillMILP (algorithm 'gapfillMILP') (default empty struct).

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 [], which is 10% of the maximum growth rate found by FBA on the merged model).

verbose : logical

true to print progress messages for the algorithms other than 'reference' (default true).

Returns

newConnected : cell

cell array with the non-reversible reactions of model that cannot carry flux in model but can in the merged model. Empty if useModelConstraints is true or if the algorithm is not 'reference'.

cannotConnect : cell

cell array with the other reactions of model that cannot carry flux in model. For the algorithms 'fastLP' and 'swiftLP' it holds the reactions that remain blocked after merging (from gapFillFastLP), and for 'topological' the ids of the metabolites that cannot be reached. Empty if useModelConstraints is true or for 'gapfillMILP'.

addedRxns : cell

cell array with the reactions 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. The RAVEN:glpkMILP warning issued with GLPK states that GLPK works for small problems and that Gurobi is much faster and more reliable for larger ones (setRavenSolver('gurobi')).

See also

ftINITFillGaps : the variant of fillGaps for ftINIT, used when the reference model already contains the reactions of the model. No merge per task is needed in that case.

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 positive result indicates that the model contains a stoichiometric leak, a futile cycle or an unconstrained exchange. The direction 'produce' corresponds to the deprecated function makeSomething and 'consume' to consumeSomething.

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 are metabolite names and false if they are IDs. Names exclude the metabolites of all compartments that carry the name, so the exact IDs are not needed. Only works if ignoreMets is a cell array (default false).
minNrFluxes : logical: true to solve the MILP problem of minimising the number of fluxes. If false, the sum of fluxes is minimised. The MILP is slower and applies when the sum gives too many fluxes (default false).
allowExcretion : logical: true to allow excretion of all metabolites. Only used when direction is 'produce' (default true).
params : struct: solver parameters, passed to solveLP when minNrFluxes is true (default a structure with relGap 0.8).
ignoreIntBounds : logical: true to ignore internal bounds (including reversibility). Exchange reactions are not affected. Finds unbalanced solutions that the default constraints do not allow (default false).

Returns

solution : double: flux vector for the solution.
metabolite : double: the index of the metabolite(s) produced or consumed. A single metabolite is reported when one suffices; some metabolites can only be exchanged in pairs or larger groups.

Examples

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

Notes

The function forces a flux of at least 1 through a pooled exchange of all metabolites, then minimises the sum of fluxes. If more than one metabolite is produced or consumed, it searches for a single metabolite that satisfies the constraint on its own.

See also

canExchange

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 tasks are gap-filled one at a time in list order, and the reactions added for one task are kept for the following tasks. The order of the tasks can therefore 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 maximises the sum of the scores of the included reactions (default is -1 for all reactions).

taskStructure : struct

structure with the tasks, as returned by parseTaskList. If supplied, inputFile is ignored.

params : struct

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

gapFillMode : char

which gap-filling function 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

true to show the MILP progression. Only read when gapFillMode is 'preMerged' (default false).

resolveTies : logical

true to fix the degenerate optimum of the gap-fill MILP to a canonical answer, which otherwise depends on the solver Seed. Only read when gapFillMode is 'preMerged'. See resolveTies of ftINITFillGaps (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 raised an error. Such tasks add no reactions, and the same holds for tasks that are feasible without gap-filling.

Examples

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

ftINITFillGaps#

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

Variant of fillGaps adapted to speed up the 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 cannot be processed further. If this
                      parameter is true, only the first cause of blocked
                      flux 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 maximises the sum of
                      the scores of the included reactions (optional,
                      default -1 for all reactions)
  params              solver parameters, passed to getMinNrFluxes (optional)
  verbose             true to show the MILP progression.
  resolveTies         true to fix the degenerate optimum of the gap-fill MILP
                      to a canonical answer, which otherwise depends on the
                      solver Seed. See resolveTies of getMinNrFluxes
                      (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 on the reference model, which already contains the
  tModel reactions. Reactions constrained to zero flux are removed
  (simplifyModel(tRefModel,false,false,true), deleteZeroInterval). 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 (for example a required biomass flux, or exchange bounds) can
  be satisfied.

  fillGaps merges the reference model with the model on every call, so it
  does not need a pre-merged reference model.

  See also: fillGaps

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 a flux of at least epsilon / 2 in the solution is "consistent" with the core. gapFillFastLP calls this function.

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 that the reaction carries flux in the L1-minimal solution. All elements are false if the LP is infeasible (the 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 or 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), for example built with getKEGGModelForOrganism with broad scope.

Name-Value Arguments

epsilon : double

minimum flux threshold for a reaction to be considered active (default 1e-4).

variant : char

the LP subroutine (default 'fast'):

'fast' : use the FASTCORE L1-norm LP (gapFillFastCore).
'swift' : use the SWIFTCORE single LP (gapFillSwiftCore), which maximises the sum of the non-core fluxes. The LP is degenerate, so results can vary between runs and solvers.

verbose : logical

true to print progress messages (default true).

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: IDs of blocked draft reactions that remain blocked after merging with the full universalModel (the database cannot rescue them).

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, including a wrong reaction direction.

The universal model contains all candidate reactions. Exchange and transport reactions are candidates only if the universal model includes them.

The function solves a single-level MILP that minimises the number of added reactions subject to the objective floor. Compared with gapFillFastLP, the MILP takes more computation time.

Parameters

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

Name-Value Arguments

minGrowth : double: minimum objective value to achieve. If empty, FBA on the merged model gives the maximum objective value and 10% of it is the floor (default []).
bigM : double: Big-M constant for the coupling constraints of universal reactions with an infinite bound. A finite bound is used as it is (default 1000).
params : struct: solver parameter struct passed to optimizeProb. Fields such as params.TimeLimit and params.intTol tune the MILP. If empty, intTol is 1e-9 and TimeLimit is 600 (default []).
verbose : logical: true to print progress messages (default true).

Returns

addedRxns : cell: reaction IDs from universalModel that are selected by the MILP.
newModel : struct: model with addedRxns incorporated.
exitFlag : double: 1 = the MILP returned a solution, -1 = no gap-filling result: the FBA or MILP is infeasible, the objective of the merged model is not positive, or the solver stops without a solution.

Notes

The MILP is solved with the solver set by setRavenSolver. RAVEN's bundled GLPK solves it for small models. The RAVEN:glpkMILP warning issued with GLPK states that GLPK works for small problems and that Gurobi is much faster and more reliable for larger ones (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 a flux of at least epsilon / 2 in the solution is "consistent" with the core.

The maximisation objective makes the LP degenerate: many optima differ only in which non-core reactions carry flux. The returned active set can therefore vary between LP solvers and between runs with slightly different numerical conditions.

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 that the reaction carries flux in the solution. All elements are false if the LP is infeasible.

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 that needs no solver. For each unreachable metabolite, finds candidate reactions in the universal database that produce it directly.

The step prunes the universal database before a 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).

Name-Value Arguments

seeds : cell: metabolite IDs available from the medium. If empty, the seeds are the metabolites of every exchange reaction whose bounds let it supply the metabolite, that is getExchangeRxns 'uptake' plus 'reverse' (default []).
targets : cell: metabolite IDs that should be produced. If empty, the targets are every metabolite of the objective (biomass) reaction, substrates and products; the first one when there are several objective reactions. Without targets and without an objective reaction the function raises RAVEN:badInput (default []).
verbose : logical: true to print a summary of the results (default true).

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 candidates 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 summarise 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 without any input or consumed without any 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, and the same results are returned as outputs.

Parameters

model : struct: a model structure.

Name-Value Arguments

templateModels : cell: a cell array of template models used 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. Such reactions, for example ATP maintenance, need bounds that allow zero flux.

Examples

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