Thermofba

addLoopLawConstraints(LPproblem, model, rxnIndex, method, reduce_vars, loopInfo)

Adds loop law constraints to LP problem or MILP problem.

USAGE:

[MILPproblem] = addLoopLawConstraints (LPproblem, model, rxnIndex)

INPUTS:

• LPproblem – Structure containing the following fields:

  • .A - LHS matrix

  • .b - RHS vector

  • .c - Objective coeff vector

  • .lb - Lower bound vector

  • .ub - Upper bound vector

  • .osense - Objective sense (-1 max, +1 min)

  • .csense - Constraint senses, a string containting the constraint sense for each row in A (‘E’, equality, ‘G’ greater than, ‘L’ less than).

  • .F - (optional) If *QP problem

  • .vartype - (optional) if MI*P problem

• model – The model for which the loops should be removed, with fields:

  • .S - m x n stoichiometric matrix

  • .lb - n x 1 lower flux bounds

  • .ub - n x 1 upper flux bounds

  • .SIntRxnBool - n x 1 boolean of internal reactions

OPTIONAL INPUT:

• rxnIndex – The index of variables in LPproblem corresponding to fluxes. Default = [1:n]

• method – Indicator which method to use: * 1 - Two variables for each reaction af, ar * 2 - One variable for each reaction af (default)

• reduce_vars – Eliminates additional integer variables. Should be faster in all cases but in practice may not be for some weird reason (default : true).

• loopInfo – Structure containing at least a field named ‘method’, for the method chosen to build loop constraints, with fields:

  • .method - method used to build loop constraints; one of: ‘original’ (original nullspace for internal reactions, Schellenberger et al., 2009), ‘fastSNP’ (minimal feasible nullspace found by Fast-SNP, Saa and Nielson, 2016), ‘LLC-NS’ (default; minimal feasible nullspace found by solving a MILP, Chan et al., 2017) or ‘LLC-EFM’ (reactions in cycles connected by EFMs for faster localized loopless constraints, Chan et al., 2017)

  • .printLevel - verbosity control

  • .N - nullspace matrix for the internal reactions

  • .isInternal - n x 1 logical, true for internal reactions

  • .rxnInLoops - reactions participating in loops

  • .conComp - connected components of the internal-reaction nullspace

  • .useRxnLink - logical, true when .rxnLink is used

  • .rxnLink - reaction linkage matrix derived from EFMs

  • .con - structure of constraint index ranges (.vU, .vL, .gU, .gL)

  • .var - structure of variable index ranges (.z, .g)

  • .rxnInLoopIds - n x 1 ids of the reactions in loops

  • .Mv - big-M value for flux constraints

  • .Mg - big-M value for energy-variable constraints

  • .BDg - default bound for energy variables

OUTPUT:

• MILPproblem – Problem structure containing the following fields describing an MILP problem:

  • A, b, c, lb, ub - same as before but longer

  • vartype - variable type of the MILP problem (‘C’, and ‘B’)

  • x0 = [] - Needed for solveMILPproblem

• loopInfo – Structure containing preprocessing data for using localized loop constraints (LLCs)

checkThermodynamicConsistency(model, q)

Check the thermodynamic consistency of the optimal flux of a model

USAGE:

v = checkThermodynamicConsistency (model, q)

INPUTS:

• model – COBRA model structure with fields:

  • .S - m x n stoichiometric matrix

  • .b - m x 1 right hand side vector

  • .c - n x 1 objective coefficient vector

  • .lb - n x 1 lower flux bounds

  • .ub - n x 1 upper flux bounds

• q – optional n x 1 weighting vector forming the diagonal weighting matrix Q; defaults to the identity

OUTPUT:

v – n x 1 vector S’ * y used to assess thermodynamic consistency

connectedRxnsByEFM(model, conComp, rxnInLoops, printLevel)

