flexibilizeEnzConcs
Flexibilize enzyme concentration of an ecModel with constrained with
proteomics data. The upper bound of the protein usage reaction is
changed, while the concentrations in ecModel.ec.concs remain unchanged.
If no (more) limiting enzyme concentrations can be found, it might be
the protein pool exchange that is limiting growth. In that case, an
attempt will be made to relax the protein pool exchange reaction, and
if the growth rate indeed increases, it is suggested to set the protein
pool exchange unconstrained (lb=-1000) before again running
flexibilizeEnzConcs. Such situations, where a proteomics integrated
ecModel is overconstrained, may occur if the ecModel should be able to
simulate maximum growth rate (from e.g. batch cultivation).
If relaxing the protein pool exchange does not increase the growth
rate, then this is not due to enzyme constraints, but rather an issue
with the metabolic network itself, or the set nutrient exchange is not
sufficient.
Input:
model an ecModel in GECKO 3 format (with ecModel.ec structure)
expGrowth estimated experimental growth rate. If not specified,
the value will be read from the model adapter
foldChange a value how much increase the enzyme concentration
(Optional, default = 2)
iterPerEnzyme the number of iterations that an enzyme can be increased.
A zero number can be defined. if zero is defined no limit
will be set, and it will increase the enzyme concentration
until reach de defined growth rate (Optional, default = 5)
modelAdapter a loaded model adapter (Optional, will otherwise use the
default model adapter)
verbose logical whether progress should be reported (Optional,
default true)
Output:
model ecModel where the constraint of measured enzyme
concentrations have been flexibilized (=relaxed) in the
S-matrix, to allow the model to reach the growth rate
defined in expGrowth, while the values in
ecModel.ec.concs have remained untouched.
flexEnz array with information about flexibilized proteins
uniprotIDs enzymes whose usage UB was flexibilized
oldConcs original concentrations, from mode.ec.concs
flexConcs flexibilized concentrations, new UB in
model
ratioIncr ratio by which the concentration increased,
the enzymes will be sorted by this field
frequence numeric how often the enzyme has been
step-wise flexibilized
Usage:
[model, flexEnz] = flexibilizeEnzConcs(model, expGrowth, foldChange, iterPerEnzyme, modelAdapter, verbose)0001 function [model, flexEnz] = flexibilizeEnzConcs(model, expGrowth, foldChange, iterPerEnzyme, modelAdapter, verbose) 0002 % flexibilizeEnzConcs 0003 % Flexibilize enzyme concentration of an ecModel with constrained with 0004 % proteomics data. The upper bound of the protein usage reaction is 0005 % changed, while the concentrations in ecModel.ec.concs remain unchanged. 0006 % 0007 % If no (more) limiting enzyme concentrations can be found, it might be 0008 % the protein pool exchange that is limiting growth. In that case, an 0009 % attempt will be made to relax the protein pool exchange reaction, and 0010 % if the growth rate indeed increases, it is suggested to set the protein 0011 % pool exchange unconstrained (lb=-1000) before again running 0012 % flexibilizeEnzConcs. Such situations, where a proteomics integrated 0013 % ecModel is overconstrained, may occur if the ecModel should be able to 0014 % simulate maximum growth rate (from e.g. batch cultivation). 0015 % 0016 % If relaxing the protein pool exchange does not increase the growth 0017 % rate, then this is not due to enzyme constraints, but rather an issue 0018 % with the metabolic network itself, or the set nutrient exchange is not 0019 % sufficient. 0020 % 0021 % Input: 0022 % model an ecModel in GECKO 3 format (with ecModel.ec structure) 0023 % expGrowth estimated experimental growth rate. If not specified, 0024 % the value will be read from the model adapter 0025 % foldChange a value how much increase the enzyme concentration 0026 % (Optional, default = 2) 0027 % iterPerEnzyme the number of iterations that an enzyme can be increased. 0028 % A zero number can be defined. if zero is defined no limit 0029 % will be set, and it will increase the enzyme concentration 0030 % until reach de defined growth rate (Optional, default = 5) 0031 % modelAdapter a loaded model adapter (Optional, will otherwise use the 0032 % default model adapter) 0033 % verbose logical whether progress should be reported (Optional, 0034 % default true) 0035 % 0036 % Output: 0037 % model ecModel where the constraint of measured enzyme 0038 % concentrations have been flexibilized (=relaxed) in the 0039 % S-matrix, to allow the model to reach the growth rate 0040 % defined in expGrowth, while the values in 0041 % ecModel.ec.concs have remained untouched. 0042 % flexEnz array with information about flexibilized proteins 0043 % uniprotIDs enzymes whose usage UB was flexibilized 0044 % oldConcs original concentrations, from mode.ec.concs 0045 % flexConcs flexibilized concentrations, new UB in 0046 % model 0047 % ratioIncr ratio by which the concentration increased, 0048 % the enzymes will be sorted by this field 0049 % frequence numeric how often the enzyme has been 0050 % step-wise flexibilized 0051 % 0052 % Usage: 0053 % [model, flexEnz] = flexibilizeEnzConcs(model, expGrowth, foldChange, iterPerEnzyme, modelAdapter, verbose) 0054 0055 if nargin < 6 || isempty(verbose) 0056 verbose = true; 0057 end 0058 if nargin < 5 || isempty(modelAdapter) 0059 modelAdapter = ModelAdapterManager.getDefault(); 0060 if isempty(modelAdapter) 0061 error('Either send in a modelAdapter or set the default model adapter in the ModelAdapterManager.') 0062 end 0063 end 0064 params = modelAdapter.getParameters(); 0065 0066 if nargin < 4 || isempty(iterPerEnzyme) 0067 iterPerEnzyme = 5; 0068 end 0069 0070 if nargin < 3 || isempty(foldChange) 0071 foldChange = 2; 0072 end 0073 0074 if nargin < 2 || isempty(expGrowth) 0075 expGrowth = modelAdapter.getParameters().gR_exp; 0076 end 0077 0078 % If a zero value is defined, not iteration limit will be set 0079 if iterPerEnzyme == 0 0080 iterPerEnzyme = inf; 0081 end 0082 0083 bioRxnIdx = getIndexes(model,params.bioRxn,'rxns'); 0084 if model.ub(bioRxnIdx) < expGrowth 0085 fprintf(['The upper bound of the biomass reaction was %s, while expGrowth '... 0086 'is %s. The upper bound has been set to expGrowth.'],... 0087 num2str(model.ub(bioRxnIdx)), num2str(expGrowth)) 0088 model.ub(bioRxnIdx) = expGrowth; 0089 end 0090 0091 % Keep track if protein pool will be flexibilized 0092 changedProtPool = false; 0093 0094 % In case the model have not been protein constrained 0095 % model = constrainEnzConcs(model); 0096 0097 [sol,hs] = solveLP(model); 0098 predGrowth = sol.f; 0099 0100 % Get those proteins with a concentration defined 0101 protConcs = find(~isnan(model.ec.concs)); 0102 0103 % Get enzymes names with a concentration value 0104 proteins = model.ec.enzymes(protConcs); 0105 0106 % Store how many time is predicted a enzyme with the highest coeff 0107 frequence = zeros(numel(proteins),1); 0108 0109 flexBreak=false; 0110 if any(protConcs) 0111 while predGrowth < expGrowth 0112 0113 % Get the control coefficients 0114 [~, controlCoeffs] = getConcControlCoeffs(model, proteins, foldChange, 0.75); 0115 0116 % Stop looking for limiting proteins all coeffs are zero 0117 if any(controlCoeffs) 0118 0119 % Find the maximum coefficient index 0120 [~,maxIdx] = max(controlCoeffs); 0121 0122 % Increase +1 the frequence 0123 frequence(maxIdx) = frequence(maxIdx) + 1; 0124 0125 % Allow to increase the UB maximum five times 0126 if frequence(maxIdx) <= iterPerEnzyme 0127 0128 % Set how much will increase the UB 0129 increase = foldChange*frequence(maxIdx); 0130 0131 % Get usage rxn for the highest coeff 0132 protUsageIdx = strcmpi(model.rxns, strcat('usage_prot_', proteins(maxIdx))); 0133 0134 % replace the UB 0135 model.lb(protUsageIdx) = -(model.ec.concs(protConcs(maxIdx)) * (1+increase)); 0136 0137 % Get the new growth rate 0138 sol = solveLP(model,0,[],hs); 0139 predGrowth = sol.f; 0140 0141 if verbose 0142 disp(['Protein ' proteins{maxIdx} ' LB adjusted. Grow: ' num2str(predGrowth)]) 0143 end 0144 else 0145 printOrange(['WARNING: Limit has been reached. Protein ' proteins{maxIdx} ' seems to be problematic. Consider changing the kcat.\n']); 0146 flexBreak=true; 0147 break 0148 end 0149 else 0150 disp('No (more) limiting enzymes have been found. Attempting to increase protein pool exchange...') 0151 protPoolIdx = getIndexes(model,'prot_pool_exchange','rxns'); 0152 oldProtPool = -model.lb(protPoolIdx); 0153 tempModel = setParam(model,'lb',protPoolIdx,-1000); 0154 tempModel = setParam(tempModel,'ub',bioRxnIdx,expGrowth); 0155 sol = solveLP(tempModel); 0156 if (sol.f-predGrowth)>1e-10 % There is improvement in growth rate 0157 % Find new protein pool constraint 0158 predGrowth = sol.f; 0159 tempModel = setParam(tempModel,'lb',bioRxnIdx,predGrowth); 0160 tempModel = setParam(tempModel,'obj',protPoolIdx,1); 0161 sol = solveLP(tempModel); 0162 newProtPool = sol.f; 0163 model.lb(protPoolIdx) = -newProtPool; 0164 fprintf(['Changing the lower bound of protein pool exchange from %s ',... 0165 'to %s enabled a\ngrowth rate of %s. It can be helpful to set ', ... 0166 'the lower bound of protein\npool exchange to -1000 before ',... 0167 'running flexibilizeEnzConcs.\n'],num2str(-oldProtPool), ... 0168 num2str(-newProtPool),num2str(predGrowth)); 0169 changedProtPool = true; 0170 else 0171 fprintf(['Protein pool exchange was also not limiting. Inability to reach growth ',... 0172 'rate is not related to\nenzyme constraints. Maximum growth rate ',... 0173 'is %s.\n'], num2str(predGrowth)) 0174 end 0175 % To return an usable model with these conditions, set the 0176 % highest growth rate reached as the experimental value 0177 expGrowth = predGrowth; 0178 break 0179 end 0180 end 0181 else 0182 error('Protein concentrations have not been defined. Please run readProteomics and constrainEnzConcs') 0183 end 0184 0185 protFlex = proteins(frequence>0); 0186 protUsageIdx = find(ismember(model.rxns, strcat('usage_prot_', protFlex))); 0187 ecProtId = find(ismember(model.ec.enzymes,protFlex)); 0188 oldConcs = model.ec.concs(ecProtId); 0189 0190 if flexBreak == false && ~isempty(protFlex) 0191 % Not all flexibilized proteins require flexibilization in the end 0192 % Test flux distribution with minimum prot_pool usage 0193 modelTemp = setParam(model,'var',params.bioRxn,expGrowth,0.5); 0194 modelTemp = setParam(modelTemp,'obj','prot_pool_exchange',1); 0195 sol = solveLP(modelTemp,1); 0196 % Check which concentrations can remain unchanged 0197 newConcs = abs(sol.x(protUsageIdx)); 0198 keepOld = newConcs<oldConcs; 0199 % Set UBs 0200 model.lb(protUsageIdx(keepOld)) = -oldConcs(keepOld); 0201 model.lb(protUsageIdx(~keepOld)) = -newConcs(~keepOld); 0202 % Output structure 0203 protFlex(keepOld) = []; 0204 oldConcs(keepOld) = []; 0205 newConcs(keepOld) = []; 0206 frequence=frequence(frequence>0); 0207 frequence(keepOld) = []; 0208 0209 flexEnz.uniprotIDs = protFlex; 0210 flexEnz.oldConcs = oldConcs; 0211 flexEnz.flexConcs = newConcs; 0212 flexEnz.frequence = frequence(frequence>0); 0213 else 0214 protFlex = proteins(frequence>0); 0215 protUsageIdx = find(ismember(model.rxns, strcat('usage_prot_', protFlex))); 0216 0217 flexEnz.uniprotIDs = protFlex; 0218 flexEnz.oldConcs = model.ec.concs(ecProtId); 0219 flexEnz.flexConcs = abs(model.lb(protUsageIdx)); 0220 flexEnz.frequence = frequence(frequence>0); 0221 end 0222 if changedProtPool 0223 flexEnz.uniprotIDs(end+1) = {'prot_pool'}; 0224 flexEnz.oldConcs(end+1) = oldProtPool; 0225 flexEnz.flexConcs(end+1) = newProtPool; 0226 flexEnz.frequence(end+1) = 1; 0227 end 0228 % Sort by biggest ratio increase 0229 flexEnz.ratioIncr = flexEnz.flexConcs./flexEnz.oldConcs; 0230 [~,b] = sort(flexEnz.ratioIncr,'descend'); 0231 flexEnz.ratioIncr = flexEnz.ratioIncr(b); 0232 flexEnz.uniprotIDs = flexEnz.uniprotIDs(b); 0233 flexEnz.oldConcs = flexEnz.oldConcs(b); 0234 flexEnz.flexConcs = flexEnz.flexConcs(b); 0235 flexEnz.frequence = flexEnz.frequence(b); 0236 end