Simulation of SBML with amici fails

[stobiblum]

Hey,

I am using amici 0.11.1 in the develop version on a WSL. I created an SBML model using COPASI. After exporting and importing the SBML file into a new COPASI session, I was able to successfully simulate the model giving the attached example plots. copasi_AhR_old_new_L.pdf copasi_AhR_old_new_AhRR.pdf (not all species are shown in every plot, allthough the labelling suggests it. Labelling is really not a strength of COPASI😸)

I wanted to create the same plots using amici. So I imported the SBML file, created a module & simulated it.

model_name = 'AhR_old_new'
sbml_file = f'AhR_PEtab/{model_name}/{model_name}.xml'
model_output_dir = 'amici_models/SBML_files/'+model_name

import importlib
import amici
import os
import sys
import numpy as np
import logging

#import sbml
sbml_importer = amici.SbmlImporter(sbml_file)
sbml_importer.sbml2amici(model_name,
                         model_output_dir,
                         verbose=logging.INFO)

#create model module
sys.path.insert(0, os.path.abspath(model_output_dir))
model_module = importlib.import_module(model_name)
model = model_module.getModel()

# set timepoints for which we want to simulate the model
model.setTimepoints(np.linspace(0, 60, 6001))

# Run simulation using default model parameters and solver options
rdata = amici.runAmiciSimulation(model)

However, amici is not capable of simulating the model, giving:

[Warning] AMICI:NaN: AMICI encountered a NaN value at index 36 of 47 in Jacobian!
[Warning] AMICI:NaN: AMICI encountered a NaN value at index 36 of 47 in Jacobian!
[Warning] AMICI:NaN: AMICI encountered a NaN value at index 36 of 47 in Jacobian!
[Warning] AMICI:NaN: AMICI encountered a NaN value at index 36 of 47 in Jacobian!
[Warning] AMICI:NaN: AMICI encountered a NaN value at index 36 of 47 in Jacobian!
[Warning] AMICI:NaN: AMICI encountered a NaN value at index 36 of 47 in Jacobian!
[Warning] AMICI:NaN: AMICI encountered a NaN value at index 36 of 47 in Jacobian!
[Warning] AMICI:NaN: AMICI encountered a NaN value at index 36 of 47 in Jacobian!
[Warning] AMICI:NaN: AMICI encountered a NaN value at index 36 of 47 in Jacobian!
[Warning] AMICI:NaN: AMICI encountered a NaN value at index 36 of 47 in Jacobian!
[Warning] AMICI:CVODES:CVode:CONV_FAILURE: AMICI ERROR: in module CVODES in function CVode : 
At t = 0 and h = 3.62543e-28, the corrector convergence test failed repeatedly or with |h| = hmin.
[Warning] AMICI:simulation: AMICI forward simulation failed at t = 0.000000:
AMICI failed to integrate the forward problem

Which puzzles me, because COPASI simulated the SBML without any notifications.

[dweindl]

Hi @stobiblum , any chance you can provide the SBML model? Otherwise it's rather hard to track it down.

[paulstapor]

Hi, agree with dweindl, that it's not possible to reproduce the problem without your model.

However, there might be another possibility to fix that (but that's just a guess...): It seems a bit like your ODE becomes too stiff to successfully integrate it with the given error tolerances. That might at least be reasonable, since COPASI uses different default error tolerances than AMICI (which uses stricter tolerances by default). Maybe try simulating with error tolerances like solver.setAbsoluteErrorTolerance(1e-10), solver.setRelativeErrorTolerance(1e-6), with solver = model.getSolver() before, and see, whether this works.

[stobiblum]

Sorry, here it is: AhR_old_new.zip It is zipped, because github doesn't support .xml files.

There is no function called setAbsoluteErrorTolerance() so I tried solver.setAbsoluteTolerance(1e-10), solver.setRelativeTolerance(1e-6) producing the same output. Are you sure that the COPASI default error tolerances are 1e-10 and 1e-6?

[dweindl]

The problem is, that the generated symbolic expressions for the Jacobian yield NaN for t=0, due to zero-division.

Similar to previous issues with Hill functions and loosely related to discussion in #1139.

Probably fixable by (not) applying certain sympy simplifications. Needs a closer look. Any volunteers? :)