Find reactions lying in internal cycles that are connected by any EFMs. Used for minimizing the number of constraints for the loopless requirement when running loopless FVA using localized loopless constraints. This function requires EFMtool (CalculateFluxModes.m) to work.

USAGE:

rxnLink = connectedRxnsByEFM (model, conComp, rxnInLoops)

INPUTS:

• model – COBRA model with field:

  • .S - m x n stoichiometric matrix

• conComp – reactions connected in the minimal nullspace for internal cycles, computed by connectedRxnsInNullSpace

• rxnInLoops – n-by-2 logical matrix where n = # of rxns in the model rxnInLoops(k, 1) = true => forward direction of the k-th reaction in internal cycles rxnInLoops(k, 2) = true => reverse direction of the k-th reaction in internal cycles Returned by findMinNull.m

OPTIONAL INPUT:

printLevel – true to show messages when the linkage matrix cannot be computed

OUTPUT:

rxnLink – n-by-n matrix. rxnLink(i, j) = 1 => reactions i and j are connected by an EFM representing an elementary internal cycle.

connectedRxnsInNullSpace(N)

Find connected components for reactions given a minimal feasible nullspace as defined in Chan et al., Bioinfo, 2017. Two reactions in different connected components imply that no EFM connecting the reactions exists and therefore EFMs can be calculated in a modular approach. If the nullspace matrix represents the nullspace for reactions in internal loops, constraints for the loopless requirement are required only for the connected components involving the reactions required to have no flux through cycles (the target set).

USAGE:

conComp = connectedRxnsInNullSpace (N)

INPUT:

N – a minimal nullspace matrix spanning the feasible flux space, having the same number row as the number of reactions Can be obtained from either findMinNull.m or fastSNP.m

OUTPUT:

conComp – connected components for any reactions connected through the nullspace E.g., conComp = [1; 0; 1; 2; 3; 2] means that the 1st and 3rd reactions are in the same connected component, 4th and 6th also in the same, 5th alone in a connected component and the 2nd reaction is not in any connected component, which means it is a blocked reaction under the condition where the nullspace is calculated.

consistentPotentials(model, printLevel)

Find a consistent set of potentials for each metabolite in a biochemical network, given the directions specified by the bounds on each reaction i.e. find y0, such that \(S^T y_0 < 0\) for a forward reaction and the opposite for reverse.

USAGE:

y0 = consistentPotentials (model, printLevel)

INPUTS:

• model – structure with fields:

  • .S - m x n stoichiometric matrix

  • .lb - n x 1 lower flux bounds

  • .ub - n x 1 upper flux bounds

  • .rxns - n x 1 cell array of reaction identifiers

  • .SIntRxnBool - n x 1 boolean of internal reactions

• printLevel – verbose level

OUTPUT:

y0 – consistent set of chemical potentials

cycleFreeFlux(V0, C, model, SConsistentRxnBool, param)

Removes stoichiometrically balanced cycles from FBA solutions when possible.

A Matlab implementation of the CycleFreeFlux algorithm from Desouki et al., 2015. Minimises the one norm of fluxes subject to bounds determined by input flux.

USAGE:

Vthermo = cycleFreeFlux (V0, C, model, SConsistentRxnBool, relaxBounds)

INPUTS:

• V0 – n x k matrix of k FBA solutions

• C – n x k matrix of k FBA objectives

• model –

  COBRA model structure with required fields:

  • .S - m x n stoichiometric matrix

  • .lb - n x 1 lower bound vector

  • .ub - n x 1 upper bound vector

  • .SIntRxnBool - n x 1 logical array, true for internal reactions

  and optional fields:

  • .b - m x 1 RHS vector (default 0’s)

  • .csense - m x 1 character vector of constraint sense (default ‘E’)

  • .C - coupling constraints

  • .d - RHS of coupling constraints

  • .dsense - coupling constraint sense

  • .SConsistentMetBool - m x 1 logical array, true for stoichiometrically consistent metabolites

  • .SConsistentRxnBool - n x 1 logical array, true for stoichiometrically consistent reactions

