Description of abc

0001 function [rmse_final,exp,simulated,growthdata,max_growth]=abc_max(ecModel,kcat_random_all,growthdata,max_growth,proc,sample_generation,j,rxn2block)
0002 % setProtPoolSize
0003 %   Internal function in BayesianSensitivityTuning. Gets the average RMSE for
0004 %   a certain set of growth data and kcats.
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 %   kcat_random_all Sampled kcat values
0011 %   growthdata      Growth data loaded from file.
0012 %   max_growth      Max growth data loaded from file.
0013 %   proc            Number of experiments with data
0014 %   sample_generation Total number of randomly generated kcat sets in this iteration
0015 %   j               Which experiment to use
0016 %   rxn2block       List of reactions to block, could be read from file.
0017 %
0018 % Output:
0019 %   rmse_final      The average RMSE
0020 %   exp             Metabolite export
0021 %   simulated       Result from simulation
0022 %   growthdata      The growth data
0023 %   max_growth      The max growth data (without input constraints)
0024 
0025 nstep = sample_generation/proc;
0026 rmse_final = zeros(1,nstep);
0027 kcat_sample = kcat_random_all(:,(j-1)*nstep+1:j*nstep);
0028 
0029 
0030 % get carbonnum for each exchange rxn to further calculation of error
0031 if ~isfield(ecModel,'excarbon')
0032     ecModel = addCarbonNum(ecModel);
0033 end
0034 
0035 for k = 1:nstep
0036     %disp(['nstep:',num2str(k),'/',num2str(nstep)])
0037     kcat_random  = kcat_sample(:,k);
0038     ecModel.ec.kcat = kcat_random;
0039     ecModel = applyKcatConstraints(ecModel);
0040     
0041     %% first search with substrate constraints
0042     objective = 'r_2111';
0043     if ~isempty(growthdata)
0044         [rmse_1,exp_1,simulated_1] = rmsecal(ecModel,growthdata,true,objective,rxn2block);
0045     else
0046         rmse_1 = [];
0047         exp_1 = [];
0048         simulated_1 = [];
0049     end
0050     %% second search for maxmial growth rate without constraints
0051     if ~isempty(max_growth)  % simulate the maximal growth rate
0052         [rmse_2,exp_2,simulated_2] = rmsecal(ecModel,max_growth,false,objective,rxn2block);
0053     else
0054         rmse_2 = [];
0055         exp_2 = [];
0056         simulated_2 = [];
0057     end
0058     exp = [exp_1;exp_2];
0059     simulated = [simulated_1;simulated_2];
0060     rmse_final(1,k) = mean([rmse_1,rmse_2],'omitnan');
0061     
0062     %% only output simulated result for one generation
0063     if nstep ~= 1 || sample_generation ~= 1
0064         simulated = [];
0065         exp = [];
0066     end
0067 end
0068 end
0069 
0070 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0071 function [rmse,exp,simulated] = rmsecal(ecModel,data,constrain,objective,rxn2block)
0072 
0073     simulated = zeros(length(data(:,1)),9);
0074     
0075     for i = 1:length(data(:,1))
0076         exp = cell2mat(data(:,3:11)); % u sub ace eth gly pyr ethyl_acetate co2 o2
0077         exp = exp.*[1,-1,1,1,1,1,1,1,-1];
0078         ex_mets = {'growth',[data{i,2},' exchange'],'acetate exchange','ethanol exchange','glycerol exchange','pyruvate exchange','ethyl acetate exchange','carbon dioxide exchange','oxygen exchange'};
0079         [~,idx] = ismember(ex_mets,ecModel.rxnNames);
0080         model_tmp = ecModel;
0081 
0082         model_tmp = changeMedia(model_tmp,'D-glucose',data(i,16));%TODO: These are currently hard-coded for yeast-GEM, add to adapter
0083         %Stop import/export of acetate and acetaldehyde
0084         model_tmp = changeRxnBounds(model_tmp,'r_1634',0,'b'); %TODO: These are currently hard-coded for yeast-GEM, add to adapter
0085         model_tmp = changeRxnBounds(model_tmp,'r_1631',0,'b');
0086 
0087         if strcmp(data(i,14),'anaerobic') ||strcmp(data(i,14),'limited') 
0088             model_tmp = anaerobicModel_GECKO(model_tmp);
0089         end
0090         if strcmp(data(i,14),'limited') 
0091              model_tmp.lb(strcmp(model_tmp.rxnNames,'oxygen exchange')) = -5;%TODO: Currently hard-coded for yeast
0092         end
0093         if ~constrain
0094             %No export of glucose
0095             model_tmp.lb(strcmp(model_tmp.rxns,'r_1714')) = 0; %TODO: Currently hard-coded for yeast
0096             model_tmp.lb(strcmp(model_tmp.rxns,ecModel.rxns(idx(2)))) = -1000; % not constrain the substrate usage
0097         else
0098             %No export of glucose
0099             model_tmp.lb(strcmp(model_tmp.rxns,'r_1714')) = 0;%TODO: Currently hard-coded for yeast
0100             if isnan(exp(i,2))
0101                 model_tmp.lb(idx(2)) = -1000;
0102             else
0103                 model_tmp.lb(idx(2)) = exp(i,2);
0104             end
0105         end
0106 
0107 
0108         model_tmp.c = double(strcmp(model_tmp.rxns, objective));
0109         sol_tmp = solveLP(model_tmp);%,objective,osenseStr,prot_cost_info,tot_prot_weight,'ibm_cplex');
0110         if checkSolution(sol_tmp)
0111             sol(:,i) = sol_tmp.x;
0112 
0113             tmp = ~isnan(exp(i,:));
0114             excarbon = ecModel.excarbon(idx);
0115             excarbon(excarbon == 0) = 1;
0116             exp_tmp = exp(i,tmp).*excarbon(tmp);
0117             simulated_tmp = sol(idx(tmp),i)'.*excarbon(tmp); % normalize the growth rate issue by factor 10
0118 
0119             exp_block = zeros(1,length(setdiff(rxn2block,model_tmp.rxns(idx(2))))); % all zeros for blocked exchange mets exchange
0120             rxnblockidx = ismember(model_tmp.rxns,setdiff(rxn2block,model_tmp.rxns(idx(2))));
0121             simulated_block = sol(rxnblockidx,i)'.* ecModel.excarbon(rxnblockidx); %
0122             exp_block = exp_block(simulated_block~=0);
0123             simulated_block = simulated_block(simulated_block~=0);
0124             if constrain
0125                 rmse_tmp(i) = sqrt(immse([exp_tmp,exp_block], [simulated_tmp,simulated_block]));
0126             else
0127                 if length(exp_tmp) >= 2
0128                     rmse_tmp(i) = sqrt(immse(exp_tmp(1:2), simulated_tmp(1:2)));
0129                 else
0130                     rmse_tmp(i) = sqrt(immse(exp_tmp(1), simulated_tmp(1)));
0131                 end
0132 
0133             end
0134             simulated(i,:) = sol(idx,i)';
0135         else
0136             simulated(i,:) = NaN;
0137             rmse_tmp(i) = NaN;
0138         end
0139     end
0140     rmse_tmp(isnan(rmse_tmp)) = []; %we just skip any case without solution, they are pretty rare, but exist
0141     rmse = sum(rmse_tmp)/length(data(:,1));
0142 end
abc_max

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
