gapFind Identifies all blocked metabolites (anything downstream of a gap) in a model. MILP algorithm that finds gaps that may be missed by simple inspection of the S matrix. To find every gap in a model, change the rxn bounds on all exchange reactions to allow uptake of every metabolite. [allGaps,rootGaps,downstreamGaps] = gapFind(model,findNCgaps,verbFlag) INPUT model a COBRA model OPTIONAL INPUTS findNCgaps find no consupmption gaps as well as no production gaps (default false) verbFlag verbose flag (default false) OUTPUTS allGaps all gaps found by GapFind rootGaps all root no production (and consumption) gaps downstreamGaps all downstream gaps based on: Kumar, V. et al. BMC Bioinformatics. 2007 Jun 20;8:212. solve problem: max ||xnp|| s.t. S(i,j)*v(j) >= e*w(i,j) S(i,j) > 0, j in IR S(i,j)*v(j) <= M*w(i,j) S(i,j) > 0, j in IR S(i,j)*v(j) >= e - M(1-w(i,j)) S(i,j) ~= 0, j in R S(i,j)*v(j) <= M*w(i,j) S(i,j) ~= 0, j in R ||w(i,j)|| >= xnp(i) lb <= v <= ub S*v >= 0 xnp(i) = {0,1} w(i,j) = {0,1} reformulated for COBRA MILP as: max sum(xnp(:)) s.t. S*v >= 0 (or = 0 if findNCgaps = true) (1) S(i,j)*v(j) - e*w(i,j) >= 0 S(i,j) > 0, j in IR (2) S(i,j)*v(j) - M*w(i,j) <= 0 S(i,j) > 0, j in IR (3) S(i,j)*v(j) - M*w(i,j) >= e-M S(i,j) ~= 0, j in R (4) S(i,j)*v(j) - M*w(i,j) <= 0 S(i,j) ~= 0, j in R (5) sum(w(i,:)) - xnp(i) >= 0 (6) lb <= v <= ub xnp and w are binary variables, v are continuous Jeff Orth 7/6/09
0001 function [allGaps,rootGaps,downstreamGaps] = gapFind(model,findNCgaps,verbFlag) 0002 %gapFind Identifies all blocked metabolites (anything downstream of a gap) 0003 %in a model. MILP algorithm that finds gaps that may be missed by simple 0004 %inspection of the S matrix. To find every gap in a model, change the rxn 0005 %bounds on all exchange reactions to allow uptake of every metabolite. 0006 % 0007 % [allGaps,rootGaps,downstreamGaps] = gapFind(model,findNCgaps,verbFlag) 0008 % 0009 %INPUT 0010 % model a COBRA model 0011 % 0012 %OPTIONAL INPUTS 0013 % findNCgaps find no consupmption gaps as well as no production gaps 0014 % (default false) 0015 % verbFlag verbose flag (default false) 0016 % 0017 %OUTPUTS 0018 % allGaps all gaps found by GapFind 0019 % rootGaps all root no production (and consumption) gaps 0020 % downstreamGaps all downstream gaps 0021 % 0022 % based on: 0023 % Kumar, V. et al. BMC Bioinformatics. 2007 Jun 20;8:212. 0024 % 0025 % solve problem: 0026 % max ||xnp|| 0027 % s.t. S(i,j)*v(j) >= e*w(i,j) S(i,j) > 0, j in IR 0028 % S(i,j)*v(j) <= M*w(i,j) S(i,j) > 0, j in IR 0029 % S(i,j)*v(j) >= e - M(1-w(i,j)) S(i,j) ~= 0, j in R 0030 % S(i,j)*v(j) <= M*w(i,j) S(i,j) ~= 0, j in R 0031 % ||w(i,j)|| >= xnp(i) 0032 % lb <= v <= ub 0033 % S*v >= 0 0034 % xnp(i) = {0,1} 0035 % w(i,j) = {0,1} 0036 % 0037 % reformulated for COBRA MILP as: 0038 % max sum(xnp(:)) 0039 % s.t. S*v >= 0 (or = 0 if findNCgaps = true) (1) 0040 % S(i,j)*v(j) - e*w(i,j) >= 0 S(i,j) > 0, j in IR (2) 0041 % S(i,j)*v(j) - M*w(i,j) <= 0 S(i,j) > 0, j in IR (3) 0042 % S(i,j)*v(j) - M*w(i,j) >= e-M S(i,j) ~= 0, j in R (4) 0043 % S(i,j)*v(j) - M*w(i,j) <= 0 S(i,j) ~= 0, j in R (5) 0044 % sum(w(i,:)) - xnp(i) >= 0 (6) 0045 % lb <= v <= ub 0046 % xnp and w are binary variables, v are continuous 0047 % 0048 % 0049 % Jeff Orth 7/6/09 0050 0051 if nargin < 2 0052 findNCgaps = false; 0053 end 0054 if nargin < 3 0055 verbFlag = false; 0056 end 0057 0058 M = length(model.rxns); %this was set to 100 in GAMS GapFind implementation 0059 N = length(model.mets); 0060 R = model.rev ~= 0; %reversible reactions 0061 R_index = find(R); 0062 IR = model.rev == 0; %irreversible reactions 0063 IR_index = find(IR); 0064 e = 0.0001; 0065 S = model.S; 0066 lb = model.lb; 0067 ub = model.ub; 0068 0069 % MILPproblem 0070 % A LHS matrix 0071 % b RHS vector 0072 % c Objective coeff vector 0073 % lb Lower bound vector 0074 % ub Upper bound vector 0075 % osense Objective sense (-1 max, +1 min) 0076 % csense Constraint senses, a string containting the constraint sense for 0077 % each row in A ('E', equality, 'G' greater than, 'L' less than). 0078 % vartype Variable types 0079 % x0 Initial solution 0080 0081 % initialize MILP fields 0082 0083 % get number of rows and cols for each constraint 0084 %rows 0085 m_c1 = N; %number of metabolites 0086 m_c2 = length(find(S(:,IR) > 0)); %number of Sij>0 in irreverisible reactions 0087 m_c3 = m_c2; 0088 m_c4 = length(find(S(:,R))); %number of Sij>0 in reversible reactions 0089 m_c5 = m_c4; 0090 m_c6 = N; %number of xnp (metabolites) 0091 %columns 0092 n_v = M; %number of reactions 0093 n_wij_IR = m_c2; 0094 n_wij_R = m_c4; 0095 n_xnp = N; 0096 0097 % LHS matrix A 0098 0099 % constraint 1 0100 A = [S sparse(m_c1,(n_wij_IR+n_wij_R+n_xnp))]; 0101 0102 % constraint 2 0103 % create Sij IR matrix and wij IR matrix 0104 Sij_IR = sparse(m_c2,n_v); 0105 wij_IR = sparse(m_c6,n_wij_IR); 0106 row = 1; 0107 for i = 1:length(IR_index) 0108 rxn_index = IR_index(i); 0109 met_index = find(S(:,rxn_index) > 0); 0110 for j = 1:length(met_index) 0111 Sij_IR(row,rxn_index) = S(met_index(j),rxn_index); 0112 wij_IR(met_index(j),row) = 1; 0113 row = row + 1; 0114 end 0115 end 0116 0117 A = [A ; Sij_IR -e*speye(m_c2,n_wij_IR) sparse(m_c2,(n_wij_R+n_xnp))]; 0118 0119 % constraint 3 0120 A = [A ; Sij_IR -M*speye(m_c3,n_wij_IR) sparse(m_c3,(n_wij_R+n_xnp))]; 0121 0122 % constraint 4 0123 % create Sij R matrix 0124 Sij_R = sparse(m_c4,n_v); 0125 wij_R = sparse(m_c6,n_wij_R); 0126 row = 1; 0127 for i = 1:length(R_index) 0128 rxn_index = R_index(i); 0129 met_index = find(S(:,rxn_index) ~= 0); 0130 for j = 1:length(met_index) 0131 Sij_R(row,rxn_index) = S(met_index(j),rxn_index); 0132 wij_R(met_index(j),row) = 1; 0133 row = row + 1; 0134 end 0135 end 0136 0137 A = [A ; Sij_R sparse(m_c4,n_wij_IR) -M*speye(m_c4,n_wij_R) sparse(m_c4,n_xnp)]; 0138 0139 % constraint 5 0140 A = [A ; Sij_R sparse(m_c5,n_wij_IR) -M*speye(m_c5,n_wij_R) sparse(m_c5,n_xnp)]; 0141 0142 % constraint 6 0143 A = [A ; sparse(m_c6,n_v) wij_IR wij_R -1*speye(m_c6,n_xnp)]; 0144 0145 % RHS vector b 0146 b = [zeros(m_c1+m_c2+m_c3,1);(e-M)*ones(m_c4,1);zeros(m_c5+m_c6,1)]; 0147 0148 % objective coefficient vector c 0149 c = [zeros(n_v+n_wij_IR+n_wij_R,1);ones(n_xnp,1)]; 0150 0151 % upper and lower bounds on variables (v,w,xnp) 0152 lb = [lb;zeros(n_wij_IR+n_wij_R+n_xnp,1)]; 0153 ub = [ub;ones(n_wij_IR+n_wij_R+n_xnp,1)]; 0154 0155 % objective sense osense 0156 osense = -1; %want to maximize objective 0157 0158 % constraint senses csense 0159 if findNCgaps 0160 csense(1:m_c1) = 'E'; 0161 else 0162 csense(1:m_c1) = 'G'; 0163 end 0164 csense((m_c1+1):(m_c1+m_c2)) = 'G'; 0165 csense((m_c1+m_c2+1):(m_c1+m_c2+m_c3)) = 'L'; 0166 csense((m_c1+m_c2+m_c3+1):(m_c1+m_c2+m_c3+m_c4)) = 'G'; 0167 csense((m_c1+m_c2+m_c3+m_c4+1):(m_c1+m_c2+m_c3+m_c4+m_c5)) = 'L'; 0168 csense((m_c1+m_c2+m_c3+m_c4+m_c5+1):(m_c1+m_c2+m_c3+m_c4+m_c5+m_c6)) = 'G'; 0169 0170 % variable types vartype 0171 vartype(1:n_v) = 'C'; 0172 vartype((n_v+1):(n_v+n_wij_IR+n_wij_R+n_xnp)) = 'B'; 0173 0174 % inital solution x0 0175 x0 = []; 0176 0177 0178 % run COBRA MILP solver 0179 gapFindMILPproblem.A = A; 0180 gapFindMILPproblem.b = b; 0181 gapFindMILPproblem.c = c; 0182 gapFindMILPproblem.lb = lb; 0183 gapFindMILPproblem.ub = ub; 0184 gapFindMILPproblem.osense = osense; 0185 gapFindMILPproblem.csense = csense; 0186 gapFindMILPproblem.vartype = vartype; 0187 gapFindMILPproblem.x0 = x0; 0188 0189 if verbFlag 0190 parameters.printLevel = 3; 0191 else 0192 parameters.printLevel = 0; 0193 end 0194 0195 solution = solveCobraMILP(gapFindMILPproblem,parameters); 0196 0197 % get the list of gaps from MILP solution 0198 metsProduced = solution.full((n_v+n_wij_IR+n_wij_R+1):(n_v+n_wij_IR+n_wij_R+n_xnp),1); 0199 allGaps = model.mets(~metsProduced); 0200 rootGaps = findRootNPmets(model,findNCgaps); %identify root gaps using findRootNPmets 0201 downstreamGaps = allGaps(~ismember(allGaps,rootGaps)); 0202 0203 0204 0205