Support RateRules on Parameters

[matthiaskoenig]

What did you expect to happen?
I tried to load one of my SBML models in AMICI. Almost all my models have RateRules on parameters.

What has happened instead?
SBML model does not load.

/home/mkoenig/envs/amici/bin/python /home/mkoenig/git/amici-git/amici_example.py
Traceback (most recent call last):
  File "/home/mkoenig/git/amici-git/amici_example.py", line 23, in <module>
    sbml_importer.sbml2amici(
  File "/home/mkoenig/envs/amici/lib/python3.8/site-packages/amici/sbml_import.py", line 318, in sbml2amici
    self._process_sbml(constant_parameters)
  File "/home/mkoenig/envs/amici/lib/python3.8/site-packages/amici/sbml_import.py", line 358, in _process_sbml
    self.check_support()
  File "/home/mkoenig/envs/amici/lib/python3.8/site-packages/amici/sbml_import.py", line 393, in check_support
    raise SBMLException('Algebraic rules are currently not supported, '
amici.sbml_import.SBMLException: Algebraic rules are currently not supported, and rate rules are only supported for species and compartments.

To Reproduce
Run code on model attached
amici_example.zip

AMICI version and system environment

• OS and version: Ubuntu 20.04
• AMICI interface: python
• AMICI version: 0.11.7

How to fix
Support RateRules on parameters.

[FFroehlich]

Hmm, supporting RateRules on parameter is a little bit tricky since it doesn't work well with the amici definition of a parameter, i.e. something with respect to which you can compute a parameter. Probably the best approach would be to create state variables for these species and use the existing machinery for rate rules.

[matthiaskoenig]

Many reasons why these rate rules cannot be encoded on parameters

• unit validation will break (because species are in amount units); having time changing parameters/states which are not in amount units is important; unit validation is very important because it allows me to automatically map and unit convert all my parameters, species and data for my parameter fitting problems and simulations
• subsequent analysis relying on the species-reaction graph will break

The problem is that AMICI is trying to do a 1:1 mapping of SBML objects on concepts in the parameter fitting problem. But what you have to do is figure out which objects in SBML correspond to which parts in your parameter fitting problem:

State variables:

• species which are constant=false and boundaryCondition=false and are not set by an AssignmentRule
• compartments which are constant=false and set by an RateRule
• parameters which are constant=false and set by an RateRule

Parameters in fitting problem

• compartments which are constant=true
• parameters which are constant=true

You have to figure out the subets of the SBML model which correspond to your parameter fitting problem.

[paulstapor]

To some degree, we already started to support parameters with assignment rules, which can be turned into splines.
So in the long run, I'd agree that it can make sense to allow also parameters to become state variables or time-dependent. The suggestion you make makes sense @matthiaskoenig ...
That would indeed mean reevaluating some concepts we use in AMICI internally, and, if changing the concepts we use, would mean quite some work... Might be worth it, though... But we should discuss that.
Anyhow, there are two problems in this:

1. How are quantities handled in AMICI internally? (We have state variables x_solver, conserved quantities tcl, parameters p, expressions w, observables y, and later we will have event observables z).
2. Which quantites are returned by AMICI and how? (We have state variables x_rdata, parameters p, and observables y, later we will have event observables z).

One could anyhow add, beyond state variables, observables, parameters, another output such as "derived quantities" to the output. For example, that could be used to output parameters with rate or assignment rules, for example, splines or whatever kind of input functions which are non-constant. The same thing we might want to do with species with assignment rules. Currently, those are internally pruned out of the state variables and returned as observables, as far as I know...

In that context, I'd vote to keep species -> x_rdata (returned state variables), as this indicates also whether a quantity way modeled as parameter or as species... and then one could return time-varying parameters as "derived quantities".
But this needs to be discussed in detail. Possibly in this thread?

[matthiaskoenig]

Sorry, I don't even understand what the problem is with AMICI and RateRules on parameters.
A rate rule on a parameter is just
dp/dt = f(p, ....)
You put p in the state variables and treat f like the right hand sides of the ODEs.
The SBML parameter p simulated via a rate rule is NOT a parameter in AMICI or parameter fitting, but a state variable.
There is not a 1:1 mapping between parameters in SBML to parameters in an optimization problem.
Why is this not working in AMICI?

[paulstapor]

In principle, there is no big deal with this, you're right. We could more or less easilyimplement it, I think...
The thing is that it touches (to some degree) the way AMICI handles things behind the scenes. And we might take this opportunity to discuss the question how things are handled internally for a moment. Maybe it makes sense change (in the long run) more than only implementing parameters with rate rules...
So, my questions were actually directed to the other developers of AMICI.

[paulstapor]

Concerning the 1:1 mapping:
Yes, we're not aiming at a 1:1 mapping between SBML parameters and optimization parameters. But I'd be happy to have a 1:1 mapping between SBML species and the state variables which AMICI returns as quantity called "x". This I still think makes sense...

[dweindl]

Sorry, I don't even understand what the problem is with AMICI and RateRules on parameters.

I'd say the biggest problem is, that it wasn't considered to be important enough so far or nobody felt like implementing it...

But I'd be happy to have a 1:1 mapping between SBML species and the state variables which AMICI returns as quantity called "x". This I still think makes sense...

Not so sure. I think it would absolutely make sense to put a wrapper around the internal result representation to map things back to model species entities. I assume most users don't really care how those things are computed/represented internally, but just want to get the value/trajectory for some given entity based on the respective ID. (Something similar would anyways be required for addressing #11, if ever).

[matthiaskoenig]

Not sure what you mean with "value/trajectory for some given entity based on the respective ID." Is this only referring to species? I have trajectories of parameters which are calculated via rate rules which I am using in parameter fitting. E.g. I have a parameter which calculates via a rate rule the time-dependent injection rate of a given substance. The parameter is a flux in [mmol/min] which cannot be encoded as a species (species are only in amount units). I want to plot it and optimize parameters in the rate rule expression to minimize the distance between experimental injection rate and model injection rate. In the field of glucose metabolism there are many tracer studies which are giving experimental rates which can be compared to reaction rates or parameters which encode these rates. I want to directly access these parameter trajectories via their id.

[dweindl]

Not sure what you mean with "value/trajectory for some given entity based on the respective ID." Is this only referring to species?

No, species/compartment/parameter/rule/... (I wrote 'species' where I didn't mean it in the previous post)

I want to directly access these parameter trajectories via their id.

Right.

[paulstapor]

I see, why you don't want to have it as a species... And I agree that this should be doable. When simply adding it to the state variables (for the moment), this should also not create too much of work...
Are there volunteers for implementing? (*hopingForSomebodyToJumpIn*...)
If not, I can have a look into this at some point this week. If it turns out to be only little work, I would give it a try...

[matthiaskoenig]

I can support you here from the SBML side and explain the different cases.
A similar issue will most likely appear in compartments depending on RateRules, but is the same solution like the RateRules on parameters. I am not using RateRules on compartments, but models of bioreactors sometimes have these.

[FFroehlich]

Will try to get this fixed over the next couple of days, shouldn't be too much work. Fixing this disconnect between SBML and AMICI definition of parameters will hopefully also resolve other issue but but may be tricky to implement this without breaking too much backwards compatibility.

[dweindl]

Fixing this disconnect between SBML and AMICI definition of parameters will hopefully also resolve other issue but but may be tricky to implement this without breaking too much backwards compatibility.

Would prioritize convenient API over backwards compatibility here.

I think it's nice to have both an internal and an SBML result interface. (The former, mostly for debugging / performance, the latter for convenience). For example, when looking at the results, I might not care whether a certain parameter was a constant in AMICI or not. With those changes we can also handle constant species better (i.e. treat as parameters, as already written by Matthias above). This will also allow changing such internals later on transparently for most users.

[matthiaskoenig]

Thanks. Great for fixing this.