applyKcatConstraints

 applyKcatConstraints
   Applies kcat-derived enzyme constraints to an ecModel. Existing enzyme
   constraints are first removed (unless updateRxns is provided), and new
   constraints are defined based on the content of model.ec.kcat.

 Input:
   model       an ecModel in GECKO 3 format (with ecModel.ec structure)
   updateRxns  if not all enzyme constraints should be updated, this can
               be given as either a logical vector of length
               model.ec.rxns, a vector of model.ec.rxns indices, or a
               (cell array of) string(s) with model.ec.rxns identifiers.
               For light models, these reactions should match model.rxns.

 Output:
   model       ecModel where reactions are constrained by enzyme usage
               if a kcat value was provided for the reaction-enzyme pair
               in model.ec.kcat

 Usage: model = applyKcatConstraints(model,updateRxns);

0001 function model = applyKcatConstraints(model,updateRxns)
0002 % applyKcatConstraints
0003 %   Applies kcat-derived enzyme constraints to an ecModel. Existing enzyme
0004 %   constraints are first removed (unless updateRxns is provided), and new
0005 %   constraints are defined based on the content of model.ec.kcat.
0006 %
0007 % Input:
0008 %   model       an ecModel in GECKO 3 format (with ecModel.ec structure)
0009 %   updateRxns  if not all enzyme constraints should be updated, this can
0010 %               be given as either a logical vector of length
0011 %               model.ec.rxns, a vector of model.ec.rxns indices, or a
0012 %               (cell array of) string(s) with model.ec.rxns identifiers.
0013 %               For light models, these reactions should match model.rxns.
0014 %
0015 % Output:
0016 %   model       ecModel where reactions are constrained by enzyme usage
0017 %               if a kcat value was provided for the reaction-enzyme pair
0018 %               in model.ec.kcat
0019 %
0020 % Usage: model = applyKcatConstraints(model,updateRxns);
0021 
0022 %these lines are for the nargin lines below only
0023 if (model.ec.geckoLight)
0024     rxns = model.rxns;
0025 else
0026     rxns = model.ec.rxns;
0027 end
0028 
0029 if nargin<2
0030     updateRxns = true(numel(rxns),1);
0031 elseif isnumeric(updateRxns)
0032     updateRxnsLog = false(numel(rxns),1);
0033     updateRxnsLog(updateRxns) = true;
0034     updateRxns = updateRxnsLog;
0035 elseif iscellstr(updateRxns) || ischar(updateRxns) || isstring(updateRxns)
0036     updateRxnsIds = convertCharArray(updateRxns);
0037     updateRxns = ismember(rxns,updateRxnsIds);
0038 end
0039  
0040 if isempty(find(updateRxns, 1)) || isempty(updateRxns)
0041      error('No reaction to update or updateRxns is logical but without any true value')
0042 end
0043 
0044 if ~isfield(model,'ec')
0045     error(['No model.ec structure could be found: the provided model is'...
0046            ' not a valid GECKO3 ecModel. First run makeEcModel(model).'])
0047 end
0048 if all(model.ec.kcat==0)
0049     printOrange('WARNING: No kcat values are provided in model.ec.kcat, model remains unchanged.\n');
0050     return
0051 end
0052 
0053 %Clear existing incorporation of enzyme usage
0054 if ~model.ec.geckoLight
0055     protMetIdx = startsWith(model.mets,'prot_') & ~strcmp(model.mets,'prot_pool');
0056     metabolRxn = unique(model.ec.rxns(updateRxns));
0057     metabolRxn = ismember(model.rxns,metabolRxn);
0058     model.S(protMetIdx,metabolRxn) = 0;
0059 end %no clearing needed for GECKO Light, will be overwritten below
0060 %For normal GECKO formulation (full model), where each enzyme is explicitly considered
0061 if ~model.ec.geckoLight 
0062     %Column 1 = rxn idx; 2 = enzyme idx; 3 = subunit copies; 4 = kcat; 5 = MW
0063     newKcats=zeros(numel(updateRxns)*10,5);
0064     updateRxns=find(updateRxns);
0065     kcatFirst=0;
0066     for i=1:numel(updateRxns)
0067         j=updateRxns(i);
0068         if model.ec.kcat(j) ~= 0
0069             enzymes   = find(model.ec.rxnEnzMat(j,:));
0070             kcatLast  = kcatFirst+numel(enzymes);
0071             kcatFirst = kcatFirst+1;
0072             newKcats(kcatFirst:kcatLast,1) = j;
0073             newKcats(kcatFirst:kcatLast,2) = enzymes;
0074             newKcats(kcatFirst:kcatLast,3) = model.ec.rxnEnzMat(j,enzymes);
0075             newKcats(kcatFirst:kcatLast,4) = model.ec.kcat(j);
0076             newKcats(kcatFirst:kcatLast,5) = model.ec.mw(enzymes);
0077             kcatFirst = kcatLast;
0078         end
0079     end
0080     if exist('kcatLast','var') % If it does not, then no kcats are found
0081         newKcats(kcatLast+1:end,:)=[];
0082 
0083         sel = newKcats(:,4) > 0; %Only apply to non-zero kcat
0084         newKcats(sel,4) = newKcats(sel,4) * 3600; %per second -> per hour
0085         newKcats(sel,4) = newKcats(sel,5) ./ newKcats(sel,4); %MW/kcat
0086         newKcats(sel,4) = newKcats(sel,3) .* newKcats(sel,4); %Multicopy subunits.
0087         newKcats(~sel,4) = 0; %Results in zero-cost
0088 
0089         %Replace rxns and enzymes with their location in model
0090         [~,newKcats(:,1)] = ismember(model.ec.rxns(newKcats(:,1)),model.rxns);
0091         [~,newKcats(:,2)] = ismember(strcat('prot_',model.ec.enzymes(newKcats(:,2))),model.mets);
0092         linearIndices     = sub2ind(size(model.S),newKcats(:,2),newKcats(:,1));
0093         model.S(linearIndices) = -newKcats(:,4); %Substrate = negative
0094     end
0095 else %GECKO light formulation, where prot_pool represents all usages
0096     prot_pool_idx = find(ismember(model.mets,'prot_pool'));
0097     %first remove the prefix of all rxns
0098     modRxns     = extractAfter(model.ec.rxns,4);
0099     % Map ec-reactions to model.rxns
0100     [hasEc,~]   = ismember(model.rxns,modRxns);
0101     hasEc       = find(hasEc & updateRxns);
0102     [~,rxnIdx]   = ismember(modRxns,model.rxns);
0103     for i = 1:numel(hasEc)
0104         % Get all isozymes per reaction
0105         ecIdx = find(rxnIdx == hasEc(i));
0106         % ecIdx = strcmp(model.rxns(hasEc(i)),modRxns);
0107         % Multiply enzymes with their MW (they are then automatically
0108         % summed per reaction), and divide by their kcat, to get a vector
0109         % of MW/kcat values.
0110         MWkcat = (model.ec.rxnEnzMat(ecIdx,:) * model.ec.mw) ./ model.ec.kcat(ecIdx);
0111         % If kcat was zero, MWkcat is Inf. If no enzyme info was found,
0112         % MWkcat is NaN. Correct both back to zero
0113         MWkcat(isinf(MWkcat) | isnan(MWkcat)) = 0;
0114         % Select the lowest MW/kcat (= most efficient), and convert to /hour
0115         model.S(prot_pool_idx, hasEc(i)) = -min(MWkcat/3600);
0116     end
0117 end
0118 end