sigmaFitter

 sigmaFitter
   Function that fits the average enzyme saturation factor in an ecModel
   according to a provided experimentally measured value for the objective
   function (i.e. growth rate at specified conditions)

 INPUTS:
   model           an ecModel in GECKO 3 format (with ecModel.ec structure)
   growthRate      growth rate that should be reached. If not
                   specified, the value will be read from the model
                   adapter.
   Ptot            Total cellular protein content in g/gDCW. If not
                   specified, the value will be read from the model
                   adapter. If not specified in model adapter, 0.5 g/gDCW
                   is assumed.
   f               Estimated fraction of enzymes in the model. If not
                   specified, the value will be read from the model
                   adapter. If not specified in model adapter, 0.5 is
                   assumed.
   makePlot        Logical whether a plot should be made. Default true.
   modelAdapter    a loaded model adapter (Optional, will otherwise use the
                   default model adapter).

 Output:
   model           ecModel with protein pool exchange upper bound adapted
                   to the optimal sigma-factor
   sigma           optimal sigma-factor

 Usage:
   [model, sigma] = sigmaFitter(model, growthRate, Ptot, f, makePlot, modelAdapter)

0001 function [model, sigma] = sigmaFitter(model, growthRate, Ptot, f, makePlot, modelAdapter)
0002 % sigmaFitter
0003 %   Function that fits the average enzyme saturation factor in an ecModel
0004 %   according to a provided experimentally measured value for the objective
0005 %   function (i.e. growth rate at specified conditions)
0006 %
0007 % INPUTS:
0008 %   model           an ecModel in GECKO 3 format (with ecModel.ec structure)
0009 %   growthRate      growth rate that should be reached. If not
0010 %                   specified, the value will be read from the model
0011 %                   adapter.
0012 %   Ptot            Total cellular protein content in g/gDCW. If not
0013 %                   specified, the value will be read from the model
0014 %                   adapter. If not specified in model adapter, 0.5 g/gDCW
0015 %                   is assumed.
0016 %   f               Estimated fraction of enzymes in the model. If not
0017 %                   specified, the value will be read from the model
0018 %                   adapter. If not specified in model adapter, 0.5 is
0019 %                   assumed.
0020 %   makePlot        Logical whether a plot should be made. Default true.
0021 %   modelAdapter    a loaded model adapter (Optional, will otherwise use the
0022 %                   default model adapter).
0023 %
0024 % Output:
0025 %   model           ecModel with protein pool exchange upper bound adapted
0026 %                   to the optimal sigma-factor
0027 %   sigma           optimal sigma-factor
0028 %
0029 % Usage:
0030 %   [model, sigma] = sigmaFitter(model, growthRate, Ptot, f, makePlot, modelAdapter)
0031 
0032 if nargin < 6 || isempty(modelAdapter)
0033     modelAdapter = ModelAdapterManager.getDefault();
0034     if isempty(modelAdapter)
0035         error('Either send in a modelAdapter or set the default model adapter in the ModelAdapterManager.')
0036     end
0037 end
0038 
0039 if nargin<5 || isempty(makePlot)
0040     makePlot = true;
0041 end
0042 if nargin<4 || isempty(f)
0043     f = modelAdapter.getParameters().f;
0044 end
0045 if nargin<3 || isempty(Ptot)
0046     Ptot = modelAdapter.getParameters().Ptot;
0047 end
0048 if nargin<2 || isempty(growthRate)
0049     growthRate = modelAdapter.getParameters().gR_exp;
0050 end
0051 
0052 objValues = [];
0053 errors    = [];
0054 sigParam  = [];
0055 objPos    = find(model.c);
0056 %Relax bounds for the objective function
0057 model.lb(objPos) = 0;
0058 model.ub(objPos) = 1000;
0059 hsSol=[];
0060 for i=1:100
0061     %Constrains the ecModel with the i-th sigma factor
0062     sigma = i/100;
0063     model = setProtPoolSize(model, Ptot, f, sigma, modelAdapter);
0064     [solution, hsSol]  = solveLP(model,0,[],hsSol);
0065     if isempty(solution.x)
0066         solution.x=zeros(length(model.rxns),1);
0067     end
0068     objValues = [objValues; solution.x(objPos)];
0069     error     = abs(((growthRate-solution.x(objPos))/growthRate)*100);
0070     errors    = [errors; error];
0071     error     = num2str(((growthRate-solution.x(objPos))/growthRate)*100);
0072     sigParam  = [sigParam; sigma];
0073 end
0074 [~, minIndx] = min(errors);
0075 sigma     = sigParam(minIndx);
0076 if makePlot
0077     figure
0078     plot(sigParam,errors,'LineWidth',5)
0079     title('Sigma fitting')
0080     xlabel('Average enzyme saturation [-]')
0081     ylabel('Absolute relative error [%]')
0082 end
0083 end