problem calculating steady state of a large model

[thomassligon]

What did you expect to happen?

steady state calculation should succeed

What has happened instead?

error messages in attached file asthma-2021-01-01.txt

To Reproduce Steps to reproduce the behavior Steady state calculation is working in Copasi. However, I need to set the resolution to 1e-6 instead of the default 1e-9.

Ideally include minimal code examples here

AMICI version and system environment

• OS and version: Windows 10
• AMICI interface: Python plus pyPESTO and PETab
• AMICI version: 0.11.7
• Additional information: Compiler name and version used, Python/Matlab version, ... Visual Studio C++ compiler

How to fix Do you know how to resolve the problem? Can you submit a pull request?

How can I improve steady state in Copasi (i.e. resolution of 1e-9)? What do I need to fix in order to get AMICI steady state working? How do I find the reactions that need attention?

[dweindl]

Hi Tom,

How can I improve steady state in Copasi (i.e. resolution of 1e-9)?

That is something you would better ask the Copasi developers.

What do I need to fix in order to get AMICI steady state working? How do I find the reactions that need attention?

From your output file:

[Warning] AMICI:NaN: AMICI encountered a NaN value at index 0/286 in fxdot!
[Warning] AMICI:NaN: AMICI encountered a NaN value at index 0/688 in p!
[Warning] AMICI:NaN: AMICI encountered a NaN value at index 55/469 in w!

This gives you some hints for which reactions to check. It's not a fun job, but you can check the AMICI-generated .cpp files to find the expressions leading to those NaNs.

[paulstapor]

What do I need to fix in order to get AMICI steady state working? How do I find the reactions that need attention?

It depends a bit on the solver settings which are made for steady state computation in AMICI. What I find particularly puzzling and what's always worth some investigation: You have a NaN value in P, so in your parameters. That's usually not a good sign an means that something has gone wrong somewhere. Try to see which parameter this is and analyze in a debugger why it might be a NaN value. With a NaN as parameter (and hence maybe in the ODE system), there are always problems... Try to find the reason, and you have at least a chance that things will work then.

[thomassligon]

Thanks Daniel and Paul! Both comments look useful. When I first saw these error messages, I thought something fundamental must be wrong. In fact, Steady state sensitvitiy computation failed due to unsuccessful factorization of RHS Jacobian might lead to all NaNs being passed from sensitivity computation to the next step, which would mean that the NaN encountered messages are irrelevant. For example, NaN value at index 0/688 in p! is at index 0, so most likely everything is a NaN here. On the other hand, NaN value at index 55/469 in w! is more likely to be a useful message, since it is somewhere inside of the structure.

[paulstapor]

Ah, now I see the trouble!

If the messages are from two consecutive steps, we probably have the following situation: step n:

• steady state was found (probably by simulating a long time)
• AMICI tries to compute state sensitivities via solving the linear steady state sensitivity system. This will go wrong if the Jacobian is singular (possibly: numerically singular)
  • for whatever reason, the optimizer doesn't reject this step and computes an update (all NaN) from the NaN-gradient/the sensitivities

step n + 1:

• the updated parameter vector is all NaN
• obviously, AMICI can't compute a steady state in this case, because the parameters are all NaN
• you get this error

What is exactly the setting of your optimization? Are you using pyPESTO or something alike? Do you have an objective function, which you have written yourself?

[paulstapor]

And: Do you really need the sensitivities, or would gradients be enough? You can run AMICI with steady-state data in adjoint mode. There, a fallback option exists, which does backward integration. This could save you from such an error.

[thomassligon]

Hi Paul, I don't know where the use of sensitivities is specified, or why, or if gradients is enough, or if adjoint mode is better. In order to answer your questions, do we need to look at the PETab .tsv files?

Here is the python code I'm using:

import amici
import pypesto
import pypesto.petab
import pypesto.optimize as optimize
import pypesto.visualize as visualize
import pypesto.engine as engine
import petab

import os
import numpy as np
import matplotlib.pyplot as plot

#folder_base = "../Models"
#model_name = "ASTHMA_V40_M08"
#petab_problem = petab.Problem.from_folder(folder_base + model_name)
petab_problem = petab.Problem.from_yaml('configuration.yaml')
testOut = petab_problem.get_optimization_to_simulation_parameter_mapping()

