Entropicfba

addExoMetToEFBA(model, exoMet, param)

generates H and h to add a min (alpha/2)(v-h)’*H*(v-h) term to an EFBA problem

USAGE

modelOut = addExoMetToEFBA (model,exoMet,param)

INPUTS

• model.S

• model.rxns

• exoMet.rxns

• exoMet.mean

• exoMet.SD

OPTIONAL INPUT

param.printLevel: param.alpha: alpha in (alpha/2)(v-h)’*H*(v-h), default = 10000; param.metabolomicWeights: String indicating the type of weights to be applied to penalise the difference

between of predicted and experimentally measured fluxes by, where ‘SD’ weights = 1/(1+exoMet.SD^2) ‘mean’ weights = 1/(1+exoMet.mean^2) (Default) ‘RSD’ weights = 1/((exoMet.SD./exoMet.mean)^2)

param.relaxBounds: True to relax bounds on reaction whose fluxes are to be fitted to exometabolomic data

OUTPUTS

modelOut

EXAMPLE:

NOTE:

Author(s):

debugInfeasibleEntropyFBA(model)

try to diagnose the reasons a model is infeasible for Entropy FBA

entropicFluxBalanceAnalysis(model, param)

Entropy maximisation of fluxes (or fluxes and concentrations) subject to mass balance, optionally coupling constraints, optionally quadratic penalisation of deviation from given fluxes.

minimize g.*vf’*(log(vf) -1) + (cf + ci)’*vf vf,vr,w,x,x0 + g.*vr’*(log(vr) -1) + (cr - ci)’*vr

• f.*x’ *(log(x) -1) + u0’*x

• f.*x0’*(log(x0) -1) + u0’*x0

• ce’*w

• (1/2)v’*Q*v

• (1/2)(v-h)’*H*(v-h)

subject to [N B]*[v w] <=> b : y_N subject to N*(vf - vr) + B*w = x - x0 = dx/dt

subject to N*(vf - vr) - x + x0 <=> by_N

C*(vf - vr) <=> dy_C

lb <= [vf - vr; w] <= ubz_v

dxl <= x - x0 <= dxu : z_dx vfl <= vf <= vfu : z_vf vrl <= vr <= vru : z_vr

xl <= x <= xu : z_x

x0l <= x0 <= x0u : z_x0

with Biochemical optimality conditions

|| N*(vf - vr) - x + x0 - b ||_inf || C*(vf - vr) - d ||_inf || g*log(vf) + ci + cf + N’*y_N + C’*y_C + z_v + z_vf ||_inf || g*log(vr) - ci + cr - N’*y_N - C’*y_C - z_v + z_vr ||_inf || f.*log(x) + u0 - y_N + z_dx - z_x ||_inf || f.*log(x0) + u0 + y_N - z_dx + z_x0 ||_inf

with Derived biochemical optimality conditions (fluxes) || g*log(vr/vf) + cr - cf - 2*(ci + N’*y_N + C’*y_C + z_vi) + z_vr - z_vf ||_inf where z_vi is the subset of z_v corresponding to internal reactions

with Derived biochemical optimality conditions (concentrations) || f.*log(x/x0) - 2*y_N + 2*z_dx + z_x - z_x0 ||_inf || f.*log(x.*x0) + 2*u0 + z_x + z_x0 ||_inf

Derived biochemical optimality conditions (fluxes and concentrations) || g*log(vf) + cf + ci + N’*(u0 + log(x) + z_dx + z_x) + C’*y_C + z_vi + z_vf ||_inf || g*log(vr) + cr - ci - N’*(u0 + log(x) + z_dx + z_x) - C’*y_C - z_vi + z_vr ||_inf

Derived biochemical optimality conditions (fluxes and concentrations, combining forward and reverse) || g*log(vr/vf) + cr - cf - 2*(ci + N’*(u0 + f*log(x) + z_dx + z_x) + C’*y_C + z_vi) - z_vf + z_vr ||_inf || g*log(vr/vf) + cr - cf - 2*(ci - N’*(u0 + f*log(x0) - z_dx + z_x0) + C’*y_C + z_vi) - z_vf + z_vr ||_inf

If (but not only if) the input data is as follows: g = 2, f = 1, cr = cf, ci = 0, C = 0, d = 0, <=> y_C = 0 vl = -inf, vu = inf <=> z_v = 0 (for internal reactions) dxl = -inf, dxu = inf, <=> z_dx = 0 ub = inf, <=> z_vf = 0 lb = - inf, <=> z_vr = 0 x0l = -inf, x0u = inf, <=> z_x0 = 0 then the above reduces to || log(vr/vf) = N’*(u0 + log(x) + z_x) ||_inf where z_x is the dual variable to the bounds on concentration x.

USAGE

