importModel

 importModel
   Import a constraint-based model from an SBML file.

 Input:
   fileName        a SBML file to import. A dialog window will open if
                   no file name is specified.
   removeExcMets   true if exchange metabolites should be removed. This is
                   needed to be able to run simulations, but it could also
                   be done using simplifyModel at a later stage (optional,
                   default true)
   removePrefix    true if identifier prefixes should be removed when
                   loading the model: G_ for genes, R_ for reactions,
                   M_ for metabolites, and C_ for compartments. These are
                   only removed if all identifiers of a certain type
                   contain the prefix. (optional, default true)
   supressWarnings true if warnings regarding the model structure should
                   be supressed (optional, default false)

 Output:
   model
       id               model ID
       name             name of model contents
       annotation       additional information about model
       rxns             reaction ids
       mets             metabolite ids
       S                stoichiometric matrix
       lb               lower bounds
       ub               upper bounds
       rev              reversibility vector
       c                objective coefficients
       b                equality constraints for the metabolite equations
       comps            compartment ids
       compNames        compartment names
       compOutside      the id (as in comps) for the compartment
                        surrounding each of the compartments
       compMiriams      structure with MIRIAM information about the
                        compartments
       rxnNames         reaction description
       rxnComps         compartments for reactions
       grRules          reaction to gene rules in text form
       rxnGeneMat       reaction-to-gene mapping in sparse matrix form
       subSystems       subsystem name for each reaction
       eccodes          EC-codes for the reactions
       rxnMiriams       structure with MIRIAM information about the reactions
       rxnNotes         reaction notes
       rxnReferences    reaction references
       rxnConfidenceScores reaction confidence scores
       genes            list of all genes
       geneComps        compartments for genes
       geneMiriams      structure with MIRIAM information about the genes
       geneShortNames   gene alternative names (e.g. ERG10)
       proteins         protein associated to each gene
       metNames         metabolite description
       metComps         compartments for metabolites
       inchis           InChI-codes for metabolites
       metFormulas      metabolite chemical formula
       metMiriams       structure with MIRIAM information about the metabolites
       metCharges       metabolite charge
       unconstrained    true if the metabolite is an exchange metabolite

 Note: A number of consistency checks are performed in order to ensure that the
 model is valid. Take these warnings seriously and modify the model
 structure to solve them.

 Usage: model = importModel(fileName, removeExcMets, removePrefix, supressWarnings)