[dweindl]

P.S.:

@stobiblum : Maybe you can have a look how in the following dwdx would have to be rearranged (using sympy) to yield a finite result in the last step:

from sympy import sympify
w = sympify('nAhR_L_ARNT**h_Tiparp/(nAhR_L_ARNT**h_Tiparp + (Kd_AhRARNT*(1 + nAhRR_ARNT)) ** h_Tiparp)')                                                                                                     

dwdx = w.diff(Symbol("nAhR_L_ARNT"))

dwdx.subs({"nAhR_L_ARNT": 0.0})

Then it should be rather easy to fix...

[stobiblum]

Thank you @dweindl

I tried to narrow the problem down, with the following code:

import sympy as sp

w1 = 'nAhR_L_ARNT**h_Tiparp/(nAhR_L_ARNT**h_Tiparp + (Kd_AhRARNT*(1 + nAhRR_ARNT)) ** h_Tiparp)'
w2 = '(Kd_AhRARNT*(1 + nAhRR_ARNT)) ** h_Tiparp'
w3 = 'nAhR_L_ARNT**h_Tiparp/(nAhR_L_ARNT**h_Tiparp)'
w4 = 'nAhR_L_ARNT**h_Tiparp'

sub1 = {'nAhR_L_ARNT': 0.0}
sub2 = {'nAhR_L_ARNT': 0.0, 'h_Tiparp': 4.0}
sub3 = {'nAhR_L_ARNT': 0.1}

e1 = sp.sympify(w1).diff(sp.Symbol('nAhR_L_ARNT')).subs(sub1)
print(e1)
e2 = sp.sympify(w2).diff(sp.Symbol('nAhR_L_ARNT')).subs(sub1)
print(e2)
e3 = sp.sympify(w3).diff(sp.Symbol('nAhR_L_ARNT')).subs(sub1)
print(e3)
e4 = sp.sympify(w4).diff(sp.Symbol('nAhR_L_ARNT')).subs(sub1)
print(e4)
e5 = sp.sympify(w4).diff(sp.Symbol('nAhR_L_ARNT')).subs(sub2)
print(e5)
e6 = sp.sympify(w4).diff(sp.Symbol('nAhR_L_ARNT')).subs(sub3)
print(e6)

e1 would be your suggestion, yielding 0.0**(2*h_Tiparp)*zoo*h_Tiparp/(0.0**h_Tiparp + (Kd_AhRARNT*(nAhRR_ARNT + 1))**h_Tiparp)**2 + 0.0**h_Tiparp*zoo*h_Tiparp/(0.0**h_Tiparp + (Kd_AhRARNT*(nAhRR_ARNT + 1))**h_Tiparp), which is not a finite solution because of zoo e2 and e3 have the finite solution 0, while e4 gives a statement including zoo: 0.0**h_Tiparp*zoo*h_Tiparp. Interestingly, there is a finite solution if h_Tiparp is defined (e5) or AhR_L_ARNT is not 0 (e6). So it seems that simplification of the differentiation nAhR_L_ARNT**h_Tiparp is problematic.

The issue might be related to https://github.com/ICB-DCM/AMICI/issues/1029

[dweindl]

So it seems that simplification of the differentiation nAhR_L_ARNT**h_Tiparp is problematic.

Hm... Not so good. Any chance to change the respective kinetic law then (i.e. the transcription function)?

[stobiblum]

Might be possible, but that would need a greater adaptation of the model.

Because there is a finite solution if h_Tiparp is defined (e5), I am wondering, why the simulation does not work for me, because h_Tiparp has a fixed value in the PEtab-files I used for the model import and should be constant throughout simulation of the SBML... I will try to remove the parametarization of h_Tiparp and hardcode instead, maybe this solves it.

[FFroehlich]

from sympy.core.parameters import evaluate
with evaluate(False):
    sp.sympify(w1).diff(sp.Symbol('nAhR_L_ARNT'))

yields

RecursionError: maximum recursion depth exceeded