[solution, modelOut] = entropicFluxBalanceAnalysis (model,param)

INPUT

model – (the following fields are required - others can be supplied)

• S - m x (n + k) Stoichiometric matrix

• c - (n + k) x 1 Linear objective coefficients

• lb - (n + k) x 1 Lower bounds on net flux

• ub - (n + k) x 1 Upper bounds on net flux

OPTIONAL INPUTS: model.osenseStr: Maximize (‘max’)/minimize (‘min’) (opt, default = ‘max’) linear part of the objective.

Nonlinear parts of the objective are always assumed to be minimised.

model.b m x 1 change in concentration with time model.csense m x 1 character array with entries in {L,E,G}

model.C: c x (n + k) Left hand side of C*v <= d model.d: c x (n + k) Right hand side of C*v <= d model.dsense c x 1 character array with entries in {L,E,G}

model.g n x 1 strictly positive weight on internal flux entropy maximisation (default 2) model.cf: n x 1 real valued linear objective coefficients on internal forward flux (default 0) model.cr: n x 1 real valued linear objective coefficients on internal reverse flux (default 0) model.vfl: n x 1 non-negative lower bound on internal forward flux (default 0) model.vfu: n x 1 non-negative upper bound on internal forward flux (default inf) model.vrl: n x 1 non-negative lower bound on internal reverse flux (default 0) model.vru: n x 1 non-negative upper bound on internal reverse flux (default 0)

model.f: m x 1 strictly positive weight on concentration entropy maximisation (default 1) model.u0: m x 1 real valued linear objective coefficients on concentrations (default 0) model.x0l: m x 1 non-negative lower bound on initial molecular concentrations model.x0u: m x 1 non-negative upper bound on initial molecular concentrations model.xl: m x 1 non-negative lower bound on final molecular concentrations model.xu: m x 1 non-negative lower bound on final molecular concentrations model.dxl: m x 1 real valued lower bound on difference between final and initial molecular concentrations dxl <= x - x0 model.dxu: m x 1 real valued upper bound on difference between final and initial initial molecular concentrations x - x0 <= dxu

model.Q (n + k) x (n + k) positive semi-definite matrix to minimise (1/2)v’*Q*v

model.H (n + k) x (n + k) positive semi-definite matrix in objective (1/2)(v-h)’*H*(v-h) to be minimised model.h (n + k) x 1 vector in objective (1/2)(v-h)’*H*(v-h) to be minimised

model.SConsistentMetBool: m x 1 boolean indicating stoichiometrically consistent metabolites model.SConsistentRxnBool: n x 1 boolean indicating stoichiometrically consistent metabolites

param.solver: {(‘pdco’),’mosek’} param.entropicFBAMethod: {(‘fluxes’),’fluxConc’)} maximise entropy of fluxes or also concentrations param.printLevel: {(0),1}

Parameters related with flux optimisation

param.maxUnidirectionalFlux: scalar real valued maximum expected value of unidirectional flux param.internalNetFluxBounds: ‘original’ (default) maintains direction and magnitude of net flux from model.lb & model.ub

‘directional’ maintains direction of net flux from model.lb & model.ub but not magnitude ‘random’ random net flux direction, replacing constraints from model.lb & model.ub

Parameters related with concentration optimisation: param.maxConc: scalar maximum permitted metabolite concentration param.externalNetFluxBounds: (‘original’) = use bounds on external reactions from model.lb and model.ub

‘dxReplacement’ = replace model.lb and model.ub with bounds on change in concentration ‘none’ = set exchange reactions to be unbounded

model.gasConstant: scalar gas constant (default 8.31446261815324 J K^-1 mol^-1) model.T: scalar temperature (default 310.15 Kelvin)

OUTPUTS: solution: solution structure with the following fields

*.v: n x 1 double net flux *.vf: n x 1 double unidirectional forward internal reaction flux *.vr: n x 1 double unidirectional reverse internal reaction flux *.vt: scalar total internal reaction flux sum(vf + vr) *.y_N: m × 1 double dual variable to steady state constraints *.y_C: z × 1 double dual variable to coupling constraints *.z_v: (n + k) x 1 double dual variable to box constraints on net flux *.z_vf: n x 1 double dual variable to box constraints on forward flux *.z_vr: n x 1 double dual variable to box constraints on reverse flux *.time: solve time *.stat: COBRA toolbox standard solution status *.origStat: solution status as provided by the solver

modelOut: solved model with optional input fields populated by defaults, if they were not provided

EXAMPLE:

NOTE:

Author(s): Ronan M.T. Fleming 2021

processConcConstraints(model, param)

USAGE

[] = processConcConstraints (model,param)

INPUTS

model – (the following fields are required - others can be supplied)

• S - m x (n + k) Stoichiometric matrix

