Relaxed Flux Balance Analysis: Toy model
Author: Ronan Fleming, Systems Biochemistry Group, University of Luxembourg.
Reviewer:
Introduction
We consider a biochemical network of m molecular species and n biochemical reactions. The biochemical network is mathematically represented by a stoichiometric matrix . In standard notation, flux balance analysis (FBA) is the linear optimisation problem where is a parameter vector that linearly combines one or more reaction fluxes to form what is termed the objective function, and where a , or , represents some fixed output, or input, of the ith molecular species. Every FBA solution must satisfy the constraints, independent of any objective chosen to optimise over the set of constraints. It may occur that the constraints on the FBA problem are not all simultaneously feasible, i.e., the system of inequalities is infeasible. This situation might be caused by an incorrectly specified reaction bound or the absence of a reaction from the stoichiometric matrix, such that a nonzero . To resolve the infeasiblility, we consider a cardinality optimisation problem that seeks to minimise the number of bounds to relax, the number of fixed outputs to relax, the number of fixed inputs to relax, or a combination of all three, in order to render the problem feasible. The cardinality optimisation problem, termed relaxed flux balance analysis, is
where denote the relaxations of the lower and upper bounds on reaction rates of the reaction rates vector v, and where denotes a relaxation of the mass balance constraint. Non-negative scalar parameters λ and can be used to trade off between relaxation of mass balance or bound constraints. A non-negative vector parameter λ can be used to prioritise relaxation of one mass balance constraint over another, e.g, to avoid relaxation of a mass balance constraint on a metabolite that is not desired to be exchanged across the boundary of the system. A non-negative vector parameter may be used to prioritise relaxation of bounds on some reactions rather than others, e.g., relaxation of bounds on exchange reactions rather than internal reactions. The optimal choice of parameters depends heavily on the biochemical context. A relaxation of the minimum number of constraints is desirable because ideally one should be able to justify the choice of bounds or choice of metabolites to be exchanged across the boundary of the system by recourse to experimental literature. This task is magnified by the number of constraints proposed to be relaxed. PROCEDURE: RelaxedFBA applied to a toy model
rxnForms = {' -> A','A -> B','B -> C', 'B -> D','D -> C','C ->'};
rxnNames = {'R1','R2','R3','R4','R5', 'R6'};
model = createModel(rxnNames, rxnNames,rxnForms);
Metabolite A[c] not in model - added to the model
R1 <=> A[c]
Metabolite B[c] not in model - added to the model
R2 A[c] <=> B[c]
Metabolite C[c] not in model - added to the model
R3 B[c] <=> C[c]
Metabolite D[c] not in model - added to the model
R4 B[c] <=> D[c]
R5 D[c] <=> C[c]
R6 C[c] <=>
Print the constraints
printConstraints(model, -1001, 1001)
MinConstraints:
R1 -1000
R2 -1000
R3 2
R4 2
R5 -1000
R6 -1000
maxConstraints:
R1 1000
R2 1000
R3 1000
R4 1000
R5 1000
R6 3
Identify the exchange reactions and biomass reaction(s) heuristically
model = findSExRxnInd(model,size(model.S,1),0);
% relaxOption: Structure containing the relaxation options:
% * 0 = do not allow to relax bounds on internal reactions
% * 1 = do not allow to relax bounds on internal reactions with finite bounds
% * 2 = allow to relax bounds on all internal reactions
relaxOption.internalRelax = 2;
% * 0 = do not allow to relax bounds on exchange reactions
% * 1 = do not allow to relax bounds on exchange reactions of the type [0,0]
% * 2 = allow to relax bounds on all exchange reactions
relaxOption.exchangeRelax = 0;
% * 0 = do not allow to relax the steady state constraint S*v = b
% * 1 = allow to relax the steady state constraint S*v = b
relaxOption.steadyStateRelax = 0;
% * toBeUnblockedReactions - n x 1 vector indicating the reactions to be unblocked (optional)
% * toBeUnblockedReactions(i) = 1 : impose v(i) to be positive
% * toBeUnblockedReactions(i) = -1 : impose v(i) to be negative
% * toBeUnblockedReactions(i) = 0 : do not add any constraint
% * excludedReactions - n x 1 bool vector indicating the reactions to be excluded from relaxation (optional)
% * excludedReactions(i) = false : allow to relax bounds on reaction i
% * excludedReactions(i) = true : do not allow to relax bounds on reaction i
% * excludedMetabolites - m x 1 bool vector indicating the metabolites to be excluded from relaxation (optional)
% * excludedMetabolites(i) = false : allow to relax steady state constraint on metabolite i
% * excludedMetabolites(i) = true : do not allow to relax steady state constraint on metabolite i
Set the tolerance to distinguish between zero and non-zero flux
%feasTol = getCobraSolverParams('LP', 'feasTol');
%relaxOption.epsilon = feasTol/100;%*100;
relaxOption.epsilon = 10e-6;
relaxOption.nbMaxIteration = 1000;
relaxOption.epsilon = 10e-6;
relaxOption.gamma0 = 10; %trade-off parameter of l0 part of v
relaxOption.gamma1 = 1; %trade-off parameter of l1 part of v
relaxOption.lambda0 = 10; %trade-off parameter of l0 part of r
relaxOption.lambda1 = 0; %trade-off parameter of l1 part of r
relaxOption.alpha0 = 0; %trade-off parameter of l0 part of p and q
relaxOption.alpha1 = 0; %trade-off parameter of l1 part of p and q
relaxOption.theta = 2; %parameter of capped l1 approximation
Check if the model is feasible
FBAsolution = optimizeCbModel(model,'max', 0, true);
disp('Model is feasible. Nothing to do.');
disp('Model is infeasible');
Call the relaxedFBA function, deal the solution, and set small values to zero
relaxOption
relaxOption =
internalRelax: 2
exchangeRelax: 0
steadyStateRelax: 0
epsilon: 1.0000e-05
solution = relaxedFBA(model,relaxOption);
[v,r,p,q] = deal(solution.v,solution.r,solution.p,solution.q);
p(p<relaxOption.epsilon) = 0;%lower bound relaxation
q(q<relaxOption.epsilon) = 0;%upper bound relaxation
r(r<relaxOption.epsilon) = 0;%steady state constraint relaxation
The output is a solution structure with a 'stat' field reporting the solver status and a set of fields matching the relaxation of constraints given in the mathematical formulation of the relaxed flux balance problem above.
% solution: Structure containing the following fields:
% * r - relaxation on steady state constraints S*v = b
% * p - relaxation on lower bound of reactions
% * q - relaxation on upper bound of reactions
Display the proposed relaxation solution
fprintf('%s\n','Relaxation of steady state constraints:')
Relaxation of steady state constraints:
fprintf('%s\n','Relaxation on lower bound of reactions:')
Relaxation on lower bound of reactions:
fprintf('%s\n','Relaxation on upper bound of reactions:')
Relaxation on upper bound of reactions:
Summarise the proposed relaxation solution
dispCutoff=relaxOption.epsilon;
fprintf('%s\n',['Relaxed flux balance analysis problem solved in ' num2str(timeTaken) ' seconds.'])
fprintf('%u%s\n',nnz(r),' steady state constraints relaxed');
fprintf('%u%s\n',nnz(abs(p)>dispCutoff & ~abs(q)>dispCutoff & model.SIntRxnBool),' internal lower bounds relaxed');
fprintf('%u%s\n',nnz(abs(q)>dispCutoff & ~abs(p)>dispCutoff & model.SIntRxnBool),' internal upper bounds relaxed');
fprintf('%u%s\n',nnz(abs(p)>dispCutoff & abs(q)>dispCutoff & model.SIntRxnBool),' internal lower and upper bounds relaxed');
fprintf('%u%s\n',nnz(abs(p)>dispCutoff & ~abs(q)>dispCutoff & ~model.SIntRxnBool),' external lower bounds relaxed');
fprintf('%u%s\n',nnz(abs(q)>dispCutoff & ~abs(p)>dispCutoff & ~model.SIntRxnBool),' external upper bounds relaxed');
fprintf('%u%s\n',nnz(abs(p)>dispCutoff & abs(q)>dispCutoff & ~model.SIntRxnBool),' external lower and upper bounds relaxed');
fprintf('%u%s\n',nnz(abs(p)>dispCutoff | abs(q)>dispCutoff & ~model.SIntRxnBool),' external lower or upper bounds relaxed');
maxUB = max(max(model.ub),-min(model.lb));
minLB = min(-max(model.ub),min(model.lb));
intRxnFiniteBound = ((model.ub < maxUB) & (model.lb > minLB));
fprintf('%u%s\n',nnz(abs(p)>dispCutoff & intRxnFiniteBound),' finite lower bounds relaxed');
fprintf('%u%s\n',nnz(abs(q)>dispCutoff & intRxnFiniteBound),' finite upper bounds relaxed');
exRxn00 = ((model.ub == 0) & (model.lb == 0));
fprintf('%u%s\n',nnz(abs(p)>dispCutoff & exRxn00),' lower bounds relaxed on fixed reactions (lb=ub=0)');
fprintf('%u%s\n',nnz(abs(q)>dispCutoff & exRxn00),' upper bounds relaxed on fixed reactions (lb=ub=0)');
disp('relaxedFBA problem infeasible, check relaxOption fields');
end
Relaxed flux balance analysis problem solved in 0.048242 seconds.
0 steady state constraints relaxed
2 internal lower bounds relaxed
0 internal upper bounds relaxed
0 internal lower and upper bounds relaxed
0 external lower bounds relaxed
0 external upper bounds relaxed
0 external lower and upper bounds relaxed
2 external lower or upper bounds relaxed
0 finite lower bounds relaxed
0 finite upper bounds relaxed
0 lower bounds relaxed on fixed reactions (lb=ub=0)
0 upper bounds relaxed on fixed reactions (lb=ub=0)
Another example
rxnNames = {'R1','R2','R3','R4','R5','R6','R7','R8','R9','R10','R11'};
model = createModel(rxnNames, rxnNames,rxnForms);
Assume all reactions are irreversible
Reaction R7 with bounds 1 <= v_7 <= 10
Print the constraints
printConstraints(model, -1001, 1001)
Identify the exchange reactions and biomass reaction(s) heuristically
model = findSExRxnInd(model,size(model.S,1),0);
% relaxOption: Structure containing the relaxation options:
% * 0 = do not allow to relax bounds on internal reactions
% * 1 = do not allow to relax bounds on internal reactions with finite bounds
% * 2 = allow to relax bounds on all internal reactions
relaxOption.internalRelax = 2;
% * 0 = do not allow to relax bounds on exchange reactions
% * 1 = do not allow to relax bounds on exchange reactions of the type [0,0]
% * 2 = allow to relax bounds on all exchange reactions
relaxOption.exchangeRelax = 0;
% * 0 = do not allow to relax the steady state constraint S*v = b
% * 1 = allow to relax the steady state constraint S*v = b
relaxOption.steadyStateRelax = 0;
% * toBeUnblockedReactions - n x 1 vector indicating the reactions to be unblocked (optional)
% * toBeUnblockedReactions(i) = 1 : impose v(i) to be positive
% * toBeUnblockedReactions(i) = -1 : impose v(i) to be negative
% * toBeUnblockedReactions(i) = 0 : do not add any constraint
% * excludedReactions - n x 1 bool vector indicating the reactions to be excluded from relaxation (optional)
% * excludedReactions(i) = false : allow to relax bounds on reaction i
% * excludedReactions(i) = true : do not allow to relax bounds on reaction i
% * excludedMetabolites - m x 1 bool vector indicating the metabolites to be excluded from relaxation (optional)
% * excludedMetabolites(i) = false : allow to relax steady state constraint on metabolite i
% * excludedMetabolites(i) = true : do not allow to relax steady state constraint on metabolite i
Set the tolerance to distinguish between zero and non-zero flux
feasTol = getCobraSolverParams('LP', 'feasTol');
relaxOption.epsilon = feasTol/100;%*100;
Call the relaxedFBA function, deal the solution, and set small values to zero
solution = relaxedFBA(model,relaxOption);
[v,r,p,q] = deal(solution.v,solution.r,solution.p,solution.q);
p(p<relaxOption.epsilon) = 0;%lower bound relaxation
q(q<relaxOption.epsilon) = 0;%upper bound relaxation
r(r<relaxOption.epsilon) = 0;%steady state constraint relaxation
The output is a solution structure with a 'stat' field reporting the solver status and a set of fields matching the relaxation of constraints given in the mathematical formulation of the relaxed flux balance problem above.
% solution: Structure containing the following fields:
% * r - relaxation on steady state constraints S*v = b
% * p - relaxation on lower bound of reactions
% * q - relaxation on upper bound of reactions
Summarise the proposed relaxation solution
dispCutoff=relaxOption.epsilon;
fprintf('%s\n',['Relaxed flux balance analysis problem solved in ' num2str(timeTaken) ' seconds.'])
fprintf('%u%s\n',nnz(r),' steady state constraints relaxed');
fprintf('%u%s\n',nnz(abs(p)>dispCutoff & ~abs(q)>dispCutoff & model.SIntRxnBool),' internal only lower bounds relaxed');
fprintf('%u%s\n',nnz(abs(q)>dispCutoff & ~abs(p)>dispCutoff & model.SIntRxnBool),' internal only upper bounds relaxed');
fprintf('%u%s\n',nnz(abs(p)>dispCutoff & abs(q)>dispCutoff & model.SIntRxnBool),' internal lower and upper bounds relaxed');
fprintf('%u%s\n',nnz(abs(p)>dispCutoff & ~abs(q)>dispCutoff & ~model.SIntRxnBool),' external only lower bounds relaxed');
fprintf('%u%s\n',nnz(abs(q)>dispCutoff & ~abs(p)>dispCutoff & ~model.SIntRxnBool),' external only upper bounds relaxed');
fprintf('%u%s\n',nnz(abs(p)>dispCutoff & abs(q)>dispCutoff & ~model.SIntRxnBool),' external lower and upper bounds relaxed');
fprintf('%u%s\n',nnz(abs(p)>dispCutoff | abs(q)>dispCutoff & ~model.SIntRxnBool),' external lower or upper bounds relaxed');
maxUB = max(max(model.ub),-min(model.lb));
minLB = min(-max(model.ub),min(model.lb));
intRxnFiniteBound = ((model.ub < maxUB) & (model.lb > minLB));
fprintf('%u%s\n',nnz(abs(p)>dispCutoff & intRxnFiniteBound),' finite lower bounds relaxed');
fprintf('%u%s\n',nnz(abs(q)>dispCutoff & intRxnFiniteBound),' finite upper bounds relaxed');
exRxn00 = ((model.ub == 0) & (model.lb == 0));
fprintf('%u%s\n',nnz(abs(p)>dispCutoff & exRxn00),' lower bounds relaxed on fixed reactions (lb=ub=0)');
fprintf('%u%s\n',nnz(abs(q)>dispCutoff & exRxn00),' upper bounds relaxed on fixed reactions (lb=ub=0)');
disp('relaxedFBA problem infeasible, check relaxOption fields');
REFERENCES
Fleming, R.M.T., et al., Cardinality optimisation in constraint-based modelling: Application to Recon 3D (submitted), 2017.
Brunk, E. et al. Recon 3D: A resource enabling a three-dimensional view of gene variation in human metabolism. (submitted) 2017.