Converting CellML

See code in GitLab.
Author: Edward Vigmond edward.vigmond@u-bordeaux.fr & Moritz Linder moritz.linder@kit.edu

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

Note that variable names may be changed since EasyML demands global uniqueness while CellML only demands section uniqueness. Thus, in CellML, there may be two a variables, but in two different channels. They will be renamed to avoid the conflict in EasyML, and variables are grouped together by section after conversion.

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 containing error messages:

   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.

4. 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

2. 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.

3. 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);

4. 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.

Example 3

Another more advanced and complex model to convert is the Severi (a.k.a. Severi-Fantini-Charawi-DiFrancesco) model (severi_fantini_charawi_difrancesco_2012.cellml). The following steps will guide you through the process of model conversion.

1. Run the EasyML generator:

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

   You are already famillar with the inconsistencies of naming for the same quantities by different authors. Thus, we need to inspect the model file in order to figure out that the transmembrane voltage is called V_ode and there is no stimulus current applied. Rerun the converter with the adapted command:

   /home/openCARP/physics/limpet/src/python/cellml_converter.py --Vm=V_ode severi_fantini_charawi_difrancesco_2012.cellml > Severi.model

   For now the warning about the stimulus current can be ignored.

2. The next step is trying to run the C translator:

   /home/openCARP/physics/limpet/src/python/limpet_fe.py Severi.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 Severi

   The following error messages are received:

   Warning, you should only use alpha/beta OR tau/infinity while formulating gates!
   
   Using the tau infinity values.
   
   Error with dL
   
   Warning, you should only use alpha/beta OR tau/infinity while formulating gates!
   
   Using the tau infinity values.
   
   Error with n
   
   Warning, you should only use alpha/beta OR tau/infinity while formulating gates!
   
   Using the tau infinity values.
   
   Error with a
   
   Error, don't assign to variables marked as differential variables
   {Var("pi")}

3. Some of these error messages are similar to the ones of the earlier examples, so let us tackle them step by step. First you need to open the (Severi.model) file in a text editor of your choice. Starting with the gating variables dL, n and a, we have to be more careful this time. Typically, tau is defined as 1/(alpha+beta). However, by taking a closer look at the definition used by Severi et al. we can observe that tau_dL is specified differently. Locate tau_dL in around line 204. So, in this particular case, it is advised to rename alpha and beta while keeping tau and infinity, otherwise the dimension of tau_dL will not be correct. This can be done, for example, by adding the prefix i_CaL_dL_gate_ to alpha_dL and beta_dL. For the n and a gate this would not be necessary, but for consistency it is advised to proceed similar.

   Since pi is marked as a differential variable related to the i_Kr_pi_gate, we also need to refactor the constant pi into e.g. PI. Similarly, the names time and do cannot be used, so just rename them into t and var_do, respectively. Afterwards, we can save the file and rerun the C translator:

   /home/openCARP/physics/limpet/src/python/limpet_fe.py Severi.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 Severi

   We receive another error:

   KeyError: Var("Kc")

   By inspecting the supplementary material, we can see that somehow a mistake as occurred in the translation of the equations containing Kc. Refactoring Kc into Ko will sort that error out. Running the C translation again leaves us with 3 warnings about the deprecated *sprintf' function which we can ignore. Nevertheless, when you run the model and plot the transmembrane voltage no action potentials arise. This means that there is still something incorrect in this model file. We now reached a good point for recaptioning some theoretical background. Severi et al. modelled a sinoatrial node cell of a rabbit. The sinoatrial node is the origin of the heart beat and therefore action potentials should arise autonomously. Now a little bit more experience is needed to overcome the remaining issues.

4. In the cellML file, all time-based variables are implemented in seconds, but openCARP uses milliseconds to calculate the next step of the simulation. As there is no flag to set or other opportunity, the conversion from seconds into milliseconds has to be executed manually. The easiest way is to directly convert the constants defined in the beginning of the model file and use the converted values in the equations. Another difficulty is based on the unit C/mol of the constant F, which in si-units is equal to As/mol. F is defined approximately in line 8. We need to check the units of all summands of the equations containing F. For the variables diff_Ca_sub and diff_Nai, it appears to be necessary to multiply F by the factor 1000. In addition, the constant F is also used in a factor called RTONF, in which the seconds are reduced. For us, this results in no further time-related model adaptions.

5. Unforunately, we are still not finished. The last error results from the definition of the current Iion, which is located around line 130. Upon closer inspection, it is obvious that the units of Iion and i_tot are unequal. To overcome the last inconsistency, the unit of Iion should be changed to nA/uF.

6. Finally we can test our model. This time we will use a different option than in the previous examples. Put your .model file in your local physics/limpet/models folder and execute the following commands in root directory of openCARP:

   cmake -S. -B_build -DUPDATE_IMPLIST=ON 
   cmake --build _build   

   This will configure CMake with updated imp_list.txt and run the CMake building process to generate the .h and .cc files for your model in physics/limpet/src/imps_src.

7. The easiest way to check the correctness is plotting some variables of a simulation. Namely, this could be the transmembrane voltage which is V_ode as well as some currents and compare them to the graphs given in the paper. Coming back to the sinoatrial node cells and the fact that they excite without any stimulus current it is the final step to tell openCARP to simulate without a stimulus current. This can be realised by adding --stim-cur 0.0. If you keep the stimulus current the cell will arise action potentials shortly after the stimulus.

   Please don't forget to add metadata to your newly generated .model file. See the shipped .model files <https://git.opencarp.org/openCARP/openCARP/-/tree/master/physics/limpet/models> for examples.

Recent questions tagged easyml, cellml