0001 function model=importModel(fileName,removeExcMets,removePrefix,supressWarnings)
0002 % importModel
0003 %   Import a constraint-based model from an SBML file.
0004 %
0005 % Input:
0006 %   fileName        a SBML file to import. A dialog window will open if
0007 %                   no file name is specified.
0008 %   removeExcMets   true if exchange metabolites should be removed. This is
0009 %                   needed to be able to run simulations, but it could also
0010 %                   be done using simplifyModel at a later stage (optional,
0011 %                   default true)
0012 %   removePrefix    true if identifier prefixes should be removed when
0013 %                   loading the model: G_ for genes, R_ for reactions,
0014 %                   M_ for metabolites, and C_ for compartments. These are
0015 %                   only removed if all identifiers of a certain type
0016 %                   contain the prefix. (optional, default true)
0017 %   supressWarnings true if warnings regarding the model structure should
0018 %                   be supressed (optional, default false)
0019 %
0020 % Output:
0021 %   model
0022 %       id               model ID
0023 %       name             name of model contents
0024 %       annotation       additional information about model
0025 %       rxns             reaction ids
0026 %       mets             metabolite ids
0027 %       S                stoichiometric matrix
0028 %       lb               lower bounds
0029 %       ub               upper bounds
0030 %       rev              reversibility vector
0031 %       c                objective coefficients
0032 %       b                equality constraints for the metabolite equations
0033 %       comps            compartment ids
0034 %       compNames        compartment names
0035 %       compOutside      the id (as in comps) for the compartment
0036 %                        surrounding each of the compartments
0037 %       compMiriams      structure with MIRIAM information about the
0038 %                        compartments
0039 %       rxnNames         reaction description
0040 %       rxnComps         compartments for reactions
0041 %       grRules          reaction to gene rules in text form
0042 %       rxnGeneMat       reaction-to-gene mapping in sparse matrix form
0043 %       subSystems       subsystem name for each reaction
0044 %       eccodes          EC-codes for the reactions
0045 %       rxnMiriams       structure with MIRIAM information about the reactions
0046 %       rxnNotes         reaction notes
0047 %       rxnReferences    reaction references
0048 %       rxnConfidenceScores reaction confidence scores
0049 %       genes            list of all genes
0050 %       geneComps        compartments for genes
0051 %       geneMiriams      structure with MIRIAM information about the genes
0052 %       geneShortNames   gene alternative names (e.g. ERG10)
0053 %       proteins         protein associated to each gene
0054 %       metNames         metabolite description
0055 %       metComps         compartments for metabolites
0056 %       inchis           InChI-codes for metabolites
0057 %       metFormulas      metabolite chemical formula
0058 %       metMiriams       structure with MIRIAM information about the metabolites
0059 %       metCharges       metabolite charge
0060 %       unconstrained    true if the metabolite is an exchange metabolite
0061 %
0062 % Note: A number of consistency checks are performed in order to ensure that the
0063 % model is valid. Take these warnings seriously and modify the model
0064 % structure to solve them.
0065 %
0066 % Usage: model = importModel(fileName, removeExcMets, removePrefix, supressWarnings)
0067 
0068 if nargin<1 || isempty(fileName)
0069     [fileName, pathName] = uigetfile({'*.xml;*.sbml'}, 'Please select the model file');
0070     if fileName == 0
0071         error('You should select a model file')
0072     else
0073         fileName = fullfile(pathName,fileName);
0074     end
0075 end
0076 fileName=char(fileName);
0077 if nargin<2 || isempty(removeExcMets)
0078     removeExcMets=true;
0079 end
0080 
0081 if nargin<3 || isempty(removePrefix)
0082     removePrefix=true;
0083 end
0084 
0085 if nargin<4
0086     supressWarnings=false;
0087 end
0088 
0089 fileName=checkFileExistence(fileName,1);
0090 % If path contains non-ASCII characters, copy file to tempdir first, as
0091 % libSBML is known to have problems with this on Windows:
0092 % https://sbml.org/software/libsbml/libsbml-docs/known-pitfalls/#matlab-on-windows-has-issues-with-unicode-filenames
0093 if ispc && any(double(fileName)>128)
0094     [~,originalFile,ext] = fileparts(fileName);
0095     tempFile = fullfile(tempdir,[originalFile ext]);
0096     copyfile(fileName,tempFile);
0097     fileName = tempFile;
0098 end
0099 
0100 %This is to match the order of the fields to those you get from importing
0101 %from Excel
0102 model=[];
0103 model.id=[];
0104 model.name=[];
0105 model.annotation=[];
0106 model.rxns={};
0107 model.mets={};
0108 model.S=[];
0109 model.lb=[];
0110 model.ub=[];
0111 model.rev=[];
0112 model.c=[];
0113 model.b=[];
0114 model.comps={};
0115 model.compNames={};
0116 model.compOutside={};
0117 model.compMiriams={};
0118 model.rxnNames={};
0119 model.rxnComps=[];
0120 model.grRules={};
0121 model.rxnGeneMat=[];
0122 model.subSystems={};
0123 model.eccodes={};
0124 model.rxnMiriams={};
0125 model.rxnNotes={};
0126 model.rxnReferences={};
0127 model.rxnConfidenceScores=[];
0128 model.genes={};
0129 model.geneComps=[];
0130 model.geneMiriams={};
0131 model.geneShortNames={};
0132 model.proteins={};
0133 model.metNames={};
0134 model.metComps=[];
0135 model.inchis={};
0136 model.metFormulas={};
0137 model.metMiriams={};
0138 model.metCharges=[];
0139 model.unconstrained=[];
0140 
0141 %Load the model using libSBML
0142 [modelSBML,errorMsg] = TranslateSBML_RAVEN(fileName,0,0,[1 1]);
0143 if exist('tempFile','var')
0144     delete(tempFile)
0145 end
0146 
0147 if isempty(modelSBML)
0148     EM=['There is a problem with the SBML file. Try using the SBML Validator at http://sbml.org/Facilities/Validator.\nlibSBML reports: ', errorMsg.message];
0149     dispEM(EM);
0150 end
0151 
0152 %Retrieve compartment names and IDs
0153 compartmentNames=cell(numel(modelSBML.compartment),1);
0154 compartmentIDs=cell(numel(modelSBML.compartment),1);
0155 compartmentOutside=cell(numel(modelSBML.compartment),1);
0156 compartmentMiriams=cell(numel(modelSBML.compartment),1);
0157 
0158 if isfield(modelSBML.compartment,'sboTerm') && isscalar(unique([modelSBML.compartment.sboTerm]))
0159     %If all the SBO terms are identical, don't add them to compMiriams
0160     modelSBML.compartment = rmfield(modelSBML.compartment,'sboTerm');
0161 end
0162 
0163 for i=1:numel(modelSBML.compartment)
0164     compartmentNames{i}=modelSBML.compartment(i).name;
0165     compartmentIDs{i}=modelSBML.compartment(i).id;
0166     if isfield(modelSBML.compartment(i),'outside')
0167         if ~isempty(modelSBML.compartment(i).outside)
0168             compartmentOutside{i}=modelSBML.compartment(i).outside;
0169         else
0170             compartmentOutside{i}='';
0171         end
0172     else
0173         compartmentOutside{i}=[];
0174     end
0175 
0176     if isfield(modelSBML.compartment(i),'annotation')
0177         compartmentMiriams{i}=parseMiriam(modelSBML.compartment(i).annotation);
0178     else
0179         compartmentMiriams{i}=[];
0180     end
0181 
0182     if isfield(modelSBML.compartment(i),'sboTerm') && ~(modelSBML.compartment(i).sboTerm==-1)
0183         compartmentMiriams{i} = addSBOtoMiriam(compartmentMiriams{i},modelSBML.compartment(i).sboTerm);
0184     end
0185 end
0186 
0187 %If there are no compartment names then use compartment id as name
0188 if all(cellfun(@isempty,compartmentNames))
0189     compartmentNames=compartmentIDs;
0190 end
0191 
0192 %Retrieve info on metabolites, genes, complexes
0193 metaboliteNames={};
0194 metaboliteIDs={};
0195 metaboliteCompartments={};
0196 metaboliteUnconstrained=[];
0197 metaboliteFormula={};
0198 metaboliteInChI={};
0199 metaboliteMiriams={};
0200 metaboliteCharges=[];
0201 
0202 geneNames={};
0203 geneIDs={};
0204 geneMiriams={};
0205 geneShortNames={};
0206 proteins={};
0207 geneCompartments={};
0208 complexIDs={};
0209 complexNames={};
0210 
0211 %If the file is not a COBRA Toolbox model. According to the format
0212 %specified in the yeast consensus model both metabolites and genes are a
0213 %type of 'species'. The metabolites have names starting with 'M_' and genes
0214 %with 'E_'
0215 geneSBOs = [];
0216 metSBOs = [];
0217 %Regex of compartment names, later to be used to remove from metabolite
0218 %names if present as suffix.
0219 regexCompNames = ['\s?\[((' strjoin({modelSBML.compartment.name},')|(') '))\]$'];
0220 for i=1:numel(modelSBML.species)
0221     if length(modelSBML.species(i).id)>=2 && strcmpi(modelSBML.species(i).id(1:2),'E_')
0222         geneNames{numel(geneNames)+1,1}=modelSBML.species(i).name;
0223 
0224         %The "E_" is included in the ID. This is because it's only used
0225         %internally in this file and it makes the matching a little
0226         %smoother
0227         geneIDs{numel(geneIDs)+1,1}=modelSBML.species(i).id;
0228         geneCompartments{numel(geneCompartments)+1,1}=modelSBML.species(i).compartment;
0229 
0230         %Get Miriam structure
0231         if isfield(modelSBML.species(i),'annotation')
0232             %Get Miriam info
0233             geneMiriam=parseMiriam(modelSBML.species(i).annotation);
0234             geneMiriams{numel(geneMiriams)+1,1}=geneMiriam;
0235         else
0236             geneMiriams{numel(geneMiriams)+1,1}=[];
0237         end
0238 
0239         %Protein short names (for example ERG10) are saved as SHORT
0240         %NAME: NAME in the notes-section of metabolites for SBML Level
0241         %2 and as PROTEIN_ASSOCIATION for each reaction in SBML Level 2
0242         %COBRA Toolbox format. For now only the SHORT NAME is loaded
0243         %and no mapping takes place
0244         if isfield(modelSBML.species(i),'notes')
0245             geneShortNames{numel(geneShortNames)+1,1}=parseNote(modelSBML.species(i).notes,'SHORT NAME');
0246         else
0247             geneShortNames{numel(geneShortNames)+1,1}='';
0248         end
0249 
0250         %Get SBO term
0251         if isfield(modelSBML.species(i),'sboTerm') && ~(modelSBML.species(i).sboTerm==-1)
0252             geneSBOs(end+1,1) = modelSBML.species(i).sboTerm;
0253         end
0254     elseif length(modelSBML.species(i).id)>=2 && strcmpi(modelSBML.species(i).id(1:3),'Cx_')
0255         %If it's a complex keep the ID and name
0256         complexIDs=[complexIDs;modelSBML.species(i).id];
0257         complexNames=[complexNames;modelSBML.species(i).name];
0258     else
0259         %If it is not gene or complex, then it must be a metabolite
0260         metaboliteNames{numel(metaboliteNames)+1,1}=modelSBML.species(i).name;
0261         metaboliteIDs{numel(metaboliteIDs)+1,1}=modelSBML.species(i).id;
0262         metaboliteCompartments{numel(metaboliteCompartments)+1,1}=modelSBML.species(i).compartment;
0263         metaboliteUnconstrained(numel(metaboliteUnconstrained)+1,1)=modelSBML.species(i).boundaryCondition;
0264 
0265         %For each metabolite retrieve the formula and the InChI code if
0266         %available First add the InChI code and the formula from the
0267         %InChI. This allows for overwriting the formula by setting the
0268         %actual formula field
0269         if ~isempty(modelSBML.species(i).annotation)
0270             %Get the formula if available
0271             startString='>InChI=';
0272             endString='</in:inchi>';
0273             formStart=strfind(modelSBML.species(i).annotation,startString);
0274             if isempty(formStart)
0275                 startString='InChI=';
0276                 endString='"/>';
0277             end
0278             formStart=strfind(modelSBML.species(i).annotation,startString);
0279             if ~isempty(formStart)
0280                 formEnd=strfind(modelSBML.species(i).annotation,endString);
0281                 formEndIndex=find(formEnd>formStart, 1 );
0282                 formula=modelSBML.species(i).annotation(formStart+numel(startString):formEnd(formEndIndex)-1);
0283                 metaboliteInChI{numel(metaboliteInChI)+1,1}=formula;
0284 
0285                 %The composition is most often present between the
0286                 %first and second "/" in the model. In some simple
0287                 %molecules, such as salts, there is no second "/". The
0288                 %formula is then assumed to be to the end of the string
0289                 compositionIndexes=strfind(formula,'/');
0290                 if numel(compositionIndexes)>1
0291                     metaboliteFormula{numel(metaboliteFormula)+1,1}=...
0292                         formula(compositionIndexes(1)+1:compositionIndexes(2)-1);
0293                 else
0294                     if isscalar(compositionIndexes)
0295                         %Probably a simple molecule which can have only
0296                         %one conformation
0297                         metaboliteFormula{numel(metaboliteFormula)+1,1}=...
0298                             formula(compositionIndexes(1)+1:numel(formula));
0299                     else
0300                         metaboliteFormula{numel(metaboliteFormula)+1,1}='';
0301                     end
0302                 end
0303             elseif isfield(modelSBML.species(i),'fbc_chemicalFormula')
0304                 metaboliteInChI{numel(metaboliteInChI)+1,1}='';
0305                 if ~isempty(modelSBML.species(i).fbc_chemicalFormula)
0306                     %Cannot extract InChi from formula, so remains
0307                     %empty
0308                     metaboliteFormula{numel(metaboliteFormula)+1,1}=modelSBML.species(i).fbc_chemicalFormula;
0309                 else
0310                     metaboliteFormula{numel(metaboliteFormula)+1,1}='';
0311                 end
0312             else
0313                 metaboliteInChI{numel(metaboliteInChI)+1,1}='';
0314                 metaboliteFormula{numel(metaboliteFormula)+1,1}='';
0315             end
0316 
0317             %Get Miriam info
0318             metMiriam=parseMiriam(modelSBML.species(i).annotation);
0319             metaboliteMiriams{numel(metaboliteMiriams)+1,1}=metMiriam;
0320         else
0321             metaboliteInChI{numel(metaboliteInChI)+1,1}='';
0322             if isfield(modelSBML.species(i),'notes')
0323                 metaboliteFormula{numel(metaboliteFormula)+1,1}=parseNote(modelSBML.species(i).notes,'FORMULA');
0324             else
0325                 metaboliteFormula{numel(metaboliteFormula)+1,1}='';
0326             end
0327             metaboliteMiriams{numel(metaboliteMiriams)+1,1}=[];
0328         end
0329         if ~isempty(modelSBML.species(i).notes)
0330             if ~isfield(modelSBML.species(i),'annotation')
0331                 metaboliteFormula{numel(metaboliteFormula)+1,1}=parseNote(modelSBML.species(i).notes,'FORMULA');
0332             end
0333         elseif ~isfield(modelSBML.species(i),'annotation')
0334             metaboliteFormula{numel(metaboliteFormula)+1,1}='';
0335         end
0336         %Get SBO term
0337         if isfield(modelSBML.species(i),'sboTerm') && ~(modelSBML.species(i).sboTerm==-1)
0338             metSBOs(end+1,1) = modelSBML.species(i).sboTerm;
0339         end
0340     end
0341 
0342     %The following lines are executed regardless isSBML2COBRA setting
0343     if isempty(modelSBML.species(i).id) || ~strcmpi(modelSBML.species(i).id(1:2),'E_')
0344         if isempty(modelSBML.species(i).id) || ~strcmpi(modelSBML.species(i).id(1:3),'Cx_')
0345             %Remove trailing [compartment] from metabolite name if present
0346             metaboliteNames{end,1}=regexprep(metaboliteNames{end,1},regexCompNames,'');
0347             metaboliteNames{end,1}=metaboliteNames{end,1};
0348             if isfield(modelSBML.species(i),'fbc_charge')
0349                 if ~isempty(modelSBML.species(i).fbc_charge) && modelSBML.species(i).isSetfbc_charge
0350                     metaboliteCharges(numel(metaboliteCharges)+1,1)=double(modelSBML.species(i).fbc_charge);
0351                 else
0352                     if isfield(modelSBML.species(i),'notes')
0353                         if strfind(modelSBML.species(i).notes,'CHARGE')
0354                             metaboliteCharges(numel(metaboliteCharges)+1,1)=str2double(parseNote(modelSBML.species(i).notes,'CHARGE'));
0355                         else
0356                             metaboliteCharges(numel(metaboliteCharges)+1,1)=NaN;
0357                         end
0358                     else
0359                         metaboliteCharges(numel(metaboliteCharges)+1,1)=NaN;
0360                     end
0361                 end
0362             elseif isfield(modelSBML.species(i),'notes')
0363                 if strfind(modelSBML.species(i).notes,'CHARGE')
0364                     metaboliteCharges(numel(metaboliteCharges)+1,1)=str2double(parseNote(modelSBML.species(i).notes,'CHARGE'));
0365                 else
0366                     metaboliteCharges(numel(metaboliteCharges)+1,1)=NaN;
0367                 end
0368             else
0369                 metaboliteCharges(numel(metaboliteCharges)+1,1)=NaN;
0370             end
0371             %Additional information from FBC format Chemical formula
0372             if isfield(modelSBML.species(i),'fbc_chemicalFormula')
0373                 if ~isempty(modelSBML.species(i).fbc_chemicalFormula)
0374                     metaboliteFormula{numel(metaboliteFormula),1}=modelSBML.species(i).fbc_chemicalFormula;
0375                 end
0376             end
0377         end
0378     end
0379 end
0380 
0381 %Add SBO terms to gene and metabolite miriam fields
0382 if numel(unique(geneSBOs)) > 1  % don't add if they're all identical
0383     for i = 1:numel(geneNames)
0384         geneMiriams{i} = addSBOtoMiriam(geneMiriams{i},geneSBOs(i));
0385     end
0386 end
0387 if numel(unique(metSBOs)) > 1
0388     for i = 1:numel(metaboliteNames)
0389         metaboliteMiriams{i} = addSBOtoMiriam(metaboliteMiriams{i},metSBOs(i));
0390     end
0391 end
0392 
0393 %Retrieve info on reactions
0394 reactionNames=cell(numel(modelSBML.reaction),1);
0395 reactionIDs=cell(numel(modelSBML.reaction),1);
0396 subsystems=cell(numel(modelSBML.reaction),1);
0397 eccodes=cell(numel(modelSBML.reaction),1);
0398 eccodes(:,:)=cellstr('');
0399 rxnconfidencescores=NaN(numel(modelSBML.reaction),1);
0400 rxnreferences=cell(numel(modelSBML.reaction),1);
0401 rxnreferences(:,:)=cellstr('');
0402 rxnnotes=cell(numel(modelSBML.reaction),1);
0403 rxnnotes(:,:)=cellstr('');
0404 grRules=cell(numel(modelSBML.reaction),1);
0405 grRules(:,:)=cellstr('');
0406 grRulesFromModifier=grRules;
0407 rxnComps=zeros(numel(modelSBML.reaction),1);
0408 rxnMiriams=cell(numel(modelSBML.reaction),1);
0409 reactionReversibility=zeros(numel(modelSBML.reaction),1);
0410 reactionUB=zeros(numel(modelSBML.reaction),1);
0411 reactionLB=zeros(numel(modelSBML.reaction),1);
0412 reactionObjective=zeros(numel(modelSBML.reaction),1);
0413 
0414 %Construct the stoichiometric matrix while the reaction info is read
0415 S=zeros(numel(metaboliteIDs),numel(modelSBML.reaction));
0416 
0417 counter=0;
0418 %If FBC, then bounds have parameter ids defined for the whole model
0419 if isfield(modelSBML,'parameter')
0420     parameter.name=cell(numel(modelSBML.parameter),1);
0421     parameter.name={modelSBML.parameter(:).id}';
0422     parameter.value={modelSBML.parameter(:).value}';
0423 end
0424 
0425 if isfield(modelSBML.reaction,'sboTerm') && isscalar(unique([modelSBML.reaction.sboTerm]))
0426     %If all the SBO terms are identical, don't add them to rxnMiriams
0427     modelSBML.reaction = rmfield(modelSBML.reaction,'sboTerm');
0428 end
0429 
0430 for i=1:numel(modelSBML.reaction)
0431 
0432     %Check that the reaction doesn't produce a complex and nothing else. If
0433     %so, then jump to the next reaction. This is because I get the genes
0434     %for complexes from the names and not from the reactions that create
0435     %them. This only applies to the non-COBRA format
0436     if isscalar(modelSBML.reaction(i).product)
0437         if length(modelSBML.reaction(i).product(1).species)>=3
0438             if strcmp(modelSBML.reaction(i).product(1).species(1:3),'Cx_')==true
0439                 continue;
0440             end
0441         end
0442     end
0443 
0444     %It didn't look like a gene complex-forming reaction
0445     counter=counter+1;
0446 
0447     reactionNames{counter}=modelSBML.reaction(i).name;
0448 
0449     reactionIDs{counter}=modelSBML.reaction(i).id;
0450     reactionReversibility(counter)=modelSBML.reaction(i).reversible;
0451 
0452     %If model is FBC, first get parameter of bound and then replace it with
0453     %the correct value. Probably faster with replace(), but this was only
0454     %introduced in Matlab R2016b
0455     if isfield(modelSBML.reaction(i),'fbc_lowerFluxBound')
0456         lb=modelSBML.reaction(i).fbc_lowerFluxBound;
0457         ub=modelSBML.reaction(i).fbc_upperFluxBound;
0458         [~,fluxBoundIdx] = ismember({lb,ub},parameter.name);
0459         reactionLB(counter) = parameter.value{fluxBoundIdx(1)};
0460         reactionUB(counter) = parameter.value{fluxBoundIdx(2)};
0461     elseif isfield(modelSBML.reaction(i).kineticLaw,'parameter')
0462         reactionLB(counter)=modelSBML.reaction(i).kineticLaw.parameter(1).value;
0463         reactionUB(counter)=modelSBML.reaction(i).kineticLaw.parameter(2).value;
0464         reactionObjective(counter)=modelSBML.reaction(i).kineticLaw.parameter(3).value;
0465     else
0466         if reactionReversibility(counter)==true
0467             reactionLB(counter)=-inf;
0468         else
0469             reactionLB(counter)=0;
0470         end
0471         reactionUB(counter)=inf;
0472         reactionObjective(counter)=0;
0473     end
0474 
0475     %Find the associated gene if available
0476     %If FBC, get gene association data from corresponding fields
0477     if isfield(modelSBML.reaction(i),'fbc_geneProductAssociation')
0478         if ~isempty(modelSBML.reaction(i).fbc_geneProductAssociation) && ~isempty(modelSBML.reaction(i).fbc_geneProductAssociation.fbc_association)
0479             grRules{counter}=modelSBML.reaction(i).fbc_geneProductAssociation.fbc_association.fbc_association;
0480         end
0481     elseif isfield(modelSBML.reaction(i),'notes')
0482         %This section was previously executed only if isSBML2COBRA is true. Now
0483         %it will be executed, if 'GENE_ASSOCIATION' is found in
0484         %modelSBML.reaction(i).notes
0485         if strfind(modelSBML.reaction(i).notes,'GENE_ASSOCIATION')
0486             geneAssociation=parseNote(modelSBML.reaction(i).notes,'GENE_ASSOCIATION');
0487         elseif strfind(modelSBML.reaction(i).notes,'GENE ASSOCIATION')
0488             geneAssociation=parseNote(modelSBML.reaction(i).notes,'GENE ASSOCIATION');
0489         else
0490             geneAssociation='';
0491         end
0492         if ~isempty(geneAssociation)
0493             %This adds the grRules. The gene list and rxnGeneMat are created
0494             %later
0495             grRules{counter}=geneAssociation;
0496         end
0497     end
0498     if isempty(grRules{counter}) && ~isempty(modelSBML.reaction(i).modifier)
0499         rules='';
0500         for j=1:numel(modelSBML.reaction(i).modifier)
0501             modifier=modelSBML.reaction(i).modifier(j).species;
0502             if ~isempty(modifier)
0503                 if strcmpi(modifier(1:2),'E_')
0504                     index=find(strcmp(modifier,geneIDs));
0505                     %This should be unique and in the geneIDs list,
0506                     %otherwise something is wrong
0507                     if numel(index)~=1
0508                         EM=['Could not get the gene association data from reaction ' reactionIDs{i}];
0509                         dispEM(EM);
0510                     end
0511                     if ~isempty(rules)
0512                         rules=[rules ' or (' geneNames{index} ')'];
0513                     else
0514                         rules=['(' geneNames{index} ')'];
0515                     end
0516                 elseif strcmp(modifier(1:2),'s_')
0517                     index=find(strcmp(modifier,metaboliteIDs));
0518                     %This should be unique and in the geneIDs list,
0519                     %otherwise something is wrong
0520                     if numel(index)~=1
0521                         EM=['Could not get the gene association data from reaction ' reactionIDs{i}];
0522                         dispEM(EM);
0523                     end
0524                     if ~isempty(rules)
0525                         rules=[rules ' or (' metaboliteIDs{index} ')'];
0526                     else
0527                         rules=['(' metaboliteIDs{index} ')'];
0528                     end
0529                 else
0530                     %It seems to be a complex. Add the corresponding
0531                     %genes from the name of the complex (not the
0532                     %reaction that creates it)
0533                     index=find(strcmp(modifier,complexIDs));
0534                     if isscalar(index)
0535                         if ~isempty(rules)
0536                             rules=[rules ' or (' strrep(complexNames{index},':',' and ') ')'];
0537                         else
0538                             rules=['(' strrep(complexNames{index},':',' and ') ')'];
0539                         end
0540                     else
0541                         %Could not find a complex
0542                         EM=['Could not get the gene association data from reaction ' reactionIDs{i}];
0543                         dispEM(EM);
0544                     end
0545                 end
0546             end
0547         end
0548         grRules{counter}=rules;
0549         grRulesFromModifier{counter}=rules;%Backup copy for grRules, useful to parse Yeast 7.6
0550     end
0551 
0552     %Add reaction compartment
0553     if isfield(modelSBML.reaction(i),'compartment')
0554         if ~isempty(modelSBML.reaction(i).compartment)
0555             rxnComp=modelSBML.reaction(i).compartment;
0556         else
0557             rxnComp='';
0558         end
0559     elseif isfield(modelSBML.reaction(i),'notes')
0560         rxnComp=parseNote(modelSBML.reaction(i).notes,'COMPARTMENT');
0561     end
0562     if ~isempty(rxnComp)
0563         %Find it in the compartment list
0564         [~, J]=ismember(rxnComp,compartmentIDs);
0565         rxnComps(counter)=J;
0566     end
0567 
0568 
0569     miriamStruct=parseMiriam(modelSBML.reaction(i).annotation);
0570     rxnMiriams{counter}=miriamStruct;
0571     if isfield(modelSBML.reaction(i),'notes')
0572         subsystems{counter,1}=cellstr(parseNote(modelSBML.reaction(i).notes,'SUBSYSTEM'));
0573         subsystems{counter,1}(cellfun('isempty',subsystems{counter,1})) = [];
0574         if strfind(modelSBML.reaction(i).notes,'Confidence Level')
0575             confScore = parseNote(modelSBML.reaction(i).notes,'Confidence Level');
0576             if isempty(confScore)
0577                 confScore = 0;
0578             end
0579             rxnconfidencescores(counter)=str2double(confScore);
0580         end
0581         rxnreferences{counter,1}=parseNote(modelSBML.reaction(i).notes,'AUTHORS');
0582         rxnnotes{counter,1}=parseNote(modelSBML.reaction(i).notes,'NOTES');
0583     end
0584 
0585     %Get SBO terms
0586     if isfield(modelSBML.reaction(i),'sboTerm') && ~(modelSBML.reaction(i).sboTerm==-1)
0587         rxnMiriams{counter} = addSBOtoMiriam(rxnMiriams{counter}, modelSBML.reaction(i).sboTerm);
0588     end
0589 
0590     %Get ec-codes
0591     eccode='';
0592     if ~isempty(modelSBML.reaction(i).annotation)
0593         if strfind(modelSBML.reaction(i).annotation,'urn:miriam:ec-code')
0594             eccode=parseAnnotation(modelSBML.reaction(i).annotation,'urn:miriam:',':','ec-code');
0595         elseif strfind(modelSBML.reaction(i).annotation,'http://identifiers.org/ec-code')
0596             eccode=parseAnnotation(modelSBML.reaction(i).annotation,'http://identifiers.org/','/','ec-code');
0597         elseif strfind(modelSBML.reaction(i).annotation,'https://identifiers.org/ec-code')
0598             eccode=parseAnnotation(modelSBML.reaction(i).annotation,'https://identifiers.org/','/','ec-code');
0599         end
0600     elseif isfield(modelSBML.reaction(i),'notes')
0601         if strfind(modelSBML.reaction(i).notes,'EC Number')
0602             eccode=[eccode parseNote(modelSBML.reaction(i).notes,'EC Number')];
0603         elseif strfind(modelSBML.reaction(i).notes,'PROTEIN_CLASS')
0604             eccode=[eccode parseNote(modelSBML.reaction(i).notes,'PROTEIN_CLASS')];
0605         end
0606     end
0607     eccodes{counter}=eccode;
0608 
0609     %Add all reactants
0610     for j=1:numel(modelSBML.reaction(i).reactant)
0611         %Get the index of the metabolite in metaboliteIDs. External
0612         %metabolites will be removed at a later stage
0613         metIndex=find(strcmp(modelSBML.reaction(i).reactant(j).species,metaboliteIDs),1);
0614         if isempty(metIndex)
0615             EM=['Could not find metabolite ' modelSBML.reaction(i).reactant(j).species ' in reaction ' reactionIDs{counter}];
0616             dispEM(EM);
0617         end
0618         S(metIndex,counter)=S(metIndex,counter)+modelSBML.reaction(i).reactant(j).stoichiometry*-1;
0619     end
0620 
0621     %Add all products
0622     for j=1:numel(modelSBML.reaction(i).product)
0623         %Get the index of the metabolite in metaboliteIDs.
0624         metIndex=find(strcmp(modelSBML.reaction(i).product(j).species,metaboliteIDs),1);
0625         if isempty(metIndex)
0626             EM=['Could not find metabolite ' modelSBML.reaction(i).product(j).species ' in reaction ' reactionIDs{counter}];
0627             dispEM(EM);
0628         end
0629         S(metIndex,counter)=S(metIndex,counter)+modelSBML.reaction(i).product(j).stoichiometry;
0630     end
0631 end
0632 
0633 %if FBC, objective function is separately defined. Multiple objective
0634 %functions can be defined, one is set as active
0635 if isfield(modelSBML, 'fbc_activeObjective')
0636     obj=modelSBML.fbc_activeObjective;
0637     for i=1:numel(modelSBML.fbc_objective)
0638         if strcmp(obj,modelSBML.fbc_objective(i).fbc_id)
0639             if ~isempty(modelSBML.fbc_objective(i).fbc_fluxObjective)
0640                 rxn=modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_reaction;
0641                 idx=ismember(reactionIDs,rxn);
0642                 reactionObjective(idx)=modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_coefficient;
0643             end
0644         end
0645     end
0646 end
0647 
0648 %subSystems can be stored as groups instead of in annotations
0649 if isfield(modelSBML,'groups_group')
0650     for i=1:numel(modelSBML.groups_group)
0651         groupreactions={modelSBML.groups_group(i).groups_member(:).groups_idRef};
0652         [~, idx] = ismember(groupreactions, reactionIDs);
0653         if any(idx)
0654             for j=1:numel(idx)
0655                 if isempty(subsystems{idx(j)}) % First subsystem
0656                     subsystems{idx(j)} = {modelSBML.groups_group(i).groups_name};
0657                 else % Consecutive subsystems: concatenate
0658                     subsystems{idx(j)} = horzcat(subsystems{idx(j)}, modelSBML.groups_group(i).groups_name);
0659                 end
0660             end
0661         end
0662     end
0663 end
0664 
0665 %Shrink the structures if complex-forming reactions had to be skipped
0666 reactionNames=reactionNames(1:counter);
0667 reactionIDs=reactionIDs(1:counter);
0668 subsystems=subsystems(1:counter);
0669 eccodes=eccodes(1:counter);
0670 rxnconfidencescores=rxnconfidencescores(1:counter);
0671 rxnreferences=rxnreferences(1:counter);
0672 rxnnotes=rxnnotes(1:counter);
0673 grRules=grRules(1:counter);
0674 rxnMiriams=rxnMiriams(1:counter);
0675 reactionReversibility=reactionReversibility(1:counter);
0676 reactionUB=reactionUB(1:counter);
0677 reactionLB=reactionLB(1:counter);
0678 reactionObjective=reactionObjective(1:counter);
0679 S=S(:,1:counter);
0680 
0681 model.name=modelSBML.name;
0682 model.id=modelSBML.id;
0683 model.rxns=reactionIDs;
0684 model.mets=metaboliteIDs;
0685 model.S=sparse(S);
0686 model.lb=reactionLB;
0687 model.ub=reactionUB;
0688 model.rev=reactionReversibility;
0689 model.c=reactionObjective;
0690 model.b=zeros(numel(metaboliteIDs),1);
0691 model.comps=compartmentIDs;
0692 model.compNames=compartmentNames;
0693 model.rxnConfidenceScores=rxnconfidencescores;
0694 model.rxnReferences=rxnreferences;
0695 model.rxnNotes=rxnnotes;
0696 
0697 %Load annotation if available. If there are several authors, only the first
0698 %author credentials are imported
0699 if isfield(modelSBML,'annotation')
0700     endString='</';
0701     I=strfind(modelSBML.annotation,endString);
0702     J=strfind(modelSBML.annotation,'<vCard:Family>');
0703     if any(J)
0704         model.annotation.familyName=modelSBML.annotation(J(1)+14:I(find(I>J(1),1))-1);
0705     end
0706     J=strfind(modelSBML.annotation,'<vCard:Given>');
0707     if any(J)
0708         model.annotation.givenName=modelSBML.annotation(J(1)+13:I(find(I>J(1),1))-1);
0709     end
0710     J=strfind(modelSBML.annotation,'<vCard:EMAIL>');
0711     if any(J)
0712         model.annotation.email=modelSBML.annotation(J(1)+13:I(find(I>J(1),1))-1);
0713     end
0714     J=strfind(modelSBML.annotation,'<vCard:Orgname>');
0715     if any(J)
0716         model.annotation.organization=modelSBML.annotation(J(1)+15:I(find(I>J(1),1))-1);
0717     end
0718     endString='"/>';
0719     I=strfind(modelSBML.annotation,endString);
0720     if strfind(modelSBML.annotation,'"urn:miriam:')
0721         J=strfind(modelSBML.annotation,'"urn:miriam:');
0722         if any(J)
0723             model.annotation.taxonomy=modelSBML.annotation(J+12:I(find(I>J,1))-1);
0724         end
0725     else
0726         J=strfind(modelSBML.annotation,'"http://identifiers.org/');
0727         if any(J)
0728             I = I(find(I>J,1))-1;
0729             J = J(find(J<I,1))+24;
0730             model.annotation.taxonomy=modelSBML.annotation(J:I);
0731         else
0732             J=strfind(modelSBML.annotation,'"https://identifiers.org/');
0733             if any(J)
0734                 I = I(find(I>J,1))-1;
0735                 J = J(find(J<I,1))+25;
0736                 model.annotation.taxonomy=modelSBML.annotation(J:I);
0737             end
0738         end
0739     end
0740 end
0741 if isfield(modelSBML,'notes')
0742     startString=strfind(modelSBML.notes,'xhtml">');
0743     endString=strfind(modelSBML.notes,'</body>');
0744     if any(startString) && any(endString)
0745         model.annotation.note=modelSBML.notes(startString(1)+7:endString-1);
0746         model.annotation.note=regexprep(model.annotation.note,'<p.*?>|</p.*?>','');
0747         model.annotation.note=strtrim(model.annotation.note);
0748         if regexp(model.annotation.note,'This file was generated using the exportModel function in RAVEN Toolbox \d\.\d and OutputSBML in libSBML')
0749             model.annotation=rmfield(model.annotation,'note'); % Default note added when running exportModel
0750         end
0751     end
0752 end
0753 
0754 if any(~cellfun(@isempty,compartmentOutside))
0755     model.compOutside=compartmentOutside;
0756 end
0757 
0758 model.rxnNames=reactionNames;
0759 model.metNames=metaboliteNames;
0760 
0761 %Match the compartments for metabolites
0762 [~, J]=ismember(metaboliteCompartments,model.comps);
0763 model.metComps=J;
0764 
0765 %If any genes have been loaded (only for the new format)
0766 if ~isempty(geneNames)
0767     %In some rare cases geneNames may not necessarily be used in grRules.
0768     %That is true for Yeast 7.6. It's therefore important to change gene
0769     %systematic names to geneIDs in sophisticated way. Gene systematic
0770     %names are not unique, since exactly the same name may be in different
0771     %compartments
0772     if all(cellfun(@isempty,strfind(grRules,geneNames{1})))
0773         geneShortNames=geneNames;
0774         %geneShortNames contain compartments as well, so these are removed
0775         geneShortNames=regexprep(geneShortNames,' \[.+$','');
0776         %grRules obtained from modifier fields contain geneNames. These are
0777         %changed into geneIDs. grRulesFromModifier is a good way to have
0778         %geneIDs and rxns association when it's important to resolve
0779         %systematic name ambiguities
0780         grRulesFromModifier=regexprep(regexprep(grRulesFromModifier,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs);
0781         grRules=regexprep(regexprep(grRules,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs);
0782 
0783         %Yeast 7.6 contains several metabolites, which were used in gene
0784         %associations. For that reason, the list of species ID is created
0785         %and we then check whether any of them have kegg.genes annotation
0786         %thereby obtaining systematic gene names
0787         geneShortNames=vertcat(geneShortNames,metaboliteNames);
0788         geneIDs=vertcat(geneIDs,metaboliteIDs);
0789         geneSystNames=extractMiriam(vertcat(geneMiriams,metaboliteMiriams),'kegg.genes');
0790         geneCompartments=vertcat(geneCompartments,metaboliteCompartments);
0791         geneMiriams=vertcat(geneMiriams,metaboliteMiriams);
0792 
0793         %Now we retain information for only these entries, which have
0794         %kegg.genes annotation
0795         geneShortNames=geneShortNames(~cellfun('isempty',geneSystNames));
0796         geneIDs=geneIDs(~cellfun('isempty',geneSystNames));
0797         geneSystNames=geneSystNames(~cellfun('isempty',geneSystNames));
0798         geneCompartments=geneCompartments(~cellfun('isempty',geneSystNames));
0799         geneMiriams=geneMiriams(~cellfun('isempty',geneSystNames));
0800         %Now we reorder geneIDs and geneSystNames by geneSystNames string
0801         %length
0802         geneNames=geneIDs;%Backuping geneIDs, since we need unsorted order for later
0803         [~, Indx] = sort(cellfun('size', geneSystNames, 2), 'descend');
0804         geneIDs = geneIDs(Indx);
0805         geneSystNames = geneSystNames(Indx);
0806         for i=1:numel(geneSystNames)
0807             for j=1:numel(grRules)
0808                 if strfind(grRules{j},geneSystNames{i})
0809                     if ~isempty(grRules{j})
0810                         if sum(ismember(geneSystNames,geneSystNames{i}))==1
0811                             grRules{j}=regexprep(grRules{j},geneSystNames{i},geneIDs{i});
0812                         elseif sum(ismember(geneSystNames,geneSystNames{i}))>1
0813                             counter=0;
0814                             ovrlpIDs=geneIDs(ismember(geneSystNames,geneSystNames{i}));
0815                             for k=1:numel(ovrlpIDs)
0816                                 if strfind(grRulesFromModifier{j},ovrlpIDs{k})
0817                                     counter=counter+1;
0818                                     grRules{j}=regexprep(grRules{j},geneSystNames{i},ovrlpIDs{k});
0819                                 end
0820                                 if counter>1
0821                                     EM=['Gene association is ambiguous for reaction ' modelSBML.reaction(j).id];
0822                                     dispEM(EM);
0823                                 end
0824                             end
0825                         end
0826                     end
0827                 end
0828             end
0829         end
0830     end
0831     model.genes=geneNames;
0832     model.grRules=grRules;
0833     [grRules,rxnGeneMat] = standardizeGrRules(model,true);
0834     model.grRules = grRules;
0835     model.rxnGeneMat = rxnGeneMat;
0836 
0837     %Match the compartments for genes
0838     [~, J]=ismember(geneCompartments,model.comps);
0839     model.geneComps=J;
0840 else
0841     if ~all(cellfun(@isempty,grRules))
0842         %If fbc_geneProduct exists, follow the specified gene order, such
0843         %that matching geneShortNames in function below will work
0844         if isfield(modelSBML,'fbc_geneProduct')
0845             genes={modelSBML.fbc_geneProduct.fbc_id};
0846 
0847             %Get gene Miriams if they were not retrieved above (this occurs
0848             %when genes are stored as fbc_geneProduct instead of species)
0849             if isempty(geneMiriams)
0850                 geneMiriams = cell(numel(genes),1);
0851                 if isfield(modelSBML.fbc_geneProduct,'sboTerm') && isscalar(unique([modelSBML.fbc_geneProduct.sboTerm]))
0852                     %If all the SBO terms are identical, don't add them to geneMiriams
0853                     modelSBML.fbc_geneProduct = rmfield(modelSBML.fbc_geneProduct,'sboTerm');
0854                 end
0855                 for i = 1:numel(genes)
0856                     geneMiriams{i}=parseMiriam(modelSBML.fbc_geneProduct(i).annotation);
0857                     if isfield(modelSBML.fbc_geneProduct(i),'sboTerm') && ~(modelSBML.fbc_geneProduct(i).sboTerm==-1)
0858                         geneMiriams{i} = addSBOtoMiriam(geneMiriams{i},modelSBML.fbc_geneProduct(i).sboTerm);
0859                     end
0860                 end
0861             end
0862             proteins={modelSBML.fbc_geneProduct.fbc_name};
0863         else
0864             genes=getGenesFromGrRules(grRules);
0865         end
0866         model.genes=genes;
0867         model.grRules=grRules;
0868         [grRules,rxnGeneMat] = standardizeGrRules(model,true);
0869         model.grRules = grRules;
0870         model.rxnGeneMat = rxnGeneMat;
0871     end
0872 end
0873 
0874 if all(cellfun(@isempty,geneShortNames))
0875     if isfield(modelSBML,'fbc_geneProduct')
0876         for i=1:numel(genes)
0877             if ~isempty(modelSBML.fbc_geneProduct(i).fbc_label)
0878                 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_label;
0879             elseif ~isempty(modelSBML.fbc_geneProduct(i).fbc_name)
0880                 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_name;
0881             else
0882                 geneShortNames{i,1}='';
0883             end
0884         end
0885     end
0886 end
0887 
0888 %If any InChIs have been loaded
0889 if any(~cellfun(@isempty,metaboliteInChI))
0890     model.inchis=metaboliteInChI;
0891 end
0892 
0893 %If any formulas have been loaded
0894 if any(~cellfun(@isempty,metaboliteFormula))
0895     model.metFormulas=metaboliteFormula;
0896 end
0897 
0898 %If any charges have been loaded
0899 if ~isempty(metaboliteCharges)
0900     model.metCharges=metaboliteCharges;
0901 end
0902 
0903 %If any gene short names have been loaded
0904 if any(~cellfun(@isempty,geneShortNames))
0905     model.geneShortNames=geneShortNames;
0906 end
0907 
0908 %If any Miriam strings for compartments have been loaded
0909 if any(~cellfun(@isempty,compartmentMiriams))
0910     model.compMiriams=compartmentMiriams;
0911 end
0912 
0913 %If any Miriam strings for metabolites have been loaded
0914 if any(~cellfun(@isempty,metaboliteMiriams))
0915     model.metMiriams=metaboliteMiriams;
0916 end
0917 
0918 %If any subsystems have been loaded
0919 if any(~cellfun(@isempty,subsystems))
0920     model.subSystems=subsystems;
0921 end
0922 if any(rxnComps)
0923     if all(rxnComps)
0924         model.rxnComps=rxnComps;
0925     else
0926         if supressWarnings==false
0927             EM='The compartments for the following reactions could not be matched. Ignoring reaction compartment information';
0928             dispEM(EM,false,model.rxns(rxnComps==0));
0929         end
0930     end
0931 end
0932 
0933 %If any ec-codes have been loaded
0934 if any(~cellfun(@isempty,eccodes))
0935     model.eccodes=eccodes;
0936 end
0937 
0938 %If any Miriam strings for reactions have been loaded
0939 if any(~cellfun(@isempty,rxnMiriams))
0940     model.rxnMiriams=rxnMiriams;
0941 end
0942 
0943 %If any Miriam strings for genes have been loaded
0944 if any(~cellfun(@isempty,geneMiriams))
0945     model.geneMiriams=geneMiriams;
0946 end
0947 
0948 %If any protein strings have been loaded
0949 if any(~cellfun(@isempty,proteins))
0950     proteins = reshape(proteins,[],1);
0951     model.proteins=proteins;
0952 end
0953 
0954 model.unconstrained=metaboliteUnconstrained;
0955 
0956 %Convert SBML IDs back into their original strings. Here we are using part
0957 %from convertSBMLID, originating from the COBRA Toolbox
0958 model.rxns=regexprep(model.rxns,'__([0-9]+)__','${char(str2num($1))}');
0959 model.mets=regexprep(model.mets,'__([0-9]+)__','${char(str2num($1))}');
0960 model.comps=regexprep(model.comps,'__([0-9]+)__','${char(str2num($1))}');
0961 model.grRules=regexprep(model.grRules,'__([0-9]+)__','${char(str2num($1))}');
0962 model.genes=regexprep(model.genes,'__([0-9]+)__','${char(str2num($1))}');
0963 model.id=regexprep(model.id,'__([0-9]+)__','${char(str2num($1))}');
0964 
0965 if removePrefix
0966     [model, hasChanged]=removeIdentifierPrefix(model);
0967     dispEM(['The following fields have prefixes removed from all entries. '...
0968     'If this is undesired, run importModel with removePrefix as false. Example: '...
0969     'importModel(''filename.xml'',[],false);'],false,hasChanged)
0970 end
0971 
0972 %Remove unused fields
0973 if isempty(model.annotation)
0974     model=rmfield(model,'annotation');
0975 end
0976 if isempty(model.compOutside)
0977     model=rmfield(model,'compOutside');
0978 end
0979 if isempty(model.compMiriams)
0980     model=rmfield(model,'compMiriams');
0981 end
0982 if isempty(model.rxnComps)
0983     model=rmfield(model,'rxnComps');
0984 end
0985 if isempty(model.grRules)
0986     model=rmfield(model,'grRules');
0987 end
0988 if isempty(model.rxnGeneMat)
0989     model=rmfield(model,'rxnGeneMat');
0990 end
0991 if isempty(model.subSystems)
0992     model=rmfield(model,'subSystems');
0993 else
0994     model.subSystems(cellfun(@isempty,subsystems))={{''}};
0995 end
0996 if isempty(model.eccodes)
0997     model=rmfield(model,'eccodes');
0998 end
0999 if isempty(model.rxnMiriams)
1000     model=rmfield(model,'rxnMiriams');
1001 end
1002 if cellfun(@isempty,model.rxnNotes)
1003     model=rmfield(model,'rxnNotes');
1004 end
1005 if cellfun(@isempty,model.rxnReferences)
1006     model=rmfield(model,'rxnReferences');
1007 end
1008 if isempty(model.rxnConfidenceScores) || all(isnan(model.rxnConfidenceScores))
1009     model=rmfield(model,'rxnConfidenceScores');
1010 end
1011 if isempty(model.genes)
1012     model=rmfield(model,'genes');
1013 elseif isrow(model.genes)
1014     model.genes=transpose(model.genes);
1015 end
1016 if isempty(model.geneComps)
1017     model=rmfield(model,'geneComps');
1018 end
1019 if isempty(model.geneMiriams)
1020     model=rmfield(model,'geneMiriams');
1021 end
1022 if isempty(model.geneShortNames)
1023     model=rmfield(model,'geneShortNames');
1024 end
1025 if isempty(model.proteins)
1026     model=rmfield(model,'proteins');
1027 end
1028 if isempty(model.inchis)
1029     model=rmfield(model,'inchis');
1030 end
1031 if isempty(model.metFormulas)
1032     model=rmfield(model,'metFormulas');
1033 end
1034 if isempty(model.metMiriams)
1035     model=rmfield(model,'metMiriams');
1036 end
1037 if ~any(model.metCharges)
1038     model=rmfield(model,'metCharges');
1039 end
1040 
1041 %This just removes the grRules if no genes have been loaded
1042 if ~isfield(model,'genes') && isfield(model,'grRules')
1043     model=rmfield(model,'grRules');
1044 end
1045 
1046 %Print warnings about bad structure
1047 if supressWarnings==false
1048     checkModelStruct(model,false);
1049 end
1050 
1051 if removeExcMets==true
1052     model=simplifyModel(model);
1053 end
1054 end
1055 
1056 function fieldContent=parseNote(searchString,fieldName)
1057 %The function obtains the particular information from 'notes' field, using
1058 %fieldName as the dummy string
1059 
1060 fieldContent='';
1061 
1062 if strfind(searchString,fieldName)
1063     [~,targetString] = regexp(searchString,['<p>' fieldName '.*?</p>'],'tokens','match');
1064     targetString=regexprep(targetString,'<p>|</p>','');
1065     targetString=regexprep(targetString,[fieldName, ':'],'');
1066     for i=1:numel(targetString)
1067         fieldContent=[fieldContent ';' strtrim(targetString{1,i})];
1068     end
1069     fieldContent=regexprep(fieldContent,'^;|;$','');
1070 else
1071     fieldContent='';
1072 end
1073 end
1074 
1075 function fieldContent=parseAnnotation(searchString,startString,midString,fieldName)
1076 
1077 fieldContent='';
1078 
1079 %Removing whitespace characters from the ending strings, which may occur in
1080 %several cases
1081 searchString=regexprep(searchString,'" />','"/>');
1082 [~,targetString] = regexp(searchString,['<rdf:li rdf:resource="' startString fieldName midString '.*?"/>'],'tokens','match');
1083 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>','');
1084 targetString=regexprep(targetString,startString,'');
1085 targetString=regexprep(targetString,[fieldName midString],'');
1086 
1087 for i=1:numel(targetString)
1088     fieldContent=[fieldContent ';' strtrim(targetString{1,i})];
1089 end
1090 
1091 fieldContent=regexprep(fieldContent,'^;|;$','');
1092 end
1093 
1094 function miriamStruct=parseMiriam(searchString)
1095 %Generates miriam structure from annotation field
1096 
1097 %Finding whether miriams are written in the old or the new way
1098 if strfind(searchString,'urn:miriam:')
1099     startString='urn:miriam:';
1100 elseif strfind(searchString,'http://identifiers.org/')
1101     startString='http://identifiers.org/';
1102 elseif strfind(searchString,'https://identifiers.org/')
1103     startString='https://identifiers.org/';
1104 else
1105     miriamStruct=[];
1106     return;
1107 end
1108 
1109 miriamStruct=[];
1110 
1111 searchString=regexprep(searchString,'" />','"/>');
1112 [~,targetString] = regexp(searchString,'<rdf:li rdf:resource=".*?"/>','tokens','match');
1113 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>','');
1114 targetString=regexprep(targetString,startString,'');
1115 
1116 fwdslash  = contains(targetString,'/');
1117 midString = cell(numel(targetString),1);
1118 midString(fwdslash) = {'/'};
1119 midString(~fwdslash) = {':'};
1120 
1121 counter=0;
1122 for i=1:numel(targetString)
1123     if isempty(regexp(targetString{1,i},'inchi|ec-code', 'once'))
1124         counter=counter+1;
1125         miriamStruct.name{counter,1} = regexprep(targetString{1,i},[midString{i} '.+'],'','once');
1126         miriamStruct.value{counter,1} = regexprep(targetString{1,i},[miriamStruct.name{counter,1} midString{i}],'','once');
1127         miriamStruct.name{counter,1} = regexprep(miriamStruct.name{counter,1},'^obo\.','');
1128     end
1129 end
1130 end
1131 
1132 function miriam = addSBOtoMiriam(miriam,sboTerm)
1133 %Appends SBO term to miriam structure
1134 
1135 sboTerm = {['SBO:' sprintf('%07u',sboTerm)]};  % convert to proper format
1136 if isempty(miriam)
1137     miriam.name = {'sbo'};
1138     miriam.value = sboTerm;
1139 elseif any(strcmp('sbo',miriam.name))
1140     currSbo = strcmp('sbo',miriam.name);
1141     miriam.value(currSbo) = sboTerm;
1142 else
1143     miriam.name(end+1) = {'sbo'};
1144     miriam.value(end+1) = sboTerm;
1145 end
1146 end