Converting CellML

See code in GitLab.
Author: Edward Vigmond edward.vigmond@u-bordeaux.fr

Converting from CellML

Intro

CellML is a markup language for describing, amongst other things, ionic models. Many models are published and available from the CellML Model Repository. However, this format can only be processed by a machine. openCARP uses its own markup language called EasyML.

The Script

Conversion of a CellML file is performed using the cellml_converter.py script:

usage: cellml_converter.py [-h] [--Vm VM] [--Istim ISTIM] [--out OUT]
                           [--trace-all] [--plugin] [--params]
                           [--keep-redundant]
                           CMLfile

Convert CellML 1.0 to model format

positional arguments:
  CMLfile           CellML file

optional arguments:
  -h, --help        show this help message and exit
  --Vm VM           transmembrane voltage variable [V]
  --Istim ISTIM     stimulus current variable [i_Stim]
  --out OUT         output file
  --trace-all       trace all variables
  --plugin          is a plug-in, not a whole cell model
  --params          parameterize all constants
  --keep-redundant  keep redundant differential expressions, eg alpha, beta
                    and diff)

The flags mean the following

The process of creating an EasyML file is not entirely automated. Manual intervention is always necessary. Generally, the process to create a model file is as follows:

1. cellml_converter.py is run on the CellML file
2. references to stimuli and extraneous variables are removed
3. redundant gating variable equations are removed
4. units are adjusted
5. integration methods are adjusted
6. lookup tables are implemented

Example 1

We will begin by converting a relatively simple model, the LuoRudy II model, which is shipped with this example in the file LuoRudy94.cellml. The python scripts in the following commands need to be called with their full path if they are not in your $PATH environment variable. They are shipped with the openCARP source code and located in the subfolder physics/limpet/src/python. To compile new ionic models, we highly recommend to build openCARP from source as described in the installation instructions on the webpage. You'll need to enable dynamic model support during the build process to follow this example. In the following, we'll assume that the openCARP source code is located in /home/openCARP but this can be easily adjusted according to your local installation.

1. run the converter to produce the EasyML model file:

   /home/openCARP/physics/limpet/src/python/cellml_converter.py LuoRudy94.cellml > LuoRudy94.model

2. run the model to C translator:

   /home/openCARP/physics/limpet/src/python/limpet_fe.py LuoRudy94.model /home/openCARP/physics/limpet/models/imp_list.txt /home/openCARP/physics/limpet/src/imps_src
   /home/openCARP/physics/limpet/src/make_dynamic_model.sh LuoRudy94

   You will see the following output which has errors:

   Syntax error at line 79, token '{'
   Warning, you should only use alpha/beta OR tau/infinity while formulating gates!
   
   Using the tau infinity values.
   
   Error with d
   
   Warning, you should only use alpha/beta OR tau/infinity while formulating gates!
   
   Using the tau infinity values.
   
   Error with f
   
   Traceback (most recent call last):
   File "/home/openCARP/physics/limpet/src/python/limpet_fe.py", line 1581, in <module>
       obj.printSource()
    File "/home/openCARP/physics/limpet/src/python/limpet_fe.py", line 786, in printSource
       declare_type = "double ",
   File "/home/openCARP/physics/limpet/src/python/limpet_fe.py", line 201, in printEquations
       for var in eqnmap.getList(target, valid):
   File "/home/openCARP/physics/limpet/src/python/im.py", line 220, in getList
    for d in self[sym].depend:
   File "/home/openCARP/physics/limpet/src/python/basic.py", line 87, in __getitem__
       return self._map[key]
   KeyError: Var("I_st")

