Issue with cvode and a Ca dyanmics mod file

[arnaudon]

Hello,

I am trying to run a model based on this old model https://senselab.med.yale.edu/ModelDB/showmodel.cshtml?model=3670&file=/NTW_NEW/capump.mod#tabs-2 (from this paper: https://journals.physiology.org/doi/pdf/10.1152/jn.1994.72.2.803) but with cvode, it crashes. I could pin point it to the michaelis-mentens term in capump.mod, which results in negative cai, then with a log(cao/cai) in IT2.mod it creates nans. I wondered why cvode is not able to integrate this term, and if there could be another way to have these mod files compatible with cvode?

Also, if I remove this term, I see a drift in cai in my simulation (and they are very different from the ones with regular timesteps and this term), so there could be some sort of conservation issue as well related to this term?

Thank you for your help on that!

[nrnhines]

When I

git clone https://github.com/ModelDBRepository/3670
cd 3670
nrnivmodl
nrngui mosinit.hoc

and select the first menu item "Fig. 5a" and press the "RunControl/Init&Run" button I see a series of action potentials on two graphs. If I then select "NEURONMainMenu/Tools/VariableStepControl" and activate "Use variable dt" the immediate error message is

./x86_64/special: GABA cannot be used with CVODE

This, is due to gabaa.mod using SOLVE release but release is a PROCEDURE. Also this mod file implements network connectivity via POINTER. So it predates the existence of both CVode and NetCon.

But you did mention that your model is "based on this old model", and your problem is associated with capump.mod using explicit Michaelis-Menten formula drive_pump = -kt * cai / (cai + kd )

My experience is that newton iterations that try to solve odes that have states in the denominator are often unstable. The work around is to formulate as a KINETIC scheme. However, CVode should be able to solve this anyway. Have you decreased the error tolerance for cai. cvode.atolscale("cai", 1e-6) may help. (i.e. if the atol for cvode is it's default of 1e-3, then the scale factor demands that at each step cai be accurate to within 1e-9 mM (or 1 picomolar).

I'm happy to help further, but would need you to send me a zip file of your model with instructions on how to reproduce the error.

[arnaudon]

Thanks a lot for your help. Indeed, I am not using the circuit code in this repository, just single cell. To make it work, I had to set

cvode.atolscale("cai", 1e-8) 
cvode.atol(1e-8)

I tried several variants to have these at maximal values, and converged to this minimalistic setting. It is satisfactory for me now, as this was to try to reproduce some old results. I may come back later with more questions if we keep using this model and this is blocking.