setExchangeBounds
Define the exchange flux bounds for a given set of metabolites.
Input:
model a model structure
mets a cell array of metabolite names (case insensitive) or
metabolite IDs, or a vector of metabolite indices
(opt, default all exchanged metabolites)
lb lower bound of exchange flux. Can be either a vector of
bounds corresponding to each of the provided metabolites,
or a single value that will be applied to all.
(opt, default to model.annotation.defaultLB if it exists,
otherwise -1000)
ub upper bound of exchange flux. Can be either a vector of
bounds corresponding to each of the provided metabolites,
or a single value that will be applied to all.
(opt, default to model.annotation.defaultUB if it exists,
otherwise 1000)
closeOthers close exchange reactions for all other exchanged
metabolites not present in the provided list. This will
prevent IMPORT of the metabolites, but their EXPORT will
not be modified.
(opt, default true)
mediaOnly only consider exchange reactions involving exchange to or
from the extracellular (media) compartment. Reactions
such as "sink" reactions that exchange metabolites
directly with an intracellular compartment will therefore
be ignored even though "getExchangeRxns" identifies such
such reactions as exchange reactions.
Note: The function will attempt to identify the
extracellular compartment by the "compNames" field, and
also requires the "metComps" field to be present,
otherwise the mediaOnly flag will be ignored.
(opt, default false)
Output:
exchModel a model structure with updated exchange flux bounds for
the provided set of metabolites
unusedMets metabolites provided by the user that were not used
because they are not involved in any exchange reactions
in the model
NOTE: Exchange reactions involving more than one metabolite will be
ignored.
Usage: exchModel = setExchangeBounds(model,mets,lb,ub,closeOthers,mediaOnly);0001 function [exchModel,unusedMets] = setExchangeBounds(model,mets,lb,ub,closeOthers,mediaOnly) 0002 % setExchangeBounds 0003 % Define the exchange flux bounds for a given set of metabolites. 0004 % 0005 % Input: 0006 % model a model structure 0007 % mets a cell array of metabolite names (case insensitive) or 0008 % metabolite IDs, or a vector of metabolite indices 0009 % (opt, default all exchanged metabolites) 0010 % lb lower bound of exchange flux. Can be either a vector of 0011 % bounds corresponding to each of the provided metabolites, 0012 % or a single value that will be applied to all. 0013 % (opt, default to model.annotation.defaultLB if it exists, 0014 % otherwise -1000) 0015 % ub upper bound of exchange flux. Can be either a vector of 0016 % bounds corresponding to each of the provided metabolites, 0017 % or a single value that will be applied to all. 0018 % (opt, default to model.annotation.defaultUB if it exists, 0019 % otherwise 1000) 0020 % closeOthers close exchange reactions for all other exchanged 0021 % metabolites not present in the provided list. This will 0022 % prevent IMPORT of the metabolites, but their EXPORT will 0023 % not be modified. 0024 % (opt, default true) 0025 % mediaOnly only consider exchange reactions involving exchange to or 0026 % from the extracellular (media) compartment. Reactions 0027 % such as "sink" reactions that exchange metabolites 0028 % directly with an intracellular compartment will therefore 0029 % be ignored even though "getExchangeRxns" identifies such 0030 % such reactions as exchange reactions. 0031 % Note: The function will attempt to identify the 0032 % extracellular compartment by the "compNames" field, and 0033 % also requires the "metComps" field to be present, 0034 % otherwise the mediaOnly flag will be ignored. 0035 % (opt, default false) 0036 % 0037 % Output: 0038 % exchModel a model structure with updated exchange flux bounds for 0039 % the provided set of metabolites 0040 % unusedMets metabolites provided by the user that were not used 0041 % because they are not involved in any exchange reactions 0042 % in the model 0043 % 0044 % NOTE: Exchange reactions involving more than one metabolite will be 0045 % ignored. 0046 % 0047 % Usage: exchModel = setExchangeBounds(model,mets,lb,ub,closeOthers,mediaOnly); 0048 0049 0050 % handle input arguments 0051 if nargin < 2 0052 mets = []; 0053 elseif ~islogical(mets) || ~isnumeric(mets) 0054 mets=convertCharArray(mets); 0055 end 0056 0057 if nargin < 3 || isempty(lb) 0058 if isfield(model,'annotation') && isfield(model.annotation,'defaultLB') 0059 lb = model.annotation.defaultLB; 0060 else 0061 lb = -1000; 0062 end 0063 end 0064 0065 if nargin < 4 || isempty(ub) 0066 if isfield(model,'annotation') && isfield(model.annotation,'defaultUB') 0067 ub = model.annotation.defaultUB; 0068 else 0069 ub = 1000; 0070 end 0071 end 0072 0073 if nargin < 5 || isempty(closeOthers) 0074 closeOthers = true; 0075 end 0076 0077 if nargin < 6 0078 mediaOnly = false; 0079 elseif mediaOnly 0080 if ~all(isfield(model,{'compNames','metComps'})) 0081 error('mediaOnly option requires the "compNames" and "metComps" model fields.'); 0082 end 0083 end 0084 0085 0086 % for models with an "unconstrained" field, generate a version of the model 0087 % where all the unconstrained metabolites have stoichiometric coeffs set to 0088 % zero, but remain in the model to avoid index changes 0089 model_temp = model; 0090 if isfield(model,'unconstrained') 0091 model_temp.S(model_temp.unconstrained == 1,:) = 0; 0092 end 0093 0094 % find exchange rxns, ignoring those involving more than one metabolite 0095 [~,exchRxnInd] = getExchangeRxns(model); 0096 exchRxnInd(sum(model_temp.S(:,exchRxnInd) ~= 0,1) > 1) = []; 0097 0098 % find all exchanged metabolites 0099 exchMetInd = find(any(model_temp.S(:,exchRxnInd) ~= 0,2)); 0100 0101 if ( mediaOnly ) 0102 % ignore exchanged metabolites in non-extracellular compartments and 0103 % any exchange reactions involving these metabolites 0104 [~,extCompInd] = ismember('extracellular',lower(model.compNames)); 0105 if extCompInd > 0 0106 ignoreMet = (model.metComps(exchMetInd) ~= extCompInd); 0107 else 0108 error('Could not find any compartments named "extracellular".'); 0109 end 0110 ignoreRxn = any(model.S(exchMetInd(ignoreMet),exchRxnInd) ~= 0,1); 0111 exchMetInd(ignoreMet) = []; 0112 exchRxnInd(ignoreRxn) = []; 0113 end 0114 0115 % Check that all exchange reactions are formulated in the same direction. 0116 % If not, this means that negative flux indicates import for some exchange 0117 % reactions, but indicates export for others. Therefore, the LB and UB 0118 % would need to be specified differently depending on the exchange reaction 0119 % direction, which is error-prone. 0120 if all(all(model.S(exchMetInd,exchRxnInd) <= 0)) 0121 importDir = 'backward'; 0122 elseif all(all(model.S(exchMetInd,exchRxnInd) >= 0)) 0123 importDir = 'forward'; 0124 else 0125 fprintf('WARNING: Some exchange reactions differ in direction, and therefore have opposite meanings of LB and UB.\n'); 0126 if closeOthers 0127 fprintf(' Therefore, the "closeOthers" option will be set to FALSE.\n'); 0128 end 0129 closeOthers = false; 0130 end 0131 0132 % prepare exchanged metabolites and bounds 0133 if ~isempty(mets) 0134 0135 % prepare lb and ub variables 0136 if numel(lb) == 1 0137 lb = lb*ones(size(mets)); 0138 elseif numel(lb) ~= numel(mets) 0139 error('lb must be a single value or a vector of equal length as mets'); 0140 end 0141 if numel(ub) == 1 0142 ub = ub*ones(size(mets)); 0143 elseif numel(ub) ~= numel(mets) 0144 error('ub must be a single value or a vector of equal length as mets'); 0145 end 0146 0147 % map provided mets to exchanged metabolites 0148 if isnumeric(mets) 0149 % mets are provided as indices 0150 [isExch,metInd] = ismember(mets,exchMetInd); 0151 elseif sum(ismember(lower(mets),lower(model.metNames))) > sum(ismember(mets,model.mets)) 0152 % assume that mets are provided as names 0153 [isExch,metInd] = ismember(lower(mets),lower(model.metNames(exchMetInd))); 0154 else 0155 % assume that mets are provided as met IDs 0156 [isExch,metInd] = ismember(mets,model.mets(exchMetInd)); 0157 end 0158 0159 % get exchanged met indices and determine unused (non-exchanged) mets 0160 unusedMets = mets(~isExch); 0161 exchMetInd = exchMetInd(metInd(isExch)); 0162 0163 % update bound vectors 0164 lb = lb(isExch); 0165 ub = ub(isExch); 0166 0167 else 0168 % if no mets were provided, use all exchanged mets 0169 if numel(lb) > 1 || numel(ub) > 1 0170 error('Only one upper and one lower bound may be provide if metabolites are not specified.'); 0171 else 0172 lb = lb*ones(size(exchMetInd)); 0173 ub = ub*ones(size(exchMetInd)); 0174 end 0175 unusedMets = {}; 0176 end 0177 0178 % check that at least one exchanged metabolite matches 0179 if isempty(exchMetInd) 0180 fprintf('None of the provided metabolites were found in any exchange reactions.\n'); 0181 exchModel = model; 0182 return 0183 end 0184 0185 % determine which metabolite is exchanged in each exchange reaction 0186 [metInd,rxnInd] = find(model_temp.S(exchMetInd,exchRxnInd) ~= 0); 0187 0188 % check for any metabolites that are exchanged in more than one reaction 0189 [tbl,i,~] = unique(metInd,'first'); 0190 repeatedInds = find(not(ismember(1:numel(tbl),i))); 0191 multiMetInd = exchMetInd(metInd(repeatedInds)); 0192 if ~isempty(multiMetInd) 0193 fprintf('WARNING: The following metabolites are involved in more than one exchange reaction:\n'); 0194 fprintf('\t%s\n',model.metNames{multiMetInd(1:min(numel(multiMetInd),10))}); 0195 if numel(multiMetInd) > 10 0196 fprintf('\t... and %u more.\n',numel(multiMetInd)-10); 0197 end 0198 end 0199 0200 % set exchange reaction bounds 0201 model = setParam(model,'lb',exchRxnInd(rxnInd),lb(metInd)); 0202 model = setParam(model,'ub',exchRxnInd(rxnInd),ub(metInd)); 0203 0204 if closeOthers 0205 % constrain import of all other exchange reactions to zero 0206 constrainInd = setdiff(exchRxnInd,exchRxnInd(rxnInd)); 0207 if strcmp(importDir,'backward') 0208 model = setParam(model,'lb',constrainInd,0); 0209 else 0210 model = setParam(model,'ub',constrainInd,0); 0211 end 0212 end 0213 0214 % assign output model 0215 exchModel = model; 0216 0217