getConcControlCoeffs

 getConcControlCoeffs
   Calculate an control coefficients of protein usage.

 Input:
   model           an ecModel in GECKO 3 format (with ecModel.ec structure)
   proteins        a list of proteins to calculate the coefficients. (Optional,
                   default = model.ec.enzymes)
   foldChange      a value how much increase the protein concentration.
                   (Optional, default = 2)
   limit           a value to determine limiting protein usage reactions.
                   Calculate as usage/concentration (Optional, default = 0)

 Output:
   enz             a logical vector of enzymes analyzed
   controlCoeffs   a vector array with the coefficients

 Usage:
    [enz, controlCoeffs] = getConcControlCoeffs(model, proteins, foldChange, limit);

0001 function [enz, controlCoeffs] = getConcControlCoeffs(model, proteins, foldChange, limit)
0002 % getConcControlCoeffs
0003 %   Calculate an control coefficients of protein usage.
0004 %
0005 % Input:
0006 %   model           an ecModel in GECKO 3 format (with ecModel.ec structure)
0007 %   proteins        a list of proteins to calculate the coefficients. (Optional,
0008 %                   default = model.ec.enzymes)
0009 %   foldChange      a value how much increase the protein concentration.
0010 %                   (Optional, default = 2)
0011 %   limit           a value to determine limiting protein usage reactions.
0012 %                   Calculate as usage/concentration (Optional, default = 0)
0013 %
0014 % Output:
0015 %   enz             a logical vector of enzymes analyzed
0016 %   controlCoeffs   a vector array with the coefficients
0017 %
0018 % Usage:
0019 %    [enz, controlCoeffs] = getConcControlCoeffs(model, proteins, foldChange, limit);
0020 
0021 if nargin < 4 || isempty(limit)
0022     limit = 0;
0023 end
0024 
0025 if nargin < 3 || isempty(foldChange)
0026     foldChange = 2;
0027 end
0028 
0029 if nargin < 2 || isempty(proteins)
0030     proteins = model.ec.enzymes;
0031 end
0032 
0033 % for now
0034 enz = false(length(proteins),1);
0035 controlCoeffs = zeros(length(proteins),1);
0036 
0037 [sol,hs] = solveLP(model);
0038 initialGrowth = sol.f;
0039 
0040 % Get enzyme index
0041 [~, protIdx] = ismember(proteins, model.ec.enzymes);
0042 
0043 % Get the protein usage reactions
0044 protUsageRxns = strcat('usage_prot_', model.ec.enzymes(protIdx));
0045 [~, protUsageRxnIdx] = ismember(protUsageRxns, model.rxns);
0046 
0047 for i = 1:numel(proteins)
0048     % Get the previous concentration
0049     prevConc = model.lb(protUsageRxnIdx(i));
0050 
0051     % Only consider those with a usage close the UB
0052     if (sol.x(protUsageRxnIdx(i))/prevConc) > limit
0053         
0054         % Update the logical vector
0055         enz(i) = 1;
0056 
0057         % Create a temporal model since coeff will be calculated one enzyme at
0058         % the time, without other change
0059         tempModel = model;
0060         % Increase the concentration by flexfactor
0061         newConc = prevConc*(foldChange);
0062         tempModel.lb(protUsageRxnIdx(i)) = newConc;
0063 
0064         % Get the new growth rate after the adjustment
0065         [tempSol,hs] = solveLP(tempModel,0,[],hs);
0066         tempGrowth = tempSol.f;
0067         
0068         % Calculate the coeff only if new growth rate is significantly
0069         % higher than initial value
0070         if (tempGrowth-initialGrowth)>1e-10
0071             controlCoeffs(i) = (tempGrowth-initialGrowth)/(prevConc-newConc);
0072         end
0073     end
0074 
0075 end
0076 end
0077