3. Open up LuoRudy94.model in the text editor of your choice. Based on the error output, look at line 79:

   dV_dt = ((I_st - (i_Na+i_Ca_L+i_K+i_K1+i_Kp+i_NaCa+i_p_Ca+i_Na_b+i_Ca_b+i_NaK+i_ns_Ca))/Cm); .units(mV/ms);
   I_st = ((({http://www.w3.org/1998/Math/MathML}rem:time,stimPeriod)<stimDuration) ? stimCurrent : 0.0); .units(uA/mm^2);

   It appears that I_st is the stimulus current which we would like to remove. We can manually remove it by deleting line 79 and eliminating it from line 78, or we can run /home/openCARP/physics/limpet/src/python/cellml_converter.py --Istim=I_st LuoRudy94.cellml > LuoRudy94.model

   Now we are left with errors with gating variables d and f. Locate them around line 106 in the LuoRudy94.model file. Note how alpha, beta, tau and _inf are defined for all. By inspecting, we see that alpha_d and alpha_f are defined using the infinity and tau values so the alphas and betas can be removed. Save the file, and rerun /home/openCARP/physics/limpet/src/python/limpet_fe.py LuoRudy94.model /home/openCARP/physics/limpet/models/imp_list.txt /home/openCARP/physics/limpet/src/imps_src which should now have no errors.

# make the linkable library and run it to be sure:

/home/openCARP/physics/limpet/src/make_dynamic_model.sh LuoRudy94
bench --load ./LuoRudy94.so

Example 2

Now, we will compile the Grandi (a.k.a. Grandi-Pasqualini-Bers) model (Grandi.cellml), which is quite complicated.

1. Run the EasyML generator:

   /home/openCARP/physics/limpet/src/python/cellml_converter.py Grandi.cellml > Grandi.model
   WARNING: stimulus current not found
   WARNING: transmembrane voltage not found

   So, not every author uses the same names for the same quantities. Inspecting the model file, we can see that transmembrane voltage is called V_m and the stimulating current is I_app. Rerun the converter:

   /home/openCARP/physics/limpet/src/python/cellml_converter.py --Vm=V_m --Istim=I_app  Grandi.cellml > Grandi.model

   Generate the C code:

   /home/openCARP/physics/limpet/src/python/limpet_fe.py Grandi.model /tmp/openCARP/physics/limpet/models/imp_list.txt /tmp/openCARP/physics/limpet/src/imps_src

   which does not complain but when we try to run the model :

   bench --load ./Grandi.so

   we get NaNs which is not good. So, at this point we need to rely on intuition and a good debugger. First we reduce the time step and hope it is just an integration issue with stiff equations:

   bench --load ./Grandi.so --dt=0.005

   which luckily works. We do not want to be restricted to a 5 microsecond timestep so we need to investigate further. The most common issues are

   1. calcium is released explosively during the upstroke so we often need to proper integrate calcium related variables as a group with a more complicated method.
   2. gating variables are not defined with alpha/beta or tau/_inf so they can go out of bounds, especially m of the sodium channel
   3. Markov submodels may also be stiff systems so need the markov_be integration method.

   In this particular case, the sodium channel, I_Na on line 170, is definitely a culprit. After inspection, we see diff_ defined as well as redundant alpha/beta/tau/_inf. This is what I finally changed it to:

   #I_Na
   I_Na = (I_Na_junc+I_Na_sl); .units(uA/uF);
   I_Na_junc = (Fjunc*GNa*(m*m*m)*h*j*(V_m - ena_junc)); .units(uA/uF);
   I_Na_sl = (Fsl*GNa*(m*m*m)*h*j*(V_m - ena_sl)); .units(uA/uF);
   ah = ((V_m>=-40.) ? 0. : (0.057*exp((-(V_m+80.)/6.8)))); .units(unitless);
   tau_h = (1./(ah+bh)); .units(ms);
   bh = ((V_m>=-40.) ? (0.77/(0.13*(1.+exp((-(V_m+10.66)/11.1))))) : ((2.7*exp((0.079*V_m)))+(3.1e5*exp((0.3485*V_m))))); .units(unitless);
   h_inf = (1./((1.+exp(((V_m+71.55)/7.43)))*(1.+exp(((V_m+71.55)/7.43))))); .units(unitless);
   tau_m = ((0.1292*exp(-(((V_m+45.79)/15.54)*((V_m+45.79)/15.54))))+(0.06487*exp(-(((V_m - 4.823)/51.12)*((V_m - 4.823)/51.12))))); .units(ms);
   m_inf = (1./((1.+exp((-(56.86+V_m)/9.03)))*(1.+exp((-(56.86+V_m)/9.03))))); .units(unitless);
   bj = ((V_m>=-40.) ? ((0.6*exp((0.057*V_m)))/(1.+exp((-0.1*(V_m+32.))))) : ((0.02424*exp((-0.01052*V_m)))/(1.+exp((-0.1378*(V_m+40.14)))))); .units(unitless);
   aj = ((V_m>=-40.) ? 0. : ((((-2.5428e4*exp((0.2444*V_m))) - (6.948e-6*exp((-0.04391*V_m))))*(V_m+37.78))/(1.+exp((0.311*(V_m+79.23)))))); .units(unitless);
   tau_j = (1./(aj+bj)); .units(ms);
   j_inf = (1./((1.+exp(((V_m+71.55)/7.43)))*(1.+exp(((V_m+71.55)/7.43))))); .units(unitless);

   We try to run it again and we still have an issue. As mentioned above, Markov chains and calcium handling can be numercially stiff. We will try more expensive integration methods at the cost of some speed. I added the following code to Grandi.model :

   group {
       Ca_i;
       Ca_j;
       Ca_sl;
       Ca_sr;
       f_Ca_Bj;
       f_Ca_Bsl;
   }.method(cvode);
   group {
       Ry_Ri;
       Ry_Ro;
       Ry_Rr;
   }.method(markov_be);

   This time, we can compile and run it. However, we are not done. Let's try to improve performance by using a lookup table for the voltage functions. We change line 2 of the model file to be :

   V_m; .external(Vm); .nodal();.lookup(-100,100,0.05);

   We generate the C file, make the linkable library and when we try to run it:

   bench --load ./Grandi.so --stim-start=100 --dt=0.005 -Oa
   
   L2 : LUT WARNING: Grandi V_m=-5 produces -nan in entry number 1!
   
   L2 : LUT WARNING: Grandi V_m=-5 produces -nan in entry number 2!
   
   Outputting the following quantities at each time:
         Time              Vm            Iion
   
   [CVODE ERROR]  CVode
     At t = 283.5 and h = 1.90735e-08, the corrector convergence test failed repeatedly or with |h| = hmin.
   
   bench: /home/openCARP/physics/limpet/src/imps_src/Grandi.cc:868: Grandi: Assertion `CVODE_flag == 0' failed.

   Hmmm. Maybe the NaN in the lookup table triggers the error. Sure enough, we see that entry for the d-gate time constant of I_Ca, tau_d, has a singularity at -5 mV and that is where the simulation failed. We need to use L'Hopital's rule for that voltage:

   tau_d = (V_m==-5.)?(-d_inf/6./0.035):((1.*d_inf*(1. - exp((-(V_m+5.)/6.))))/(0.035*(V_m+5.))); .units(ms);
   d_inf = (1./(1.+exp((-(V_m+5.)/6.)))); .units(unitless);

   This model file </downloads/Grandi.good.model> now works!!!! For further practice, determine which subset of Ca variables can be removed from the group and/or if another integration method can be used.

   Also don't forget to add metadata to your newly generated .model file. See the shipped .model files. for examples.