importer = pypesto.petab.PetabImporter(petab_problem)

model = importer.create_model()
print(model.getParameterScale())
print("Model parameters:", list(model.getParameterIds()), '\n')
print("Model const parameters:", list(model.getFixedParameterIds()), '\n')
print("Model outputs:   ", list(model.getObservableIds()), '\n')
print("Model states:    ", list(model.getStateIds()), '\n')

import libsbml
converter_config = libsbml.SBMLLocalParameterConverter().getDefaultProperties()
petab_problem.sbml_document.convert(converter_config)

obj = importer.create_objective(guess_steadystate=False)
# setting guess_steadystate=True might be interesting if all conservation laws are removed from the model.
# for some models, hyperparamters need to be adjusted
#obj.amici_solver.setMaxSteps(10000)
#obj.amici_solver.setRelativeTolerance(1e-7)
#obj.amici_solver.setAbsoluteTolerance(1e-7)

problem = importer.create_problem(obj)

#optimizer = pypesto.ScipyOptimizer()
optimizer = optimize.ScipyOptimizer()

# engine = pypesto.SingleCoreEngine()
#engine = pypesto.MultiProcessEngine()
engine = engine.SingleCoreEngine()

# do the optimization
#result = pypesto.minimize(problem=problem, optimizer=optimizer,
#                          n_starts=10, engine=engine)
result = optimize.minimize(problem=problem, optimizer=optimizer, n_starts=10, engine=engine)
result.optimize_result.get_for_key('fval')
ref = visualize.create_references(
    x=petab_problem.x_nominal_scaled, fval=obj(petab_problem.x_nominal_scaled))

visualize.waterfall(result, reference=ref, scale_y='lin')
visualize.parameters(result, reference=ref)

print('end of asthmaPE')

[paulstapor]

Hi,

sorry, just realized that this one is still open: just set obj.amici_solver.setSensitivityMethod(2) before optimizing, then you'll use adjoints.

[thomassligon]

Thanks! I'll definitely try that. But, first I need to finish testing in Copasi, updating the initial concentrations in SBML, and moving that to the parameters table in PETab.

[thomassligon]

I have spent a lot of time improving the model and testing in Copasi, but I am still getting exactly the same errors in AMICI. The errors say that we are getting NaNs in the first element of fxdot and p, and in the 57th element of w. Here are 2 lines of w:

#define s3035 w[56]
#define flux_r0 w[57]

This shows that the NaNs appear when w changes from species to fluxes. I don't have any NaNs in the SBML. It looks like AMICI is failing to compute one step, and then hands the NaNs to the next step, which fails with these errors (see original post above).

One other observation is that Copasi succeeds in finding what it wants after integrating to 10e+8 (seconds). When integrating manually in Copasi, the plots look good before 1e+5. I have attached the Copasi protocol file. Copasi-Steady-State-Protocol-2021-02-26.txt

[thomassligon]

@paulstapor I am now running petablint -y configuration.yaml -v The output is

Checking SBML model...
libSBML Warning (General SBML conformance): LibSBML expected to read the annotation into a ModelHistory object. Unfortunately, some attributes were not present or correct and the resulting ModelHistory object will not correctly produce the annotation.  This functionality will be improved in later versions of libSBML.
Reference: L2V4 Section 6.3
 An invalid ModelHistory element has been stored.

Checking measurement table...
Checking condition table...
Checking observable table...
Checking parameter table...
Traceback (most recent call last):
  File "c:\program files (x86)\microsoft visual studio\shared\python37_64\lib\runpy.py", line 193, in _run_module_as_main
    "__main__", mod_spec)
  File "c:\program files (x86)\microsoft visual studio\shared\python37_64\lib\runpy.py", line 85, in _run_code
    exec(code, run_globals)
  File "C:\Program Files (x86)\Microsoft Visual Studio\Shared\Python37_64\Scripts\petablint.exe\__main__.py", line 7, in <module>
  File "c:\program files (x86)\microsoft visual studio\shared\python37_64\lib\site-packages\petab\petablint.py", line 164, in main
    ret = petab.lint.lint_problem(problem)
  File "c:\program files (x86)\microsoft visual studio\shared\python37_64\lib\site-packages\petab\lint.py", line 720, in lint_problem
    problem.measurement_df, problem.condition_df)
  File "c:\program files (x86)\microsoft visual studio\shared\python37_64\lib\site-packages\petab\lint.py", line 208, in check_parameter_df
    assert_parameter_estimate_is_boolean(df)
  File "c:\program files (x86)\microsoft visual studio\shared\python37_64\lib\site-packages\petab\lint.py", line 524, in assert_parameter_estimate_is_boolean
    if int(estimate) not in [True, False]:
