0001 function flux = geometricFBA(model,varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033 param = struct('epsilon',1e-6,'flexRel',0,'printLevel',1);
0034 field = fieldnames(param);
0035 if mod(nargin,2) ~= 1
0036 error('incorrect number of input parameters')
0037 else
0038 for k = 1:2:(nargin-1)
0039 param.(field{strcmp(varargin{k},field)}) = varargin{k+1};
0040 end
0041 end
0042 param.flexTol = param.flexRel * param.epsilon;
0043
0044
0045 FBAsolution = optimizeCbModel(model);
0046 ind = find(model.c);
0047 if length(ind) == 1
0048 model.lb(ind) = FBAsolution.f;
0049 end
0050
0051 A = model.S;
0052 b = model.b;
0053 L = model.lb;
0054 U = model.ub;
0055
0056
0057 b = b(:); L = L(:); U = U(:);
0058
0059
0060 J = any(A,2); A = A(J,:); b = b(J);
0061
0062
0063 v = nan(size(L));
0064 J = (U-L < param.epsilon);
0065 v(J) = (L(J)+U(J))/2;
0066 J = find(isnan(v));
0067
0068 if param.printLevel
0069 fprintf('%s\t%g\n\n%s\t@%s\n','# reactions:',length(v),'iteration #0',datestr(now,16));
0070 end
0071
0072 L0 = L; U0 = U;
0073 for k = J(:)'
0074 f = zeros(length(v),1); f(k) = -1;
0075 [dummy,opt,conv] = easyLP(f,A,b,L0,U0);
0076 if conv
0077 vL = max(-opt,L(k));
0078 else
0079 vL = L(k);
0080 end
0081 [dummy,opt,conv] = easyLP(-f,A,b,L0,U0);
0082 if conv
0083 vU = min(opt,U(k));
0084 else vU = U(k);
0085 end
0086 if abs(vL) < param.epsilon
0087 vL = 0;
0088 end
0089 if abs(vU) < param.epsilon
0090 vU = 0;
0091 end
0092 vM = (vL + vU)/2;
0093 if abs(vM) < param.epsilon
0094 vM = 0;
0095 end
0096 if abs(vU - vL) < param.epsilon
0097 vL = (1-sign(vM)* param.flexTol)*vM;
0098 vU = (1+sign(vM)* param.flexTol)*vM;
0099 end
0100 L(k) = vL; U(k) = vU;
0101 end
0102
0103 v = nan(size(L));
0104 J = (U-L < param.epsilon);
0105 v(J) = (L(J)+U(J))/2; v = v.*(abs(v) > param.epsilon);
0106
0107 if param.printLevel
0108 fprintf('%s\t\t%g\n%s\t\t%g\n\n','fixed:',sum(J),'@ zero:',sum(v==0));
0109 end
0110
0111
0112 J = find(U-L >= param.epsilon);
0113 n = 1;
0114 mu = [];
0115 Z = [];
0116
0117 while ~isempty(J)
0118
0119 if param.printLevel
0120 fprintf('%s #%g\t@%s\n','iteration',n,datestr(now,16));
0121 end
0122
0123 if n == 1
0124 M = zeros(size(L));
0125 else
0126 M = (L+U)/2;
0127 end
0128
0129 mu(:,n) = M;
0130 allL = L; allU = U; allA = A; allB = b;
0131 [a1,a2] = size(A);
0132
0133
0134 for k = 1:(n-1)
0135 [b1,b2] = size(allA);
0136 f = sparse(b2+2*a2,1); f((b2+1):end) = -1;
0137 opt = -Z(k);
0138 allA = [allA,sparse(b1,2*a2);
0139 speye(a2,a2),sparse(a2,b2-a2),-speye(a2),speye(a2);
0140 f(:)'];
0141 allB = [allB;mu(:,k);opt];
0142 allL = [allL;zeros(2*a2,1)];
0143 allU = [allU;inf*ones(2*a2,1)];
0144 end
0145
0146 [b1,b2] = size(allA);
0147 f = zeros(b2+2*a2,1); f((b2+1):end) = -1;
0148 allA = [allA,sparse(b1,2*a2);
0149 speye(a2,a2),sparse(a2,b2-a2),-speye(a2),speye(a2)];
0150 allB = [allB;M];
0151 allL = [allL;zeros(2*a2,1)];
0152 allU = [allU;inf*ones(2*a2,1)];
0153
0154 [v,opt,conv] = easyLP(f,allA,allB,allL,allU);
0155 if ~conv, disp('error: no convergence'); flux = (L+U)/2; return; end
0156
0157 opt = ceil(-opt/eps)*eps;
0158 Z(n) = opt;
0159 allA = [allA; sparse(f(:)')];
0160 allB = [allB; -opt];
0161
0162 for k = J(:)'
0163 f = zeros(length(allL),1); f(k) = -1;
0164 [dummy,opt,conv] = easyLP(f,allA,allB,allL,allU);
0165 if conv
0166 vL = max(-opt,L(k));
0167 else
0168 vL = L(k);
0169 end
0170 [dummy,opt,conv] = easyLP(-f,allA,allB,allL,allU);
0171 if conv
0172 vU = min(opt,U(k));
0173 else
0174 vU = U(k);
0175 end
0176 if abs(vL) < param.epsilon
0177 vL = 0;
0178 end
0179 if abs(vU) < param.epsilon
0180 vU = 0;
0181 end
0182 vM = (vL + vU)/2;
0183 if abs(vM) < param.epsilon
0184 vM = 0;
0185 end
0186 if abs(vU - vL) < param.epsilon
0187 vL = (1-sign(vM)* param.flexTol)*vM;
0188 vU = (1+sign(vM)* param.flexTol)*vM;
0189 end
0190 L(k) = vL;
0191 U(k) = vU;
0192 end
0193
0194 v = nan(size(L));
0195 J = (U-L < param.epsilon);
0196 v(J) = (L(J)+U(J))/2; v = v.*(abs(v) > param.epsilon);
0197
0198 if param.printLevel
0199 fprintf('%s\t\t%g\n%s\t\t%g\n\n','fixed:',sum(J),'@ zero:',sum(v==0));
0200 end
0201
0202 n = n+1;
0203 J = find(U-L >= param.epsilon);
0204
0205 flux = v;
0206 end
0207
0208 function [v,fOpt,conv] = easyLP(c,A,b,lb,ub)
0209
0210
0211
0212
0213
0214
0215
0216
0217
0218
0219
0220
0221
0222
0223
0224
0225
0226
0227
0228
0229
0230
0231
0232 csense(1:length(b)) = 'E';
0233 model = struct('A',A,'b',b,'c',full(c),'lb',lb,'ub',ub,'osense',-1,'csense',csense);
0234 solution = solveCobraLP(model);
0235 v = solution.full;
0236 fOpt = solution.obj;
0237 conv = solution.stat == 1;