0001 function [improvedRxns,intermediateSlns] = analyzeGCdesign(modelRed,selectedRxns,target,deletions,maxKOs,objFunction,delPenalty,intermediateSlns)
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
0034
0035 if (nargin < 5)
0036 maxKOs = 10;
0037 end
0038 if (nargin < 6)
0039 objFunction = 1;
0040 end
0041 if (nargin < 7)
0042 delPenalty = .99;
0043 end
0044 if (nargin < 8)
0045 intermediateSlns = {deletions};
0046 end
0047
0048
0049 switch objFunction
0050 case 1
0051 objectiveFunction = 'maxRate';
0052 hasSlope = false;
0053 case 2
0054 objectiveFunction = 'growth*maxRate';
0055 hasSlope = false;
0056 case 3
0057 objectiveFunction = 'maxRate*(delPenalty^numDels)';
0058 hasSlope = false;
0059 case 4
0060 objectiveFunction = 'growth*maxRate*(delPenalty^numDels)';
0061 hasSlope = false;
0062 case 5
0063 objectiveFunction = 'maxRate*(slope^(-1))';
0064 hasSlope = true;
0065 case 6
0066 objectiveFunction = 'growth*maxRate*(slope^(-1))';
0067 hasSlope = true;
0068 case 7
0069 objectiveFunction = 'maxRate*(delPenalty^numDels)*(slope^(-1))';
0070 hasSlope = true;
0071 case 8
0072 objectiveFunction = 'growth*maxRate*(delPenalty^numDels)*(slope^(-1))';
0073 hasSlope = true;
0074 end
0075
0076 if isempty(deletions)
0077 error('no knockout reactions defined')
0078 end
0079
0080 delArraySize = size(deletions);
0081 if delArraySize(1) > 1
0082 rxns = deletions';
0083 else
0084 rxns = deletions;
0085 end
0086
0087 BOF = modelRed.rxns(modelRed.c==1);
0088
0089
0090 modelKO = changeRxnBounds(modelRed,rxns,0,'b');
0091 FBAsol1 = optimizeCbModel(modelKO,'max',true);
0092 if FBAsol1.stat>0
0093 modelKOfixed = changeRxnBounds(modelKO,BOF,FBAsol1.f-1e-6,'l');
0094 modelKOfixed = changeObjective(modelKOfixed,target);
0095 FBAsol2 = optimizeCbModel(modelKOfixed,'min',true);
0096 growth = FBAsol1.f;
0097 maxRate = FBAsol2.f;
0098 numDels = length(rxns);
0099
0100 if hasSlope
0101 modelTarget = changeObjective(modelKO,target);
0102 FBAsol4 = optimizeCbModel(modelTarget,'min',true);
0103 modelTargetFixed = changeRxnBounds(modelKO,target,FBAsol4.f,'b');
0104 FBAsol5 = optimizeCbModel(modelTargetFixed,'max',true);
0105 minProdRate = FBAsol4.f;
0106 maxGrowthMinRate = FBAsol5.f;
0107
0108 if growth ~= maxGrowthMinRate
0109 slope = (maxRate-minProdRate)/(growth-maxGrowthMinRate);
0110 else
0111 slope = 1;
0112 end
0113 end
0114
0115 objective = eval(objectiveFunction);
0116
0117 bestObjective = objective
0118 bestRxns = rxns;
0119
0120 else
0121 bestObjective = 0
0122 bestRxns = rxns;
0123 end
0124
0125
0126
0127 h = waitbar(0,'improving knockout design');
0128 for i = 1:length(rxns)+1
0129 bestObjective2 = bestObjective;
0130 bestRxns2 = bestRxns;
0131 for j = 1:length(selectedRxns)+1
0132 waitbar((j+(i-1)*length(selectedRxns))/((length(rxns)+1)*(length(selectedRxns)+1)),h);
0133 newRxns = rxns;
0134 if (i==length(rxns)+1)&&(j==length(selectedRxns)+1)
0135
0136 elseif j ~= length(selectedRxns)+1
0137 newRxns{i} = selectedRxns{j};
0138 elseif i == 1
0139 newRxns = rxns(2:length(rxns));
0140 elseif i == length(rxns)
0141 newRxns = rxns(1:length(rxns)-1);
0142 else
0143 newRxns = cat(2,rxns(1:i-1),rxns(i+1:length(rxns)));
0144 end
0145
0146 if length(newRxns) <= maxKOs
0147
0148 modelKO = changeRxnBounds(modelRed,newRxns,0,'b');
0149 FBAsol1 = optimizeCbModel(modelKO,'max',true);
0150 if FBAsol1.stat>0
0151 modelKOfixed = changeRxnBounds(modelKO,BOF,FBAsol1.f-1e-6,'l');
0152 modelKOfixed = changeObjective(modelKOfixed,target);
0153 FBAsol2 = optimizeCbModel(modelKOfixed,'min',true);
0154 FBAsol3 = optimizeCbModel(modelKOfixed,'max',true);
0155 growth = FBAsol1.f;
0156 maxRate = FBAsol2.f;
0157 numDels = length(newRxns);
0158
0159 if hasSlope
0160 modelTarget = changeObjective(modelKO,target);
0161 FBAsol4 = optimizeCbModel(modelTarget,'min',true);
0162 modelTargetFixed = changeRxnBounds(modelKO,target,FBAsol4.f,'b');
0163 FBAsol5 = optimizeCbModel(modelTargetFixed,'max',true);
0164 minProdRate = FBAsol4.f;
0165 maxGrowthMinRate = FBAsol5.f;
0166
0167 if growth ~= maxGrowthMinRate
0168 slope = (maxRate-minProdRate)/(growth-maxGrowthMinRate);
0169 else
0170 slope = 1;
0171 end
0172 end
0173
0174 newObjective = eval(objectiveFunction);
0175
0176
0177 if newObjective > bestObjective2
0178 bestObjective2 = newObjective
0179 bestRxns2 = newRxns
0180 intermediateSlns{length(intermediateSlns)+1} = bestRxns2;
0181 end
0182 end
0183 end
0184 end
0185
0186 if bestObjective2 > bestObjective
0187 bestObjective = bestObjective2
0188 bestRxns = bestRxns2
0189 end
0190 end
0191 close(h);
0192
0193 bestObjective
0194 bestRxns
0195
0196
0197 if length(bestRxns) ~= length(rxns)
0198 [bestRxns,intermediateSlns] = analyzeGCdesign(modelRed,selectedRxns,target,bestRxns,maxKOs,objFunction,delPenalty,intermediateSlns);
0199 elseif length(find(strcmp(bestRxns,rxns)))~=length(rxns)
0200 [bestRxns,intermediateSlns] = analyzeGCdesign(modelRed,selectedRxns,target,bestRxns,maxKOs,objFunction,delPenalty,intermediateSlns);
0201 end
0202
0203
0204 improvedRxns = sort(bestRxns)
0205
0206
0207
0208