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=cellfun(@(x) cell(0,0),filler,'UniformOutput',false);
0231 largeFiller=cell(nOldRxns,1);
0232 largeFiller(:)={''};
0233 celllargefiller=cellfun(@(x) cell(0,0),largeFiller,'UniformOutput',false);
0234 
0235 %***Add everything to the model except for the equations.
0236 if numel(rxnsToAdd.equations)~=nRxns
0237     EM='rxnsToAdd.equations must have the same number of elements as rxnsToAdd.rxns';
0238     dispEM(EM);
0239 end
0240 
0241 %Parse the equations. This is done at this early stage since I need the
0242 %reversibility info
0243 [S, mets, badRxns, reversible]=constructS(rxnsToAdd.equations);
0244 EM='The following equations have one or more metabolites both as substrate and product. Only the net equations will be added:';
0245 dispEM(EM,false,rxnsToAdd.rxns(badRxns));
0246 
0247 newModel.rev=[newModel.rev;reversible];
0248 newModel.rxns=[newModel.rxns;rxnsToAdd.rxns(:)];
0249 
0250 if isfield(rxnsToAdd,'rxnNames')
0251     rxnsToAdd.rxnNames=convertCharArray(rxnsToAdd.rxnNames);
0252     if numel(rxnsToAdd.rxnNames)~=nRxns
0253         EM='rxnsToAdd.rxnNames must have the same number of elements as rxnsToAdd.rxns';
0254         dispEM(EM);
0255     end
0256     %Fill with standard if it doesn't exist
0257     if ~isfield(newModel,'rxnNames')
0258         newModel.rxnNames=largeFiller;
0259     end
0260     newModel.rxnNames=[newModel.rxnNames;rxnsToAdd.rxnNames(:)];
0261 else
0262     %Fill with standard if it doesn't exist
0263     if isfield(newModel,'rxnNames')
0264         newModel.rxnNames=[newModel.rxnNames;filler];
0265     end
0266 end
0267 
0268 if isfield(newModel,'annotation') & isfield(newModel.annotation,'defaultLB')
0269     newLb=newModel.annotation.defaultLB;
0270 else
0271     newLb=-inf;
0272 end
0273 
0274 if isfield(rxnsToAdd,'lb')
0275     if numel(rxnsToAdd.lb)~=nRxns
0276         EM='rxnsToAdd.lb must have the same number of elements as rxnsToAdd.rxns';
0277         dispEM(EM);
0278     end
0279     %Fill with standard if it doesn't exist
0280     if ~isfield(newModel,'lb')
0281         newModel.lb=zeros(nOldRxns,1);
0282         newModel.lb(newModel.rev~=0)=newLb;
0283     end
0284     newModel.lb=[newModel.lb;rxnsToAdd.lb(:)];
0285 else
0286     %Fill with standard if it doesn't exist
0287     if isfield(newModel,'lb')
0288         I=zeros(nRxns,1);
0289         I(reversible~=0)=newLb;
0290         newModel.lb=[newModel.lb;I];
0291     end
0292 end
0293 
0294 if isfield(newModel,'annotation') & isfield(newModel.annotation,'defaultUB')
0295     newUb=newModel.annotation.defaultUB;
0296 else
0297     newUb=inf;
0298 end
0299 
0300 if isfield(rxnsToAdd,'ub')
0301     if numel(rxnsToAdd.ub)~=nRxns
0302         EM='rxnsToAdd.ub must have the same number of elements as rxnsToAdd.rxns';
0303         dispEM(EM);
0304     end
0305     %Fill with standard if it doesn't exist
0306     if ~isfield(newModel,'ub')
0307         newModel.ub=repmat(newUb,nOldRxns,1);
0308     end
0309     newModel.ub=[newModel.ub;rxnsToAdd.ub(:)];
0310 else
0311     %Fill with standard if it doesn't exist
0312     if isfield(newModel,'ub')
0313         newModel.ub=[newModel.ub;repmat(newUb,nRxns,1)];
0314     end
0315 end
0316 
0317 if isfield(rxnsToAdd,'c')
0318     if numel(rxnsToAdd.c)~=nRxns
0319         EM='rxnsToAdd.c must have the same number of elements as rxnsToAdd.rxns';
0320         dispEM(EM);
0321     end
0322     %Fill with standard if it doesn't exist
0323     if ~isfield(newModel,'c')
0324         newModel.c=zeros(nOldRxns,1);
0325     end
0326     newModel.c=[newModel.c;rxnsToAdd.c(:)];
0327 else
0328     %Fill with standard if it doesn't exist
0329     if isfield(newModel,'c')
0330         newModel.c=[newModel.c;zeros(nRxns,1)];
0331     end
0332 end
0333 
0334 if isfield(rxnsToAdd,'eccodes')
0335     rxnsToAdd.eccodes=convertCharArray(rxnsToAdd.eccodes);
0336     if numel(rxnsToAdd.eccodes)~=nRxns
0337         EM='rxnsToAdd.eccodes must have the same number of elements as rxnsToAdd.rxns';
0338         dispEM(EM);
0339     end
0340     %Fill with standard if it doesn't exist
0341     if ~isfield(newModel,'eccodes')
0342         newModel.eccodes=largeFiller;
0343     end
0344     newModel.eccodes=[newModel.eccodes;rxnsToAdd.eccodes(:)];
0345 else
0346     %Fill with standard if it doesn't exist
0347     if isfield(newModel,'eccodes')
0348         newModel.eccodes=[newModel.eccodes;filler];
0349     end
0350 end
0351 
0352 if isfield(rxnsToAdd,'subSystems')
0353     if numel(rxnsToAdd.subSystems)~=nRxns
0354         EM='rxnsToAdd.subSystems must have the same number of elements as rxnsToAdd.rxns';
0355         dispEM(EM);
0356     end
0357     for i=1:numel(rxnsToAdd.subSystems)
0358         if ischar(rxnsToAdd.subSystems{i})
0359             rxnsToAdd.subSystems{i}=rxnsToAdd.subSystems(i);
0360         end
0361     end
0362     %Fill with standard if it doesn't exist
0363     if ~isfield(newModel,'subSystems')
0364         newModel.subSystems=celllargefiller;
0365     end
0366     newModel.subSystems=[newModel.subSystems;rxnsToAdd.subSystems];
0367 else
0368     %Fill with standard if it doesn't exist
0369     if isfield(newModel,'subSystems')
0370         newModel.subSystems=[newModel.subSystems;cellfiller];
0371     end
0372 end
0373 if isfield(rxnsToAdd,'rxnMiriams')
0374     if numel(rxnsToAdd.rxnMiriams)~=nRxns
0375         EM='rxnsToAdd.rxnMiriams must have the same number of elements as rxnsToAdd.rxns';
0376         dispEM(EM);
0377     end
0378     %Fill with standard if it doesn't exist
0379     if ~isfield(newModel,'rxnMiriams')
0380         newModel.rxnMiriams=cell(nOldRxns,1);
0381     end
0382     newModel.rxnMiriams=[newModel.rxnMiriams;rxnsToAdd.rxnMiriams(:)];
0383 else
0384     %Fill with standard if it doesn't exist
0385     if isfield(newModel,'rxnMiriams')
0386         newModel.rxnMiriams=[newModel.rxnMiriams;cell(nRxns,1)];
0387     end
0388 end
0389 if isfield(rxnsToAdd,'rxnComps')
0390     if numel(rxnsToAdd.rxnComps)~=nRxns
0391         EM='rxnsToAdd.rxnComps must have the same number of elements as rxnsToAdd.rxns';
0392         dispEM(EM);
0393     end
0394     %Fill with standard if it doesn't exist
0395     if ~isfield(newModel,'rxnComps')
0396         newModel.rxnComps=ones(nOldRxns,1);
0397         EM='Adding reactions with compartment information to a model without such information. All existing reactions will be assigned to the first compartment';
0398         dispEM(EM,false);
0399     end
0400     newModel.rxnComps=[newModel.rxnComps;rxnsToAdd.rxnComps(:)];
0401 else
0402     %Fill with standard if it doesn't exist
0403     if isfield(newModel,'rxnComps')
0404         newModel.rxnComps=[newModel.rxnComps;ones(nRxns,1)];
0405         %fprintf('NOTE: The added reactions will be assigned to the first
0406         %compartment\n');
0407     end
0408 end
0409 if isfield(rxnsToAdd,'grRules')
0410     rxnsToAdd.grRules=convertCharArray(rxnsToAdd.grRules);
0411     if numel(rxnsToAdd.grRules)~=nRxns
0412         EM='rxnsToAdd.grRules must have the same number of elements as rxnsToAdd.rxns';
0413         dispEM(EM);
0414     end
0415     %Fill with standard if it doesn't exist
0416     if ~isfield(newModel,'grRules')
0417         newModel.grRules=largeFiller;
0418     end
0419     newModel.grRules=[newModel.grRules;rxnsToAdd.grRules(:)];
0420 else
0421     %Fill with standard if it doesn't exist
0422     if isfield(newModel,'grRules')
0423         newModel.grRules=[newModel.grRules;filler];
0424     end
0425 end
0426 
0427 if isfield(rxnsToAdd,'rxnFrom')
0428     rxnsToAdd.rxnFrom=convertCharArray(rxnsToAdd.rxnFrom);
0429     if numel(rxnsToAdd.rxnFrom)~=nRxns
0430         EM='rxnsToAdd.rxnFrom must have the same number of elements as rxnsToAdd.rxns';
0431         dispEM(EM);
0432     end
0433     %Fill with standard if it doesn't exist
0434     if ~isfield(newModel,'rxnFrom')
0435         newModel.rxnFrom=largeFiller;
0436     end
0437     newModel.rxnFrom=[newModel.rxnFrom;rxnsToAdd.rxnFrom(:)];
0438 else
0439     %Fill with standard if it doesn't exist
0440     if isfield(newModel,'rxnFrom')
0441         newModel.rxnFrom=[newModel.rxnFrom;filler];
0442     end
0443 end
0444 
0445 if isfield(rxnsToAdd,'rxnNotes')
0446     rxnsToAdd.rxnNotes=convertCharArray(rxnsToAdd.rxnNotes);
0447     if numel(rxnsToAdd.rxnNotes)~=nRxns
0448         EM='rxnsToAdd.rxnNotes must have the same number of elements as rxnsToAdd.rxns';
0449         dispEM(EM);
0450     end
0451     %Fill with standard if it doesn't exist
0452     if ~isfield(newModel,'rxnNotes')
0453         newModel.rxnNotes=largeFiller;
0454     end
0455     newModel.rxnNotes=[newModel.rxnNotes;rxnsToAdd.rxnNotes(:)];
0456 else
0457     %Fill with standard if it doesn't exist
0458     if isfield(newModel,'rxnNotes')
0459         newModel.rxnNotes=[newModel.rxnNotes;filler];
0460     end
0461 end
0462 
0463 if isfield(rxnsToAdd,'rxnReferences')
0464     rxnsToAdd.rxnReferences=convertCharArray(rxnsToAdd.rxnReferences);
0465     if numel(rxnsToAdd.rxnReferences)~=nRxns
0466         EM='rxnsToAdd.rxnReferences must have the same number of elements as rxnsToAdd.rxns';
0467         dispEM(EM);
0468     end
0469     %Fill with standard if it doesn't exist
0470     if ~isfield(newModel,'rxnReferences')
0471         newModel.rxnReferences=largeFiller;
0472     end
0473     newModel.rxnReferences=[newModel.rxnReferences;rxnsToAdd.rxnReferences(:)];
0474 else
0475     %Fill with standard if it doesn't exist
0476     if isfield(newModel,'rxnReferences')
0477         newModel.rxnReferences=[newModel.rxnReferences;filler];
0478     end
0479 end
0480 
0481 if isfield(rxnsToAdd,'pwys')
0482     rxnsToAdd.pwys=convertCharArray(rxnsToAdd.pwys);
0483     if numel(rxnsToAdd.pwys)~=nRxns
0484         EM='rxnsToAdd.pwys must have the same number of elements as rxnsToAdd.rxns';
0485         dispEM(EM);
0486     end
0487     %Fill with standard if it doesn't exist
0488     if ~isfield(newModel,'pwys')
0489         newModel.pwys=largeFiller;
0490     end
0491     newModel.pwys=[newModel.pwys;rxnsToAdd.pwys(:)];
0492 else
0493     %Fill with standard if it doesn't exist
0494     if isfield(newModel,'pwys')
0495         newModel.pwys=[newModel.pwys;filler];
0496     end
0497 end
0498 
0499 if isfield(rxnsToAdd,'rxnConfidenceScores')
0500     if numel(rxnsToAdd.rxnConfidenceScores)~=nRxns
0501         EM='rxnsToAdd.rxnConfidenceScores must have the same number of elements as rxnsToAdd.rxns';
0502         dispEM(EM);
0503     end
0504     %Fill with standard if it doesn't exist
0505     if ~isfield(newModel,'rxnConfidenceScores')
0506         newModel.rxnConfidenceScores=NaN(nOldRxns,1);
0507     end
0508     newModel.rxnConfidenceScores=[newModel.rxnConfidenceScores;rxnsToAdd.rxnConfidenceScores(:)];
0509 else
0510     %Fill with standard if it doesn't exist
0511     if isfield(newModel,'rxnConfidenceScores')
0512         newModel.rxnConfidenceScores=[newModel.rxnConfidenceScores;NaN(nRxns,1)];
0513     end
0514 end
0515 
0516 if isfield(rxnsToAdd,'rxnDeltaG')
0517     if numel(rxnsToAdd.rxnDeltaG)~=nRxns
0518         EM='rxnsToAdd.rxnDeltaG must have the same number of elements as rxnsToAdd.rxns';
0519         dispEM(EM);
0520     end
0521     %Fill with standard if it doesn't exist
0522     if ~isfield(newModel,'rxnDeltaG')
0523         newModel.rxnDeltaG=NaN(nOldRxns,1);
0524     end
0525     newModel.rxnDeltaG=[newModel.rxnDeltaG;rxnsToAdd.rxnDeltaG(:)];
0526 else
0527     %Fill with standard if it doesn't exist
0528     if isfield(newModel,'rxnDeltaG')
0529         newModel.rxnDeltaG=[newModel.rxnDeltaG;NaN(nRxns,1)];
0530     end
0531 end
0532 
0533 
0534 %***Start parsing the equations and adding the info to the S matrix The
0535 %mets are matched to model.mets
0536 if eqnType==1
0537     [I, J]=ismember(mets,model.mets);
0538     if ~all(I)
0539         if allowNewMets==true || ischar(allowNewMets)
0540             %Add the new mets
0541             metsToAdd.mets=mets(~I);
0542             metsToAdd.metNames=metsToAdd.mets;
0543             metsToAdd.compartments=compartment;
0544             if ischar(allowNewMets)
0545                 newModel=addMets(newModel,metsToAdd,true,allowNewMets);
0546             else
0547                 newModel=addMets(newModel,metsToAdd,true);
0548             end
0549         else
0550             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?';
0551             dispEM(EM);
0552         end
0553     end
0554     %Calculate the indexes of the metabolites and add the info
0555     metIndexes=J;
0556     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0557 end
0558 
0559 %Do some stuff that is the same for eqnType=2 and eqnType=3
0560 if eqnType==2 || eqnType==3
0561     %For later..
0562     t2=strcat(model.metNames,'***',model.comps(model.metComps));
0563 end
0564 
0565 %The mets are matched to model.metNames and assigned to "compartment"
0566 if eqnType==2
0567     %%Check that the metabolite names aren't present in the same
0568     %%compartment.
0569     %Not the neatest way maybe..
0570     t1=strcat(mets,'***',compartment);
0571     [I, J]=ismember(t1,t2);
0572     
0573     if ~all(I)
0574         if allowNewMets==true || ischar(allowNewMets)
0575             %Add the new mets
0576             metsToAdd.metNames=mets(~I);
0577             metsToAdd.compartments=compartment;
0578             if ischar(allowNewMets)
0579                 newModel=addMets(newModel,metsToAdd,true,allowNewMets);
0580             else
0581                 newModel=addMets(newModel,metsToAdd,true);
0582             end
0583         else
0584             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';
0585             dispEM(EM);
0586         end
0587     end
0588     
0589     %Calculate the indexes of the metabolites
0590     metIndexes=J;
0591     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0592 end
0593 
0594 %The equations are on the form metNames[compName]
0595 if eqnType==3
0596     %Parse the metabolite names
0597     metNames=cell(numel(mets),1);
0598     compartments=metNames;
0599     for i=1:numel(mets)
0600         starts=max(strfind(mets{i},'['));
0601         ends=max(strfind(mets{i},']'));
0602         
0603         %Check that the formatting is correct
0604         if isempty(starts) || isempty(ends) || ends<numel(mets{i})
0605             EM=['The metabolite ' mets{i} ' is not correctly formatted for eqnType=3'];
0606             dispEM(EM);
0607         end
0608         
0609         %Check that the compartment is correct
0610         compartments{i}=mets{i}(starts+1:ends-1);
0611         I=ismember(compartments(i),newModel.comps);
0612         if ~I
0613             EM=['The metabolite ' mets{i} ' has a compartment that is not in model.comps'];
0614             dispEM(EM);
0615         end
0616         metNames{i}=mets{i}(1:starts-1);
0617     end
0618     
0619     %Check if the metabolite exists already
0620     t1=strcat(metNames,'***',compartments);
0621     [I, J]=ismember(t1,t2);
0622     
0623     if ~all(I)
0624         if allowNewMets==true | ischar(allowNewMets)
0625             %Add the new mets
0626             metsToAdd.metNames=metNames(~I);
0627             metsToAdd.compartments=compartments(~I);
0628             if ischar(allowNewMets)
0629                 newModel=addMets(newModel,metsToAdd,true,allowNewMets);
0630             else
0631                 newModel=addMets(newModel,metsToAdd,true);
0632             end
0633         else
0634             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';
0635             dispEM(EM);
0636         end
0637     end
0638     
0639     %Calculate the indexes of the metabolites
0640     metIndexes=J;
0641     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
0642 end
0643 
0644 %Add the info to the stoichiometric matrix
0645 newModel.S=[newModel.S sparse(size(newModel.S,1),nRxns)];
0646 
0647 for i=1:nRxns
0648     newModel.S(metIndexes,nOldRxns+i)=S(:,i);
0649     %Parse the grRules and check whether all genes in grRules appear in
0650     %genes
0651     if isfield(newModel,'grRules')
0652         rule=newModel.grRules{nOldRxns+i};
0653         rule=strrep(rule,'(','');
0654         rule=strrep(rule,')','');
0655         rule=strrep(rule,' or ',' ');
0656         rule=strrep(rule,' and ',' ');
0657         genes=regexp(rule,' ','split');
0658         [I, J]=ismember(genes,newModel.genes);
0659         if ~all(I) && any(rule)
0660             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'];
0661             dispEM(EM);
0662         end
0663     end
0664 end
0665 
0666 %Make temporary minimal model structure with only new rxns, to parse to
0667 %standardizeGrRules
0668 newRxnsModel.genes=newModel.genes;
0669 newRxnsModel.grRules=newModel.grRules(length(model.rxns)+1:end);
0670 newRxnsModel.rxns=newModel.rxns(length(model.rxns)+1:end);
0671 
0672 %Fix grRules and reconstruct rxnGeneMat
0673 [grRules,rxnGeneMat] = standardizeGrRules(newRxnsModel,true);
0674 newModel.rxnGeneMat = [newModel.rxnGeneMat; rxnGeneMat];
0675 newModel.grRules = [newModel.grRules(1:nOldRxns); grRules];
0676 end