addRxns

 addRxns
   Adds reactions to a model

 Input:
   model            a model structure
   rxnsToAdd        the reaction structure can have the following fields:
       rxns                cell array with unique strings that identifies
                           each reaction
       equations           cell array with equation strings. Decimal
                           coefficients are expressed as "1.2".
                           Reversibility is indicated by "<=>" or "=>"
       mets                (alternative to equations) cell array with the
                           metabolites involved in each reaction as nested
                           arrays. E.g.: {{'met1','met2'},{'met1','met3','met4'}}
                           In the case of one single reaction added, it
                           can be a string array: {'met1','met2'}
       stoichCoeffs        (alternative to equations) cell array with the
                           corresponding stoichiometries as nested vectors
                           E.g.: {[-1,+2],[-1,-1,+1]}. In the case of one
                           single reaction added, it can be a vector: [-1,+2]
       rxnNames            cell array with the names of each reaction
                           (optional, default '')
       lb                  vector with the lower bounds (optional, default
                           model.annotations.defaultLB or -inf for
                           reversible reactions and 0 for irreversible
                           when "equations" is used. When "mets" and
                           "stoichCoeffs" are ,used it defaults for all
                           to model.annotations.defaultLB or -inf)
       ub                  vector with the upper bounds (optional, default
                           model.annotations.defaultUB or inf)
       c                   vector with the objective function coefficients
                           (optional, default 0)
       eccodes             cell array with the EC-numbers for each
                           reactions. Delimit several EC-numbers with ";"
                           (optional, default '')
       subSystems          cell array with the subsystems for each
                           reaction (optional, default '')
       grRules             cell array with the gene-reaction relationship
                           for each reaction. E.g. "(A and B) or (C)"
                           means that the reaction could be catalyzed by a
                           complex between A & B or by C on its own. All
                           the genes have to be present in model.genes.
                           Add genes with addGenesRaven before calling
                           this function if needed (optional, default '')
       rxnMiriams          cell array with Miriam structures (optional,
                           default [])
       rxnComps            cell array with compartments (as in
                           model.comps) (optional, default {})
       rxnNotes            cell array with reaction notes (optional,
                           default '')
       rxnDeltaG           Gibbs free energy at biochemical standard
                           condition in kJ/mole (optional, default NaN)
       rxnReferences       cell array with reaction references (optional,
                           default '')
       rxnConfidenceScores vector with reaction confidence scores
                           (optional, default NaN)
   eqnType          double describing how the equation string should be
                    interpreted
                    1 - The metabolites are matched to model.mets. New
                        metabolites (if allowed) are added to
                        "compartment" (default)
                    2 - The metabolites are matched to model.metNames and
                        all metabolites are assigned to "compartment". Any
                        new metabolites that are added will be assigned
                        IDs "m1", "m2"... If IDs on the same form are
                        already used in the model then the numbering will
                        start from the highest used integer+1
                    3 - The metabolites are written as
                        "metNames[comps]". Only compartments in
                        model.comps are allowed. Any
                        new metabolites that are added will be assigned
                        IDs "m1", "m2"... If IDs on the same form are
                        already used in the model then the numbering will
                        start from the highest used integer+1
   compartment      a string with the compartment the metabolites should
                    be placed in when using eqnType=2. Must match
                    model.comps (optional when eqnType=1 or eqnType=3)
   allowNewMets     true if the function is allowed to add new
                    metabolites. Can also be a string, which will be used
                    as prefix for the new metabolite IDs. It is highly
                    recommended to first add any new metabolites with
                    addMets rather than automatically through this
                    function. addMets supports more annotation of
                    metabolites, allows for the use of exchange
                    metabolites, and using it reduces the risk of parsing
                    errors (optional, default false)
   allowNewGenes    true if the functions is allowed to add new genes
                    (optional, default false)

 Output:
   newModel         an updated model structure

   This function does not make extensive checks about formatting of
   gene-reaction rules.

   When adding metabolites to a compartment where they previously do not
   the function will copy any available information from the metabolite in
   another compartment.

 Usage: newModel = addRxns(model, rxnsToAdd, eqnType, compartment,...
                       allowNewMets, allowNewGenes)

