setConstraintsIrrevModel

PURPOSE ^

setConstraintsIrrevModel Set constraints for a subset of rxns while

SYNOPSIS ^

function constrOptIrrev = setConstraintsIrrevModel(constrOpt,model,modelIrrev,rev2irrev)

DESCRIPTION ^

setConstraintsIrrevModel Set constraints for a subset of rxns while
converting reversible to irreversible reaction names and handling the
constraint directions correctly

 [constrIndIrrev,constrValIrrev,constrSenseIrrev] = ...
           setConstraintsIrrevModel(rxnNameList,constrValue,constrSense,model,modelIrrev)
 
INPUTS
 constrOpt    Constraint options
   rxnList     Reaction selection cell array (for reversible
               representation)
   values      Constraint values
   sense       Constraint senses ordered as rxnNameList

 model         Model in reversible format
 modelIrrev    Model in irreversible format
 rev2irrev     Reversible to irreversible reaction index conversion
               obtained from convertToIrreversible

OUTPUT
 constrOptIrrev   Constraint options in irrev model
   rxnList         Reaction selection cell array
     rxnInd          Selection index for constraints in irreversible model (e.g. [2 4 5 9 10])
   values          Correctly ordered constraint values
     sense           Correctly ordered constraint senses

 6/9/05 Changed this so that it allows multiple occurences of the same
 rxn

 Markus Herrgard 10/14/03

 Completely rewritten 1/22/07

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function constrOptIrrev = setConstraintsIrrevModel(constrOpt,model,modelIrrev,rev2irrev)
0002 %setConstraintsIrrevModel Set constraints for a subset of rxns while
0003 %converting reversible to irreversible reaction names and handling the
0004 %constraint directions correctly
0005 %
0006 % [constrIndIrrev,constrValIrrev,constrSenseIrrev] = ...
0007 %           setConstraintsIrrevModel(rxnNameList,constrValue,constrSense,model,modelIrrev)
0008 %
0009 %INPUTS
0010 % constrOpt    Constraint options
0011 %   rxnList     Reaction selection cell array (for reversible
0012 %               representation)
0013 %   values      Constraint values
0014 %   sense       Constraint senses ordered as rxnNameList
0015 %
0016 % model         Model in reversible format
0017 % modelIrrev    Model in irreversible format
0018 % rev2irrev     Reversible to irreversible reaction index conversion
0019 %               obtained from convertToIrreversible
0020 %
0021 %OUTPUT
0022 % constrOptIrrev   Constraint options in irrev model
0023 %   rxnList         Reaction selection cell array
0024 %     rxnInd          Selection index for constraints in irreversible model (e.g. [2 4 5 9 10])
0025 %   values          Correctly ordered constraint values
0026 %     sense           Correctly ordered constraint senses
0027 %
0028 % 6/9/05 Changed this so that it allows multiple occurences of the same
0029 % rxn
0030 %
0031 % Markus Herrgard 10/14/03
0032 %
0033 % Completely rewritten 1/22/07
0034 
0035 constrIndIrrev = [];
0036 constrValIrrev = [];
0037 constrSenseIrrev = '';
0038 
0039 % Bit vector describing the processing status of this constraint
0040 constrProcessed = false*ones(length(constrOpt.rxnList),1);
0041 
0042 for i = 1:length(constrOpt.rxnList)
0043     
0044     % Has this rxn already been processed?
0045     if (~constrProcessed(i))
0046         
0047         % Get reaction name
0048         rxnName = constrOpt.rxnList{i};
0049         % Find reaction index in reversible reaction set
0050         rxnID = find(ismember(model.rxns,rxnName));
0051         
0052         % Find out if there are more than one constraint for this rxn
0053         thisConstrInd = find(strcmp(constrOpt.rxnList,rxnName));
0054         % How many matching rxns
0055         nThisConstr = length(thisConstrInd);
0056         % Mark as processed
0057         constrProcessed(thisConstrInd) = true;
0058         
0059         % Get constraint values
0060         thisConstrValue = constrOpt.values(thisConstrInd);
0061         % Get constraint directions
0062         thisConstrSense = constrOpt.sense(thisConstrInd);
0063         
0064         if (~isempty(rxnID))
0065             % Find reaction index in irreversible reaction set
0066             irrevRxnID = rev2irrev{rxnID};
0067             
0068             % Add reaction indices and constraint values into the list
0069             if (length(irrevRxnID) == 1) % Irreversible
0070                 irrevRxnName = modelIrrev.rxns{irrevRxnID};
0071                 if length(rxnName)>2 && ~strcmp(irrevRxnName(end-1:end),'_r')
0072                     % Value of map is directly the index of the reaction
0073                     for j = 1:nThisConstr
0074                         constrIndIrrev(end+1) = irrevRxnID;
0075                         constrValIrrev(end+1) = thisConstrValue(j);
0076                         constrSenseIrrev(end+1) = thisConstrSense(j);
0077                     end
0078                 else %Reaction is reversed in model. Flip sign of value
0079                     for j = 1:nThisConstr
0080                         constrIndIrrev(end+1) = irrevRxnID;
0081                         constrValIrrev(end+1) = -thisConstrValue(j);
0082                         constrSenseIrrev(end+1) = thisConstrSense(j);
0083                     end
0084                 end
0085             else % Reversible
0086                 % Only one constraint or an equality constraint represented
0087                 % through two inequality constraints
0088                 if (nThisConstr == 1 | (thisConstrValue(1) == thisConstrValue(2)))
0089                     % This would be an equality constraint
0090                     if (nThisConstr == 2)
0091                         thisConstrSense = 'E';
0092                     end
0093                     % Map contains both forward and reverse reaction indices
0094                     constrIndIrrev(end+1) = irrevRxnID(1);
0095                     constrIndIrrev(end+1) = irrevRxnID(2);
0096                     if (thisConstrValue(1) >= 0)
0097                         constrValIrrev(end+1) = thisConstrValue(1);
0098                         constrValIrrev(end+1) = 0;
0099                         % Deal with different directions of constraints
0100                         switch thisConstrSense
0101                             case 'E'
0102                                 constrSenseIrrev(end+1) = 'E';
0103                                 constrSenseIrrev(end+1) = 'E';
0104                             case 'G'
0105                                 constrSenseIrrev(end+1) = 'G';
0106                                 constrSenseIrrev(end+1) = 'E';
0107                             case 'L'
0108                                 constrSenseIrrev(end+1) = 'L';
0109                                 constrSenseIrrev(end+1) = 'G';
0110                         end
0111                     else
0112                         constrValIrrev(end+1) = 0;
0113                         constrValIrrev(end+1) = -thisConstrValue(1);
0114                         switch thisConstrSense
0115                             case 'E'
0116                                 constrSenseIrrev(end+1) = 'E';
0117                                 constrSenseIrrev(end+1) = 'E';
0118                             case 'G'
0119                                 constrSenseIrrev(end+1) = 'G';
0120                                 constrSenseIrrev(end+1) = 'L';
0121                             case 'L'
0122                                 constrSenseIrrev(end+1) = 'E';
0123                                 constrSenseIrrev(end+1) = 'G';
0124                         end
0125                     end
0126                 else % More than one constraint (the only case that makes sense is a <= v <= b)
0127                     lowestConstrVal = min(thisConstrValue);        
0128                     highestConstrVal = max(thisConstrValue);
0129                     
0130                     if ((lowestConstrVal > 0) & (highestConstrVal > 0))  % Both positive
0131                         constrIndIrrev(end+1:end+3) = irrevRxnID([1 1 2]);
0132                         constrValIrrev(end+1:end+3) = [lowestConstrVal highestConstrVal 0];
0133                         constrSenseIrrev(end+1:end+3) = 'GLE';
0134                     elseif ((lowestConstrVal < 0) & (highestConstrVal < 0)) % Both negative
0135                         constrIndIrrev(end+1:end+3) = irrevRxnID([1 2 2]);
0136                         constrValIrrev(end+1:end+3) = [0 -highestConstrVal -lowestConstrVal];
0137                         constrSenseIrrev(end+1:end+3) = 'EGL';
0138                     else % Low positive, hi negative
0139                         constrIndIrrev(end+1:end+2) = irrevRxnID([1 2]);
0140                         constrValIrrev(end+1:end+2) = [highestConstrVal -lowestConstrVal];
0141                         constrSenseIrrev(end+1:end+2) = 'LL';
0142                     end
0143                 end
0144             end
0145         end
0146     end
0147 end
0148 
0149 constrOptIrrev.rxnList = columnVector(modelIrrev.rxns(constrIndIrrev));
0150 constrOptIrrev.rxnInd = columnVector(constrIndIrrev);
0151 constrOptIrrev.values = columnVector(constrValIrrev);
0152 constrOptIrrev.sense = columnVector(constrSenseIrrev);

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003