Intracellular Na and K concentrations do not reach steady state with the O'Hara model

Question

Hi,

First, I implemented the O'Hara model myself in Python (adapted from the original Matlab code). When I run it for endocardial cells with bcl=1000 the intracellular Na and K concentrations reach steady state within 500 beats (they stabilize around 144.6 mM / 7.33 mM).

However, when I use the O'Hara model in bench the intracellular Na and K concentrations do not reach steady state but still decrease even within 1000 beats. It also happens when I use the exact same initial conditions as in Python.

I believe that it might be due to the timestep. I tried to decrease dt in bench but when I try to decrease it below dt=1e-3 the simulation crashes. Do you have an idea what I could do to reach steady state?

Thanks in advance!

Best,

Sophia

Comments

Hi,

I think I found the problem!

When I converted the .model file to a .mmt file and did the simulation with Myokit I observed the same unstable behavior of the Na and K concentrations as in OpenCARP.
However, when I manually added the stimulus current to the equation for diff_Ki the concentrations stabilized.

So I think the problem is that in bench the stimulus current is not added properly to the diff_Ki equation.

I didn't expect this, since in my command file I used ['--stim-assign', 'on',  '--stim-species', 'Ki'] so I though the stimulus would affect the change of Ki (and I checked it, the variable name is indeed Ki).

Do you have an idea why it's still not working?

Best,
Sophia

---

This seems to be a bug. Let's move the discussion to the issue tracker: https://git.opencarp.org/openCARP/openCARP/-/issues/338#note_6452

---

Hi,

First of all thank you a lot for addressing this problem in the latest commit!
Unfortunately, even when loading the latest version of opencarp and re-converting the model using make_dynamic_model the ion concentrations still don't reach steady state for me.
I didn't specify the surface_to_volume command and just used the default value. Do you think I should use another value?

Best,
Sophia

---

There is virtually no way of figuring out the cell surface to cell volume ratio automatically. The new "--surface-to-volume" option allows you to set it manually though.
The default is 0.14, to get the same behavior as adding it to difF_Ki directly, you should probably set it to 0.0594 if I see correctly.

---

I set it to 0.0594 now but unfortunately the behavior still does not change..

---

Probably best to look at the values computed here https://git.opencarp.org/openCARP/openCARP/-/blob/master/physics/limpet/src/MULTI_ION_IF.cc?ref_type=heads#L1763
in a debugger or with additional print statements and compare them to the ones you would add directly to diff_Ki

---

Hi :) so when we recalculated the value we got a result of 0.326, but with this value it still did not work. I did it brute force now trying different values and somehow, a value around 635 seems to yield the most stable results (Nai stays stable over 10,000 s but Ki starts to slightly increase again after 1,000 s). For values <=600 Ki decreases all the time and for values >= 635 Ki first decreases but then increases again. Do you have an idea why the surface to volume value seems to be in that range?

Answer 1

Hi,

Although there aren't any limitations stated in the manual about dt, there are some about spacedt, timedt, and several others. For instance spacedt can only be as small as 'dt/1000'. Perhaps, that is the issue. I recommend double checking these values and their restrictions. Kindly refer to 24.19 in openCARP manual for them.

Best,

Kocak

Comments

Okay thanks! And do you have other ideas except for dt that I could try to bring the intracellular ion concentrations to steady state?

---

Unfortunately, no. I bet more experienced members will have an opinion, though.

Answer 2

Interesting phenomenon indeed.

dt is given in microseconds with a default value of 5. Decreasing it to 1e-3 would mean to integrate the ODEs with a time step of 1ns, which seems unnecessarily small to me.

To dig further into this issue, you could

- compare the currents between your Python implementation and the openCARP one to identify differences

- check for the numerical resolution of the state variables. Are all implemented as double or are floats used in any of the implementations?

- check with a third tool (e.g. MyoKit or OpenCor as they can load the CellML files directly)

Comments

Another thing that comes to my mind:
As mentioned in the comments in the metadata of the openCARP O'Hara et al. model (https://git.opencarp.org/openCARP/openCARP/-/blob/master/physics/limpet/models/OHara.model#L7), the sodium current was replaced with the formulation proposed by ten Tusscher et al. as done in Dutta et al. (https://doi.org/10.1016/j.pbiomolbio.2017.02.007) and suggested by O'Hara et al. to make it suitable for tissue level simulations (https://journals.plos.org/ploscompbiol/article/comment?id=10.1371/annotation/c9345189-f5f9-46c0-a4cf-d774cd858917) . This might also affect the sodium / potassium balance and thus the existence of a stable limit cycle.

---

Thanks a lot for your answer!

To your first comment: dt=1e-3ms means integration of the ODEs with a timestep of 1us not ns, or did I misunderstand something? Because in my python code I get problems (division by zero error) when I use a minimal timestep >10ns, so I still think that this might be the problem..

To your last comment: In the OHara.model code it says

TT2INa = 0;
ORdINa = 1;
INa_Type = ORdINa; .flag(TT2INa, ORdINa);

So as far as I understand the OHara formulation is the default, right? Because then the formulation is identical to my python code.

I will compare the currents and state variables now and check out other tools as you suggested! :)

---

Yes, you are right!

For whatever reason dt in in openCARP is defined in microseconds while it's defined in milliseconds on bench

---

So I compared the currents and mainly INa, Ito and IK differ.
However, the formulations are exactly the same so I think I need to try another method to solve the ODEs.
I also did the simulation with openCOR and here the intracellular Na and K concentrations also reach steady state within 500 beats. They use the cvode solver with a maximum timestep of 0.1ms.
But when I try to use this method for Nai and Ki I get the warning:

WARNING: Model OHara_test: MLIR code generation does not support CVODE integration method. Using original codegen...

And if I try to use the model anyways then the simulation does not start.. Do you know why cvode does not seem to be supported? (I'm using opencarp version 12.0)

---

I guess you compiled openCARP yourself from source. If you disable MLIR code generation (setting ENABLE_MLIR_CODEGEN to false in CMake) and have CVODE on your system, it should work

---

Hi Axel, we compiled without MLIR. However, `make_dynamic_model.sh` (seemingly) relies on `my_switches.def` for the CVODE path and this file is not generated/installed when using CMake. This is the case on v12.0, as in the newer versions limpet_fe doesn't work for us anyway...

EDIT: Generating the file manually does not fix it though ... must be a deeper rooted problem, but I did not find anything so far. sundials/cvode is installed (with headers) in standard paths. MLIR is disabled. my_switches.def has CVODE=1.

---

I've never used CVODE for solving the IMP myself. Might take some time until someone can look into it as many of us are still on vacation

---

Work in progress: https://git.opencarp.org/openCARP/openCARP/-/issues/330