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