ValueError: cannot convert float NaN to integer

This is telling me that there is a problem in the parameters.tsv file, but which of the 573 parameters is it? It would be nice if the -v parameter were really verbose...

[paulstapor]

That shouldn't happen at all. That's not an actual output of petablint, but an error... Bot sure whether command line version and current python version do exactly the same, but would hope so... In case, retry with petab.lint_problem(petab_problem). Ideally check your PEtab version and update beforehand... simple pip install petab should do the trick.

[dweindl]

Bot sure whether command line version and current python version do exactly the same, but would hope so...

They do.

File "c:\program files (x86)\microsoft visual studio\shared\python37_64\lib\site-packages\petab\lint.py", line 524, in assert_parameter_estimate_is_boolean if int(estimate) not in [True, False]: ValueError: cannot convert float NaN to integer

Seems things are not converted to the right type. Would consider that a bug there. But easy for you to evade by just putting 0 or 1 for estimate in all rows instead of leaving some empty.

[thomassligon]

Very good. The last time I changed the model, I added some parameters, and the estimate value was left empty. Setting 1 explicitly fixes that. Now, petablint runs successfully, producing only a warning:

libSBML Warning (General SBML conformance): LibSBML expected to read the annotation into a ModelHistory object. Unfortunately, some attributes were not present or correct and the resulting ModelHistory object will not correctly produce the annotation.  This functionality will be improved in later versions of libSBML.
Reference: L2V4 Section 6.3
 An invalid ModelHistory element has been stored.

But I am still getting NaNs:

...
[Warning] AMICI:NaN: AMICI encountered a NaN value at index 0/286 in fxdot!
[Warning] AMICI:NaN: AMICI encountered a NaN value at index 0/696 in p!
[Warning] AMICI:NaN: AMICI encountered a NaN value at index 57/515 in w!
[Warning] AMICI:CVODES:CVode:OTHER: AMICI ERROR: in module CVODES in function CVode : The right-hand side routine failed at the first call.
[Warning] AMICI:simulation: AMICI simulation failed:
Steady state computation failed. First run of Newton solver failed: No convergence was achieved. Simulation to steady state failed. Second run of Newton solver failed: No convergence was achieved.
Error occurred in:
stacktrace not available on windows platforms
[Warning] AMICI:CVODES:CVode:OTHER: AMICI ERROR: in module CVODES in function CVode : At t = 2.19835e+199, the right-hand side routine failed in an unrecoverable manner. 
[Warning] AMICI:simulation: AMICI simulation failed:
Steady state computation failed. First run of Newton solver failed: No convergence was achieved. Simulation to steady state failed. Second run of Newton solver failed: No convergence was achieved.
Error occurred in:
stacktrace not available on windows platforms

It might be worth noting that Copasi needs multiple runs in order to succeed in calculating a steady state.

[paulstapor]

Could you just upload or send your current PEtab files including some code, that imports the model and tries to evaluate it?

[thomassligon]

Yes, look at https://github.com/LeonardSchmiester/AsthmaMapModel You will find AsthmaMapModel/PE.py, which contains the parameter estimation code, and AsthmaMapModel/configuration.yaml, which contains pointers to everything else that is needed.

[thomassligon]

Here are some observations of the problem in AMICI.

The first and second calls to AMICI run to a very large time (1.9e+199) and fail with the right-hand side routine failed in an unrecoverable manner The third call to AMICI has NaNs in the model. This means that the model itself is presumably free of NaNs; they are produced by AMICI.

<<call to amici.py, line 6070 return _amici.runAmiciSimulations(solver, edatas, model, failfast, num_threads)>>

