0001 function [x, population, scores, optGeneSol] = optGene(model, targetRxn, substrateRxn, generxnList, MaxKOs, population)
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 global HTABLE
0028 HTABLE = java.util.Hashtable;
0029 global MaxKnockOuts
0030
0031 ngenes = length(generxnList);
0032
0033
0034 rxnok = 1;
0035 geneok = 1;
0036 for i = 1:length(generxnList)
0037 if(~ ismember(generxnList{i}, model.rxns)), rxnok = 0; end
0038 if(~ ismember(generxnList{i}, model.genes)),geneok = 0; end
0039 end
0040 if geneok
0041 display('assuming list is genes');
0042 elseif rxnok
0043 display('assuming list is reactions');
0044 else
0045 display('list appears to be neither genes nor reactions: aborting');
0046 return;
0047 end
0048
0049
0050
0051
0052 if nargin < 5
0053 MaxKnockOuts = 10;
0054 else
0055 MaxKnockOuts = MaxKOs;
0056 end
0057 mutationRate = 1/ngenes;
0058 crossovermutationRate = mutationRate*.2;
0059 CrossoverFraction = .80;
0060 PopulationSize = [125 125 125 125];
0061 Generations = 10000;
0062 TimeLimit = 3600*24*2;
0063 StallTimeLimit = 3600*24*1;
0064 StallGenLimit = 10000;
0065 PlotFcns = {@gaplotscores, @gaplotbestf, @gaplotscorediversity, @gaplotstopping, @gaplotmutationdiversity};
0066 crossfun = @(a,b,c,d,e,f) crossoverCustom(a,b,c,d,e,f,crossovermutationRate);
0067 MigrationFraction = .1;
0068 MigrationInterval = 100;
0069
0070
0071
0072
0073
0074 InitialPopulation = [];
0075 if nargin > 5
0076 InitialPopulation = double(population);
0077 end
0078 options = gaoptimset( ...
0079 'PopulationType', 'bitstring', ...
0080 'CreationFcn', @lowmutationcreation, ...
0081 'MutationFcn', {@mutationUniformEqual, mutationRate}, ...
0082 'PopulationSize', PopulationSize, ...
0083 'StallTimeLimit', StallTimeLimit, ...
0084 'TimeLimit', TimeLimit, ...
0085 'PlotFcns', PlotFcns, ...
0086 'InitialPopulation', InitialPopulation, ...
0087 'CrossoverFraction', CrossoverFraction, ...
0088 'CrossoverFcn', crossfun, ...
0089 'StallGenLimit', StallGenLimit, ...
0090 'Generations', Generations, ...
0091 'TolFun', 1e-10, ...
0092 'Vectorize', 'on', ...
0093 'MigrationFraction', MigrationFraction, ...
0094 'MigrationInterval', MigrationInterval ...
0095 ...
0096 );
0097
0098
0099
0100
0101
0102
0103 FitnessFunction = @(x) optGeneFitnessTilt(x,model,targetRxn, generxnList, geneok);
0104
0105 gap.fitnessfcn = FitnessFunction;
0106 gap.nvars = ngenes;
0107 gap.options = options;
0108
0109 [x,FVAL,REASON,OUTPUT,population, scores] = ga(gap);
0110
0111
0112 [optGeneSol] = GetOptGeneSol(model, targetRxn, substrateRxn, generxnList, population, x, scores, geneok);
0113
0114 return;
0115
0116
0117
0118
0119
0120 function [Population] = lowmutationcreation(GenomeLength,FitnessFcn,options)
0121 totalPopulation = sum(options.PopulationSize);
0122 initPopProvided = size(options.InitialPopulation,1);
0123 individualsToCreate = totalPopulation - initPopProvided;
0124
0125
0126 Population = true(totalPopulation,GenomeLength);
0127
0128 if initPopProvided > 0
0129 Population(1:initPopProvided,:) = options.InitialPopulation;
0130 end
0131
0132 Population(initPopProvided+1:end,:) = logical(1/GenomeLength > rand(individualsToCreate,GenomeLength));
0133 return;
0134
0135
0136
0137
0138
0139 function mutationChildren = mutationUniformEqual(parents,options,GenomeLength,FitnessFcn,state,thisScore,thisPopulation,mutationRate)
0140 global MaxKnockOuts
0141 if(nargin < 8)
0142 mutationRate = 0.01;
0143 end
0144 mutationChildren = zeros(length(parents),GenomeLength);
0145 for i=1:length(parents)
0146 child = thisPopulation(parents(i),:);
0147 kos = sum(child);
0148 mutationPoints = find(rand(1,length(child)) < mutationRate);
0149 child(mutationPoints) = ~child(mutationPoints);
0150
0151 if MaxKnockOuts > 0
0152 while(sum(child(:))> MaxKnockOuts)
0153 ind2 = find(child);
0154 removeindex = ind2(randint(1,1,length(ind2))+1);
0155 child(removeindex) = 0;
0156 end
0157 end
0158
0159
0160
0161 if rand > .5 && kos > 1
0162 while(sum(child(:))>= kos)
0163 ind2 = find(child);
0164 removeindex = ind2(randint(1,1,length(ind2))+1);
0165 child(removeindex) = 0;
0166 end
0167 end
0168 mutationChildren(i,:) = child;
0169 end
0170 return;
0171
0172
0173
0174
0175
0176 function xoverKids = crossoverCustom(parents,options,GenomeLength,FitnessFcn,unused,thisPopulation, mutationRate)
0177 nKids = length(parents)/2;
0178
0179
0180 xoverKids = zeros(nKids,GenomeLength);
0181
0182 global MaxKnockOuts;
0183
0184
0185
0186 index = 1;
0187
0188
0189 for i=1:nKids
0190
0191 r1 = parents(index);
0192 index = index + 1;
0193 r2 = parents(index);
0194 index = index + 1;
0195
0196
0197
0198
0199 for j = 1:GenomeLength
0200 if(rand > 0.5)
0201 xoverKids(i,j) = thisPopulation(r1,j);
0202 else
0203 xoverKids(i,j) = thisPopulation(r2,j);
0204 end
0205 end
0206 if MaxKnockOuts>0
0207 while(sum(xoverKids(i,:))> MaxKnockOuts)
0208 ind2 = find(xoverKids(i,:));
0209 removeindex = ind2(randint(1,1,length(ind2))+1);
0210 xoverKids(i,removeindex) = 0;
0211 end
0212 end
0213 end
0214
0215 xoverKids = mutationUniformEqual(1:size(xoverKids,1) ,[],GenomeLength,[],[],[],xoverKids,mutationRate);
0216 return;
0217
0218
0219 function state = gaplotmutationdiversity(options,state,flag,p1)
0220
0221
0222
0223
0224
0225
0226
0227
0228
0229
0230
0231 global MaxKnockOuts
0232 if nargin < 4
0233 p1 = 10;
0234 end
0235
0236 p1 = MaxKnockOuts+1;
0237 p1 = [0:(MaxKnockOuts)];
0238 switch flag
0239 case 'init'
0240 title('Mutation Histogram','interp','none')
0241 xlabel('number of mutations');
0242 ylabel('Number of individuals');
0243 case 'iter'
0244
0245 if size(state.Score,2) > 1 && isfield(state,'Rank')
0246 index = (state.Rank == 1);
0247
0248
0249 if nnz(index) > 1
0250 set(gca,'ylimmode','auto');
0251 hist(sum(state.Population(index,:)),p1);
0252 else
0253 set(gca,'ylim',[0 1]);
0254 hist([sum(state.Population(index,:)); sum(state.Population(index,:))],p1);
0255 end
0256
0257 nObj = size(state.Score,2);
0258 fminval = min(state.Score(index,:),[],1);
0259 fmaxval = max(state.Score(index,:),[],1);
0260 legendText = cell(1,nObj);
0261 for i = 1:nObj
0262 legendText{i} = ['fun',num2str(i),' [',sprintf('%g ',fminval(i)), ...
0263 sprintf('%g',fmaxval(i)),']'];
0264 end
0265 legend(legendText);
0266 else
0267 hist(sum(state.Population,2),p1);
0268 end
0269 end
0270