addRxns

 addRxns
   Adds reactions to a model

   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
                               (opt, default '')
            lb                 vector with the lower bounds (opt, default
                               to 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
                               reactions to model.annotations.defaultLB or
                               -inf)
            ub                 vector with the upper bounds (opt, default
                               to model.annotations.defaultUB or inf)
            c                  vector with the objective function
                               coefficients (opt, default 0)
            eccodes            cell array with the EC-numbers for each
                               reactions. Delimit several EC-numbers with
                               ";" (opt, default '')
            subSystems         cell array with the subsystems for each
                               reaction (opt, default '')
            grRules            cell array with the gene-reaction
                               relationship for each reaction. For example
                               "(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 (opt, default '')
            rxnMiriams         cell array with Miriam structures (opt,
                               default [])
            rxnComps           cell array with compartments (as in
                               model.comps) (opt, default {})
            rxnNotes           cell array with reaction notes (opt,
                               default '')
            rxnDeltaG          Gibbs free energy at biochemical standard
                               condition in kJ/mole (opt, default NaN)
            rxnReferences      cell array with reaction references (opt,
                               default '')
            rxnConfidenceScores   vector with reaction confidence scores
                               (opt, 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 (opt 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 (opt, default false)
   allowNewGenes    true if the functions is allowed to add new genes
                    (opt, default false)

   newModel         an updated model structure

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

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