[Warning] AMICI:CVODES:CVode:OTHER: AMICI ERROR: in module CVODES in function CVode : At t = 1.91183e+199, the right-hand side routine failed in an unrecoverable manner. 
[Warning] AMICI:simulation: AMICI simulation failed:
Steady state computation failed. First run of Newton solver failed: No convergence was achieved. Simulation to steady state failed. Second run of Newton solver failed: No convergence was achieved.
Error occurred in:
stacktrace not available on windows platforms

<<call to amici.py, line 6070 return _amici.runAmiciSimulations(solver, edatas, model, failfast, num_threads)>>

[Warning] AMICI:CVODES:CVode:OTHER: AMICI ERROR: in module CVODES in function CVode : At t = 1.91183e+199, the right-hand side routine failed in an unrecoverable manner. 
[Warning] AMICI:simulation: AMICI simulation failed:
Steady state computation failed. First run of Newton solver failed: No convergence was achieved. Simulation to steady state failed. Second run of Newton solver failed: No convergence was achieved.
Error occurred in:
stacktrace not available on windows platforms

<<call to amici.py, line 6070 return _amici.runAmiciSimulations(solver, edatas, model, failfast, num_threads)>>

[Warning] AMICI:NaN: AMICI encountered a NaN value at index 0/286 in fxdot!
[Warning] AMICI:NaN: AMICI encountered a NaN value at index 0/696 in p!
[Warning] AMICI:NaN: AMICI encountered a NaN value at index 57/515 in w!
[Warning] AMICI:CVODES:CVode:OTHER: AMICI ERROR: in module CVODES in function CVode : The right-hand side routine failed at the first call.
[Warning] AMICI:simulation: AMICI simulation failed:
Steady state computation failed. First run of Newton solver failed: No convergence was achieved. Simulation to steady state failed. Second run of Newton solver failed: No convergence was achieved.
Error occurred in:
stacktrace not available on windows platforms

[thomassligon]

After Paul helped me look at the problem, I tried some more debugging steps, and it looks to me like there might be a problem with get_fval and get_grad. Here are the details:

I can see One call to PyPESTO minimize (optimizer.py line 398) which calls scipy.optimize.minimize and x0 looks good. One call to ScalarFunction (scipy base.py line 261) where x0 looks good. Three or more calls to

• call_unprocessed (pyPESTO base.py line 146) where x_full is NaN on the third call
• amici.parameter_mapping.fill_in_parameters (amici_calculator.py line 85 and 94) where x_dot is NaN on the third call
• _amici.runAmiciSimulations (amici.py line 6070)

Multiple calls to get_fval (pyPESTO base.py line 254) which always returns inf Multiple calls to get_grad (pyPESTO base.py line 261) which always returns an array of NaNs

The call to get_grad that returns NaN occurs before the call to call_unprocessed with x_full = NaNs.

[thomassligon]

Here is the current status of this problem. In a discussion on the Copasi forum, Pedro Mendes found that steady state calculation in Copasi works well for this model if Target Criterion is set to Rate instead of Rate and Distance. In this case, Resolution can be left at 1e-9. In a discussion on pyPESTO, we discovered that AMICI forward integration fails after 1e+199 seconds. With that information, I tried a long integration with Copasi and found that it fails after 1e+19. Next, I can see from the successful Copasi steady state that two species require a long time to reach steady state (1e+19 and 1e+20 seconds). This looks like it might be a good hint for finding the integration error. However, these species are definitely not typical unbounded species, and the plots of particle numbers or concentrations do not reveal any problems. In summary, Copasi is able to ignore the failed integration and find a steady state, but AMICI produces NaNs, preventing the steady state calculation from succeeding. The failure in a long forward simulation points to a potential problem, but I haven't been able to find the cause for it yet.

[paulstapor]

Hi Tom, If AMICI fails after integrating until t = 1e199, it just means that the right hand side did not get sufficiently small until that point and that AMICI abandoned all hope that this could become better. Roughly 1e200 is the encoded somewhat final timepoint in AMICI, i.e. as "if things didn't converge until now, they just won't do at all". Alongsie with this error, you should get a warning like either "exceeded maximum number of steps" or "simulated to late time point without convergence of RHS".

My guess would be, that there's a species which indeeds grows unboundedly in a linear fashion.... (e.g., a species creation without a decay for all "children" of this species in a reaction network)

[dweindl]

Closing due to inactivity.