OPTIONAL INPUTS:

• SConsistentRxnBool – n x 1 logical array, true for stoichiometrically consistent reactions. Defaults to model.SConsistentRxnBool or model.SIntRxnBool when present.

• param – structure of optional parameters:

  • .printLevel - print level

  • .eta - minimum change in flux considered nonzero. Default feasTol*10

  • .relaxBounds - relax bounds that don’t include zero. Default is false

  • .parallelize - turn parfor use on or off. Default is true if k > 12

  • .enforceCoupling - {(0), 1} where 1 = enforce coupling constraint

  • .approach - {‘lp’, (‘regularised’)} formulation of the cycleFreeFlux problem; ‘lp’ is a linear optimisation based formulation, prone to infeasibility due to numerical issues; ‘regularised’ is a quadratically regularised version, less sensitive to numerical issues but slower

  • .debug - {(0), 1} where 1 = extra debugging steps

  • .removeFixedBool - {(0), 1} where 1 = move variables with equal upper and lower bounds to rhs

  • .steadyStateRelax - relax the steady-state constraint

  • .internalRelax - relax bounds on internal reactions

  • .exchangeRelax - relax bounds on exchange reactions

OUTPUTS:

• Vthermo – n x k matrix of cycle free flux vectors

• thermoConsistentFluxBool – n x 1 logical array, true for thermodynamically consistent flux

Example

% Remove cycles from a single flux vector solution = optimizeCbModel(model); Vthermo = cycleFreeFlux(solution.v, model.c, model);

% Remove cycles from multiple flux vectors [minFlux, maxFlux, Vmin, Vmax] = fluxVariability(model, 0, ‘max’, model.rxns, 0, 1, ‘FBA’); V0 = [Vmin, Vmax]; n = size(model.S, 2); C = [eye(n), eye(n)]; Vthermo = cycleFreeFlux(V0, C, model);

fastSNP(model, varargin)

Generate a minimal feasible basis for all internal cycles using the method Fast-SNP introduced in Saa and Nielson, Bioinformatics, 2016.

USAGE:

N = fastSNP (model, ‘name’, ‘value’, …)

INPUTS:

model – COBRA model with field:

• .S - m x n stoichiometric matrix

OPTIONAL INPUTS:

parameters – solver-specific parameter structure or name-value pair argument for solverCobraLP

OUTPUT:

• N – a minimal feasible basis generating all internal cycles

• resInfo – structure containing the following used parameters/information:

  • .iter - number of iterations

  • .iterTime - time for each iteration

  • .weight - the random weight vector used for finding new basis vector

  • .M - the bound for minimum/maximum flux

  • .feasTol - feasibility tolerance for checking solution feasibility

  • .tol0 - tolerance for zeros in the basis vector

  • .epsilon - tolerance for a new basis vector not

    lying in the projection of the current null-space, i.e. w’(I - P)v >= epsilon or w’(I - P)v <= -epsilon

Siu Hung Joshua Chan 2017 Nov

findMinNull(model, formulation, varargin)

Find a minimal null-space for all internal cycles by solving a MILP, as proposed in Chan et al., 2017.

USAGE:

[rxnInLoops, N, loopInfo] = findMinNull (model, formulation, parameters)

INPUT:

model – COBRA model structure with fields:

• .S - m x n stoichiometric matrix

• .rxns - n x 1 cell array of reaction identifiers

• .lb - n x 1 lower flux bounds

• .ub - n x 1 upper flux bounds

OPTIONAL INPUTS:

• formulation –

  1: solve the MILP problem for a minimal null-space basis

  in an interlaced fashion by presolving some relaxed LPs to simultaneously determine the directions of reactions that participate in internal cycles (default, the quickest way as of 2017 Nov)

  2: directly solve the MILP. Then determine the directions

  of reactions that participate in internal cycles

• parameters – solver-specific parameter structure or name-value pair argument for solverCobraMILP

