fitParameters Fits parameters such as maintenance ATP by quadratic programming model a model structure xRxns cell array with the IDs of the reactions that will be fixed for each data point xValues matrix with the corresponding values for each xRxns (columns are reactions) rxnsToFit cell array with the IDs of reactions that will be fitted to valuesToFit matrix with the corresponding values for each rxnsToFit (columns are reactions) parameterPositions stucture that determines where the parameters are in the stoichiometric matrix. Contains the fields: position cell array of vectors where each element contains the positions in the S-matrix for that parameter isNegative cell array of vectors where the elements are true if that position should be the negative of the fitted value (to differentiate between production/consumption) fitToRatio if the ratio of simulated to measured values should be fitted instead of the absolute value. Used to prevent large fluxes from having too large impact (opt, default true) initialGuess initial guess of the parameters (opt) plotFitting true if the resulting fitting should be plotted (opt, default false) parameters fitted parameters in the same order as in parameterPositions fitnessScore the correponding residual sum of squares newModel updated model structure with the fitted parameters Usage: [parameters fitnessScore exitFlag newModel]=fitParameters(model,... xRxns,xValues,rxnsToFit,valuesToFit,parameterPositions,fitToRatio,... initialGuess,plotFitting) Rasmus Agren, 2013-08-01
0001 function [parameters fitnessScore exitFlag newModel]=fitParameters(model,xRxns,xValues,rxnsToFit,valuesToFit,parameterPositions,fitToRatio,initialGuess,plotFitting) 0002 % fitParameters 0003 % Fits parameters such as maintenance ATP by quadratic programming 0004 % 0005 % model a model structure 0006 % xRxns cell array with the IDs of the reactions that will be fixed for each data point 0007 % xValues matrix with the corresponding values for each 0008 % xRxns (columns are reactions) 0009 % rxnsToFit cell array with the IDs of reactions that will be fitted to 0010 % valuesToFit matrix with the corresponding values for each 0011 % rxnsToFit (columns are reactions) 0012 % parameterPositions stucture that determines where the parameters are in the 0013 % stoichiometric matrix. Contains the fields: 0014 % position cell array of vectors where each element contains 0015 % the positions in the S-matrix for that parameter 0016 % isNegative cell array of vectors where the elements are true 0017 % if that position should be the negative of the 0018 % fitted value (to differentiate between 0019 % production/consumption) 0020 % fitToRatio if the ratio of simulated to measured values should 0021 % be fitted instead of the absolute value. Used to prevent 0022 % large fluxes from having too large impact (opt, 0023 % default true) 0024 % initialGuess initial guess of the parameters (opt) 0025 % plotFitting true if the resulting fitting should be plotted 0026 % (opt, default false) 0027 % 0028 % parameters fitted parameters in the same order as in 0029 % parameterPositions 0030 % fitnessScore the correponding residual sum of squares 0031 % newModel updated model structure with the fitted parameters 0032 % 0033 % Usage: [parameters fitnessScore exitFlag newModel]=fitParameters(model,... 0034 % xRxns,xValues,rxnsToFit,valuesToFit,parameterPositions,fitToRatio,... 0035 % initialGuess,plotFitting) 0036 % 0037 % Rasmus Agren, 2013-08-01 0038 % 0039 0040 if nargin<7 0041 fitToRatio=true; 0042 end 0043 if nargin<8 0044 initialGuess=ones(numel(parameterPositions.position),1); 0045 end 0046 if isempty(initialGuess) 0047 initialGuess=ones(numel(parameterPositions.position),1); 0048 end 0049 if nargin<9 0050 plotFitting=false; 0051 end 0052 0053 %Find the indexes of reactions that will be fitted 0054 [I rxnsToFitIndexes]=ismember(rxnsToFit,model.rxns); 0055 0056 if ~all(I) 0057 dispEM('Could not find all reactions in rxnsToFit'); 0058 end 0059 0060 %Find the indexes of reactions that will be used for constraints. 0061 [I xRxnsIndexes]=ismember(xRxns,model.rxns); 0062 0063 if ~all(I) 0064 dispEM('Could not find all reactions in xRxns'); 0065 end 0066 0067 [parameters fitnessScore exitFlag]=fminsearch(@(parameters) getRSS(parameters,model,xRxnsIndexes,xValues,rxnsToFitIndexes,valuesToFit,parameterPositions,fitToRatio),initialGuess); 0068 0069 parameters=abs(parameters); 0070 0071 if plotFitting==true 0072 %Set the resulting parameters 0073 [rss resultingFluxes newModel]=getRSS(parameters,model,xRxnsIndexes,xValues,rxnsToFitIndexes,valuesToFit,parameterPositions,true); 0074 plot(xValues,valuesToFit,'o',xValues,resultingFluxes,'-*'); 0075 end 0076 end 0077 0078 function [rss resultingFluxes newModel]=getRSS(parameters,model,xRxnsIndexes,xValues,rxnsToFitIndexes,valuesToFit,parameterPositions,fitToRatio) 0079 parameters=abs(parameters); 0080 0081 %Set the parameters at the positions specified in parameterPositions 0082 for i=1:numel(parameterPositions.position) 0083 %Set positive 0084 model.S(parameterPositions.position{i}(parameterPositions.isNegative{i}==false))=parameters(i); 0085 0086 %Set negative 0087 model.S(parameterPositions.position{i}(parameterPositions.isNegative{i}==true))=parameters(i)*-1; 0088 end 0089 0090 %Also return an updated model 0091 newModel=model; 0092 0093 %Loop through each data point, set xRxns to xValues and calculate the sum 0094 %of squares for the rxnsToFit 0095 rss=0; 0096 resultingFluxes=[]; 0097 for i=1:size(xValues,1) 0098 %Fix for more xRxns! 0099 model.lb(xRxnsIndexes)=xValues(i,:); 0100 model.ub(xRxnsIndexes)=xValues(i); 0101 0102 sol=solveLP(model); 0103 0104 %Calculate the rss 0105 if fitToRatio==false 0106 rs=sol.x(rxnsToFitIndexes)'-valuesToFit(i,:); 0107 else 0108 rs=sol.x(rxnsToFitIndexes)'./valuesToFit(i,:)-ones(1,size(valuesToFit,2)); 0109 end 0110 rss=rss+rs*rs'; 0111 resultingFluxes=[resultingFluxes sol.x(rxnsToFitIndexes)]; 0112 end 0113 end