• c - (n + k) x 1 Linear objective coefficients

• lb - (n + k) x 1 Lower bounds on net flux

• ub - (n + k) x 1 Upper bounds on net flux

OPTIONAL INPUTS model.SConsistentMetBool: m x 1 boolean indicating stoichiometrically consistent metabolites model.SConsistentRxnBool: n x 1 boolean indicating stoichiometrically consistent metabolites model.rxns:

model.f: m x 1 strictly positive weight on concentration entropy maximisation (default 1) model.u0: m x 1 standard transformed Gibbs energy of formation (default 0) model.c0l: m x 1 non-negative lower bound on initial molecular concentrations model.c0u: m x 1 non-negative upper bound on initial molecular concentrations model.cl: m x 1 non-negative lower bound on final molecular concentrations model.cu: m x 1 non-negative lower bound on final molecular concentrations model.dcl: m x 1 real valued lower bound on difference between final and initial molecular concentrations (default -inf) model.dcu: m x 1 real valued upper bound on difference between final and initial initial molecular concentrations (default inf) model.gasConstant: scalar gas constant (default 8.31446261815324 J K^-1 mol^-1) model.temperature: scalar temperature (default 310.15 Kelvin) param.maxConc: (1e4) maximim micromolar concentration allowed param.maxConc: (1e-4) minimum micromolar concentration allowed param.externalNetFluxBounds: (‘original’) =

‘dxReplacement’ = when model.dcl or model.dcu is provided then they set the upper and lower bounds on metabolite exchange

param.printLevel:

OUTPUTS: f: m x 1 strictly positive weight on concentration entropy maximisation (default 1) u0: m x 1 standard transformed Gibbs energy of formation, divided by RT (default 0) c0l: m x 1 non-negative lower bound on initial molecular concentrations c0u: m x 1 non-negative upper bound on initial molecular concentrations cl: m x 1 non-negative lower bound on final molecular concentrations cu: m x 1 non-negative lower bound on final molecular concentrations dcl: m x 1 real valued lower bound on difference between final and initial molecular concentrations dcu: m x 1 real valued upper bound on difference between final and initial initial molecular concentrations

wl: k x 1 lower bound on external net flux wu: k x 1 upper bound on external net flux

B: m x k External stoichiometric matrix b: m x 1 RHS of S*v = b

EXAMPLE:

NOTE:

Author(s): Ronan Fleming

processFluxConstraints(model, param)

USAGE

• processFluxConstraints (model)

• processFluxConstraints (model,param)

INPUTS

model – (the following fields are required - others can be supplied)

• S - m x (n + k) Stoichiometric matrix

• c - (n + k) x 1 Linear objective coefficients

• lb - (n + k) x 1 Lower bounds on net flux

• ub - (n + k) x 1 Upper bounds on net flux

OPTIONAL INPUTS

model.osenseStr: (‘max’) or ‘min’. The default linear objective sense is maximisation, i.e. ‘max’ model.cf: n x 1 real valued linear objective coefficients on internal forward flux (default 0) model.cr: n x 1 real valued linear objective coefficients on internal reverse flux (default 0) model.g n x 1 strictly positive weight on internal flux entropy maximisation (default 2) model.SConsistentRxnBool: n x 1 boolean indicating stoichiometrically consistent metabolites model.vfl: n x 1 non-negative lower bound on internal forward flux (default 0) model.vfu: n x 1 non-negative upper bound on internal forward flux (default inf) model.vrl: n x 1 non-negative lower bound on internal reverse flux (default 0) model.vru: n x 1 non-negative upper bound on internal reverse flux (default 0)

param.printLevel: param.solver: {‘pdco’,(‘mosek’)} param.debug: {(0),1} 1 = run in debug mode param.entropicFBAMethod: {(‘fluxes’),’fluxesConcentrations’} maximise entropy of fluxes (default) or also concentrations param.maxUnidirectionalFlux: maximum unidirectional flux (1e5 by default) param.minUnidirectionalFlux: minimum unidirectional flux (zero by default) param.internalNetFluxBounds: (‘original’) = use model.lb and model.ub to set the direction and magnitude of internal net flux bounds

‘directional’ = use model.lb and model.ub to set the direction of net flux bounds (ignoring magnitude) ‘none’ = ignore model.lb and model.ub and allow all net flues to be reversible

OUTPUTS: vl: n x 1 lower bound on internal net flux vu: n x 1 upper bound on internal net flux vel: k x 1 lower bound on external net flux veu: k x 1 upper bound on external net flux vfl: n x 1 non-negative lower bound on internal forward flux vfu: n x 1 non-negative upper bound on internal forward flux vrl: n x 1 non-negative lower bound on internal reverse flux vru: n x 1 non-negative upper bound on internal reverse flux ci: n x 1 linear objective coefficients corresponding to internal net fluxes vel: k x 1 linear objective coefficients corresponding to external net fluxes cf: n x 1 real valued linear objective coefficients on internal forward flux cr: n x 1 real valued linear objective coefficients on internal reverse flux g n x 1 strictly positive weight on internal flux entropy maximisation

