Skip to content

Commit 33f4db5

Browse files
rmtflemingrmtfleming
authored
Merge pull request #2100 from opencobra/develop
Develop
2 parents 44a3e1e + 9919575 commit 33f4db5

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

44 files changed

+15644
-2
lines changed

‎papers

Submodule papers updated 17 files

‎src/base/solvers/entropicFBA/addExoMetToEFBA.m

Lines changed: 112 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,112 @@
1+
function modelOut = addExoMetToEFBA(model,exoMet,param)
2+
% generates H and h to add a min (alpha/2)(v-h)'*H*(v-h) term to an EFBA problem
3+
%
4+
% USAGE:
5+
% modelOut = addExoMetToEFBA(model,exoMet,param)
6+
%
7+
% INPUTS:
8+
% model.S:
9+
% model.rxns:
10+
%
11+
% exoMet.rxns:
12+
% exoMet.mean:
13+
% exoMet.SD:
14+
%
15+
% OPTIONAL INPUT
16+
% param.printLevel:
17+
% param.alpha: alpha in (alpha/2)(v-h)'*H*(v-h), default = 10000;
18+
% param.metabolomicWeights: String indicating the type of weights to be applied to penalise the difference
19+
% between of predicted and experimentally measured fluxes by, where
20+
% 'SD' weights = 1/(1+exoMet.SD^2)
21+
% 'mean' weights = 1/(1+exoMet.mean^2) (Default)
22+
% 'RSD' weights = 1/((exoMet.SD./exoMet.mean)^2)
23+
% param.relaxBounds: True to relax bounds on reaction whose fluxes are to be fitted to exometabolomic data
24+
%
25+
% OUTPUTS:
26+
% modelOut:
27+
%
28+
% EXAMPLE:
29+
%
30+
% NOTE:
31+
%
32+
% Author(s):
33+
34+
if ~exist('param','var')
35+
param=struct;
36+
end
37+
if ~isfield(param,'printLevel')
38+
param.printLevel=0;
39+
end
40+
if ~isfield(param,'alpha')
41+
param.alpha=10000;
42+
end
43+
if ~isfield(param,'relaxBounds')
44+
param.relaxBounds=1;
45+
end
46+
if ~isfield(param,'metabolomicWeights')
47+
param.metabolomicWeights = 'mean';
48+
end
49+
50+
[bool, locb] = ismember(exoMet.rxns, model.rxns);
51+
if any(locb == 0)
52+
if printLevel > 0
53+
fprintf('%s\n','The following exometabolomic exchange reactions were not found in the model:')
54+
disp(exoMet.rxns(locb == 0))
55+
end
56+
end
57+
58+
if length(unique(exoMet.rxns)) ~= length(exoMet.rxns) && printLevel > 0
59+
disp('There are duplicate rxnID entries in the metabolomic data! Only using data corresponding to first occurance')
60+
end
61+
62+
% Assume mean model flux is equal to the mean experimental reaction flux.
63+
[~,nRxn] = size(model.S);
64+
vExp = NaN * ones(nRxn,1);
65+
vExp(locb(bool)) = exoMet.mean(bool);
66+
67+
vSD = NaN * ones(nRxn,1);
68+
vSD(locb(bool)) = exoMet.SD(bool);
69+
70+
% Set the weight on the Euclidean distance of the predicted steady state
71+
% flux from the experimental steady state flux. Uniform at first.
72+
weightExp = ones(nRxn,1);
73+
% The weight penalty on relaxation from the experimental flux vector should be greater than that
74+
% on the lower and upper bounds, if the experimental data is considered more reliable
75+
% than the lower and upper bounds of the model.
76+
% Assumes the lower and upper bound of standard deviation of
77+
% experimental reaction flux are separated by two standard
78+
% deviations.
79+
80+
% Penalise the relaxation from the experimental flux value by 1/(1+weights^2)
81+
switch param.metabolomicWeights
82+
case 'SD'
83+
weightExp(locb(bool)) = 1 ./ (1 + (vSD(locb(bool))).^2);
84+
case 'mean'
85+
weightExp(locb(bool)) = 1 ./ (1 + (vExp(locb(bool))).^2);
86+
case 'RSD'
87+
weightExp(locb(bool)) = 1 ./ ((vSD(locb(bool))./vExp(locb(bool))).^2);
88+
otherwise
89+
weightExp(locb(bool)) = 2;
90+
end
91+
% Weights are ignored on the reactions without experimental data, i.e. vExp(n) == NaN.
92+
weightExp(~isfinite(vExp)) = 0;
93+
94+
%add
95+
model.H=diag(sparse(weightExp*param.alpha));
96+
model.h=vExp;
97+
if param.printLevel>1
98+
figure;
99+
histogram(weightExp)
100+
xlabel('weights on the diagonal of model.H')
101+
ylabel('#reactions')
102+
end
103+
104+
if param.relaxBounds
105+
bool=ismember(model.rxns,exoMet.rxns);
106+
model.lb(bool) = model.lb_preconstrainRxns(bool);
107+
model.ub(bool) = model.ub_preconstrainRxns(bool);
108+
end
109+
110+
modelOut = model;
111+
end
112+

‎src/base/solvers/entropicFBA/debugInfeasibleEntropyFBA.m

