Entropicfba

addExoMetToEFBA(model, exoMet, param)[source]

Generates .H and .h to add a min (alpha/2)(v-h)’*H*(v-h) term to an entropic FBA (EFBA) problem, fitting predicted fluxes to exometabolomic data

USAGE:

modelOut = addExoMetToEFBA (model, exoMet, param)

INPUTS:

  • model – COBRA model structure with fields:

    • .S - m x n stoichiometric matrix

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

  • exoMet – Structure of exometabolomic flux data with fields:

    • .rxns - k x 1 cell array of reaction identifiers, matched against model.rxns

    • .mean - k x 1 mean experimentally measured flux for each reaction in exoMet.rxns

    • .SD - k x 1 standard deviation of the experimentally measured flux for each reaction in exoMet.rxns

OPTIONAL INPUT:

param – Structure with fields:

  • .printLevel - verbose level (default: 0)

  • .alpha - alpha in (alpha/2)(v-h)’*H*(v-h) (default: 10000)

  • .metabolomicWeights - string indicating the type of weights to be applied to penalise the difference between 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)

  • .relaxBounds - true to relax the bounds on the reactions whose fluxes are to be fitted to exometabolomic data (default: true)

OUTPUTS:

modelOut – model, with .H (quadratic penalty matrix) and .h (linear penalty target) added, and with .lb and .ub relaxed on the fitted reactions when param.relaxBounds is true

debugInfeasibleEntropyFBA(model)[source]

Tries to diagnose the reasons a model is infeasible for entropic flux balance analysis (EFBA) by progressively relaxing coupling constraints, internal flux bounds, and exchange bounds until a feasible EFBA solution is found

USAGE:

[sol, modelOut] = debugInfeasibleEntropyFBA (model)

INPUT:

model – COBRA model structure, compatible with entropicFluxBalanceAnalysis. Diagnosis progressively removes/relaxes its optional coupling constraint fields (.C, .d) and rescales/relaxes .lb, .ub, and .b on the reactions not in .SConsistentRxnBool if needed

OUTPUTS:

  • sol – Solution structure returned by the last entropicFluxBalanceAnalysis call attempted (see entropicFluxBalanceAnalysis for fields), with .stat equal to 1 if that attempt was feasible

  • modelOut – The (possibly relaxed) model corresponding to sol; only assigned once a relaxation is found that yields a feasible solution

entropicFluxBalanceAnalysis(model, param)[source]

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.

Optional enzyme-constrained (GECKO/ecModel) extension

When the enzyme-usage column variables are present (model.E, and optionally model.D, with bounds model.evarlb/evarub and linear objective coefficients model.evarc; see the field documentation below), nEvar additional variables e are appended AFTER all flux/concentration columns, matching the buildOptProblemFromModel [S E; C D] block ordering: model.E occupies the metabolite (N) rows and model.D the coupling (C) rows. The problem above then extends to

minimize (all terms as above) vf,vr,w,x,x0,e + evarc’*e

  • ge.*e’*(log(e) -1)

where ge = param.enzymeEntropyWeight (default 0 => the enzyme columns are LINEAR additional variables and the entropy term drops out; a positive ge adds an experimental entropy term on the enzyme columns).

subject to N*(vf - vr) - x + x0 + E*e <=> by_N
C*(vf - vr) + D*e <=> dy_C

evarlb <= e <= evarub : z_e

with the additional enzyme Biochemical optimality condition

|| evarc + E’*y_N + D’*y_C + z_e ||_inf

where z_e is the reduced cost of the enzyme-usage variables (dual to the evarlb <= e <= evarub bounds), returned as solution.z_e, and e is returned as solution.e. The +z_e sign matches the reduced-cost convention of both the mosek and pdco backends (the same sign as the external-reaction condition || ce + B’*y_N + z_ve ||_inf). When the enzyme fields are absent (model.E empty) the formulation reduces exactly to the one above.

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, split into

    internal and external as follows:

    ci: n x 1 linear objective coefficients corresponding to internal net fluxes ce: k x 1 linear objective coefficients corresponding to internal net fluxes

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

set H(j,j) = finite quadratic penalty, for reaction j OR set h(j) = 0 if no quadratic penalty, for reaction j

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

set h(j) = finite desired net flux value, for reaction j OR set h(j) = NaN if no desired net flux value exists for reaction j

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