💩

Also interesting why we end up with complex infinity zoo here, and not regular infinity. Regular real calculus should give a finite solution that should evaluate to finite values. Maybe we aren't setting the suitable set of assumptions?

[FFroehlich]

Indeed:

import sympy as sp
nAhR_L_ARNT = sp.Symbol('x', real=True, positive=True)
nAhRR_ARNT = sp.Symbol('y', real=True, positive=True)
h_Tiparp = sp.Symbol('h', real=True, positive=True)
Kd_AhRARNT = sp.Symbol('K_d', real=True, positive=True)

w1 = nAhR_L_ARNT**h_Tiparp/(nAhR_L_ARNT**h_Tiparp + (Kd_AhRARNT*(1 + nAhRR_ARNT)) ** h_Tiparp)

from sympy.core.parameters import evaluate
sp.simplify(sp.sympify(w1).diff(nAhR_L_ARNT))

yields

h*xh(K_dy + K_d)h/(x*(x*h + (K_dy + K_d)h)2)

which is finite unless nAhR_L_ARNT = x == 0 since powers of x appears as product both in the nominator and the denominator. Thats a removable singularity though, not sure why sympy is not further simplifying that. Having x**(h-1) in the nominator should remove the singularity and enable simulation.

As far as I remember both real and positive are assumed by default in amici. You can try to apply any of the functions from https://docs.sympy.org/latest/modules/simplify/simplify.html to get sympy to produce a denominator that won't evaluate to 0. Looks like its currently not possible to directly pass simplification functions to sbml2amici, but we should fix that and by passing the simplifications you found, amici should hopefully be able to simulate the model you provided without any issues.

[dilpath]

Looks like powsimp might be sufficient

>>> import sympy as sp
>>> sp.sympify('a^b').diff('a').subs({'a':0})
0**b*zoo*b
>>> sp.sympify('a^b').diff('a').powsimp().subs({'a':0})
0**(b - 1)*b

[dweindl]

yields

RecursionError: maximum recursion depth exceeded

hankey

Been there...

Also interesting why we end up with complex infinity zoo here, and not regular infinity. Regular real calculus should give a finite solution that should evaluate to finite values. Maybe we aren't setting the suitable set of assumptions?

That was just my example, was too lazy to set to real, but may make a difference in the simplifications.

As far as I remember both real and positive are assumed by default in amici.

Only real, not positive.

May be interesting to have as additional option.

Looks like its currently not possible to directly pass simplification functions to sbml2amici

Yeah, also noticed that. Was sure I had that in there, but apparently not. Should be changed. Will do.

Looks like powsimp might be sufficient

That's already being used. Added it because of trouble with some Hill function some time ago. There it did the job. In this case unfortunately not.

[dilpath]

Looks like powsimp might be sufficient

That's already being used. Added it because of trouble with some Hill function some time ago. There it did the job. In this case unfortunately not.

Ahh why not?

>>> import sympy as sp
>>> w = 'nAhR_L_ARNT**h_Tiparp/(nAhR_L_ARNT**h_Tiparp + (Kd_AhRARNT*(1 + nAhRR_ARNT)) ** h_Tiparp)'
>>> dwdx = sp.sympify(w).diff('nAhR_L_ARNT')
>>> dwdx.subs({'nAhR_L_ARNT': 0.0})
0.0**(2*h_Tiparp)*zoo*h_Tiparp/(0.0**h_Tiparp + (Kd_AhRARNT*(nAhRR_ARNT + 1))**h_Tiparp)**2 + 0.0**h_Tiparp*zoo*h_Tiparp/(0.0**h_Tiparp + (Kd_AhRARNT*(nAhRR_ARNT + 1))**h_Tiparp)

Above is the zoo issue. Below is the "fix" with powsimp.

>>> dwdx.powsimp().subs({'nAhR_L_ARNT': 0.0})
0.0**(h_Tiparp - 1)*h_Tiparp/(0.0**h_Tiparp + (Kd_AhRARNT*(nAhRR_ARNT + 1))**h_Tiparp) - 0.0**(2*h_Tiparp - 1)*h_Tiparp/(0.0**h_Tiparp + (Kd_AhRARNT*(nAhRR_ARNT + 1))**h_Tiparp)**2

