Optimal solver setting and other ways to improve numerical stability of the models

[PavelBal]

@msevestre @Yuri05 @Huan-Yang @StephanSchaller @abdelr @VanessaBa

I am not sure where this discussion belongs to, so I start it here and we can transfer it later.

We should have some kind of a wiki for best practices how to parametrize the solver settings available in PK-Sim/MoBi. Furthermore, we can collect ideas which improvements we can implement, and how to make models more robust. Feel free to edit this post.

Solver settings

AbsTol

Default value is 1e-10

• What it does - ???
• General rule - decreasing the value improves the stability but increases computational time
• Experience: for many "simple" models, values between 1e-7 and 1e-9 are sufficient
• For complex models (e.g. QSP, DDI), decreasing the value up to 1e-12 may be required
• Values lower than 1e-12 often lead to less stability

RelTol

Default value is 1e-5

• What it does - ???
• General rule - decreasing the value improves the stability but increases computational time
• Experience: I sometimes decrease this value to 1e-7, but I do not really have a rule of thumb what to decrease - AbsTol or RelTol
• I heard that AbsTol and RelTol should not have values too far apart - but I do not really know what "too far apart" means and how to set these values

H0

Default value is 1e-10

• Initial step size (in minutes?)
• Experience: Never changed it

HMin

Default value is 0

• Minimal allowed step size
• Experience: Never changed it

HMax

Default value is 60

• Maximal allowed step size (minutes)
• Reducing this value might improve the stability of the model but increases the computational time. The value should be selected with regard to the expected kinetics of the model. E.g. for the diabetes model, a value of 5 minutes makes more sense, as we do not expect any simulation scenarios where the RHS remain unchanged for 60 minutes.

MxStep

Default value is 100000

• The default value is probably too high and leads sometimes to long computational times before the solver actually breaks and returns an error.

UseJacobian

Default is TRUE

• Are there some scenarios where setting this parameter to FALSE makes sense?

Modeling

Events

• In event start conditions, when using TIME as condition, use >= resp <= istead of = (or was it the other way around?)
• Avoid using conditions in equations of parameters or reactions. Instead, use events. Example: do not use eqauation Time > 1000 ? X : Y. Create an event at Time = 1000 and chage the value of the parameter from X to Y.

[msevestre]

Great initiative
I would suggest to put in in the docs once finalized

[PavelBal]

Solver implementation

• Add sparse solver implementation #1 - what is this issue about?
• @Huan-Yang you were talking about CVODES providing two different solvers - for stiff and non-stiff models, and that we should allow selection of the solver... Could you elaborate a bit more?
• Maybe this can be of some help: https://computing.llnl.gov/projects/sundials/cvode

Scale factors

• 2DO: Benchmark current implementation in MoBi
• Add implementation in R?

[Huan-Yang]

Hi @PavelBal ,
Regarding the CVODES/CVODE, they are multiple solvers. Two often used ones are bdf (for stiff) and adams (for non-stiff). Stiff case is hard to integrate than non-stiff case. But bdf takes a bit longer integration time than adams. In practice, for a new set of ODEs, I would start with bdf solver first. If the simulation looks fine, check the simulation from adams.
Somehow, the choice of solvers would depend on the ultimate goal of the ODE-model: if it is only for one-time simulation. I prefer to choose bdf. If it is related to optimization/inference, some accuracy-efficiency trade-off needs to consider. I believe that can vary in a case-by-case manner and depends on available approaches for optimization/inference.

Besides the concern of numerical solver, I believe the failure of 'ideal' solver can result from purely mathematical expression of the right-hand side (RHS).
Please check this. Indeed wrong choice of a RHS can cause later numerically unstable simulation, which is independent of the choice of numerical solvers.

[Yuri05]

https://computing.llnl.gov/projects/sundials/faq#cvode

https://sundials.readthedocs.io/en/latest/cvodes/Mathematics_link.html

[wilbertdew]

Hi All,

I think it would be good to have some more details on this in the documentation and some suggestions for optimization.
Especially if you want to speed up a batch run of simulations, this can also make a big difference for computation time in my experience.
I found this explanation helpful.
It also points to the importance of the relative and absolute tolerance ratio, as this decides at which point the absolute tolerance becomes dominant.
In general, I think reducing the relative tolerance is mainly relevant for models which are highly nonlinear and contain rapidly changing values, e.g., target binding at a high and saturating dose.
Reducing the absolute tolerance mainly makes sense if you simulate values below 1e-10 (in base units I assume) which are still relevant and need to be precise, or if you reduce the relative tolerance.
Reducing Hmax might be mainly relevant for problems that include tipping points rather than a gradual approach to steady-state. If there is a long phase (relative to 60 minutes) with rather linear behavior, increasing Hmax can speed up the simulations significantly without losing precision.