Optional enzyme-constrained (GECKO/ecModel) column variables (feature 010-gecko-entropic-fba). When these fields are present, the [S E; C D] block and the enzyme-usage bounds/objective are folded into the entropic problem so that enzyme-constrained (GECKO/ecModel) models can be solved by entropic FBA under param.solver ‘mosek’ or ‘pdco’. When absent, the function behaves exactly as it does without them (the enzyme fields are the same ones understood by buildOptProblemFromModel). model.E: m x nEvar enzyme-usage columns in the metabolite (S) rows (default: absent) model.D: c x nEvar enzyme-usage columns in the coupling (C) rows model.evarlb: nEvar x 1 lower bounds on the enzyme-usage variables model.evarub: nEvar x 1 upper bounds on the enzyme-usage variables model.evarc: nEvar x 1 linear objective coefficients on the enzyme-usage variables model.evars: nEvar x 1 cell array of enzyme-usage variable names

param: optional parameter structure controlling the solver and

algorithm, with fields:

param.solver: {(‘pdco’),’mosek’} param.entropicFBAMethod: {(‘fluxes’),’fluxConc’)} maximise entropy of fluxes or also concentrations param.printLevel: {(0),1} param.debug: (default false) boolean flag reserved for internal debugging branches param.entropicMethod: legacy input name for param.entropicFBAMethod, retained for backward compatibility param.enzymeEntropyWeight: (default 0) scalar or nEvar x 1 weight on the entropy of the

enzyme-usage variables (the model.E/model.D columns). 0 => the enzyme columns are LINEAR additional variables (default); a positive value adds an experimental entropy term on the enzyme columns (its enzyme-dual correctness is not yet validated - use with caution).

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 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 x 1 double dual variable to steady state constraints

    • .y_C - z x 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

    • .e - nEvar x 1 double enzyme-usage variable values (only when model.E is present)

    • .z_e - nEvar x 1 double reduced cost of the enzyme-usage variables (only when model.E is present)

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

EXAMPLE:

NOTE:

mosekParamSetEFBA(param)[source]

Creates a structure of pertinent user-defined options for the MOSEK solver, adding/overriding the MOSEK-specific integer and double parameters used by entropic flux balance analysis

USAGE:

param = mosekParamSetEFBA (param)

INPUT:

param – Parameter structure to be passed to the MOSEK solver; the following MOSEK-specific fields are set (any other fields already present in param are preserved):

  • .MSK_IPAR_LOG_PRESOLVE - amount of output printed by the presolve procedure; a higher value logs more information (MOSEK default: 1; set here to 1)

  • .MSK_IPAR_INTPNT_SCALING - how the problem is scaled before the interior-point optimizer is used (MSK_SCALING_NONE, MSK_SCALING_MODERATE, MSK_SCALING_AGGRESSIVE, or MSK_SCALING_FREE; MOSEK default is MSK_SCALING_FREE; set here to MSK_SCALING_FREE)

  • .MSK_IPAR_INTPNT_REGULARIZATION_USE - whether regularization is allowed (MSK_ON/MSK_OFF; MOSEK default is MSK_ON; set here to MSK_OFF)

  • .MSK_DPAR_INTPNT_CO_TOL_DFEAS - dual feasibility tolerance used by the interior-point optimizer for conic problems (MOSEK default: 1.0e-8; set here to 1.0e-11)

  • .MSK_DPAR_INTPNT_CO_TOL_PFEAS - primal feasibility tolerance used by the interior-point optimizer for conic problems (MOSEK default: 1.0e-8; set here to 1.0e-11)

  • .MSK_DPAR_INTPNT_CO_TOL_REL_GAP - relative gap termination tolerance used by the interior-point optimizer for conic problems (MOSEK default: 1.0e-8; set here to 1.0e-11; useful for ensuring dual feasibility is as good as primal)

  • .MSK_IPAR_INTPNT_MAX_ITERATIONS - maximum number of iterations allowed in the interior-point optimizer (MOSEK default: 400; set here to 400)

  • .MSK_IPAR_BI_IGNORE_MAX_ITER - if MSK_IPAR_INTPNT_BASIS is MSK_BI_NO_ERROR and the interior-point optimizer terminated on the maximum number of iterations, whether basis identification is still performed (MSK_ON/MSK_OFF; MOSEK default is MSK_OFF; set here to MSK_OFF)

  • .MSK_IPAR_INTPNT_SOLVE_FORM - whether the primal or the dual problem is solved (MSK_SOLVE_PRIMAL, MSK_SOLVE_DUAL, or MSK_SOLVE_FREE; MOSEK default is MSK_SOLVE_FREE; set here to MSK_SOLVE_FREE)

  • .MSK_DPAR_INTPNT_TOL_INFEAS - controls when the optimizer declares the model primal or dual infeasible; a smaller value makes the optimizer more conservative about declaring infeasibility (MOSEK default: 1.0e-8; set here to 1e-8)

  • .MSK_IPAR_LOG_INTPNT - amount of output printed by the interior-point optimizer; a higher value logs more information (MOSEK default: 4; set here to 5)

  • .MSK_IPAR_INFEAS_REPORT_AUTO - amount of information presented in an automatically generated infeasibility report (MSK_ON/MSK_OFF; MOSEK default is MSK_OFF; set here to MSK_OFF)

  • .MSK_IPAR_INFEAS_REPORT_LEVEL - amount of information presented in an infeasibility report; higher values give more information (MOSEK default: 1; set here to 100)

