Best way to get started with SBML and AMICI

[winkmal]

For my research in simulation of anaerobic digesters, it seemed AMICI had some very useful features for my usecase. I have mostly used Matlab/Simulink, and sometimes Python-based tools like SciPy. However, since I am new to SBML standard, tools and notation, I cannot really find a good way to get started. My models are stiff nonlinear ODE systems of about 12-34 differential states, and several algebraic output variables computed from those states, which are needed to compare them to the (scarcely available) measurements, and eventually to construct a paramter estimation routine. What I tried so far:

• [x] Code and test a very basic model in Matlab SymBiology Toolbox --> worked til a certain point
• [x] Export the working Matlab example to SBML
• [x] Import the SBML file to AMICI --> Gives an error about different SBML standards
• [x] Ask a popular chatbot 🤖 to write/translate the code for me --> Gives syntactically correct, but non-functional Python code

All the examples in the AMICI doc are based on importing an already existing SBML file. Could you provide some hints on what is the best way to convert a set of equations/Matlab/SciPy-Code into an Amici ODE model?

[stephanmg]

Dear @winkmal,

could you provide details (log) of the import error when importing your SBML file into AMICI (... and the actual SBML file)?

[FFroehlich]

hi @winkmal, you probably want to check out https://github.com/yaml2sbml-dev/yaml2sbml, which allows you to directly specify model equations and is probably what you want.

Nevertheless, I would encourage you to look at https://simplesbml.readthedocs.io/en/latest/ or https://pysb.readthedocs.io/en/stable/ . 12-34 differential states is quite a bit and, in my experience, errors are easy to make but hard to notice, so adding some abstraction and improving readability will go a long way.

[winkmal]

could you provide details (log) of the import error when importing your SBML file into AMICI (... and the actual SBML file)?

Sure, please find the Matlab code, the SBML file, the Python import code snippet and the error log in this public gist.

[stephanmg]

Thanks for providing the details:

SBMLException: Algebraic rules are currently not supported, and rate rules are only supported for species, compartments, and parameters.

Seems like DAEs are simply not fully supported in SBML at the moment? (see the SBML Level 2 specification)

[FFroehlich]

Thanks for providing the details:

SBMLException: Algebraic rules are currently not supported, and rate rules are only supported for species, compartments, and parameters.

Seems like DAEs are simply not fully supported in SBML at the moment? (see the SBML Level 2 specification)

DAEs can be specified in SBML, but are currently not supported in AMICI.

[stephanmg]

Thanks for providing the details: SBMLException: Algebraic rules are currently not supported, and rate rules are only supported for species, compartments, and parameters. Seems like DAEs are simply not fully supported in SBML at the moment? (see the SBML Level 2 specification)

DAEs can be specified in SBML, but are currently not supported in AMICI.

Thanks for clarifying! Found contradicting statements by a quick web search.

[winkmal]

DAEs can be specified in SBML, but are currently not supported in AMICI.

Well, this is actually not a DAE. These are just two ODEs (which could be solved algebraically). Eq. 3-4 are only algebraic output variables depending on the differential state variables.

\begin{align}
\frac{\textnormal{d} m_\mathrm{DVS1}}{\textnormal{d} t} &= -k_1 \cdot m_\mathrm{DVS1} \\
\frac{\textnormal{d} m_\mathrm{DVS2}}{\textnormal{d} t} &= -k_2 \cdot m_\mathrm{DVS2} \\
m_\mathrm{VS}       &= m_\mathrm{DVS1} + m_\mathrm{DVS2} + m_\mathrm{NDVS} \\
m_\mathrm{NDVS}     &= (1 - f_1 - f_2)*m_\mathrm{VS,0} \\ 
m_\mathrm{DVS1}(0)  &= f_1 \cdot m_\mathrm{VS,0} \\
m_\mathrm{DVS2}(0)  &= f_2 \cdot m_\mathrm{VS,0} \\
\end{align}

Eq. 5-6 give the initial conditions. The reason to write it that way is that only $m_\mathrm{VS}$ can be determined experimentally, while the other three cannot.

[FFroehlich]

In your gist you specify

% Algebraic Rule: Must sum up to 0
addrule(mass2frak,'m_VS - (m_DVS1 + m_DVS2 + m_NDVS)', 'RuleType','algebraic');

AMICI in general does not support algebraic SBML rules, whether they define a DAE or not. Nevertheless, I would argue that what you specify pretty much looks like a (trivial) DAE where the temporal evolution of m_VS only defined by an algebraic component and no differential component. Should be straightforward to solve that algebraic equation in closed form and set that as assignment rule.

[winkmal]

AMICI in general does not support algebraic SBML rules, whether they define a DAE or not. Nevertheless, I would argue that what you specify pretty much looks like a (trivial) DAE where the temporal evolution of m_VS only defined by an algebraic component and no differential component. Should be straightforward to solve that algebraic equation in closed form and set that as assignment rule.

Well, if I understand correctly, you are saying that I need to change the RuleType? According to the documentation of SimBiology®, there are four types of rules:

initialAssignment — Lets you specify the initial value of a parameter, species, or compartment capacity, as a function of other model component values in the model.

repeatedAssignment — Lets you specify a value that holds at all times during simulation, and is a function of other model component values in the model.

algebraic — Lets you specify mathematical constraints on one or more parameters, species, or compartments that must hold during a simulation.

rate — Lets you specify the time derivative of a parameter value, species amount, or compartment capacity.

So instead of algebraic, I should just use repeatedAssignment, and it will work?

[dweindl]

So instead of algebraic, I should just use repeatedAssignment, and it will work?

It will work, in the sense of you will be able to import and simulate the model. m_VS will take the value of m_DVS1 + m_DVS2 + m_NDVS.

It will not work in the sense of enforcing m_VS - (m_DVS1 + m_DVS2 + m_NDVS) == 0. But what I understand from your problem description is that this is not what you wanted anyways, and only a misunderstanding of the different types of rules.

[winkmal]

With some help, I was able to replace the algebraic rule with repeatedAssignment, and regenerate the SBML file.

Now, AMICI complains about the unit liter. I updated the gist with the new *.m script, SBML file and importError.log. It says:

amici.sbml_import - ERROR - libSBML SBML component consistency (Error): The value of a <species>'s 'units' attribute is restricted. 
Reference: L2V4 Section 4.8.5
 The value of a <species>'s 'substanceUnits' attribute can only be one of the following: 'substance', 'mole', 'item', 'gram', 'kilogram', 'dimensionless', or the identifier of a <unitDefinition> derived from 'mole' (with an 'exponent' of '1'), 'item' (with an 'exponent' of '1'), 'gram' (with an 'exponent' of '1'), 'kilogram' (with an 'exponent' of '1'), or 'dimensionless'.  The current value ('MWBUILTINUNIT_liter') is not allowed.

[dweindl]

Now, AMICI complains about the unit liter.

I think changing V_ch4 from species to parameter should do the job.

(It is libsbml complaining about an invalid SBML model. AMICI itself does not care about units.)

[winkmal]

I think changing V_ch4 from species to parameter should do the job.

Well, I introduced m_ch4 in kilogram, and added V_ch4 as an Observable, so import is finally working! I updated the SBML files in the Gist. (Not the *.m script, may do that eventually.) Thanks for all your help.