EXAMPLE:

NOTE:

Author(s): Ronan Fleming

solveCobraEP(EPproblem, varargin)

Solves the following optimisation problem:

minimize osense*(c.*d)’x + d.*x’(log(x) -1) + (1/2)*x’*Q*x

x

subject to A*x <=> by

lb <= x <= ub : z

or subject to blc <= A*x <= buc : y

lb <= x <= ub : z

However, when EPproblem.P is present, the following optimisation problem is solved:

minimize osense*(c.*d)’*x + (d.*x)’*log(x./q) + (1/2)*x’*Q*x = osense*(c.*d)’*x + (d.*x)’*log(x) - (d.*x)’*log(q) + (1/2)*x’*Q*x

x,q

subject to A*x <=> by

P*x - q = 0 : r lb <= x <= ub : z

or subject to blc <= A*x <= buc : y

P*x -q = 0 : r lb <= x <= ub : z

USAGE

sol = solveCobraEP (EPproblem, varargin)

INPUT

EPproblem –

Structure containing the following fields describing the EP problem to be solved

• .A - m x n Linear constraint matrix

• .c - n x 1 Linear objective coeff vector

• .lb - n x 1 Lower bound vector

• .ub - n x 1 Upper bound vector

• .d - n x 1 Non-negative vector indicating the non-negative

  variables whose entropy is maximised. If d(i)==0 then there is only a linear objective on x(i).

• .osense - Linear objective sense (-1 means maximise, 1 means minimise)

With either the following fields

• .b - m x 1 right hand side vector i.e. A <=> b

• .csense - m x 1 string containting the constraint sense for

  each row in A (‘E’, equality, ‘G’ greater than, ‘L’ less than).

Or with the following fields

• .blc - m x 1 left hand side vector i.e. blc <= A*x

• .buc - m x 1 right hand side vector i.e. A*x <= buc

OPTIONAL INPUTS

• EPproblem – Structure containing the following fields describing the EP problem to be solved

  • .P - p x n matrix with entries {0,1}, such that q_i := P_{i,:}*x is the sums

    of the x corresponding to nonzero columns of the ith row of P, i.e. P_{i,:}. Used for normalised entropy maximisation.

  • .Q - positive semidefinite matrix for quadratic part of objective (see above)

  • varargin –

    Additional parameters either as parameter struct, or as

    parameter/value pairs. A combination is possible, if the parameter struct is either at the beginning or the end of the optional input. All fields of the struct which are not COBRA parameters (see getCobraSolverParamsOptionsForType) for this problem type will be passed on to the solver in a solver specific manner. Some optional parameters which can be passed to the function as parameter value pairs, or as part of the options struct are listed below:

    ‘verify’ ‘printLevel’ ‘debug’ ‘feasTol’ ‘optTol’ ‘solver’

  • printLevel – Printing level

    • 0 - Silent (Default)

    • 1 - Warnings and Errors

    • 2 - Summary information

    • 3 - More detailed information

    • > 10 - Pause statements, and maximal printing (debug mode)

solver: Optimisation solver used, {(‘mosek’),’pdco’} feasTol: Feasibility tolerance optTol: Optimality tolerance

OUTPUT

sol – Structure containing the following fields describing a LP sol: *.obj: Objective value *.objLinear osense*c’*x; *.objEntropy d.*x’*(log(x) -1); *.objQuadratic (1/2)*x’*Q*x; *.v: n+k ×1 double *.vf: n × 1 double *.vr: n × 1 double *.vt: 1’*vt + 1’*vr *.y_N: m x 1 double dual sol to constraints :math: A*x (‘E’ | ‘G’ | ‘L’) b *.z_dx: 0 *.z_vf: n × 1 double dual sol to \(lb <= vr <= ub\) *.z_vr: n × 1 double dual sol to \(lb <= vf <= ub\) *.z_vi: n × 1 double dual sol to \(lb <= v <= ub\) *.z_v: n + k × 1 double dual sol to \(lb <= w <= ub\) * .stat: Solver status in standardized form

• 0 - Infeasible problem

• 1 - Optimal sol

• 2 - Unbounded sol

• 3 - Almost optimal sol

• -1 - Some other problem (timelimit, numerical problem etc)

• .origStat: Original status returned by the specific solver

• .origStatText: Original status text returned by the specific solver

• .time: Solve time in seconds

• .solver: Solver used to solve EP problem

• .epmethod: solver method used e.g. ‘CONIC’