OUTPUT:

param – param, with the MOSEK-specific fields listed under INPUT set/overridden

prepareEnzymeConstrainedEP(EPproblem, model, m, n, nCoupling, enzymeEntropyWeight)[source]

Fold optional enzyme-constrained (GECKO) column variables into an assembled entropic problem, for entropicFluxBalanceAnalysis.

USAGE:

[EPproblem, nEvar] = prepareEnzymeConstrainedEP (EPproblem, model, m, n, nCoupling, enzymeEntropyWeight)

INPUTS:

  • EPproblem – partially-assembled entropic problem with fields .A (rows ordered metabolite(m), internal-net(n), coupling(nCoupling), other), .c, .lb, .ub, .d

  • model – COBRA model carrying the optional column-variable fields model.E (enzyme columns in the metabolite rows), model.D (enzyme columns in the coupling rows), model.evarlb, model.evarub, model.evarc

  • m – number of metabolite (S) rows

  • n – number of internal (SConsistent) reactions

  • nCoupling – number of coupling (model.C) rows (0 if none)

OPTIONAL INPUTS:

enzymeEntropyWeight – entropy weight applied to the enzyme-usage columns: 0 or [] (default) => linear-only enzyme variables; a positive scalar or nEvar-vector => enzyme columns carry an entropy term (added to EPproblem.d)

OUTPUTS:

  • EPproblem – the problem with nEvar enzyme columns appended to .A/.c/.lb/.ub/.d

  • nEvar – number of enzyme-usage variables folded in

Note

The enzyme columns are appended AFTER all existing columns (matching the buildOptProblemFromModel [S E; C D] ordering): model.E occupies the metabolite rows, model.D the coupling rows, zeros elsewhere. When enzymeEntropyWeight is 0/[] the enzyme variables are linear-only (their .d entries are 0), so nnz(.d) and the exponential-cone count are unchanged.

processConcConstraints(model, param)[source]

Derives the concentration-related bounds, weights, and external stoichiometry used to build an entropic flux balance analysis (EFBA) concentration subproblem from a COBRA model and EFBA parameter structure

USAGE:

[f, u0, c0l, c0u, l_c, u_c, dcl, dcu, l_w, u_w, B, b, l_r, u_r, paramOut] = processConcConstraints (model, param)

