See code in GitLab.
Author: Edward Vigmond edward.vigmond@ubordeaux-1.fr
This tutorial explains the basics of using the code generation tool limpet_fe.py
for generating ODE solver code for models of cellular dynamics. The code needs to be compiled, thus you need
to have a working developer environment in place (PETSc and other required libraries). The best
way to make sure this is the case and to avoid compatibility issues is to compile openCARP
from source as outlined here.
If you want to use additional models without recompiling openCARP itself each time,
make sure to enable dynamic library support when building openCARP. If you use CMake, this is done by
passing -DDLOPEN=ON
when calling CMake. If you use the makefile, set DLOPEN=1
in your
my_switches.def
file.
limpet_fe.py
is the EasyML to C translator for writing ionic models.
The easiest way to see how to write a model is to examine a simpler example.
To work on this tutorial do
cd ${TUTORIALS}/01_EP_single_cell/04_limpet_fe
Let's look at the Beeler-Reuter model (in your current folder: my_MBRDR.model) as modified by Drouhard and Roberge:
Iion; .nodal(); .external();
V; .nodal(); .external(Vm);
V; .lookup(-800, 800, 0.05);
Ca_i; .lookup(0.001, 30, 0.001); .units(uM);
V_init = -86.926861;
# sodium current
GNa = 15;
ENa = 40.0;
I_Na = GNa*m*m*m*h*(V-ENa);
I_Na *= sv->j;
a_m = ((V < 100)
? 0.9*(V+42.65)/(1.-exp(-0.22*(V+42.65)))
: 890.94379*exp(.0486479163*(V-100.))/
(1.+5.93962526*exp(.0486479163*(V-100.)))
);
b_m = ((V < -85)
? 1.437*exp(-.085*(V+39.75))
: 100./(1.+.48640816*exp(.2597503577*(V+85.)))
);
a_h = ((V>-90.)
? 0.1*exp(-.193*(V+79.65))
: .737097507-.1422598189*(V+90.)
);
b_h = 1.7/(1.+exp(-.095*(V+20.5)));
a_d = APDshorten*(0.095*exp(-0.01*(V-5.)))/
(exp(-0.072*(V-5.))+1.);
b_d = APDshorten*(0.07*exp(-0.017*(V+44.)))/
(exp(0.05*(V+44.))+1.) ;
a_f = APDshorten*(0.012*exp(-0.008*(V+28.)))/
(exp(0.15*(V+28.))+1.);
b_f = APDshorten*(0.0065*exp(-0.02*(V+30.)))/
(exp(-0.2*(V+30.))+1.);
a_X = ((V<400.)
? (0.0005*exp(0.083*(V+50.)))/
(exp(0.057*(V+50.))+1.)
: 151.7994692*exp(.06546786198*(V-400.))/
(1.+1.517994692*exp(.06546786198*(V-400.)))
);
b_X = (0.0013*exp(-0.06*(V+20.)))/(exp(-0.04*(V+20.))+1.);
xti = 0.8*(exp(0.04*(V+77.))-1.)/exp(0.04*(V+35.));
I_K = (( V != -23. )
? 0.35*(4.*(exp(0.04*(V+85.))-1.)/(exp(0.08*(V+53.))+
exp(0.04*(V+53.)))-
0.2*(V+23.)/expm1(-0.04*(V+23.)))
:
0.35*(4.*(exp(0.04*(V+85.))-1.)/(exp(0.08*(V+53.))+
exp(0.04*(V+53.))) + 0.2/0.04 )
);
# slow inward
Gsi = 0.09;
Esi = -82.3-13.0287*log(Ca_i/1.e6);
I_si = Gsi*d*f*(V-Esi);
I_X = X*xti;
Iion= I_Na+I_si+I_X+I_K;
Ca_i_init = 3.e-1;
diff_Ca_i = ((V<200.)
? (-1.e-1*I_si+0.07*1.e6*(1.e-7-Ca_i/1.e6))
: 0
);
group {
GNa;
Gsi;
APDshorten = 1;
} .param();
group {
I_Na;
I_si;
I_X;
I_K;
} .trace();
Generally, it looks like C code with a few extra commands. Let's break it down a bit:
Iion; .nodal(); .external();
V; .nodal(); .external(Vm);
These lines are very important. .external()
tells us that these are global variables and through these variables
we will interact with the simulation. The argument to external() is the actual name of the global variable if we want to call it something else locally. Here, the local variable V
is actually the global variable Vm
while Iion
is the
global name which we use locally. To alter the transmembrane voltage, Iion
must be assigned a value. openCARP then uses this value of Iion to adjust the transmembrane voltage.
You do not specify the change of Vm.
Note
Ionic models are fully functioning models. Plugins are additional components which must be added to an ionic model. They modify or expand behaviour, e.g., a new channel or a force genration description. For an ionic model, the ionic current is set, e.g., Iion = 2.3, while plug-ins add an additional current to an existing current, i.e., Iion = Iion + 2.3
The next two lines :
V; .lookup(-800, 800, 0.05);
Ca_i; .lookup(0.001, 30, 0.001); .units(uM);
tell the translator to use lookup tables for variables which are functions of these variables in order to avoid expensive function calls. The table range for V is -100 to 100 in steps of 0.05. This should be done after the model is verified to be working correctly.
There are several Hodgkin-Huxley type gating variables. Here is the f gate :
a_f = APDshorten*(0.012*exp(-0.008*(V+28.)))/
(exp(0.15*(V+28.))+1.);
b_f = APDshorten*(0.0065*exp(-0.02*(V+30.)))/
(exp(-0.2*(V+30.))+1.);
Variables of the form a_XXX and b_XXX are automatically recognized as the \(\alpha\) and \(\beta\) coefficients. The differential equation for the variable XXX will be automatically generated as well as an intial condition. Alternatively, one can assign tau_XXX and XXX_inf, which correspond to \(\tau_{XXX}\) and \(XXX_{\inf}\).
Note
if a variable is a differential variable but is NOT a recognized gating variable, its initial condition init_XXX must be explicitly set. That is why you will find Ca_i_init in the .model file.
The block:
group {
GNa;
Gsi;
APDshorten = 1;
} .param();
declares variables with default values. However, their values can be changed at the start of a simulation with the parameter changing syntax.
The section:
group {
I_Na;
I_si;
I_X;
I_K;
} .trace();
declares a group of trace variables. If a node is selected as being traced, these quantities will be output as functions of time at the node.
Units are a never ending source of grief. All sorts are used in various papers. openCARP expects certain standards so that interactions are possible between components. The table below gives the recommended units
quantity | units |
---|---|
membrane current | \(\mu A/cm^2 = pA/pF\) |
\(dX/dt\) | \(ms^{-1}\) |
\([Ca]_i\) | \(\mu M\) |
\([X]\) | \(mM\) |
voltage/potential | \(mV\) |
membrane conductance | \(mS/cm^2\) |
\(\tau_X\) | \(ms\) |
\(\alpha,\beta\) | \(ms^{-1}\) |
Internally, you can use what you like but you need to convert any variables passed between openCARP and your model, i.e., the total membrane current, concentrations, derivatives, and variables associated with gates.
Now, we need to make a runtime shared object library:
make_dynamic_model.sh my_MBRDR
make_dynamic_model.sh
(located in physics/limpet/src
) calls limpet_fe.py
(located in physics/limpet/src/python
) to convert the model file to C++ code for compilation.
Provided you have no errors, this will produce the files, my_MBRDR.cc and my_MBRDR.h.
make_dynamic_model.sh
will then call make_dynamic_model.py
,
which will produce my_MBRDR.so. This shared library can then be loaded by openCARP at runtime
and provide an ionic model with the name my_MBRDR
.
NOTE: The best way to make ensure a working development environment including thirdparty dependencies like PETSc and to avoid compatibility issues between your openCARP binary and the newly compiled model library is to compile openCARP from source as outlined here. To create a dynamic model, you need to have compiled openCARP from its source. The script needs the path of the source code as well as compiled libraries from openCARP.
The program bench is used to run single cell experiments. See bench for more details running this program. In this tutorial, run.py is just a wrapper for bench:
./run.py --load-module ./my_MBRDR.so
NOTE: For the --load option, an absolute path must be specified. That is why we used the "./" in "./my_MBRDR.so". If the --imp option is not specified, it will attempt to load the ionic model with the same name as the loaded library.
Remembering that we set APDshorten
as a parameter, we can see its effect:
./run.py --load-module ./my_MBRDR.so --imp-par='APDshorten=3'
As an exercise, try adding Ca_i
as a trace variable, or add a j gate to the sodium channel described by the following equations:
\(\alpha_j = 0.055\frac{exp(-0.25*(V+78))}{exp(-(V+78)/5)+1}\) and
\(\beta_j = \frac{0.3}{(exp(-0.1*(V+32))+1}\)
To use the ionic model in a tissue simulation, it is specified in the parameter file as :
num_external_imp = 1
external_imp = ./my_MBRDR.so
imp_region[0].im = my_MBRDR
Obviously, multiple IMPs can be added by increasing num_external_imp
. Make sure that the shared library was built in the same
environment (operating system, compiler etc.) that openCARP was built in.
© Copyright 2020 openCARP project Supported by DFG and EuroHPC Contact Imprint and data protection