Lines changed: 114 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,114 @@
1+
function [sol,modelOut] = debugInfeasibleEntropyFBA(model)
2+
%try to diagnose the reasons a model is infeasible for Entropy FBA
3+
4+
sol = optimizeCbModel(model);
5+
if sol.stat~=1
6+
warning('Model does not admit a flux balance analysis solution')
7+
return
8+
end
9+
10+
param.printLevel = 1;
11+
12+
%test with default parameters
13+
[sol, ~] = entropicFluxBalanceAnalysis(model,param);
14+
15+
if sol.stat ==1
16+
fprintf('%s\n',['EFBA feasible with default solver: ' sol.solver '.'])
17+
else
18+
fprintf('%s\n',['EFBA infeasible with default solver: ' sol.solver '.'])
19+
end
20+
21+
switch sol.solver
22+
case 'mosek'
23+
solMOSEK = sol;
24+
25+
%test with pdco
26+
param.solver = 'pdco';
27+
[solPDCO, ~] = entropicFluxBalanceAnalysis(model,param);
28+
29+
if solPDCO.stat ==1
30+
fprintf('%s\n',['EFBA feasible with ' solPDCO.solver '.'])
31+
else
32+
fprintf('%s\n',['EFBA infeasible with ' solPDCO.solver '.'])
33+
end
34+
case 'pdco'
35+
solPDCO = sol;
36+
37+
%test with mosek
38+
param.solver ='mosek';
39+
[solMOSEK, ~] = entropicFluxBalanceAnalysis(model,param);
40+
if solMOSEK.stat ==1
41+
fprintf('%s\n',['EFBA feasible with ' solMOSEK.solver '.'])
42+
else
43+
fprintf('%s\n',['EFBA infeasible with ' solMOSEK.solver '.'])
44+
end
45+
end
46+
47+
48+
if solPDCO.stat ~= solMOSEK.stat
49+
warning('pdco and mosek sol.stat are inconsistent')
50+
end
51+
52+
53+
if solMOSEK.stat==1
54+
fprintf('%s\n','EFBA feasible with mosek')
55+
end
56+
fprintf('%s\n','EFBA infeasible with mosek')
57+
58+
%test without coupling constraints
59+
modelTmp = rmfield(model,'C');
60+
modelTmp = rmfield(modelTmp,'d');
61+
[sol, ~] = entropicFluxBalanceAnalysis(modelTmp,param);
62+
if sol.stat==1
63+
fprintf('%s\n','Coupling constraints are causing thermodynamic infeasibility, removed.')
64+
modelOut = modelTmp;
65+
return
66+
else
67+
param.internalNetFluxBounds = 'directional';
68+
[sol, modelOut] = entropicFluxBalanceAnalysis(modelTmp,param);
69+
if sol.stat==1
70+
fprintf('%s\n','Internal finite flux bounds are causing thermodynamic infeasibility, removed.')
71+
else
72+
73+
param.internalNetFluxBounds = 'max';
74+
[sol, modelOut] = entropicFluxBalanceAnalysis(modelTmp,param);
75+
if sol.stat==1
76+
fprintf('%s\n','Small finite innternal directional flux bounds are causing thermodynamic infeasibility, removed.')
77+
else
78+
param.internalNetFluxBounds = 'none';
79+
[sol, modelOut] = entropicFluxBalanceAnalysis(modelTmp,param);
80+
if sol.stat==1
81+
fprintf('%s\n','Internal directional flux bounds are causing thermodynamic infeasibility, removed.')
82+
else
83+
%try rescaling finite model constraints, i.e. rhs, lb, ub
84+
param.internalNetFluxBounds = 'none';
85+
scaleFactor = 1e-2;
86+
modelTmp2 = modelTmp;
87+
modelTmp2.lb(~modelTmp.SConsistentRxnBool) = modelTmp.lb(~modelTmp.SConsistentRxnBool)*scaleFactor;
88+
modelTmp2.ub(~modelTmp.SConsistentRxnBool) = modelTmp.ub(~modelTmp.SConsistentRxnBool)*scaleFactor;
89+
modelTmp2.b = modelTmp.b*scaleFactor;
90+
if isfield(modelTmp,'d')
91+
modelTmp2.d = modelTmp.d*scaleFactor;
92+
end
93+
[sol, modelOut] = entropicFluxBalanceAnalysis(modelTmp2,param);
94+
if sol.stat==1
95+
fprintf('%s\n','Multiscale exchange bounds are causing thermodynamic infeasibility, rescaled.')
96+
else
97+
%try relaxing the exchange bounds
98+
modelTmp2 = modelTmp;
99+
modelTmp2.lb(~modelTmp.SConsistentRxnBool) = modelTmp.lb(~modelTmp.SConsistentRxnBool) - 1000;
100+
modelTmp2.ub(~modelTmp.SConsistentRxnBool) = modelTmp.ub(~modelTmp.SConsistentRxnBool) + 1000;
101+
param.internalNetFluxBounds = 'none';
102+
[sol, modelOut] = entropicFluxBalanceAnalysis(modelTmp2,param);
103+
if sol.stat==1
104+
fprintf('%s\n','Exchange bounds are too tight and causing thermodynamic infeasibility, relaxed.')
105+
else
106+
107+
end
108+
end
109+
end
110+
end
111+
end
112+
end
113+
end
114+

0 commit comments

Comments
 (0)