bayesianSensitivityTuning

 bayesianSensitivityTuning
   Modifies kcats to better fit experimental data. Currently only works with 
   light models.

 Input:
   ecModel         ecModel that was generated by makeEcModel, or loaded from
                   an earlier run. Not compatible with ecModels generated by
                   earlier GECKO versions (pre 3.0).
   modelAdapter    a loaded model adapter (Optional, will otherwise use the
                   default model adapter).
   maxIterations   The maximum number of iterations to run (Optional, default 150)

 Output:
   model           ecModel with new kcats, not applied.

0001 function ecModel = bayesianSensitivityTuning(ecModel,modelAdapter,maxIterations)
0002 % bayesianSensitivityTuning
0003 %   Modifies kcats to better fit experimental data. Currently only works with
0004 %   light models.
0005 %
0006 % Input:
0007 %   ecModel         ecModel that was generated by makeEcModel, or loaded from
0008 %                   an earlier run. Not compatible with ecModels generated by
0009 %                   earlier GECKO versions (pre 3.0).
0010 %   modelAdapter    a loaded model adapter (Optional, will otherwise use the
0011 %                   default model adapter).
0012 %   maxIterations   The maximum number of iterations to run (Optional, default 150)
0013 %
0014 % Output:
0015 %   model           ecModel with new kcats, not applied.
0016 %
0017 
0018 if nargin < 3
0019     maxIterations = 150;
0020 end
0021 
0022 if nargin < 2 || isempty(modelAdapter)
0023     modelAdapter = ModelAdapterManager.getDefault();
0024     if isempty(modelAdapter)
0025         error('Either send in a modelAdapter or set the default ecModel adapter in the ModelAdapterManager.')
0026     end
0027 end
0028 
0029 basePath = modelAdapter.getParameters().path;
0030 
0031 growthdata = table2cell(readtable(fullfile(basePath,'data','BayesianGrowthRates.csv')));
0032 max_growth = table2cell(readtable(fullfile(basePath,'data','BayesianMaxGrowth.csv')));
0033 rxn2block = table2cell(readtable(fullfile(basePath,'data','BayesianRxn2Block.csv'), 'Delimiter', ','));
0034 
0035 proc = 18;
0036 numPerGeneration = 126;
0037 rejectnum = 0.2;
0038 
0039 %these will be updated in the loop below
0040 kcats = ecModel.ec.kcat;
0041 kcat_var = ones(length(ecModel.ec.kcat),1);
0042 
0043 D = abc_max(ecModel,kcats,growthdata,max_growth,1,1,1,rxn2block);
0044 D_100 = D;
0045 theta_100 = [];
0046 kcat_100 = [];
0047 sampledgeneration = 1;
0048 while D > rejectnum && D_100 > 0.5 % range for D_1 to D_100 is within 0.3 and D_100 should < 0.5
0049     if sampledgeneration <= maxIterations
0050         disp(['Running ' num2str(sampledgeneration) ' of ' num2str(maxIterations) ': D = ' num2str(D) ' D_100 = ' num2str(D_100)])
0051         old = theta_100;
0052         kcat_old_100 = kcat_100;
0053 
0054         if sampledgeneration == 1
0055             sample_generation = 144;
0056         else
0057             sample_generation = numPerGeneration;
0058         end
0059         % generate one maxIterations sample of kcats
0060         kcat_random_all = arrayfun(@getrSample, kcats, kcat_var,repmat(sample_generation,length(kcats),1),'UniformOutput',false);
0061         kcat_random_all = cell2mat(kcat_random_all);
0062         %Simulate and measure RMSE
0063         parfor j = 1:proc
0064             rmse_final = abc_max(ecModel,kcat_random_all,growthdata,max_growth,proc,sample_generation,j,rxn2block)
0065             new_tmp{j} = rmse_final;
0066         end
0067         
0068         new = cell2mat(new_tmp);
0069         theta = [new,old];
0070         kcat = [kcat_random_all,kcat_old_100];
0071 
0072         % Initialize an empty set to store the best 100  after each step
0073         [~,D_idx]= sort(theta,'ascend');
0074         theta_100 = theta(D_idx(1:100));
0075         D = abs(theta_100(100)-theta_100(1)); % the largest one theta - smallest theta
0076         D_100 = theta_100(100);
0077         kcat_100 = kcat(:,D_idx(1:100));
0078         %writematrix([kcat_100;repmat(tot_prot_weight,1,length(theta_100));theta_100],['kcat_genra',num2str(sampledgeneration),'.txt'])
0079 
0080         % recalculate the sigma and mu
0081         ss = num2cell(kcat_100',1);
0082         [a,b] = arrayfun(@updateprior,ss);
0083         kcats = a';
0084         kcat_var = b';
0085         sampledgeneration = sampledgeneration + 1;
0086     else
0087         D = rejectnum;
0088         D_100 = D;
0089     end
0090 end
0091 
0092 ecModel.ec.kcat = kcats;
0093 
0094 end
0095 %code for generating data files from the .mat files in dlkcat
0096 %x = load('C:/Code/Components/GECKO/Gecko3/GECKO/tutorials/tutorial_yeast-GEM/data/Saccharomyces_cerevisiae_dl.mat')
0097 %growthrates = x.growthrates;
0098 %growthrates(strcmp(growthrates(:,1),'Saccharomyces_cerevisiae'),:)
0099 %size(growthrates)
0100 %T = cell2table(growthrates(strcmp(growthrates(:,1),'Saccharomyces_cerevisiae'),:),'VariableNames',{'Species','Substrate','Uptake','sub','ace','eth','gly','pyr','ethyl_acetate','co2','o2','Unused1','Unused2','OxAvail','Unused3','Media'});
0101 %writetable(T,'C:/Code/Components/GECKO/Gecko3/GECKO/tutorials/tutorial_yeast-GEM/data/BayesianGrowthRates.csv');
0102 %maxGrowth = cell2table(x.max_growth, 'VariableNames', {'Species','Substrate','Uptake','sub','ace','eth','gly','pyr','ethyl_acetate','co2','o2','Unused1','Unused2','OxAvail','Unused3','Media'});
0103 %writetable(maxGrowth, 'C:/Code/Components/GECKO/Gecko3/GECKO/tutorials/tutorial_yeast-GEM/data/BayesianMaxGrowth.csv')
0104 %rxn2block = x.rxn2block;
0105 %T = cell2table(rxn2block,'VariableNames',{'Rxns'});
0106 %writetable(T,'C:/Code/Components/GECKO/Gecko3/GECKO/tutorials/tutorial_yeast-GEM/data/BayesianRxn2Block.csv');
0107 
0108