Assuming nAhRR_ARNT>=0, Kd_AhRARNT>0, and h_Tiparp>=1, this should evaluate to zero? There are issues if h_Tiparp<1 but h_Tiparp==4 in the SBML?

If I include the initial values from the SBML, without powsimp:

>>> dwdx.subs({'nAhR_L_ARNT': 0.0}).subs({'h_Tiparp': 4})
nan
>>> dwdx.subs({'nAhR_L_ARNT': 0.0}).subs({'h_Tiparp': 4, 'Kd_AhRARNT': 0.1, 'nAhRR_ARNT': 0})
nan

With powsimp:

>>> dwdx.powsimp().subs({'nAhR_L_ARNT': 0.0}).subs({'h_Tiparp': 4})
0
>>> dwdx.powsimp().subs({'nAhR_L_ARNT': 0.0}).subs({'h_Tiparp': 4, 'Kd_AhRARNT': 0.1, 'nAhRR_ARNT': 0})
0

[dweindl]

Assuming nAhRR_ARNT>=0, Kd_AhRARNT>0

That may be the missing ingredient. Currently we only assume real, but no non-negativity.

[dweindl]

>>> dwdx.powsimp().subs({'nAhR_L_ARNT': 0.0}).subs({'h_Tiparp': 4, 'Kd_AhRARNT': 0.1, 'nAhRR_ARNT': 0})
0

@dilpath : The difference there is that with this substitution the expression can still be rearranged for evaluation. In AMICI we generate static expressions that are then evaluated as is. Another problem related to that is that even 0*nan = nan (not 0) .

[FFroehlich]

import sympy as sp
nAhR_L_ARNT = sp.Symbol('x', real=True)
nAhRR_ARNT = sp.Symbol('y', real=True)
h_Tiparp = sp.Symbol('h', real=True)
Kd_AhRARNT = sp.Symbol('K_d', real=True)

w1 = nAhR_L_ARNT**h_Tiparp/(nAhR_L_ARNT**h_Tiparp + (Kd_AhRARNT*(1 + nAhRR_ARNT)) ** h_Tiparp)

print(sp.powsimp(sp.sympify(w1).diff(nAhR_L_ARNT)))

Yields h*x**(h - 1)/(x**h + (K_d*(y + 1))**h) - h*x**(2*h - 1)/(x**h + (K_d*(y + 1))**h)**2, which, as far as I can tell, should evaluate fine, so the positivity assumption does not seem necessary. Will have to check whether thats actually what is happening in AMICI and whether thats is indeed whats causing the problems.

[dweindl]

I hope I copied the right one :grimacing: . But it was not the only issue, there were a couple more NaNs in J. But all related to that transcription kinetic law afaik.

[FFroehlich]

Looks like we have to have a serious discussion with sympy:

sp.powsimp applied to the whole matrix doesnt to anything. sp.powsimp applied to the individual entries does the job ...

This may have changed with the 1.6 update but we didn't notice so far.

💩

[dweindl]

:hankey: :hankey: :hankey: :hankey: :hankey:

Okay, that explains that...

[FFroehlich]

Okay looks like powsimp applied to the whole matrix was interpreted to simplify the respective matrix-valued expression (which is just the matrix, so no simplification performed), which actually makes sense. Some deprecation warning would have been nice though ...

[lcontento]

I think we should use powsimp with the argument deep=True. It looks like it broadcasts inside matrices (and also inside function arguments, which I think is something we want).

[FFroehlich]

I think we should use powsimp with the argument deep=True. It looks like it broadcasts inside matrices (and also inside function arguments, which I think is something we want).

Hmm documentation states it broadcasts inside functions, can't find anything on matrices, but that may be a side effect. I fixed this. It's already fixed in #1158, by using applyfunc, which is probably the cleanest option. Specifying deep=True nevertheless looks like a good default, so we should change that.

[FFroehlich]

@stobiblum this issue should be resolved in the new 0.11.2 release that should be available from pip shortly, if not, please let us know.