https://github.com/Data2Dynamics/d2d
Tip revision: 83d6eedd5f6c33ae94dbce6c4e620d2f1bb57b1f authored by Joep Vanlier on 09 April 2015, 08:41:55 UTC
Forgot to update C version number
Forgot to update C version number
Tip revision: 83d6eed
arExportSBML.m
% export model to SBML using libSBML
%
% function arExportSBML(m, c, copasi)
%
% m: model index
% c: condition index
% copasi: copasi compatibility mode
function arExportSBML(m, c, copasi)
global ar
if(~exist([cd '/SBML' ], 'dir'))
mkdir([cd '/SBML' ])
end
if(~exist('copasi','var'))
copasi = false;
end
M = TranslateSBML(which('empty.xml'));
F = TranslateSBML(which('filled.xml'));
M.id = ar.model(m).name;
M.notes = ar.model(m).description{1};
%% compartements
if(~isempty(ar.model(m).c))
for jc = 1:length(ar.model(m).c)
M.compartment(jc).typecode = 'SBML_COMPARTMENT';
M.compartment(jc).metaid = '';
M.compartment(jc).notes = '';
M.compartment(jc).annotation = '';
M.compartment(jc).sboTerm = -1;
M.compartment(jc).name = '';
M.compartment(jc).id = ar.model(m).c{jc};
M.compartment(jc).compartmentType = '';
M.compartment(jc).spatialDimensions = 3;
M.compartment(jc).constant = 1;
M.compartment(jc).units = '';
M.compartment(jc).outside = '';
M.compartment(jc).isSetSize = 1;
M.compartment(jc).isSetVolume = 1;
M.compartment(jc).level = 2;
M.compartment(jc).version = 4;
qp = ismember(ar.pLabel, ar.model(m).pc{jc}); %R2013a compatible
if(sum(qp)==1)
pvalue = ar.p(qp);
if(ar.qLog10(qp))
pvalue = 10^pvalue;
end
M.compartment(jc).size = pvalue;
elseif(sum(qp)==0)
qp = ismember(ar.model(m).condition(c).pold, ar.model(m).pc{jc}); %R2013a compatible
if(sum(qp)==1)
pvalue = ar.model(m).condition(c).fp{qp};
M.compartment(jc).size = str2double(pvalue);
else
pvalue = str2double(ar.model(m).pc{jc});
if(~isnan(pvalue))
M.compartment(jc).size = pvalue;
else
error('%s not found', ar.model(m).pc{jc});
end
end
else
error('%s not found', ar.model(m).pc{jc});
end
end
else
M.compartment(1).typecode = 'SBML_COMPARTMENT';
M.compartment(1).metaid = '';
M.compartment(1).notes = '';
M.compartment(1).annotation = '';
M.compartment(1).sboTerm = -1;
M.compartment(1).name = '';
M.compartment(1).id = 'default';
M.compartment(1).compartmentType = '';
M.compartment(1).spatialDimensions = 3;
M.compartment(1).constant = 1;
M.compartment(1).units = '';
M.compartment(1).outside = '';
M.compartment(1).isSetSize = 1;
M.compartment(1).isSetVolume = 1;
M.compartment(1).level = 2;
M.compartment(1).version = 4;
M.compartment(1).size = 1;
end
%% species
Crules = {};
for jx = 1:length(ar.model(m).x)
M.species(jx).typecode = 'SBML_SPECIES';
M.species(jx).metaid = '';
M.species(jx).notes = '';
M.species(jx).annotation = '';
M.species(jx).sboTerm = -1;
M.species(jx).name = '';
M.species(jx).id = ar.model(m).x{jx};
M.species(jx).speciesType = '';
if(~isempty(ar.model(m).cLink))
M.species(jx).compartment = ar.model(m).c{ar.model(m).cLink(jx)};
else
M.species(jx).compartment = 'default';
end
M.species(jx).initialAmount = NaN;
M.species(jx).substanceUnits = '';
M.species(jx).hasOnlySubstanceUnits = 0;
M.species(jx).boundaryCondition = 0;
M.species(jx).charge = 0;
M.species(jx).constant = 0;
M.species(jx).isSetInitialAmount = 0;
M.species(jx).isSetInitialConcentration = 1;
M.species(jx).isSetCharge = 0;
M.species(jx).level = 2;
M.species(jx).version = 4;
qp = ismember(ar.pLabel, ar.model(m).px0{jx}); %R2013a compatible
if(sum(qp)==1)
M.species(jx).initialConcentration = 1;
Crules{end+1,1} = ar.model(m).x{jx}; %#ok<AGROW>
Crules{end,2} = ar.pLabel{qp}; %#ok<AGROW>
elseif(sum(qp)==0)
qp = ismember(ar.model(m).condition(c).pold, ar.model(m).px0{jx}); %R2013a compatible
if(sum(qp)==1)
pvalue = char(sym(ar.model(m).condition(c).fp{qp}));
if(~isnan(str2double(pvalue)))
pvalue = str2double(pvalue);
M.species(jx).initialConcentration = pvalue;
else
Crules{end+1,1} = ar.model(m).x{jx}; %#ok<AGROW>
Crules{end,2} = pvalue; %#ok<AGROW>
M.species(jx).initialConcentration = 1;
end
else
error('%s not found', ar.model(m).pc{jc});
end
else
error('%s not found', ar.model(m).pc{jc});
end
end
%% parameters
for jp = 1:length(ar.model(m).condition(c).p)
M.parameter(jp).typecode = 'SBML_PARAMETER';
M.parameter(jp).metaid = '';
M.parameter(jp).notes = '';
M.parameter(jp).annotation = '';
M.parameter(jp).sboTerm = -1;
M.parameter(jp).name = '';
M.parameter(jp).id = ar.model(m).condition(c).p{jp};
M.parameter(jp).units = '';
M.parameter(jp).constant = 1;
M.parameter(jp).isSetValue = 1;
M.parameter(jp).level = 2;
M.parameter(jp).version = 4;
qp = ismember(ar.pLabel, ar.model(m).condition(c).p{jp}); %R2013a compatible
if(sum(qp)==1)
pvalue = ar.p(qp);
if(ar.qLog10(qp) == 1)
pvalue = 10^pvalue;
end
else
pvalue = 1;
end
M.parameter(jp).value = pvalue;
end
%% rules
for jr = 1:size(Crules,1)
M.initialAssignment(jr).typecode = 'SBML_INITIAL_ASSIGNMENT';
M.initialAssignment(jr).metaid = '';
M.initialAssignment(jr).notes = '';
M.initialAssignment(jr).annotation = '';
M.initialAssignment(jr).sboTerm = -1;
M.initialAssignment(jr).symbol = Crules{jr,1};
M.initialAssignment(jr).math = Crules{jr,2};
M.initialAssignment(jr).level = 2;
M.initialAssignment(jr).version = 4;
end
%% reactions
fv = ar.model(m).fv;
fv = sym(fv);
fv = subs(fv, ar.model(m).u, ar.model(m).fu);
fv = subs(fv, ar.model(m).condition(c).pold, ar.model(m).condition(c).fp');
vcount = 1;
arWaitbar(0);
for jv = 1:length(ar.model(m).fv)
arWaitbar(jv, length(ar.model(m).fv));
ratetemplate = fv(jv);
if(ratetemplate~=0)
M.reaction(vcount).typecode = 'SBML_REACTION';
M.reaction(vcount).metaid = '';
M.reaction(vcount).notes = '';
M.reaction(vcount).annotation = '';
M.reaction(vcount).sboTerm = -1;
M.reaction(vcount).name = '';
M.reaction(vcount).reversible = 0;
M.reaction(vcount).fast = -1;
M.reaction(vcount).isSetFast = 0;
M.reaction(vcount).level = 2;
M.reaction(vcount).version = 4;
if(isfield(ar.model(m),'v') && ~isempty(ar.model(m).v{jv}))
M.reaction(vcount).id = ar.model(m).v{jv};
else
M.reaction(vcount).id = sprintf('reaction%i', jv);
end
scount = 1;
scomp = [];
for jsource = find(ar.model(m).N(:,jv)<0)'
M.reaction(vcount).reactant(scount).typecode = 'SBML_SPECIES_REFERENCE';
M.reaction(vcount).reactant(scount).metaid = '';
M.reaction(vcount).reactant(scount).notes = '';
M.reaction(vcount).reactant(scount).annotation = '';
M.reaction(vcount).reactant(scount).sboTerm = -1;
M.reaction(vcount).reactant(scount).species = ar.model(m).x{jsource};
M.reaction(vcount).reactant(scount).id = '';
M.reaction(vcount).reactant(scount).name = '';
M.reaction(vcount).reactant(scount).stoichiometry = abs(ar.model(m).N(jsource,jv));
M.reaction(vcount).reactant(scount).stoichiometryMath = F.reaction.reactant.stoichiometryMath;
M.reaction(vcount).reactant(scount).level = 2;
M.reaction(vcount).reactant(scount).version = 4;
scount = scount + 1;
if(~isempty(scomp) && scomp~=ar.model.cLink(jsource))
error('influx from different compartments in reaction %i', jv);
end
if(~isempty(ar.model(m).cLink))
scomp = ar.model.cLink(jsource);
end
end
scount = 1;
tcomp = [];
for jsource = find(ar.model(m).N(:,jv)>0)'
M.reaction(vcount).product(scount).typecode = 'SBML_SPECIES_REFERENCE';
M.reaction(vcount).product(scount).metaid = '';
M.reaction(vcount).product(scount).notes = '';
M.reaction(vcount).product(scount).annotation = '';
M.reaction(vcount).product(scount).sboTerm = -1;
M.reaction(vcount).product(scount).species = ar.model(m).x{jsource};
M.reaction(vcount).product(scount).id = '';
M.reaction(vcount).product(scount).name = '';
M.reaction(vcount).product(scount).stoichiometry = abs(ar.model(m).N(jsource,jv));
M.reaction(vcount).product(scount).stoichiometryMath = F.reaction.reactant.stoichiometryMath;
M.reaction(vcount).product(scount).level = 2;
M.reaction(vcount).product(scount).version = 4;
scount = scount + 1;
if(~isempty(tcomp) && tcomp~=ar.model.cLink(jsource))
error('efflux to different compartments in reaction %i', jv);
end
if(~isempty(ar.model(m).cLink))
tcomp = ar.model.cLink(jsource);
end
end
vars = symvar(ratetemplate);
vars = setdiff(vars, sym(ar.model(m).x(ar.model(m).N(:,jv)<0))); %R2013a compatible
vars = setdiff(vars, sym(ar.model(m).condition(c).p)); %R2013a compatible
if(~isempty(vars))
for jmod = 1:length(vars);
M.reaction(vcount).modifier(jmod).typecode = 'SBML_MODIFIER_SPECIES_REFERENCE';
M.reaction(vcount).modifier(jmod).metaid = '';
M.reaction(vcount).modifier(jmod).notes = '';
M.reaction(vcount).modifier(jmod).annotation = '';
M.reaction(vcount).modifier(jmod).sboTerm = -1;
M.reaction(vcount).modifier(jmod).species = char(vars(jmod));
M.reaction(vcount).modifier(jmod).id = '';
M.reaction(vcount).modifier(jmod).name = '';
M.reaction(vcount).modifier(jmod).level = 2;
M.reaction(vcount).modifier(jmod).version = 4;
end
end
M.reaction(vcount).kineticLaw.typecode = 'SBML_KINETIC_LAW';
M.reaction(vcount).kineticLaw.metaid = '';
M.reaction(vcount).kineticLaw.notes = '';
M.reaction(vcount).kineticLaw.annotation = '';
M.reaction(vcount).kineticLaw.sboTerm = -1;
if(~isempty(ar.model(m).cLink))
if(~copasi)
M.reaction(vcount).kineticLaw.formula = char(ratetemplate);
elseif(~isempty(scomp) && ~isempty(tcomp) && scomp~=tcomp) % multi-compartment reaction
M.reaction(vcount).kineticLaw.formula = [ar.model(m).c{scomp} ' * (' char(ratetemplate) ')'];
else
if(~isempty(scomp))
M.reaction(vcount).kineticLaw.formula = [ar.model(m).c{scomp} ' * (' char(ratetemplate) ')'];
elseif(~isempty(tcomp))
M.reaction(vcount).kineticLaw.formula = [ar.model(m).c{tcomp} ' * (' char(ratetemplate) ')'];
else
error('scomp and tcomp empty');
end
end
else
M.reaction(vcount).kineticLaw.formula = char(ratetemplate);
end
M.reaction(vcount).kineticLaw.math = '';
M.reaction(vcount).kineticLaw.parameter = F.reaction.kineticLaw.parameter;
M.reaction(vcount).kineticLaw.level = 2;
M.reaction(vcount).kineticLaw.version = 4;
vcount = vcount + 1;
end
end
arWaitbar(-1);
[a,b] = isSBML_Model(M);
if(a == 1)
if(~copasi)
OutputSBML(M, ['SBML/' ar.model(m).name '_l2v4.xml']);
else
OutputSBML(M, ['SBML/' ar.model(m).name '_copasi_l2v4.xml']);
end
else
error('%s', b);
end