sensitivityTuning

 sensitivityTuning
    Function that relaxes the most limiting kcats until a certain growth rate
    is reached. The function will update kcats in model.ec.kcat.

 Input:
   model              an ecModel in GECKO 3 format (with ecModel.ec structure)
   desiredGrowthRate  kcats will be relaxed until this growth rate is reached
   modelAdapter       a loaded model adapter (Optional, will otherwise use the
                      default model adapter).
   foldChange         kcat values will be increased by this fold-change.
                      (Opt, default 10)
   protToIgnore       vector of protein ids to be ignore in tuned kcats.
                      e.g. {'P38122', 'Q99271'} (Optional, default = [])
   verbose            logical whether progress should be reported (Optional,
                      default true)

 Output:
   model              ecModel with updated model.ec.kcat
   tunedKcats         structure with information on tuned kcat values
                      rxns     identifiers of reactions with tuned kcat
                               values
                      rxnNames names of the reactions in tunedKcats.rxns
                      enzymes  enzymes that catalyze the reactions in
                               tunedKcats.rxns, whose kcat value has been
                               tuned.
                      oldKcat  kcat values in the input model
                      newKcat  kcat values in the output model, after tuning

 Note: The model.ec.notes field will contain the original kcat value and
 source, unless the kcat has previously been set by sensitivityTuning, in
 which case the notes field remains unchanged.

 Usage:
   [model, tunedKcats] = sensitivityTuning(model, desiredGrowthRate, modelAdapter, foldChange, protToIgnore, verbose)

