exportModel

 exportModel
   Exports a constraint-based model to an SBML file (L3V1 FBCv2)

 Input:
   model               a model structure
   fileName            filename to export the model to. A dialog window
                       will open if no file name is specified.
   neverPrefix         true if prefixes are never added to identifiers,
                       even if start with e.g. digits. This might result
                       in invalid SBML files (optional, default false)
   supressWarnings     true if warnings should be supressed. This might
                       results in invalid SBML files, as no checks are
                       performed (optional, default false)
   sortIds             logical whether metabolites, reactions and genes
                       should be sorted alphabetically by their
                       identifiers (optional, default false)

 Usage: exportModel(model,fileName,neverPrefix,supressWarnings,sortIds)

0001 function exportModel(model,fileName,neverPrefix,supressWarnings,sortIds)
0002 % exportModel
0003 %   Exports a constraint-based model to an SBML file (L3V1 FBCv2)
0004 %
0005 % Input:
0006 %   model               a model structure
0007 %   fileName            filename to export the model to. A dialog window
0008 %                       will open if no file name is specified.
0009 %   neverPrefix         true if prefixes are never added to identifiers,
0010 %                       even if start with e.g. digits. This might result
0011 %                       in invalid SBML files (optional, default false)
0012 %   supressWarnings     true if warnings should be supressed. This might
0013 %                       results in invalid SBML files, as no checks are
0014 %                       performed (optional, default false)
0015 %   sortIds             logical whether metabolites, reactions and genes
0016 %                       should be sorted alphabetically by their
0017 %                       identifiers (optional, default false)
0018 %
0019 % Usage: exportModel(model,fileName,neverPrefix,supressWarnings,sortIds)
0020 
0021 if nargin<2 || isempty(fileName)
0022     [fileName, pathName] = uiputfile({'*.xml;*.sbml'}, 'Select file for model export',[model.id '.xml']);
0023     if fileName == 0
0024         error('You should provide a file location')
0025     else
0026         fileName = fullfile(pathName,fileName);
0027     end
0028 end
0029 fileName=char(fileName);
0030 if nargin<3 || isempty(neverPrefix)
0031     neverPrefix=false;
0032 end
0033 if nargin<4 || isempty(supressWarnings)
0034     supressWarnings=false;
0035 end
0036 if nargin<5 || isempty(sortIds)
0037     sortIds=false;
0038 end
0039 if sortIds==true
0040     model=sortIdentifiers(model);
0041 end
0042 
0043 if isfield(model,'ec')
0044     warning("exportModel does not store information from the 'model.ec' structure. Use 'writeYAMLmodel(model)' to export all content from a GECKO model.")
0045 end
0046 
0047 %If no subSystems are defined, then no need to use groups package
0048 if isfield(model,'subSystems')
0049     modelHasSubsystems=true;
0050 else
0051     modelHasSubsystems=false;
0052 end
0053 
0054 %The default SBML format settings, which are used as input for appropriate
0055 %libSBML functions to generate the blank SBML model structure before using
0056 %exporting in with OutputSBML to xml file
0057 sbmlLevel=3;
0058 sbmlVersion=1;
0059 sbmlPackages={'fbc'};
0060 sbmlPackageVersions=2;
0061 if modelHasSubsystems
0062     sbmlPackages={sbmlPackages,'groups'};
0063     sbmlPackageVersions=[sbmlPackageVersions,1];
0064 end
0065 
0066 %Check if the "unconstrained" field is still present. This shows if
0067 %exchange metabolites have been removed
0068 if ~isfield(model,'unconstrained')
0069     model.unconstrained=zeros(numel(model.mets),1);
0070 end
0071 
0072 %If model id and name do not exist, make sure that default
0073 %strings are included
0074 if ~isfield(model,'id')
0075     fprintf('WARNING: The model is missing the "id" field. Uses "blankID". \n');
0076     model.id='blankID';
0077 end
0078 if ~isfield(model,'name')
0079     fprintf('WARNING: The model is missing the "name" field. Uses "blankName". \n');
0080     model.name='blankName';
0081 end
0082 
0083 % Add prefixes if required
0084 if ~neverPrefix
0085     [model,hasChanged] = addIdentifierPrefix(model);
0086     dispEM(['The following fields have one or more entries that do not start '...
0087         'with a letter or _ (conflicting with SBML specifications). Prefixes '...
0088         'are added to all entries in those fields:'],false,hasChanged)
0089 end
0090 
0091 %Check the model structure
0092 if supressWarnings==false
0093     checkModelStruct(model);
0094 end
0095 
0096 %Add several blank fields, if they do not exist already. This is to reduce
0097 %the number of conditions below
0098 if ~isfield(model,'compMiriams')
0099     model.compMiriams=cell(numel(model.comps),1);
0100 end
0101 if ~isfield(model,'inchis')
0102     model.inchis=cell(numel(model.mets),1);
0103 end
0104 if ~isfield(model,'metFormulas')
0105     model.metFormulas=cell(numel(model.mets),1);
0106 end
0107 if ~isfield(model,'metMiriams')
0108     model.metMiriams=cell(numel(model.mets),1);
0109 end
0110 if ~isfield(model,'geneMiriams') && isfield(model,'genes')
0111     model.geneMiriams=cell(numel(model.genes),1);
0112 end
0113 if ~isfield(model,'geneShortNames') && isfield(model,'genes')
0114     model.geneShortNames=cell(numel(model.genes),1);
0115 end
0116 if ~isfield(model,'proteins') && isfield(model,'genes')
0117     model.proteins=cell(numel(model.genes),1);
0118 end
0119 if ~isfield(model,'subSystems')
0120     model.subSystems=cell(numel(model.rxns),1);
0121 end
0122 if ~isfield(model,'eccodes')
0123     model.eccodes=cell(numel(model.rxns),1);
0124 end
0125 if ~isfield(model,'rxnReferences')
0126     model.rxnReferences=cell(numel(model.rxns),1);
0127 end
0128 if ~isfield(model,'rxnConfidenceScores')
0129     model.rxnConfidenceScores=NaN(numel(model.rxns),1);
0130 end
0131 if ~isfield(model,'rxnNotes')
0132     model.rxnNotes=cell(numel(model.rxns),1);
0133 end
0134 if ~isfield(model,'rxnMiriams')
0135     model.rxnMiriams=cell(numel(model.rxns),1);
0136 end
0137 
0138 if sbmlLevel<3
0139     %Check if genes have associated compartments
0140     if ~isfield(model,'geneComps') && isfield(model,'genes')
0141         if supressWarnings==false
0142             EM='There are no compartments specified for genes. All genes will be assigned to the first compartment. This is because the SBML structure requires all elements to be assigned to a compartment';
0143             dispEM(EM,false);
0144         end
0145         model.geneComps=ones(numel(model.genes),1);
0146     end
0147 end
0148 
0149 %Convert ids to SBML-convenient format. This is to avoid the data loss when
0150 %unsupported characters are included in ids. Here we are using part from
0151 %convertSBMLID, originating from the COBRA Toolbox
0152 model.rxns=regexprep(model.rxns,'([^0-9_a-zA-Z])','__${num2str($1+0)}__');
0153 model.mets=regexprep(model.mets,'([^0-9_a-zA-Z])','__${num2str($1+0)}__');
0154 model.comps=regexprep(model.comps,'([^0-9_a-zA-Z])','__${num2str($1+0)}__');
0155 if isfield(model,'genes')
0156     problemGenes=find(~cellfun('isempty',regexp(model.genes,'([^0-9_a-zA-Z])')));
0157     originalGenes=model.genes(problemGenes);
0158     replacedGenes=regexprep(model.genes(problemGenes),'([^0-9_a-zA-Z])','__${num2str($1+0)}__');
0159     model.genes(problemGenes)=replacedGenes;
0160     for i=1:numel(problemGenes)
0161         model.grRules = regexprep(model.grRules, ['(^|\s|\()' originalGenes{i} '($|\s|\))'], ['$1' replacedGenes{i} '$2']);
0162     end
0163 end
0164 
0165 %Generate an empty SBML structure
0166 modelSBML=getSBMLStructure(sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0167 modelSBML.metaid=model.id;
0168 modelSBML.id=regexprep(model.id,'([^0-9_a-zA-Z])','__${num2str($1+0)}__');
0169 modelSBML.name=model.name;
0170 
0171 if isfield(model,'annotation')
0172     if isfield(model.annotation,'note')
0173         modelSBML.notes=['<notes><body xmlns="http://www.w3.org/1999/xhtml"><p>',regexprep(model.annotation.note,'<p>|</p>',''),'</p></body></notes>'];
0174     end
0175 else
0176     modelSBML.notes='<notes><body xmlns="http://www.w3.org/1999/xhtml"><p>This file was generated using the exportModel function in RAVEN Toolbox 2 and OutputSBML in libSBML </p></body></notes>';
0177 end
0178 
0179 if isfield(model,'annotation')
0180     nameString='';
0181     if isfield(model.annotation,'familyName')
0182         if ~isempty(model.annotation.familyName)
0183             nameString=['<vCard:Family>' model.annotation.familyName '</vCard:Family>'];
0184         end
0185     end
0186     if isfield(model.annotation,'givenName')
0187         if ~isempty(model.annotation.givenName)
0188             nameString=[nameString '<vCard:Given>' model.annotation.givenName '</vCard:Given>'];
0189         end
0190     end
0191     email='';
0192     if isfield(model.annotation,'email')
0193         if ~isempty(model.annotation.email)
0194             email=['<vCard:EMAIL>' model.annotation.email '</vCard:EMAIL>'];
0195         end
0196     end
0197     org='';
0198     if isfield(model.annotation,'organization')
0199         if ~isempty(model.annotation.organization)
0200             org=['<vCard:ORG rdf:parseType="Resource"><vCard:Orgname>' model.annotation.organization '</vCard:Orgname></vCard:ORG>'];
0201         end
0202     end
0203     if ~isempty(nameString) || ~isempty(email) || ~isempty(org) % Only fill .annotation if ownership data is provided
0204         modelSBML.annotation=['<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#meta_' model.id '">'];
0205         modelSBML.annotation=[modelSBML.annotation '<dc:creator><rdf:Bag><rdf:li rdf:parseType="Resource">'];
0206         if ~isempty(nameString)
0207             modelSBML.annotation=[modelSBML.annotation '<vCard:N rdf:parseType="Resource">' nameString '</vCard:N>'];
0208         end
0209         modelSBML.annotation=[modelSBML.annotation email org '</rdf:li></rdf:Bag></dc:creator>'];
0210         modelSBML.annotation=[modelSBML.annotation '<dcterms:created rdf:parseType="Resource">'...
0211             '<dcterms:W3CDTF>' datestr(now,'yyyy-mm-ddTHH:MM:SSZ') '</dcterms:W3CDTF></dcterms:created><dcterms:modified rdf:parseType="Resource">'...
0212             '<dcterms:W3CDTF>' datestr(now,'yyyy-mm-ddTHH:MM:SSZ') '</dcterms:W3CDTF></dcterms:modified>'];
0213         if isfield(model.annotation,'taxonomy')
0214             modelSBML.annotation=[modelSBML.annotation '<bqbiol:is><rdf:Bag><rdf:li rdf:resource="https://identifiers.org/taxonomy/' regexprep(model.annotation.taxonomy,'taxonomy/','') '"/></rdf:Bag></bqbiol:is>'];
0215         end
0216         modelSBML.annotation=[modelSBML.annotation '</rdf:Description></rdf:RDF></annotation>'];
0217     end
0218 end
0219 
0220 %Prepare compartments
0221 for i=1:numel(model.comps)
0222     %Add the default values, as these will be the same in all entries
0223     if i==1
0224         if isfield(modelSBML.compartment, 'sboTerm')
0225             modelSBML.compartment(i).sboTerm=290;
0226         end
0227         if isfield(modelSBML.compartment, 'spatialDimensions')
0228             modelSBML.compartment(i).spatialDimensions=3;
0229         end
0230         if isfield(modelSBML.compartment, 'size')
0231             modelSBML.compartment(i).size=1;
0232         end
0233         if isfield(modelSBML.compartment, 'constant')
0234             modelSBML.compartment(i).constant=1;
0235         end
0236         if isfield(modelSBML.compartment, 'isSetSize')
0237             modelSBML.compartment(i).isSetSize=1;
0238         end
0239         if isfield(modelSBML.compartment, 'isSetSpatialDimensions')
0240             modelSBML.compartment(i).isSetSpatialDimensions=1;
0241         end
0242     end
0243     %Copy the default values to the next entry as long as it is not the
0244     %last one
0245     if i<numel(model.comps)
0246         modelSBML.compartment(i+1)=modelSBML.compartment(i);
0247     end
0248 
0249     if isfield(modelSBML.compartment,'metaid')
0250         modelSBML.compartment(i).metaid=model.comps{i};
0251     end
0252     %Prepare Miriam strings
0253     if ~isempty(model.compMiriams{i})
0254         [~,sbo_ind] = ismember('sbo',model.compMiriams{i}.name);
0255         if sbo_ind > 0
0256             modelSBML.compartment(i).sboTerm=str2double(regexprep(model.compMiriams{i}.value{sbo_ind},'SBO:','','ignorecase'));
0257             % remove the SBO term from compMiriams so the information is
0258             % not duplicated in the "annotation" field later on
0259             model.compMiriams{i}.name(sbo_ind) = [];
0260             model.compMiriams{i}.value(sbo_ind) = [];
0261         end
0262     end
0263     if ~isempty(model.compMiriams{i}) && isfield(modelSBML.compartment(i),'annotation')
0264         modelSBML.compartment(i).annotation=['<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#meta_' model.comps{i} '">'];
0265         modelSBML.compartment(i).annotation=[modelSBML.compartment(i).annotation '<bqbiol:is><rdf:Bag>'];
0266         modelSBML.compartment(i).annotation=[modelSBML.compartment(i).annotation getMiriam(model.compMiriams{i}) '</rdf:Bag></bqbiol:is></rdf:Description></rdf:RDF></annotation>'];
0267     end
0268     if isfield(modelSBML.compartment, 'name')
0269         modelSBML.compartment(i).name=model.compNames{i};
0270     end
0271     if isfield(modelSBML.compartment, 'id')
0272         modelSBML.compartment(i).id=model.comps{i};
0273     end
0274 
0275 end
0276 
0277 %Begin writing species
0278 for i=1:numel(model.mets)
0279     %Add the default values, as these will be the same in all entries
0280     if i==1
0281         if isfield(modelSBML.species, 'sboTerm')
0282             modelSBML.species(i).sboTerm=247;
0283         end
0284         if isfield(modelSBML.species, 'initialAmount')
0285             modelSBML.species(i).initialAmount=1;
0286         end
0287         if isfield(modelSBML.species, 'initialConcentration')
0288             modelSBML.species(i).initialConcentration=0;
0289         end
0290         if isfield(modelSBML.species, 'isSetInitialAmount')
0291             modelSBML.species(i).isSetInitialAmount=1;
0292         end
0293         if isfield(modelSBML.species, 'isSetInitialConcentration')
0294             modelSBML.species(i).isSetInitialConcentration=1;
0295         end
0296     end
0297     %Copy the default values to the next entry as long as it is not the
0298     %last one
0299     if i<numel(model.mets)
0300         modelSBML.species(i+1)=modelSBML.species(i);
0301     end
0302 
0303     if isfield(modelSBML.species,'metaid')
0304         modelSBML.species(i).metaid=model.mets{i};
0305     end
0306     if isfield(modelSBML.species, 'name')
0307         modelSBML.species(i).name=model.metNames{i};
0308     end
0309     if isfield(modelSBML.species, 'id')
0310         modelSBML.species(i).id=model.mets{i};
0311     end
0312     if isfield(modelSBML.species, 'compartment')
0313         modelSBML.species(i).compartment=model.comps{model.metComps(i)};
0314     end
0315     if isfield(model,'unconstrained')
0316         if model.unconstrained(i)
0317             modelSBML.species(i).boundaryCondition=1;
0318         end
0319     end
0320     if isfield(modelSBML.species, 'fbc_charge') && isfield(model,'metCharges')
0321         if ~isnan(model.metCharges(i))
0322             modelSBML.species(i).fbc_charge=model.metCharges(i);
0323             modelSBML.species(i).isSetfbc_charge=1;
0324         else
0325             modelSBML.species(i).isSetfbc_charge=0;
0326         end
0327     end
0328     if ~isempty(model.metMiriams{i})
0329         [~,sbo_ind] = ismember('sbo',model.metMiriams{i}.name);
0330         if sbo_ind > 0
0331             modelSBML.species(i).sboTerm=str2double(regexprep(model.metMiriams{i}.value{sbo_ind},'SBO:','','ignorecase'));
0332             % remove the SBO term from metMiriams so the information is
0333             % not duplicated in the "annotation" field later on
0334             model.metMiriams{i}.name(sbo_ind) = [];
0335             model.metMiriams{i}.value(sbo_ind) = [];
0336         end
0337     end
0338     if isfield(modelSBML.species,'annotation')
0339         if ~isempty(model.metMiriams{i}) || ~isempty(model.metFormulas{i})
0340             hasInchi=false;
0341             if ~isempty(model.metFormulas{i})
0342                 %Only export formula if there is no InChI. This is because
0343                 %the metFormulas field is populated by InChIs if available
0344                 if ~isempty(model.inchis{i})
0345                     hasInchi=true;
0346                 end
0347                 if hasInchi==false
0348                     modelSBML.species(i).fbc_chemicalFormula=model.metFormulas{i};
0349                 end
0350             end
0351             if ~isempty(model.metMiriams{i}) || hasInchi==true
0352                 modelSBML.species(i).annotation=['<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#meta_' model.mets{i} '">'];
0353                 modelSBML.species(i).annotation=[modelSBML.species(i).annotation '<bqbiol:is><rdf:Bag>'];
0354                 if ~isempty(model.metMiriams{i})
0355                     modelSBML.species(i).annotation=[modelSBML.species(i).annotation getMiriam(model.metMiriams{i})];
0356                 end
0357                 if hasInchi==true
0358                     modelSBML.species(i).annotation=[modelSBML.species(i).annotation '<rdf:li rdf:resource="https://identifiers.org/inchi/InChI=' regexprep(model.inchis{i},'^InChI=','') '"/>'];
0359                     modelSBML.species(i).fbc_chemicalFormula=char(regexp(model.inchis{i}, '/(\w+)/', 'tokens', 'once'));
0360                 end
0361                 modelSBML.species(i).annotation=[modelSBML.species(i).annotation '</rdf:Bag></bqbiol:is></rdf:Description></rdf:RDF></annotation>'];
0362             end
0363         end
0364     end
0365 end
0366 
0367 if isfield(model,'genes')
0368     for i=1:numel(model.genes)
0369         %Add the default values, as these will be the same in all entries
0370         if i==1
0371             if isfield(modelSBML.fbc_geneProduct, 'sboTerm')
0372                 modelSBML.fbc_geneProduct(i).sboTerm=243;
0373             end
0374         end
0375         %Copy the default values to the next index as long as it is not the
0376         %last one
0377         if i<numel(model.genes)
0378             modelSBML.fbc_geneProduct(i+1)=modelSBML.fbc_geneProduct(i);
0379         end
0380 
0381         if isfield(modelSBML.fbc_geneProduct,'metaid')
0382             modelSBML.fbc_geneProduct(i).metaid=model.genes{i};
0383         end
0384         if ~isempty(model.geneMiriams{i})
0385             [~,sbo_ind] = ismember('sbo',model.geneMiriams{i}.name);
0386             if sbo_ind > 0
0387                 modelSBML.fbc_geneProduct(i).sboTerm=str2double(regexprep(model.geneMiriams{i}.value{sbo_ind},'SBO:','','ignorecase'));
0388                 % remove the SBO term from compMiriams so the information is
0389                 % not duplicated in the "annotation" field later on
0390                 model.geneMiriams{i}.name(sbo_ind) = [];
0391                 model.geneMiriams{i}.value(sbo_ind) = [];
0392             end
0393         end
0394         if ~isempty(model.geneMiriams{i}) && isfield(modelSBML.fbc_geneProduct(i),'annotation')
0395             modelSBML.fbc_geneProduct(i).annotation=['<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#meta_' model.genes{i} '">'];
0396             modelSBML.fbc_geneProduct(i).annotation=[modelSBML.fbc_geneProduct(i).annotation '<bqbiol:is><rdf:Bag>'];
0397             modelSBML.fbc_geneProduct(i).annotation=[modelSBML.fbc_geneProduct(i).annotation getMiriam(model.geneMiriams{i}) '</rdf:Bag></bqbiol:is></rdf:Description></rdf:RDF></annotation>'];
0398         end
0399         if isfield(modelSBML.fbc_geneProduct, 'fbc_id')
0400             modelSBML.fbc_geneProduct(i).fbc_id=model.genes{i};
0401         end
0402         if isfield(modelSBML.fbc_geneProduct, 'fbc_label') && isfield(model,'geneShortNames')
0403             if isempty(model.geneShortNames{i})
0404                 modelSBML.fbc_geneProduct(i).fbc_label=model.genes{i};
0405             else
0406                 modelSBML.fbc_geneProduct(i).fbc_label=model.geneShortNames{i};
0407             end
0408         end
0409         if isfield(modelSBML.fbc_geneProduct, 'fbc_name') && isfield(model,'proteins')
0410             if ~isempty(model.proteins{i})
0411                 modelSBML.fbc_geneProduct(i).fbc_name=model.proteins{i};
0412             end
0413         end
0414     end
0415 end
0416 
0417 %Generate a list of unique fbc_bound names
0418 totalValues=[model.lb; model.ub];
0419 totalNames=cell(size(totalValues,1),1);
0420 
0421 listUniqueValues=unique(totalValues);
0422 
0423 for i=1:length(listUniqueValues)
0424     listUniqueNames{i,1}=['FB',num2str(i),'N',num2str(abs(round(listUniqueValues(i))))]; % create unique flux bound IDs.
0425     ind=find(ismember(totalValues,listUniqueValues(i)));
0426     totalNames(ind)=listUniqueNames(i,1);
0427 end
0428 
0429 for i=1:length(listUniqueNames)
0430     %Add the default values, as these will be the same in all entries
0431     if i==1
0432         if isfield(modelSBML.parameter, 'constant')
0433             modelSBML.parameter(i).constant=1;
0434         end
0435         if isfield(modelSBML.parameter, 'isSetValue')
0436             modelSBML.parameter(i).isSetValue=1;
0437         end
0438     end
0439     %Copy the default values to the next index as long as it is not the
0440     %last one
0441     if i<numel(listUniqueNames)
0442         modelSBML.parameter(i+1)=modelSBML.parameter(i);
0443     end
0444     modelSBML.parameter(i).id=listUniqueNames{i};
0445     modelSBML.parameter(i).value=listUniqueValues(i);
0446 end
0447 
0448 for i=1:numel(model.rxns)
0449     %Add the default values, as these will be the same in all entries
0450     if i==1
0451         if isfield(modelSBML.reaction, 'sboTerm')
0452             modelSBML.reaction(i).sboTerm=176;
0453         end
0454         if isfield(modelSBML.reaction, 'isSetFast')
0455             modelSBML.reaction(i).isSetFast=1;
0456         end
0457     end
0458     %Copy the default values to the next index as long as it is not the
0459     %last one
0460     if i<numel(model.rxns)
0461         modelSBML.reaction(i+1)=modelSBML.reaction(i);
0462     end
0463 
0464     if isfield(modelSBML.reaction,'metaid')
0465         modelSBML.reaction(i).metaid=model.rxns{i};
0466     end
0467 
0468     %Export notes information
0469     if (~isnan(model.rxnConfidenceScores(i)) || ~isempty(model.rxnReferences{i}) || ~isempty(model.rxnNotes{i}))
0470         modelSBML.reaction(i).notes='<notes><body xmlns="http://www.w3.org/1999/xhtml">';
0471         if ~isnan(model.rxnConfidenceScores(i))
0472             modelSBML.reaction(i).notes=[modelSBML.reaction(i).notes '<p>Confidence Level: ' num2str(model.rxnConfidenceScores(i)) '</p>'];
0473         end
0474         if ~isempty(model.rxnReferences{i})
0475             modelSBML.reaction(i).notes=[modelSBML.reaction(i).notes '<p>AUTHORS: ' model.rxnReferences{i} '</p>'];
0476         end
0477         if ~isempty(model.rxnNotes{i})
0478             modelSBML.reaction(i).notes=[modelSBML.reaction(i).notes '<p>NOTES: ' model.rxnNotes{i} '</p>'];
0479         end
0480         modelSBML.reaction(i).notes=[modelSBML.reaction(i).notes '</body></notes>'];
0481     end
0482 
0483     % Export SBO terms from rxnMiriams
0484     if ~isempty(model.rxnMiriams{i})
0485         [~,sbo_ind] = ismember('sbo',model.rxnMiriams{i}.name);
0486         if sbo_ind > 0
0487             modelSBML.reaction(i).sboTerm=str2double(regexprep(model.rxnMiriams{i}.value{sbo_ind},'SBO:','','ignorecase'));
0488             % remove the SBO term from rxnMiriams so the information is not
0489             % duplicated in the "annotation" field later on
0490             model.rxnMiriams{i}.name(sbo_ind) = [];
0491             model.rxnMiriams{i}.value(sbo_ind) = [];
0492         end
0493     end
0494 
0495     %Export annotation information from rxnMiriams
0496     if (~isempty(model.rxnMiriams{i}) && isfield(modelSBML.reaction(i),'annotation')) || ~isempty(model.eccodes{i})
0497         modelSBML.reaction(i).annotation=['<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#meta_' model.rxns{i} '">'];
0498         modelSBML.reaction(i).annotation=[modelSBML.reaction(i).annotation '<bqbiol:is><rdf:Bag>'];
0499         if ~isempty(model.eccodes{i})
0500             eccodes=regexp(model.eccodes{i},';','split');
0501             for j=1:numel(eccodes)
0502                 modelSBML.reaction(i).annotation=[modelSBML.reaction(i).annotation  '<rdf:li rdf:resource="https://identifiers.org/ec-code/' regexprep(eccodes{j},'ec-code/|EC','') '"/>'];
0503             end
0504         end
0505         modelSBML.reaction(i).annotation=[modelSBML.reaction(i).annotation getMiriam(model.rxnMiriams{i}) '</rdf:Bag></bqbiol:is></rdf:Description></rdf:RDF></annotation>'];
0506     end
0507 
0508     if isfield(modelSBML.reaction, 'name')
0509         modelSBML.reaction(i).name=model.rxnNames{i};
0510     end
0511     if isfield(modelSBML.reaction, 'id')
0512         modelSBML.reaction(i).id=model.rxns{i};
0513     end
0514 
0515     %Add the information about reactants and products
0516     involvedMets=addReactantsProducts(model,modelSBML,i);
0517     for j=1:numel(involvedMets.reactant)
0518         if j<numel(involvedMets.reactant)
0519             modelSBML.reaction(i).reactant(j+1)=modelSBML.reaction(i).reactant(j);
0520         end
0521         modelSBML.reaction(i).reactant(j).species=involvedMets.reactant(j).species;
0522         modelSBML.reaction(i).reactant(j).stoichiometry=involvedMets.reactant(j).stoichiometry;
0523         modelSBML.reaction(i).reactant(j).isSetStoichiometry=involvedMets.reactant(j).isSetStoichiometry;
0524         modelSBML.reaction(i).reactant(j).constant=involvedMets.reactant(j).constant;
0525     end
0526     if numel(involvedMets.reactant)==0
0527         modelSBML.reaction(i).reactant='';
0528     end
0529     for j=1:numel(involvedMets.product)
0530         if j<numel(involvedMets.product)
0531             modelSBML.reaction(i).product(j+1)=modelSBML.reaction(i).product(j);
0532         end
0533         modelSBML.reaction(i).product(j).species=involvedMets.product(j).species;
0534         modelSBML.reaction(i).product(j).stoichiometry=involvedMets.product(j).stoichiometry;
0535         modelSBML.reaction(i).product(j).isSetStoichiometry=involvedMets.product(j).isSetStoichiometry;
0536         modelSBML.reaction(i).product(j).constant=involvedMets.product(j).constant;
0537     end
0538     if numel(involvedMets.product)==0
0539         modelSBML.reaction(i).product='';
0540     end
0541     %Export reversibility information. Reactions are irreversible by
0542     %default
0543     if model.rev(i)==1
0544         modelSBML.reaction(i).reversible=1;
0545     end
0546     if isfield(model, 'rxnComps')
0547         modelSBML.reaction(i).compartment=model.comps{model.rxnComps(i)};
0548     end
0549     if isfield(model, 'grRules')
0550         modelSBML.reaction(i).fbc_geneProductAssociation.fbc_association.fbc_association=model.grRules{i};
0551     end
0552     modelSBML.reaction(i).fbc_lowerFluxBound=totalNames{i};
0553     modelSBML.reaction(i).fbc_upperFluxBound=totalNames{length(model.lb)+i};
0554 end
0555 
0556 %Prepare subSystems
0557 if modelHasSubsystems
0558     modelSBML.groups_group.groups_kind = 'partonomy';
0559     modelSBML.groups_group.sboTerm     = 633;
0560     grpTemplate = modelSBML.groups_group;
0561 
0562     % === 1) Normalize: make every entry a cell array of chars (vectorized) ===
0563     isChar = cellfun(@ischar, model.subSystems);
0564     model.subSystems(isChar)  = cellfun(@(s){s}, model.subSystems(isChar),  'UniformOutput', false);
0565     model.subSystems(cellfun(@isempty, model.subSystems)) = {{}};
0566 
0567     % If some entries contain string scalars by mistake, coerce them:
0568     model.subSystems = cellfun(@(c) cellfun(@char, c, 'UniformOutput', false), model.subSystems, 'UniformOutput', false);
0569 
0570     % === 2) Flatten once: names and their reaction indices (vectorized) ===
0571     flatNames = vertcat(model.subSystems{:});                    % 1×M cellstr of all subsystem labels
0572     if isempty(flatNames)
0573         % Nothing to do: no subsystems present
0574         return
0575     end
0576 
0577     counts   = cellfun(@numel, model.subSystems);         % reactions -> how many subsystems
0578     % For each reaction r, repeat r exactly counts(r) times
0579     flatIdx  = arrayfun(@(r,c) repmat(r, c, 1), (1:numel(model.subSystems)).', counts, 'UniformOutput', false);
0580     flatIdx  = vertcat(flatIdx{:});           % M×1 vector of reaction indices
0581 
0582     % === 3) Group in one shot: unique subsystems + members per subsystem ===
0583     [subSystems, ~, g] = unique(flatNames, 'stable');          % stable preserves first appearance
0584     membersIdx = accumarray(g, flatIdx, [], @(v){v});          % cell: one vector of r-idx per group
0585     nSubs = numel(subSystems);
0586 
0587     % === 4) Preallocate and build SBML groups (single simple loop) ===
0588     modelSBML.groups_group(1:nSubs) = grpTemplate;
0589     groupIDs = "group" + (1:nSubs);
0590 
0591     if isfield(grpTemplate,'groups_member')
0592         memTemplate = grpTemplate.groups_member;
0593     else
0594         memTemplate = struct('groups_idRef','');
0595     end
0596 
0597     for i = 1:nSubs
0598         rIdx       = membersIdx{i};
0599         groupRXNs  = model.rxns(rIdx);
0600         cgroup     = grpTemplate;
0601 
0602         % Preallocate and fill members
0603         nM = numel(groupRXNs);
0604         if nM > 0
0605             cgroup.groups_member(1:nM) = memTemplate;
0606             [cgroup.groups_member.groups_idRef] = deal(groupRXNs{:});
0607         else
0608             cgroup.groups_member = repmat(memTemplate, 0, 1);
0609         end
0610 
0611         cgroup.groups_id   = char(groupIDs(i));     % keep as char for SBML compatibility
0612         cgroup.groups_name = subSystems{i};
0613         modelSBML.groups_group(i) = cgroup;
0614     end
0615 end
0616 
0617 
0618 %Prepare fbc_objective subfield
0619 
0620 modelSBML.fbc_objective.fbc_type='maximize';
0621 modelSBML.fbc_objective.fbc_id='obj';
0622 
0623 ind=find(model.c);
0624 
0625 if isempty(ind)
0626     modelSBML.fbc_objective.fbc_fluxObjective.fbc_coefficient=0;
0627 else
0628     for i=1:length(ind)
0629         %Copy the default values to the next index as long as it is not the
0630         %last one
0631         if i<numel(ind)
0632             modelSBML.reaction(i+1)=modelSBML.reaction(i);
0633         end
0634         values=model.c(model.c~=0);
0635         modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_reaction=modelSBML.reaction(ind(i)).id;
0636         modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_coefficient=values(i);
0637         modelSBML.fbc_objective(i).fbc_fluxObjective.isSetfbc_coefficient=1;
0638     end
0639 end
0640 
0641 modelSBML.fbc_activeObjective=modelSBML.fbc_objective.fbc_id;
0642 
0643 fbcStr=['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/fbc/version',num2str(sbmlPackageVersions(1))];
0644 if modelHasSubsystems
0645     groupStr=['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/groups/version',num2str(sbmlPackageVersions(2))];
0646     modelSBML.namespaces=struct('prefix',{'','fbc','groups'},...
0647         'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],...
0648         fbcStr,groupStr});
0649 else
0650     modelSBML.namespaces=struct('prefix',{'','fbc'},...
0651         'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],...
0652         fbcStr});
0653 end
0654 
0655 if sbmlPackageVersions(1) == 2
0656     modelSBML.fbc_strict=1;
0657     modelSBML.isSetfbc_strict = 1;
0658 end
0659 
0660 modelSBML.rule=[];
0661 modelSBML.constraint=[];
0662 
0663 [ravenDir,prevDir]=findRAVENroot();
0664 fileName=checkFileExistence(fileName,1,true,false);
0665 
0666 OutputSBML_RAVEN(modelSBML,fileName,1,0,[1,0]);
0667 end
0668 
0669 
0670 function modelSBML=getSBMLStructure(sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions)
0671 %Returns the blank SBML model structure by using appropriate libSBML
0672 %functions. This creates structure by considering three levels
0673 
0674 sbmlFieldNames=getStructureFieldnames('model',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0675 sbmlDefaultValues=getDefaultValues('model',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0676 
0677 for i=1:numel(sbmlFieldNames)
0678     modelSBML.(sbmlFieldNames{1,i})=sbmlDefaultValues{1,i};
0679     sbmlSubfieldNames=getStructureFieldnames(sbmlFieldNames{1,i},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0680     sbmlSubfieldValues=getDefaultValues(sbmlFieldNames{1,i},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0681     if ~strcmp(sbmlFieldNames{1,i},'event') && ~strcmp(sbmlFieldNames{1,i},'functionDefinition') && ~strcmp(sbmlFieldNames{1,i},'initialAssignment')
0682         for j=1:numel(sbmlSubfieldNames)
0683             modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j})=sbmlSubfieldValues{1,j};
0684             sbmlSubsubfieldNames=getStructureFieldnames(sbmlSubfieldNames{1,j},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0685             sbmlSubsubfieldValues=getDefaultValues(sbmlSubfieldNames{1,j},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0686             if ~strcmp(sbmlSubfieldNames{1,j},'modifier') && ~strcmp(sbmlSubfieldNames{1,j},'kineticLaw')
0687                 for k=1:numel(sbmlSubsubfieldNames)
0688                     %'compartment' and 'species' fields are not supposed to
0689                     %have their standalone structures if they are subfields
0690                     %or subsubfields
0691                     if ~strcmp(sbmlSubfieldNames{1,j},'compartment') && ~strcmp(sbmlSubfieldNames{1,j},'species')
0692                         modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j}).(sbmlSubsubfieldNames{1,k})=sbmlSubsubfieldValues{1,k};
0693                     end
0694                     %If it is fbc_association in the third level, we need
0695                     %to establish the fourth level, since libSBML requires
0696                     %it
0697                     if strcmp(sbmlSubsubfieldNames{1,k},'fbc_association')
0698                         fbc_associationFieldNames=getStructureFieldnames('fbc_association',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0699                         fbc_associationFieldValues=getDefaultValues('fbc_association',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0700                         for l=1:numel(fbc_associationFieldNames)
0701                             modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j}).(sbmlSubsubfieldNames{1,k}).(fbc_associationFieldNames{1,l})=fbc_associationFieldValues{1,l};
0702                         end
0703                     end
0704                 end
0705             end
0706         end
0707     end
0708     if ~isstruct(modelSBML.(sbmlFieldNames{1,i}))
0709         modelSBML.(sbmlFieldNames{1,i})=sbmlDefaultValues{1,i};
0710     end
0711 end
0712 
0713 modelSBML.unitDefinition.id='mmol_per_gDW_per_hr';
0714 
0715 unitFieldNames=getStructureFieldnames('unit',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0716 unitDefaultValues=getDefaultValues('unit',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0717 
0718 kinds={'mole','gram','second'};
0719 exponents=[1 -1 -1];
0720 scales=[-3 0 0];
0721 multipliers=[1 1 1*60*60];
0722 
0723 for i=1:numel(unitFieldNames)
0724     modelSBML.unitDefinition.unit(1).(unitFieldNames{1,i})=unitDefaultValues{1,i};
0725     for j=1:3
0726         modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=unitDefaultValues{1,i};
0727         if strcmp(unitFieldNames{1,i},'kind')
0728             modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=kinds{j};
0729         elseif strcmp(unitFieldNames{1,i},'exponent')
0730             modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=exponents(j);
0731         elseif strcmp(unitFieldNames{1,i},'scale')
0732             modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=scales(j);
0733         elseif strcmp(unitFieldNames{1,i},'multiplier')
0734             modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=multipliers(j);
0735         end
0736     end
0737 end
0738 end
0739 
0740 function miriamString=getMiriam(miriamStruct)
0741 %Returns a string with list elements for a miriam structure ('<rdf:li
0742 %rdf:resource="https://identifiers.org/go/GO:0005739"/>' for example). This
0743 %is just to speed up things since this is done many times during the
0744 %exporting
0745 
0746 miriamString='';
0747 if isfield(miriamStruct,'name')
0748     for i=1:numel(miriamStruct.name)
0749         miriamString=[miriamString '<rdf:li rdf:resource="https://identifiers.org/' miriamStruct.name{i} '/' miriamStruct.value{i} '"/>'];
0750     end
0751 end
0752 end
0753 
0754 function [tmp_Rxn]=addReactantsProducts(model,sbmlModel,i)
0755 %This function provides reactants and products for particular reaction. The
0756 %function was 'borrowed' from writeSBML in COBRA toolbox, lines 663-679
0757 
0758 met_idx = find(model.S(:, i));
0759 tmp_Rxn.product=[];
0760 tmp_Rxn.reactant=[];
0761 for j_met=1:size(met_idx,1)
0762     tmp_idx = met_idx(j_met,1);
0763     sbml_tmp_species_ref.species = sbmlModel.species(tmp_idx).id;
0764     met_stoich = model.S(tmp_idx, i);
0765     sbml_tmp_species_ref.stoichiometry = abs(met_stoich);
0766     sbml_tmp_species_ref.isSetStoichiometry=1;
0767     sbml_tmp_species_ref.constant=1;
0768     if (met_stoich > 0)
0769         tmp_Rxn.product = [ tmp_Rxn.product, sbml_tmp_species_ref ];
0770     else
0771         tmp_Rxn.reactant = [ tmp_Rxn.reactant, sbml_tmp_species_ref];
0772     end
0773 end
0774 end
0775 
0776 function vecT = columnVector(vec)
0777 % Code below taken from COBRA Toolbox under GNU General Public License v3.0
0778 % license file in readme/GPL.MD.
0779 %
0780 % Converts a vector to a column vector
0781 %
0782 % USAGE:
0783 %
0784 %   vecT = columnVector(vec)
0785 %
0786 % INPUT:
0787 %   vec:     a vector
0788 %
0789 % OUTPUT:
0790 %   vecT:    a column vector
0791 
0792 [n, m] = size(vec);
0793 
0794 if n < m
0795     vecT = vec';
0796 else
0797     vecT = vec;
0798 end
0799 end