OUTPUTS:

• rxnInLoops – #rxns-by-2 matrix. rxnInLoops(j, 1) = true => reverse direction of rxn j in loops rxnInLoops(j, 2) = true => forward direction of rxn j in loops

• N – Minimal feasible null-space matrix for internal cycles

• loopInfo – structure containing the following parameters/information:

  • .M - the bound for minimum/maximum flux

  • .minFlux - minimum flux required for a reaction to be active

  • .ignoreRxns - rxns with small coefficients that are prechecked before solving the MILP to avoid numerical issues

  • .nsTime - wall time for finding the null-space

  • .nsCPU - CPU time for finding the null-space

  • .loopPreprocessCPU - CPU time for finding the directions of reactions participating in loops

  • .loopPreprocessTime - wall time for finding the directions of reactions participating in loops

metaboliteMassBalancePlot(model, metAbbr, solution, N)

Plots the top N reactions producing and consuming a metabolite in a flux solution

USAGE:

metaboliteMassBalancePlot (model, metAbbr, solution, N)

INPUTS:

• model – COBRA model structure with fields:

  • .S - m x n stoichiometric matrix

  • .mets - m x 1 cell array of metabolite identifiers

  • .rxns - n x 1 cell array of reaction identifiers

• metAbbr – metabolite abbreviation

• solution – solveCobraLP output of a solution to FBA problem, with field:

  • .full - full flux vector of the solution

• N – Number of reactions to include for production/consumption

printDirectionalityFromBounds(model, lb, ub)

Prints the directionality for each reaction depending on the bounds for each reaction. Defaults to using model.lb & model.ub if none provided

USAGE:

directionality = printDirectionalityFromBounds (model, lb, ub)

INPUT:

model – COBRA model structure with field:

• .S - m x n stoichiometric matrix

OPTIONAL INPUTS:

• lb – flux lower bounds

• ub – flux upper bounds

OUTPUT:

directionality – n x 1 cell array of strings with directionality for each reaction

processingLLCs(process, varargin)

This function contains two different processes for manipulating the localized loop constraints (LLCs) to formulate a model-specific and objective-function-specific MILP for finding loopless flux distributions with a minimal number of binary variables

USAGE:

• 1. Preprocess loop information after calling `addLoopLawConstraints` – [solveLP, MILPproblem, loopInfo] = processingLLCs(‘preprocess’, loopInfo, LPproblem, model, nRxns, osenseStr, MILPproblem)

• 2. Update the loop constraints for a specific objective vector – [solveLP, MILPproblem] = processingLLCs(‘update’, loopInfo, osenseStr, MILPproblem, objVector)

INPUTS:

• process – operation to perform: ‘preprocess’ or ‘update’

• loopInfo – structure containing info about the loops, initially outputed from addLoopLawConstraints and updated by the current function

• LPproblem – original COBRA LP problem structure

• model – the COBRA model from which LPproblem is constructed

• nRxns – number of reactions in the model (default size(model.S, 2))

• osenseStr – optimization sense of the current problem (e.g., FBA max v_biomass, or FVA min v_PYK) ‘max’ (defaulted) or ‘min’

• MILPproblem – COBRA MILP problem generated from LPproblem using addLoopLawConstraints (default generated from LPproblem if not given in the preprocessing call)

• objVector – nRxns-by-1 objective vector for the current optimization problem for updating LLCs (e.g., FVA min v_PYK given fixed v_biomass)

OUTPUTS:

• varargout – comma-separated outputs depending on process; the ‘preprocess’ call returns {solveLP, MILPproblem, loopInfo} and the ‘update’ call returns {solveLP, MILPproblem}.

• solveLP – true if solving LP is sufficient to gaurantee the objective function value is the same as solving MILP with loop constraints

• MILPproblem – updated MILP problem with LLCs

• loopInfo – updated loopInfo with info about the loops (outputed only for the preprocess call)

testCheckThermodynamicConsistency

load textbook.mat