0001 function newModel=addRxns(model,rxnsToAdd,eqnType,compartment,allowNewMets,allowNewGenes)
0002 % addRxns
0003 %   Adds reactions to a model
0004 %
0005 % Input:
0006 %   model            a model structure
0007 %   rxnsToAdd        the reaction structure can have the following fields:
0008 %       rxns                cell array with unique strings that identifies
0009 %                           each reaction
0010 %       equations           cell array with equation strings. Decimal
0011 %                           coefficients are expressed as "1.2".
0012 %                           Reversibility is indicated by "<=>" or "=>"
0013 %       mets                (alternative to equations) cell array with the
0014 %                           metabolites involved in each reaction as nested
0015 %                           arrays. E.g.: {{'met1','met2'},{'met1','met3','met4'}}
0016 %                           In the case of one single reaction added, it
0017 %                           can be a string array: {'met1','met2'}
0018 %       stoichCoeffs        (alternative to equations) cell array with the
0019 %                           corresponding stoichiometries as nested vectors
0020 %                           E.g.: {[-1,+2],[-1,-1,+1]}. In the case of one
0021 %                           single reaction added, it can be a vector: [-1,+2]
0022 %       rxnNames            cell array with the names of each reaction
0023 %                           (optional, default '')
0024 %       lb                  vector with the lower bounds (optional, default
0025 %                           model.annotations.defaultLB or -inf for
0026 %                           reversible reactions and 0 for irreversible
0027 %                           when "equations" is used. When "mets" and
0028 %                           "stoichCoeffs" are ,used it defaults for all
0029 %                           to model.annotations.defaultLB or -inf)
0030 %       ub                  vector with the upper bounds (optional, default
0031 %                           model.annotations.defaultUB or inf)
0032 %       c                   vector with the objective function coefficients
0033 %                           (optional, default 0)
0034 %       eccodes             cell array with the EC-numbers for each
0035 %                           reactions. Delimit several EC-numbers with ";"
0036 %                           (optional, default '')
0037 %       subSystems          cell array with the subsystems for each
0038 %                           reaction (optional, default '')
0039 %       grRules             cell array with the gene-reaction relationship
0040 %                           for each reaction. E.g. "(A and B) or (C)"
0041 %                           means that the reaction could be catalyzed by a
0042 %                           complex between A & B or by C on its own. All
0043 %                           the genes have to be present in model.genes.
0044 %                           Add genes with addGenesRaven before calling
0045 %                           this function if needed (optional, default '')
0046 %       rxnMiriams          cell array with Miriam structures (optional,
0047 %                           default [])
0048 %       rxnComps            cell array with compartments (as in
0049 %                           model.comps) (optional, default {})
0050 %       rxnNotes            cell array with reaction notes (optional,
0051 %                           default '')
0052 %       rxnDeltaG           Gibbs free energy at biochemical standard
0053 %                           condition in kJ/mole (optional, default NaN)
0054 %       rxnReferences       cell array with reaction references (optional,
0055 %                           default '')
0056 %       rxnConfidenceScores vector with reaction confidence scores
0057 %                           (optional, default NaN)
0058 %   eqnType          double describing how the equation string should be
0059 %                    interpreted
0060 %                    1 - The metabolites are matched to model.mets. New
0061 %                        metabolites (if allowed) are added to
0062 %                        "compartment" (default)
0063 %                    2 - The metabolites are matched to model.metNames and
0064 %                        all metabolites are assigned to "compartment". Any
0065 %                        new metabolites that are added will be assigned
0066 %                        IDs "m1", "m2"... If IDs on the same form are
0067 %                        already used in the model then the numbering will
0068 %                        start from the highest used integer+1
0069 %                    3 - The metabolites are written as
0070 %                        "metNames[comps]". Only compartments in
0071 %                        model.comps are allowed. Any
0072 %                        new metabolites that are added will be assigned
0073 %                        IDs "m1", "m2"... If IDs on the same form are
0074 %                        already used in the model then the numbering will
0075 %                        start from the highest used integer+1
0076 %   compartment      a string with the compartment the metabolites should
0077 %                    be placed in when using eqnType=2. Must match
0078 %                    model.comps (optional when eqnType=1 or eqnType=3)
0079 %   allowNewMets     true if the function is allowed to add new
0080 %                    metabolites. Can also be a string, which will be used
0081 %                    as prefix for the new metabolite IDs. It is highly
0082 %                    recommended to first add any new metabolites with
0083 %                    addMets rather than automatically through this
0084 %                    function. addMets supports more annotation of
0085 %                    metabolites, allows for the use of exchange
0086 %                    metabolites, and using it reduces the risk of parsing
0087 %                    errors (optional, default false)
0088 %   allowNewGenes    true if the functions is allowed to add new genes
0089 %                    (optional, default false)
0090 %
0091 % Output:
0092 %   newModel         an updated model structure
0093 %
0094 %   This function does not make extensive checks about formatting of
0095 %   gene-reaction rules.
0096 %
0097 %   When adding metabolites to a compartment where they previously do not
0098 %   the function will copy any available information from the metabolite in
0099 %   another compartment.
0100 %
0101 % Usage: newModel = addRxns(model, rxnsToAdd, eqnType, compartment,...
0102 %                       allowNewMets, allowNewGenes)
0103 
0104 if nargin<3
0105     eqnType=1;
0106 elseif ~isnumeric(eqnType)
0107     EM='eqnType must be numeric';
0108     dispEM(EM);
0109 elseif ~ismember(eqnType,[1 2 3])
0110     EM='eqnType must be 1, 2, or 3';
0111     dispEM(EM);
0112 end    
0113 
0114 if nargin<4
0115     compartment=[];
0116 else
0117     compartment=char(compartment);
0118 end
0119 if nargin<5
0120     allowNewMets=false;
0121 elseif ~islogical(allowNewMets)
0122     allowNewMets=char(allowNewMets);
0123 end
0124 if nargin<6
0125     allowNewGenes=false;
0126 end
0127 
0128 if allowNewGenes & isfield(rxnsToAdd,'grRules')
0129     genesToAdd.genes = strjoin(convertCharArray(rxnsToAdd.grRules));
0130     genesToAdd.genes = regexp(genesToAdd.genes,' |)|(|and|or','split'); % Remove all grRule punctuation
0131     genesToAdd.genes = genesToAdd.genes(~cellfun(@isempty,genesToAdd.genes));  % Remove spaces and empty genes
0132     genesToAdd.genes = setdiff(unique(genesToAdd.genes),model.genes); % Only keep new genes
0133     if isfield(model,'geneComps')
0134         genesToAdd.geneComps(1:numel(genesToAdd.genes)) = repmat(11,numel(genesToAdd.genes),1);
0135     end
0136     if ~isempty(genesToAdd.genes)
0137         fprintf('\nNew genes added to the model:\n')
0138         fprintf([strjoin(genesToAdd.genes,'\n') '\n'])
0139         newModel=addGenesRaven(model,genesToAdd);
0140     else
0141         newModel=model;
0142     end
0143 else
0144     newModel=model;
0145 end
0146 
0147 %If no reactions should be added
0148 if isempty(rxnsToAdd)
0149     return;
0150 end
0151 
0152 if eqnType==2 || (eqnType==1 && allowNewMets==true)
0153     compartment=char(compartment);
0154     if ~ismember(compartment,model.comps)
0155         EM='compartment must match one of the compartments in model.comps';
0156         dispEM(EM);
0157     end
0158 end
0159 
0160 if ~isfield(rxnsToAdd,'rxns')
0161     EM='rxns is a required field in rxnsToAdd';
0162     dispEM(EM);
0163 else
0164     rxnsToAdd.rxns=convertCharArray(rxnsToAdd.rxns);
0165     %To fit with some later printing
0166     rxnsToAdd.rxns=rxnsToAdd.rxns(:);
0167 end
0168 if ~(isfield(rxnsToAdd,'equations') || (isfield(rxnsToAdd,'mets') && isfield(rxnsToAdd,'stoichCoeffs')))
0169     EM='Either "equations" or "mets"+"stoichCoeffs" are required fields in rxnsToAdd';
0170     dispEM(EM);
0171 end
0172 
0173 if any(ismember(rxnsToAdd.rxns,model.rxns))
0174     EM='One or more reaction id was already present in the model. Use changeRxns or remove the existing reactions before adding new ones';
0175     dispEM(EM);
0176 end
0177 
0178 %Normal case: equations provided
0179 if isfield(rxnsToAdd,'equations')
0180     rxnsToAdd.equations=convertCharArray(rxnsToAdd.equations);
0181 
0182     %Alternative case: mets+stoichiometry provided
0183 else
0184     %In the case of 1 rxn added (array of strings + vector), transform to
0185     %cells of length=1:
0186     if iscellstr(rxnsToAdd.mets)
0187         rxnsToAdd.mets={rxnsToAdd.mets};
0188     end
0189     if isnumeric(rxnsToAdd.stoichCoeffs)
0190         rxnsToAdd.stoichCoeffs = {rxnsToAdd.stoichCoeffs};
0191     end
0192     %Now both rxnsToAdd.mets and rxnsToAdd.stoichCoeffs should be cell
0193     %arrays & of the same size:
0194     if ~iscell(rxnsToAdd.mets) || ~iscell(rxnsToAdd.stoichCoeffs)
0195         EM='rxnsToAdd.mets & rxnsToAdd.stoichCoeffs must be cell arrays';
0196         dispEM(EM);
0197     elseif length(rxnsToAdd.stoichCoeffs) ~= length(rxnsToAdd.mets)
0198         EM = 'rxnsToAdd.stoichCoeffs must have the same number of elements as rxnsToAdd.mets';
0199         dispEM(EM);
0200     end
0201     %In this case we need lb to decide if the reaction is reversible or not:
0202     if ~isfield(rxnsToAdd,'lb')
0203         %Fill with standard if it doesn't exist
0204         rxnsToAdd.lb=-inf(size(rxnsToAdd.mets));
0205     elseif ~isnumeric(rxnsToAdd.lb)
0206         EM = 'rxnsToAdd.lb must be a vector';
0207         dispEM(EM);
0208     elseif length(rxnsToAdd.lb) ~= length(rxnsToAdd.mets)
0209         EM = 'rxnsToAdd.lb must have the same number of elements as rxnsToAdd.mets';
0210         dispEM(EM);
0211     end
0212     %Now we construct equations, to comply with the rest of the script:
0213     rxnsToAdd.equations = cell(size(rxnsToAdd.mets));
0214     for i = 1:length(rxnsToAdd.mets)
0215         mets         = rxnsToAdd.mets{i};
0216         stoichCoeffs = rxnsToAdd.stoichCoeffs{i};
0217         if isfield(rxnsToAdd,'ub')
0218             isrev        = rxnsToAdd.lb(i) < 0 && rxnsToAdd.ub(i) > 0;
0219         else
0220             isrev        = rxnsToAdd.lb(i) < 0;
0221         end
0222         rxnsToAdd.equations{i} = buildEquation(mets,stoichCoeffs,isrev);
0223     end
0224 end
0225 
0226 nRxns=numel(rxnsToAdd.rxns);
0227 nOldRxns=numel(model.rxns);
0228 filler=cell(nRxns,1);
0229 filler(:)={''};
0230 cellfiller=cell(nRxns,1);
0231 cellfiller(:)={{''}};
0232 largeFiller=cell(nOldRxns,1);
0233 largeFiller(:)={''};
0234 largeCellFiller=cell(nOldRxns,1);
0235 largeCellFiller(:)={{''}};
0236 
0237 %***Add everything to the model except for the equations.
0238 if numel(rxnsToAdd.equations)~=nRxns
0239     EM='rxnsToAdd.equations must have the same number of elements as rxnsToAdd.rxns';
0240     dispEM(EM);
0241 end
0242 
0243 %Parse the equations. This is done at this early stage since I need the
0244 %reversibility info
0245 [S, mets, badRxns, reversible]=constructS(rxnsToAdd.equations);
0246 EM='The following equations have one or more metabolites both as substrate and product. Only the net equations will be added:';
0247 dispEM(EM,false,rxnsToAdd.rxns(badRxns));
0248 
0249 newModel.rev=[newModel.rev;reversible];
0250 newModel.rxns=[newModel.rxns;rxnsToAdd.rxns(:)];
0251 
0252 if isfield(rxnsToAdd,'rxnNames')
0253     rxnsToAdd.rxnNames=convertCharArray(rxnsToAdd.rxnNames);
0254     if numel(rxnsToAdd.rxnNames)~=nRxns
0255         EM='rxnsToAdd.rxnNames must have the same number of elements as rxnsToAdd.rxns';
0256         dispEM(EM);
0257     end
0258     %Fill with standard if it doesn't exist
0259     if ~isfield(newModel,'rxnNames')
0260         newModel.rxnNames=largeFiller;
0261     end
0262     newModel.rxnNames=[newModel.rxnNames;rxnsToAdd.rxnNames(:)];
0263 else
0264     %Fill with standard if it doesn't exist
0265     if isfield(newModel,'rxnNames')
0266         newModel.rxnNames=[newModel.rxnNames;filler];
0267     end
0268 end
0269 
0270 if isfield(newModel,'annotation') & isfield(newModel.annotation,'defaultLB')
0271     newLb=newModel.annotation.defaultLB;
0272 else
0273     newLb=-inf;
0274 end
0275 
0276 if isfield(rxnsToAdd,'lb')
0277     if numel(rxnsToAdd.lb)~=nRxns
0278         EM='rxnsToAdd.lb must have the same number of elements as rxnsToAdd.rxns';
0279         dispEM(EM);
0280     end
0281     %Fill with standard if it doesn't exist
0282     if ~isfield(newModel,'lb')
0283         newModel.lb=zeros(nOldRxns,1);
0284         newModel.lb(newModel.rev~=0)=newLb;
0285     end
0286     newModel.lb=[newModel.lb;rxnsToAdd.lb(:)];
0287 else
0288     %Fill with standard if it doesn't exist
0289     if isfield(newModel,'lb')
0290         I=zeros(nRxns,1);
0291         I(reversible~=0)=newLb;
0292         newModel.lb=[newModel.lb;I];
0293     end
0294 end
0295 
0296 if isfield(newModel,'annotation') & isfield(newModel.annotation,'defaultUB')
0297     newUb=newModel.annotation.defaultUB;
0298 else
0299     newUb=inf;
0300 end
0301 
0302 if isfield(rxnsToAdd,'ub')
0303     if numel(rxnsToAdd.ub)~=nRxns
0304         EM='rxnsToAdd.ub must have the same number of elements as rxnsToAdd.rxns';
0305         dispEM(EM);
0306     end
0307     %Fill with standard if it doesn't exist
0308     if ~isfield(newModel,'ub')
0309         newModel.ub=repmat(newUb,nOldRxns,1);
0310     end
0311     newModel.ub=[newModel.ub;rxnsToAdd.ub(:)];
0312 else
0313     %Fill with standard if it doesn't exist
0314     if isfield(newModel,'ub')
0315         newModel.ub=[newModel.ub;repmat(newUb,nRxns,1)];
0316     end
0317 end
0318 
0319 if isfield(rxnsToAdd,'c')
0320     if numel(rxnsToAdd.c)~=nRxns
0321         EM='rxnsToAdd.c must have the same number of elements as rxnsToAdd.rxns';
0322         dispEM(EM);
0323     end
0324     %Fill with standard if it doesn't exist
0325     if ~isfield(newModel,'c')
0326         newModel.c=zeros(nOldRxns,1);
0327     end
0328     newModel.c=[newModel.c;rxnsToAdd.c(:)];
0329 else
0330     %Fill with standard if it doesn't exist
0331     if isfield(newModel,'c')
0332         newModel.c=[newModel.c;zeros(nRxns,1)];
0333     end
0334 end
0335 
0336 if isfield(rxnsToAdd,'eccodes')
0337     rxnsToAdd.eccodes=convertCharArray(rxnsToAdd.eccodes);
0338     if numel(rxnsToAdd.eccodes)~=nRxns
0339         EM='rxnsToAdd.eccodes must have the same number of elements as rxnsToAdd.rxns';
0340         dispEM(EM);
0341     end
0342     %Fill with standard if it doesn't exist
0343     if ~isfield(newModel,'eccodes')
0344         newModel.eccodes=largeFiller;
0345     end
0346     newModel.eccodes=[newModel.eccodes;rxnsToAdd.eccodes(:)];
0347 else
0348     %Fill with standard if it doesn't exist
0349     if isfield(newModel,'eccodes')
0350         newModel.eccodes=[newModel.eccodes;filler];
0351     end
0352 end
0353 
0354 if isfield(rxnsToAdd,'subSystems')
0355     % Has to be cell array
0356     if ischar(rxnsToAdd.subSystems)
0357         rxnsToAdd.subSystems = {rxnsToAdd.subSystems};
0358     end
0359     % If all nested cells are 1x1, then unnest
0360     if all(cellfun(@(x) iscell(x) && isscalar(x), rxnsToAdd.subSystems))
0361         rxnsToAdd.subSystems = transpose([rxnsToAdd.subSystems{:}]);
0362     end
0363     % Cell array should now be as simple as possible. Check if it is nested
0364     subSysRxnsNested = any(cellfun(@(x) iscell(x), rxnsToAdd.subSystems));
0365     if isfield(newModel,'subSystems')
0366         subSysModelNested = any(cellfun(@(x) iscell(x), newModel.subSystems));
0367     else
0368         subSysModelNested = subSysRxnsNested;
0369         if subSysRxnsNested
0370             newModel.subSystems=largeCellFiller;
0371         else
0372             newModel.subSystems=largeFiller;
0373         end
0374     end
0375     if subSysRxnsNested && ~subSysModelNested
0376         % Make all existing subSystems nested
0377         newModel.subSystems = cellfun(@(x) {x}, newModel.subSystems, 'uni', 0);
0378     elseif ~subSysRxnsNested && subSysModelNested
0379         rxnsToAdd.subSystems = cellfun(@(x) {x}, rxnsToAdd.subSystems, 'uni', 0);
0380     end
0381     if numel(rxnsToAdd.subSystems)~=nRxns
0382         EM='rxnsToAdd.subSystems must have the same number of elements as rxnsToAdd.rxns';
0383         dispEM(EM);
0384     end
0385     newModel.subSystems=[newModel.subSystems;rxnsToAdd.subSystems(:)];
0386 else
0387     %Fill with standard if it does not exist
0388     if isfield(newModel,'subSystems')
0389         if any(cellfun(@(x) iscell(x), newModel.subSystems))
0390             newModel.subSystems=[newModel.subSystems;cellfiller];
0391         else
0392             newModel.subSystems=[newModel.subSystems;filler];
0393         end
0394     end
0395 end
0396 if isfield(rxnsToAdd,'rxnMiriams')
0397     if numel(rxnsToAdd.rxnMiriams)~=nRxns
0398         EM='rxnsToAdd.rxnMiriams must have the same number of elements as rxnsToAdd.rxns';
0399         dispEM(EM);
0400     end
0401     %Fill with standard if it doesn't exist
0402     if ~isfield(newModel,'rxnMiriams')
0403         newModel.rxnMiriams=cell(nOldRxns,1);
0404     end
0405     newModel.rxnMiriams=[newModel.rxnMiriams;rxnsToAdd.rxnMiriams(:)];
0406 else
0407     %Fill with standard if it doesn't exist
0408     if isfield(newModel,'rxnMiriams')
0409         newModel.rxnMiriams=[newModel.rxnMiriams;cell(nRxns,1)];
0410     end
0411 end
0412 if isfield(rxnsToAdd,'rxnComps')
0413     if numel(rxnsToAdd.rxnComps)~=nRxns
0414         EM='rxnsToAdd.rxnComps must have the same number of elements as rxnsToAdd.rxns';
0415         dispEM(EM);
0416     end
0417     %Fill with standard if it doesn't exist
0418     if ~isfield(newModel,'rxnComps')
0419         newModel.rxnComps=ones(nOldRxns,1);
0420         EM='Adding reactions with compartment information to a model without such information. All existing reactions will be assigned to the first compartment';
0421         dispEM(EM,false);
0422     end
0423     newModel.rxnComps=[newModel.rxnComps;rxnsToAdd.rxnComps(:)];
0424 else
0425     %Fill with standard if it doesn't exist
0426     if isfield(newModel,'rxnComps')
0427         newModel.rxnComps=[newModel.rxnComps;ones(nRxns,1)];
0428         %fprintf('NOTE: The added reactions will be assigned to the first
0429         %compartment\n');
0430     end
0431 end
0432 if isfield(rxnsToAdd,'grRules')
0433     rxnsToAdd.grRules=convertCharArray(rxnsToAdd.grRules);
0434     if numel(rxnsToAdd.grRules)~=nRxns
0435         EM='rxnsToAdd.grRules must have the same number of elements as rxnsToAdd.rxns';
0436         dispEM(EM);
0437     end
0438     %Fill with standard if it doesn't exist
0439     if ~isfield(newModel,'grRules')
0440         newModel.grRules=largeFiller;
0441     end
0442     newModel.grRules=[newModel.grRules;rxnsToAdd.grRules(:)];
0443 else
0444     %Fill with standard if it doesn't exist
0445     if isfield(newModel,'grRules')
0446         newModel.grRules=[newModel.grRules;filler];
0447     end
0448 end
0449 
0450 if isfield(rxnsToAdd,'rxnFrom')
0451     rxnsToAdd.rxnFrom=convertCharArray(rxnsToAdd.rxnFrom);
0452     if numel(rxnsToAdd.rxnFrom)~=nRxns
0453         EM='rxnsToAdd.rxnFrom must have the same number of elements as rxnsToAdd.rxns';
0454         dispEM(EM);
0455     end
0456     %Fill with standard if it doesn't exist
0457     if ~isfield(newModel,'rxnFrom')
0458         newModel.rxnFrom=largeFiller;
0459     end
0460     newModel.rxnFrom=[newModel.rxnFrom;rxnsToAdd.rxnFrom(:)];
0461 else
0462     %Fill with standard if it doesn't exist
0463     if isfield(newModel,'rxnFrom')
0464         newModel.rxnFrom=[newModel.rxnFrom;filler];
0465     end
0466 end
0467 
0468 if isfield(rxnsToAdd,'rxnNotes')
0469     rxnsToAdd.rxnNotes=convertCharArray(rxnsToAdd.rxnNotes);
0470     if numel(rxnsToAdd.rxnNotes)~=nRxns
0471         EM='rxnsToAdd.rxnNotes must have the same number of elements as rxnsToAdd.rxns';
0472         dispEM(EM);
0473     end
0474     %Fill with standard if it doesn't exist
0475     if ~isfield(newModel,'rxnNotes')
0476         newModel.rxnNotes=largeFiller;
0477     end
0478     newModel.rxnNotes=[newModel.rxnNotes;rxnsToAdd.rxnNotes(:)];
0479 else
0480     %Fill with standard if it doesn't exist
0481     if isfield(newModel,'rxnNotes')
0482         newModel.rxnNotes=[newModel.rxnNotes;filler];
0483     end
0484 end
0485 
0486 if isfield(rxnsToAdd,'rxnReferences')
0487     rxnsToAdd.rxnReferences=convertCharArray(rxnsToAdd.rxnReferences);
0488     if numel(rxnsToAdd.rxnReferences)~=nRxns
0489         EM='rxnsToAdd.rxnReferences must have the same number of elements as rxnsToAdd.rxns';
0490         dispEM(EM);
0491     end
0492     %Fill with standard if it doesn't exist
0493     if ~isfield(newModel,'rxnReferences')
0494         newModel.rxnReferences=largeFiller;
0495     end
0496     newModel.rxnReferences=[newModel.rxnReferences;rxnsToAdd.rxnReferences(:)];
0497 else
0498     %Fill with standard if it doesn't exist
0499     if isfield(newModel,'rxnReferences')
0500         newModel.rxnReferences=[newModel.rxnReferences;filler];
0501     end
0502 end
0503 
0504 if isfield(rxnsToAdd,'pwys')
0505     rxnsToAdd.pwys=convertCharArray(rxnsToAdd.pwys);
0506     if numel(rxnsToAdd.pwys)~=nRxns
0507         EM='rxnsToAdd.pwys must have the same number of elements as rxnsToAdd.rxns';
0508         dispEM(EM);
0509     end
0510     %Fill with standard if it doesn't exist
0511     if ~isfield(newModel,'pwys')
0512         newModel.pwys=largeFiller;
0513     end
0514     newModel.pwys=[newModel.pwys;rxnsToAdd.pwys(:)];
0515 else
0516     %Fill with standard if it doesn't exist
0517     if isfield(newModel,'pwys')
0518         newModel.pwys=[newModel.pwys;filler];
0519     end
0520 end
0521 
0522 if isfield(rxnsToAdd,'rxnConfidenceScores')
0523     if numel(rxnsToAdd.rxnConfidenceScores)~=nRxns
0524         EM='rxnsToAdd.rxnConfidenceScores must have the same number of elements as rxnsToAdd.rxns';
0525         dispEM(EM);
0526     end
0527     %Fill with standard if it doesn't exist
0528     if ~isfield(newModel,'rxnConfidenceScores')
0529         newModel.rxnConfidenceScores=NaN(nOldRxns,1);
0530     end
0531     newModel.rxnConfidenceScores=[newModel.rxnConfidenceScores;rxnsToAdd.rxnConfidenceScores(:)];
0532 else
0533     %Fill with standard if it doesn't exist
0534     if isfield(newModel,'rxnConfidenceScores')
0535         newModel.rxnConfidenceScores=[newModel.rxnConfidenceScores;NaN(nRxns,1)];
0536     end
0537 end
0538 
0539 if isfield(rxnsToAdd,'rxnDeltaG')
0540     if numel(rxnsToAdd.rxnDeltaG)~=nRxns
0541         EM='rxnsToAdd.rxnDeltaG must have the same number of elements as rxnsToAdd.rxns';
0542         dispEM(EM);
0543     end
0544     %Fill with standard if it doesn't exist
0545     if ~isfield(newModel,'rxnDeltaG')
0546         newModel.rxnDeltaG=NaN(nOldRxns,1);
0547     end
0548     newModel.rxnDeltaG=[newModel.rxnDeltaG;rxnsToAdd.rxnDeltaG(:)];
0549 else
0550     %Fill with standard if it doesn't exist
0551     if isfield(newModel,'rxnDeltaG')
0552         newModel.rxnDeltaG=[newModel.rxnDeltaG;NaN(nRxns,1)];
0553     end
0554 end
0555 
0556 
0557 %***Start parsing the equations and adding the info to the S matrix The
0558 %mets are matched to model.mets
0559 if eqnType==1
0560     [I, J]=ismember(mets,model.mets);
0561     if ~all(I)
0562         if allowNewMets==true || ischar(allowNewMets)
0563             %Add the new mets
0564             metsToAdd.mets=mets(~I);
0565             metsToAdd.metNames=metsToAdd.mets;
0566             metsToAdd.compartments=compartment;
0567             if ischar(allowNewMets)
0568                 newModel=addMets(newModel,metsToAdd,true,allowNewMets);
0569             else
0570                 newModel=addMets(newModel,metsToAdd,true);
0571             end
0572         else
0573             EM='One or more equations contain metabolites that are not in model.mets. Set allowNewMets to true to allow this function to add metabolites or use addMets to add them before calling this function. Are you sure that eqnType=1?';
0574             dispEM(EM);
0575         end
0576     end
0577     %Calculate the indexes of the metabolites and add the info
0578     metIndexes=J;
0579     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0580 end
0581 
0582 %Do some stuff that is the same for eqnType=2 and eqnType=3
0583 if eqnType==2 || eqnType==3
0584     %For later..
0585     t2=strcat(model.metNames,'***',model.comps(model.metComps));
0586 end
0587 
0588 %The mets are matched to model.metNames and assigned to "compartment"
0589 if eqnType==2
0590     %%Check that the metabolite names aren't present in the same
0591     %%compartment.
0592     %Not the neatest way maybe..
0593     t1=strcat(mets,'***',compartment);
0594     [I, J]=ismember(t1,t2);
0595     
0596     if ~all(I)
0597         if allowNewMets==true || ischar(allowNewMets)
0598             %Add the new mets
0599             metsToAdd.metNames=mets(~I);
0600             metsToAdd.compartments=compartment;
0601             if ischar(allowNewMets)
0602                 newModel=addMets(newModel,metsToAdd,true,allowNewMets);
0603             else
0604                 newModel=addMets(newModel,metsToAdd,true);
0605             end
0606         else
0607             EM='One or more equations contain metabolites that are not in model.mets. Set allowNewMets to true to allow this function to add metabolites or use addMets to add them before calling this function';
0608             dispEM(EM);
0609         end
0610     end
0611     
0612     %Calculate the indexes of the metabolites
0613     metIndexes=J;
0614     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0615 end
0616 
0617 %The equations are on the form metNames[compName]
0618 if eqnType==3
0619     %Parse the metabolite names
0620     metNames=cell(numel(mets),1);
0621     compartments=metNames;
0622     for i=1:numel(mets)
0623         starts=max(strfind(mets{i},'['));
0624         ends=max(strfind(mets{i},']'));
0625         
0626         %Check that the formatting is correct
0627         if isempty(starts) || isempty(ends) || ends<numel(mets{i})
0628             EM=['The metabolite ' mets{i} ' is not correctly formatted for eqnType=3'];
0629             dispEM(EM);
0630         end
0631         
0632         %Check that the compartment is correct
0633         compartments{i}=mets{i}(starts+1:ends-1);
0634         I=ismember(compartments(i),newModel.comps);
0635         if ~I
0636             EM=['The metabolite ' mets{i} ' has a compartment that is not in model.comps'];
0637             dispEM(EM);
0638         end
0639         metNames{i}=mets{i}(1:starts-1);
0640     end
0641     
0642     %Check if the metabolite exists already
0643     t1=strcat(metNames,'***',compartments);
0644     [I, J]=ismember(t1,t2);
0645     
0646     if ~all(I)
0647         if allowNewMets==true | ischar(allowNewMets)
0648             %Add the new mets
0649             metsToAdd.metNames=metNames(~I);
0650             metsToAdd.compartments=compartments(~I);
0651             if ischar(allowNewMets)
0652                 newModel=addMets(newModel,metsToAdd,true,allowNewMets);
0653             else
0654                 newModel=addMets(newModel,metsToAdd,true);
0655             end
0656         else
0657             EM='One or more equations contain metabolites that are not in model.metNames. Set allowNewMets to true to allow this function to add metabolites or use addMets to add them before calling this function';
0658             dispEM(EM);
0659         end
0660     end
0661     
0662     %Calculate the indexes of the metabolites
0663     metIndexes=J;
0664     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0665 end
0666 
0667 %Add the info to the stoichiometric matrix
0668 newModel.S=[newModel.S sparse(size(newModel.S,1),nRxns)];
0669 
0670 for i=1:nRxns
0671     newModel.S(metIndexes,nOldRxns+i)=S(:,i);
0672     %Parse the grRules and check whether all genes in grRules appear in
0673     %genes
0674     if isfield(newModel,'grRules')
0675         rule=newModel.grRules{nOldRxns+i};
0676         rule=strrep(rule,'(','');
0677         rule=strrep(rule,')','');
0678         rule=strrep(rule,' or ',' ');
0679         rule=strrep(rule,' and ',' ');
0680         genes=regexp(rule,' ','split');
0681         [I, J]=ismember(genes,newModel.genes);
0682         if ~all(I) && any(rule)
0683             EM=['Not all genes for reaction ' rxnsToAdd.rxns{i} ' were found in model.genes. If needed, add genes with addGenesRaven before calling this function, or set the ''allowNewGenes'' flag to true'];
0684             dispEM(EM);
0685         end
0686     end
0687 end
0688 
0689 %Make temporary minimal model structure with only new rxns, to parse to
0690 %standardizeGrRules
0691 newRxnsModel.genes=newModel.genes;
0692 newRxnsModel.grRules=newModel.grRules(length(model.rxns)+1:end);
0693 newRxnsModel.rxns=newModel.rxns(length(model.rxns)+1:end);
0694 
0695 %Fix grRules and reconstruct rxnGeneMat
0696 [grRules,rxnGeneMat] = standardizeGrRules(newRxnsModel,true);
0697 newModel.rxnGeneMat = [newModel.rxnGeneMat; rxnGeneMat];
0698 newModel.grRules = [newModel.grRules(1:nOldRxns); grRules];
0699 end