0001 function [model, tunedKcats] = sensitivityTuning(model, desiredGrowthRate, modelAdapter, foldChange, protToIgnore, verbose)
0002 % sensitivityTuning
0003 %    Function that relaxes the most limiting kcats until a certain growth rate
0004 %    is reached. The function will update kcats in model.ec.kcat.
0005 %
0006 % Input:
0007 %   model              an ecModel in GECKO 3 format (with ecModel.ec structure)
0008 %   desiredGrowthRate  kcats will be relaxed until this growth rate is reached
0009 %   modelAdapter       a loaded model adapter (Optional, will otherwise use the
0010 %                      default model adapter).
0011 %   foldChange         kcat values will be increased by this fold-change.
0012 %                      (Opt, default 10)
0013 %   protToIgnore       vector of protein ids to be ignore in tuned kcats.
0014 %                      e.g. {'P38122', 'Q99271'} (Optional, default = [])
0015 %   verbose            logical whether progress should be reported (Optional,
0016 %                      default true)
0017 %
0018 % Output:
0019 %   model              ecModel with updated model.ec.kcat
0020 %   tunedKcats         structure with information on tuned kcat values
0021 %                      rxns     identifiers of reactions with tuned kcat
0022 %                               values
0023 %                      rxnNames names of the reactions in tunedKcats.rxns
0024 %                      enzymes  enzymes that catalyze the reactions in
0025 %                               tunedKcats.rxns, whose kcat value has been
0026 %                               tuned.
0027 %                      oldKcat  kcat values in the input model
0028 %                      newKcat  kcat values in the output model, after tuning
0029 %
0030 % Note: The model.ec.notes field will contain the original kcat value and
0031 % source, unless the kcat has previously been set by sensitivityTuning, in
0032 % which case the notes field remains unchanged.
0033 %
0034 % Usage:
0035 %   [model, tunedKcats] = sensitivityTuning(model, desiredGrowthRate, modelAdapter, foldChange, protToIgnore, verbose)
0036 
0037 if nargin < 6 || isempty(verbose)
0038     verbose = true;
0039 end
0040 if nargin < 5 || isempty(protToIgnore)
0041     protToIgnore = {};
0042 end
0043 if nargin < 4 || isempty(foldChange)
0044     foldChange = 10;
0045 end
0046 if nargin < 3 || isempty(modelAdapter)
0047     modelAdapter = ModelAdapterManager.getDefault();
0048     if isempty(modelAdapter)
0049         error('Either send in a modelAdapter or set the default model adapter in the ModelAdapterManager.')
0050     end
0051 end
0052 params = modelAdapter.params;
0053 if nargin < 2 || isempty(desiredGrowthRate)
0054     desiredGrowthRate = params.gR_exp;
0055 end
0056 
0057 kcatList = [];
0058 m = model;
0059 m.c = double(strcmp(m.rxns, params.bioRxn));% Make sure that growth is maximized
0060 
0061 [res,hs] = solveLP(m);
0062 if isempty(res.x)
0063     error('FBA of input model gives no valid result. Reduce protein pool constraint with setProtPoolSize and/or check if exchange constraints are correctly defined.')
0064 end
0065 lastGrowth = 0;
0066 if ~m.ec.geckoLight
0067     %for the full model, we first find the draw reaction with the most flux
0068     drawRxns = startsWith(m.rxns, 'usage_prot_');
0069     idxToIgnore = cellfun(@(x) find(strcmpi(model.rxns, ['usage_prot_' x])), protToIgnore);
0070     iteration = 1;
0071     while true
0072         [res,hs] = solveLP(m,0,[],hs); %skip parsimonius, only takes time
0073         if (lastGrowth == res.f)
0074             printOrange('WARNING: No growth increase from increased kcats - check if the constraints on the uptake reactions are too tight.\n');
0075             break;
0076         end
0077         lastGrowth = res.f;
0078         if verbose; disp(['Iteration ' num2str(iteration) ': Growth: ' num2str(lastGrowth)]); end
0079         if (lastGrowth >= desiredGrowthRate)
0080             break;
0081         end
0082         %If you get an error here, it is likely due to numerical issues in the solver
0083         %The trick where we don't allow low kcats is to fix that, but maybe
0084         %it is not enough.
0085         iteration            = iteration + 1;
0086         %find the highest draw_prot rxn flux
0087         drawFluxes           = zeros(length(drawRxns),1);
0088         drawFluxes(drawRxns) = res.x(drawRxns);
0089         % Remove from the list user defined proteins
0090         drawFluxes(idxToIgnore) = 0;
0091         [~,sel]              = min(drawFluxes); % since bounds -1000 to 0
0092         %Now get the metabolite
0093         metSel               = m.S(:,sel) < 0; % negative coeff
0094         %now find the reaction with the largest consumption of this protein
0095         protFluxes           = m.S(metSel,:).' .* res.x; %negative
0096         [~,rxnSel]           = min(protFluxes);
0097         kcatList             = [kcatList, rxnSel];
0098         rxn                  = m.rxns(rxnSel);
0099         targetSubRxn         = strcmp(m.ec.rxns, rxn);
0100         if ~strcmp(m.ec.source(targetSubRxn),'sensitivityTuning')
0101             oldNote          = m.ec.notes{targetSubRxn};
0102             newNote          = ['preTuneKcat=' num2str(m.ec.kcat(targetSubRxn)) ' | source:' m.ec.source{targetSubRxn}];
0103             if ~isempty(oldNote)
0104                 newNote      = [oldNote '; ' newNote];
0105             end
0106             m.ec.notes{targetSubRxn}    = newNote;
0107         end
0108         m.ec.kcat(targetSubRxn)     = m.ec.kcat(targetSubRxn) .* foldChange;
0109         m.ec.source(targetSubRxn)   = {'sensitivityTuning'};
0110         m                    = applyKcatConstraints(m,targetSubRxn);
0111     end
0112 
0113 else
0114     origRxns = extractAfter(m.ec.rxns,4);
0115     %find the reactions involved in proteins to be ignored
0116     idxToIgnore = cellfun(@(x) find(m.ec.rxnEnzMat(:, strcmpi(m.ec.enzymes, x))), protToIgnore, 'UniformOutput', false);
0117     %create an unique vector
0118     idxToIgnore = unique(cat(1, idxToIgnore{:}));
0119     %get the correct idx in model.rxns
0120     idxToIgnore = cellfun(@(x) find(strcmpi(m.rxns, x)), origRxns(idxToIgnore));
0121     iteration = 1;
0122     while true
0123         res = solveLP(m,0); %skip parsimonius, only takes time
0124         if (lastGrowth == res.f)
0125             printOrange('No growth increase from increased kcats - check if the constraints on the uptake reactions are too tight.\n');
0126             break;
0127         end
0128         lastGrowth = res.f;
0129         if (lastGrowth >= desiredGrowthRate)
0130             break;
0131         end
0132         %If you get an error here, it is likely due to numerical issues in the solver
0133         %The trick where we don't allow low kcats is to fix that, but maybe
0134         %it is not enough.
0135         if verbose; disp(['Iteration ' num2str(iteration) ': Growth: ' num2str(lastGrowth)]); end
0136         iteration       = iteration + 1;
0137         %find the highest protein usage flux
0138         protPoolStoich  = m.S(strcmp(m.mets, 'prot_pool'),:).';
0139         protPoolStoich(idxToIgnore) = 0;
0140         [~,sel]         = min(res.x .* protPoolStoich); %max consumption
0141         kcatList        = [kcatList, sel];
0142         rxn             = m.rxns(sel.');
0143         targetSubRxns   = strcmp(origRxns, rxn);
0144         m.ec.kcat(targetSubRxns) = m.ec.kcat(targetSubRxns) .* foldChange;
0145         m               = applyKcatConstraints(m,rxn);
0146     end
0147 end
0148 
0149 kcatList        = unique(kcatList);
0150 tunedKcats.rxns = m.rxns(kcatList);
0151 tunedKcats.rxnNames = m.rxnNames(kcatList);
0152 if ~model.ec.geckoLight
0153     [~, rxnIdx]     = ismember(tunedKcats.rxns,m.ec.rxns);
0154 else
0155     [~, rxnIdx]     = ismember(tunedKcats.rxns,origRxns);
0156 end
0157 tunedKcats.enzymes  = cell(numel(kcatList),1);
0158 for i=1:numel(rxnIdx)
0159     [~, metIdx]     = find(m.ec.rxnEnzMat(rxnIdx(i),:));
0160     tunedKcats.enzymes{i}   = strjoin(m.ec.enzymes(metIdx),';');
0161 end
0162 tunedKcats.oldKcat  = model.ec.kcat(rxnIdx);
0163 tunedKcats.newKcat  = m.ec.kcat(rxnIdx);
0164 tunedKcats.source   = model.ec.source(rxnIdx);
0165 
0166 model = m;
0167 end