INPUTS:

  • model – COBRA model structure with the following required fields:

    • .S - m x (n + k) stoichiometric matrix, where n is the number of stoichiometrically consistent (internal) reactions and k is the number of stoichiometrically inconsistent (external) reactions

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

    • .SConsistentRxnBool - (n + k) x 1 boolean indicating the stoichiometrically consistent (internal) reactions

    and the following optional fields:

    • .rxns - (n + k) x 1 cell array of reaction identifiers, used only to report inconsistent reactions

    • .pp - m x 1 boolean vector indicating the independent rows of [N, b] (default: all rows); also used to restrict b, B, l_r, and u_r to the independent rows

    • .b - m x 1 right hand side of N*v = b (default: 0)

    • .f - m x 1 strictly positive weight on concentration entropy maximisation, or ‘rand’, ‘one’, or ‘two’ (default ‘one’)

    • .u0 - m x 1 standard transformed Gibbs energy of formation, or ‘rand’, ‘one’, or ‘zero’ (default ‘zero’)

    • .c0l - m x 1 non-negative lower bound on initial molecular concentrations (default 0)

    • .c0u - m x 1 non-negative upper bound on initial molecular concentrations (default param.maxConc)

    • .l_c - m x 1 non-negative lower bound on final molecular concentrations

    • .u_c - m x 1 non-negative upper bound on final molecular concentrations

    • .dcl - m x 1 real valued lower bound on the difference between final and initial molecular concentrations (default -inf)

    • .dcu - m x 1 real valued upper bound on the difference between final and initial molecular concentrations (default inf)

    • .gasConstant - gas constant in kJ/(K*mol) (default 8.3144621e-3)

    • .T - scalar temperature; takes precedence over .temperature if both are present

    • .temperature - scalar temperature in Kelvin (default 310.15)

  • param – Structure with the following optional fields:

    • .concUnit - concentration unit conversion factor (default 10-3, i.e. assumes concentrations are in mMol)

    • .concentrationBounds - {(‘none’)} whether/how to set bounds on final concentration: ‘none’ (l_c = 0, u_c = inf), ‘setToGiven’ (use model.l_c/model.u_c, selected automatically when both are present), or ‘maximimumFiniteRange’ (use param.minConc/param.maxConc)

    • .maxConc - (1e4) maximum micromolar concentration allowed

    • .minConc - (1e-4) minimum micromolar concentration allowed

    • .externalNetFluxBounds - (‘original’) how to set the bounds on external net flux and the initial/final concentration difference: ‘original’ (use model.lb/model.ub unmodified), ‘identities’, ‘bReplacingB’, ‘none’, or ‘dxReplacement’ (when model.dcl or model.dcu is provided, use them to set the exchange bounds); defaults to ‘dxReplacement’ if model.dcl or model.dcu is present, else ‘original’

    • .printLevel - verbose level, used to print diagnostic messages about the option chosen for .externalNetFluxBounds

    • .qpMassBalance - backward-compatible alias; copied to .massBalancePenalty if .massBalancePenalty is absent

    • .strictMassBalance - if present, its logical negation overwrites .qpMassBalance

    • .massBalancePenalty - (‘none’) penalty applied to the mass-balance regularisation term: ‘quadratic’ (l_r/u_r unbounded) or ‘none’ (l_r = u_r = 0)

OUTPUTS:

  • f – m x 1 strictly positive weight on concentration entropy maximisation

  • u0 – m x 1 standard transformed Gibbs energy of formation, divided by gasConstant*temperature if both are given

  • c0l – m x 1 non-negative lower bound on initial molecular concentrations

  • c0u – m x 1 non-negative upper bound on initial molecular concentrations

  • l_c – m x 1 non-negative lower bound on final molecular concentrations

  • u_c – m x 1 non-negative upper bound on final molecular concentrations

  • dcl – m x 1 real valued lower bound on the difference between final and initial molecular concentrations

  • dcu – m x 1 real valued upper bound on the difference between final and initial molecular concentrations

  • l_w – k x 1 lower bound on external net flux

  • u_w – k x 1 upper bound on external net flux

  • B – m x k external stoichiometric matrix (rows restricted to model.pp when present)

  • b – right hand side of N*v = b (rows restricted to model.pp when present)

  • l_r – lower bound on the regularisation term in N*v + r = b (default -inf)

  • u_r – upper bound on the regularisation term in N*v + r = b (default inf)

  • paramOut – param structure, returned to capture any fields set internally (e.g. defaulted fields)

processFluxConstraints(model, param)[source]

Derives the internal and external net-flux bounds, unidirectional flux bounds, and objective/entropy weights used to build an entropic flux balance analysis (EFBA) flux subproblem from a COBRA model and EFBA parameter structure

USAGE:

[vl, vu, l_w, u_w, vfl, vfu, vrl, vru, ci, ce, cf, cr, g, paramOut] = processFluxConstraints (model, param)

INPUTS:

model – COBRA model structure with the following required fields:

  • .S - m x (n + k) stoichiometric matrix, where n is the number of stoichiometrically consistent (internal) reactions and k is the number of stoichiometrically inconsistent (external) reactions

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

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

