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

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 whose value can be: * ‘original’: use the original nullspace for internal reactions (Schellenberger et al., 2009) * ‘fastSNP’ : use the minimal feasible nullspace found by Fast-SNP (Saa and Nielson, 2016) * ‘LLC-NS’ : (default): use the minimal feasible nullspace found by solving a MILP (Chan et al., 2017) * ‘LLC-EFM’ : find whether reactions in cycles are connected by EFMs or not

for faster localized loopless constraints (Chan et al., 2017)

Can contain other fields for LLC preprocessing which might be updated during this function

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)

USAGE:

v = checkThermodynamicConsistency(model, q)

INPUTS:

model: q:

OUTPUT:

v:

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 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

• .lb

• .ub

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

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

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

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 metAbbr: metabolite abbreviation solution: solveCobraLP output of a solution to FBA problem 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

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:

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:

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