and the following optional fields:

  • .SConsistentRxnBool - (n + k) x 1 boolean indicating the stoichiometrically consistent (internal) reactions; together with .SConsistentMetBool, computed via findStoichConsistentSubset if either is absent

  • .SConsistentMetBool - m x 1 boolean indicating the stoichiometrically consistent metabolites (see .SConsistentRxnBool)

  • .c - (n + k) x 1 linear objective coefficients (default: all zero, giving ci = ce = 0)

  • .osenseStr - (‘max’) or ‘min’, the linear objective sense

  • .rxns - (n + k) x 1 cell array of reaction identifiers, used only to report reactions ignored in the ‘identities’ external-bound option

  • .dcl - m x 1 real valued lower bound on the difference between final and initial molecular concentrations, read when param.externalNetFluxBounds is ‘original’ or ‘dxReplacement’

  • .dcu - m x 1 real valued upper bound on the difference between final and initial molecular concentrations (see .dcl)

  • .vfl - n x 1 non-negative lower bound on internal forward flux (default: derived from param.minUnidirectionalFlux and vl)

  • .vfu - n x 1 non-negative upper bound on internal forward flux (default: param.maxUnidirectionalFlux)

  • .vrl - n x 1 non-negative lower bound on internal reverse flux (default: derived from param.minUnidirectionalFlux and vu)

  • .vru - n x 1 non-negative upper bound on internal reverse flux (default: param.maxUnidirectionalFlux)

  • .cf - n x 1 real valued linear objective coefficients on internal forward flux, or ‘rand’, ‘one’, or ‘zero’ (default ‘zero’; if absent, taken from .c_vf)

  • .cr - n x 1 real valued linear objective coefficients on internal reverse flux, or ‘rand’, ‘one’, or ‘zero’ (default ‘zero’; if absent, taken from .c_vr)

  • .c_vf - fallback source for .cf if .cf is absent

  • .c_vr - fallback source for .cr if .cr is absent

  • .g - n x 1 strictly positive weight on internal flux entropy maximisation, or ‘zero’, ‘rand’, ‘one’, or ‘two’ (default ‘one’ if param.entropicFBAMethod is ‘fluxes’, else ‘two’)

OPTIONAL INPUT:

param – Structure with the following optional fields:

  • .printLevel - verbose level

  • .debug - {(0), 1} 1 = run in debug mode, printing and validating problematic vfl > vfu bounds using .feasTol

  • .feasTol - feasibility tolerance; read only when .debug is true, to set vfu = vfl + feasTol where vfl > vfu

  • .entropicFBAMethod - {(‘fluxes’), ‘fluxesConcentrations’} maximise entropy of fluxes (default) or also concentrations; if absent, taken from the legacy .method field

  • .method - legacy alias for .entropicFBAMethod; consumed and removed from param if .entropicFBAMethod is absent and .method contains ‘flux’

  • .internalBounds - deprecated; if present, an error is raised directing the caller to .internalNetFluxBounds/ .externalNetFluxBounds instead

  • .maxUnidirectionalFlux - maximum unidirectional flux (default: min(1e5, max(abs(model.ub))))

  • .minUnidirectionalFlux - minimum unidirectional flux (default: 0)

  • .internalNetFluxBounds - (‘original’) how to set the internal net flux bounds: ‘original’ (use model.lb/model.ub unmodified), ‘directional’ (direction only, magnitude set to param.maxUnidirectionalFlux), ‘max’ (+-param.maxUnidirectionalFlux on all internal reactions), ‘none’ (+-inf), ‘random’, ‘rangeNt’, or ‘expRangeNt’

  • .externalNetFluxBounds - (‘original’) how to set the bounds on external net flux and the initial/final concentration difference: ‘original’ (use model.lb/model.ub unmodified), ‘none’ (+-inf), ‘identities’, ‘bReplacingB’, or ‘dxReplacement’

OUTPUTS:

  • vl – n x 1 lower bound on internal net flux

  • vu – n x 1 upper bound on internal net flux

  • l_w – k x 1 lower bound on external net flux

  • u_w – 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

  • ce – 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

  • paramOut – param structure, returned to capture any fields set internally (e.g. defaulted fields)

solveCobraEP(EPproblem, varargin)[source]

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)

    • .F - backward-compatible alias for .Q (the naming used by

      solveCobraQP); copied to .Q and removed if .Q is absent

    • .sumFluxes - upper bound(s) on the auxiliary sum-of-flux conic

      variable(s) passed to mosek (default: [], i.e. no bound)

    • .sumConc - upper bound(s) on the auxiliary sum-of-concentration

      conic variable(s) passed to mosek (default: [])

    • .sumConc0 - upper bound(s) on the auxiliary sum-of-initial-

      concentration conic variable(s) passed to mosek (default: [])

    • 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 x1 double * .vf: n x 1 double * .vr: n x 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 x 1 double dual sol to \(lb <= vr <= ub\) * .z_vr: n x 1 double dual sol to \(lb <= vf <= ub\) * .z_vi: n x 1 double dual sol to \(lb <= v <= ub\) * .z_v: n + k x 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’