CEPS
Using CEPS: Input files

Once CEPS executables are generated, you can run them from command line:

$> myCepsApplication myInputFile

CEPS uses a text file as primary input of the program. This file sets most of the simulation parameters. Input files commonly have the .in extension, but it can be any or none.

After a word about the syntax of this file, you will find information on parameters that are used for cardiac simulations. Several input files can be found in the example directory. There should be one for each application.




Syntax of main input file

Parameters are given by lines containing the colon ':' character:

a keyword : a value

Keywords are used by the application to identify parameters. If any mandatory keyword is missing, the program prints an error message about the missing key and stops. If the parameter has a built-in default value, a warning message is displayed when run with option -w.

Keywords are not case sensitive, and all white spaces are removed from them.

Values can be given on several lines if the backslash '\' character is used:

another keyword : a value \
on \
several lines

Several entries can be given for the same keyword. In this case a single parameter is created by concatenating all the entries related to the same keyword, using commas ',' as separators of the new parameter string:

key : set of params 1, set of params 2
# is equivalent to
key : set of params 1
key : set of params 2




Geometry

3d mesh: <meshName1> [,<meshName2>,...]
2d mesh: <meshName1> [,<meshName2>,...]
1d mesh: <meshName1> [,<meshName2>,...]

Specify meshes that constitute the computational domain with these entries, respectively for volumic, surfacic and linear elements of the geometry.

Names can be given with or without extensions. If no extension is provided, CEPS will then check for files

  • meshName.mesh (medit,gmsh,...) in 3-dimension format
  • meshName.ele (Tetgen)
  • meshName.pvtu (parallel XML VTK file, if installed)
  • meshName.vtk (legacy VTK file, if installed) We strongly suggest to use extensions for all file names.

The expected unit for coordinates is centimeter. See the GeometryScale parameter for easy unit change.

A note about fibers: in order to pass fibers orientation to CEPS cardiac applications, it is possible to do it in two ways:

  • if the mesh is provided in a vtk file, fibers can be given in the same file, in 3 vectorial cell-data arrays named "fibers" for longitudinal direction, and "fibers_transverse" and "fibers_normal" for other directions.
  • for other mesh formats, in a separate .msh file (not .mesh) where orientations are given element per element (one line per element) eight lines after the keyword "$ElementData". Each line must have at least 3 scalars for main direction, then optionally 3 for transverse and 3 more for normal direction. An example of these files can be found in the example directory of the project.
coupled nodes: <fileName>

Give a list of coupled nodes to make links between two meshes. This is particularly useful for the bilayer model. To work properly, the file, given with its full name, must contain the following pattern for each pair of coupled meshes:

COUPLED_NODES <meshA> <meshB> <N>
<nodeA1> <nodeB1>
<nodeA2> <nodeB2>
...
<nodeAN> <nodeBN>

where meshA and meshB are the base names of the two linked meshes, as given in the main input file, N is the number of couples to be read and nodeABX are the respective IDs of the coupled nodes in those meshes.

geometry scale: <scalar>

Multiplies all the mesh points coordinates by this scalar.

partitioning method: <method>

How to distribute the elements amongst processors for parallel runs. If CEPS is installed with ParMetis, method can be

  • ParmetisElem: element based partitioning,
  • ParmetisElemGeom: faster element based partitioning,

If CEPS is installed with PTScotch, method can also be

  • PTScotchNode: node partitioning using primal graph




Tissue parameters

Am: <options>
Cm: <options>

$A_\mathrm{m}$ is the surface of cell membrane per $\si{\centi\metre\cubed}$, expressed in $\si{\per\centi\metre}$. $C_\mathrm{m}$ is the average membrane capacitance, homogeneous through the whole tissue, in $\si{\micro\farad\per\centi\metre\squared}$.

The options string must start with either:

  • CONSTANT , then a single scalar must be given. For example: CONSTANT 1.0 sets the value to 1. everywhere in the domain.
  • PIECEWISE , then pairs of attributes and scalar values must be given. The attribute -1 can be used to give default values. For example, if the mesh has regions of attributes 3, 4, 5 and 6: PIECEWISE 5 2.0 6 1.0 -1 0.5 will set the coefficient to 2. in region with attribute 5, 1. in region with attribute 6, and the value will be 0.5 in regions with attributes 3 or 4.
  • VTK , then the name of a vtk file containing the data must be given, as well as the name of the array in the file. With this option, you can pass a coefficient in the same file as the mesh. Data must be given as "point data", not cell data.

See the units page for more information about physical quantities. Note to developpers: the syntax of these coefficients is that of ScalarPhysicalCoefficient, and can be used for others coefficients as well.

extracellular conductivity: <region1> <g_l1> <g_t1> <g_n1> \
[, <region2> <g_l2> <g_t2> <g_n2>, ...]
intracellular conductivity: ...
# or
extracellular conductivity: <region1> <g_l1> <g_t1> <g_n1>
extracellular conductivity: <region2> <g_l2> <g_t2> <g_n2>
intracellular conductivity: ...
# ...

Specify the tensorial conductivities of tissue regions, in $\si{\milli\siemens\per\centi\metre}$ (cf units). Several conductivities can be given using comma-separated substrings, or multiple entries.

Regions are specified by their mesh attribute. If conductivity is to apply everywhere, use -1 as region. In that case, it is still possible to give a region specific value with another attribute (no order in attributes required).

The scalar values are the terms of the diagonalized conductivity tensors, locally oriented along the fibers direction.

For example the line

extracellular conductivity: 1 0.5 0.1 0.1, -1 1 0.5 0.5

sets $\sigma_{\mathrm{ext}}(x) = \mathrm{diag}(1,0.5,0.5)_{(l,t,n)(x)}$ everywhere except in region 1 where $\sigma_{\mathrm{ext}} = \mathrm{diag}(0.5,0.1,0.1)_{(l,t,n)(x)}$.

Important note** Even when using the monodomain model, where both conductivities can be combined in the case of equal anisotropy ratio, both conductivities must be precised. The harmonic average of the conductivities at each node is then used.

volume fraction file: <file>
volume fracion array name: <name> (default: volumeFraction)

If conductivity is altered by volume fraction, the value of this volume fraction can be given as data defined on each mesh cell. Input file can be .msh or .vtk. For vtk files, the data array named after the above parameter can be in the same file as all other data.

volume fraction maps intracellular: attr1 file1 component1, attr2 file2 component1, attr1 file3 component2, ...
volume fraction maps extracellular: ...

If conductivities are altered by volume fraction, interpolation maps must be provided in order to evaluate the modified conductivities at a given point with a given volume fraction. Each parameter substring (between commas) consists in the region attribute where the interpolation map is applied (-1 for everywhere), the file name where the map is stored and finally the component on which the conductivity applies.

Interpolation maps are stored in separate files. Each file is made of 7 columns: the volume fraction, then the 6 values of the modified symmetrical conductivity tensor ( $g_{11}$, $g_{12}$, $g_{13}$, $g_{22}$, $g_{23}$, $g_{33}$).




Ionic model

ionic model: <name> <region> <component> <ODEsolver> \
[, <name2> <region2> <component2> <ODEsolver2>, ...]

Ionic model and how to numerically solve it.

  • name must be one of the following (old CEPS keywords are still valid):
    • TTP06::Endo
    • TTP06::MidMyo
    • TTP06::Epi
    • MS03
    • CRN98
    • BR77
  • region is the integer mesh attribute of where to apply the model, several ionic models can be given for different regions.
  • component is the index of the problem unknown concerned by the ionic model (e.g. extra or intra cellular potential for the bidomain problem).

  • ODEsolver must be one of the following:

    • EXPLICIT_EULER (or just EULER)
    • EULER_TAU_WINF (or abbr. ETW)
    • IMPLICIT_EULER (or FBE)
    • RUSH_LARSEN 1-4 (or RL 1-4)
    • EXP_ADAMS_BASHFORTH 1-4 (or EAB 1-4)
    • BACKWARD_DIFFERENTIATION 1-4 (or SBDF 1-4)
    • RUNGE_KUTTA opt (or RK opt) where opt can be an integer between 1 and 5 or Heun2 Heun3 Ralston3 SSP3 38 Ralston4 Nystrom5 Fehlberg5.

    Refer to the numerical methods page for details about the solvers.

ode time step: <small scalar>

The time step used by ODE solvers, in milliseconds. It must be lower than the PDE time step.

Options that are specific to a given ionic model**

MS normalization voltage min: <scalar> (default -90 mV)
MS normalization voltage max: <scalar> (default 50 mV)
MS tau in : <scalar> (default 0.3 ms)
MS tau out : <scalar> (default 6 ms)
MS tau open : <scalar> (default 130 ms)
MS tau close : <scalar> (default 150 ms)
MS vGate : <scalar> (default 0.13 )

Parameters of the Mitchell-Schaeffer 2003 model can be tuned with these options. The normalization voltages can be used to keep the voltage unit in mV or adimensionned (when the values are 0 and 1).

CRN cell surface : <scalar> (default 5.439e-5)
CRN gNaFactor : <scalar> (default 1.)
CRN gtoFactor : <scalar> (default 1.)
CRN gKrFactor : <scalar> (default 1.)
CRN gCaLFactor : <scalar> (default 1.)

Tuning parameters for the Courtemanche Ramirez Nattel 1998 model, which is normalized by surface. The cell surface parameter allows to match CEPS units.




Stimulation

Sets how the tissue is excited. You can defined two global type of stimulations :

  • volumic stimulation to perfom stimulation protocol inside a volume (or surface in 2D)
  • boundary stimulation to perform stimulation protocol on the boundary of the domain
  • pacer stimulation to perform stimulation protocol on boundary with 2 electrodes and a total intensity

All stimulation controls are defined with three type of input file entry

time function : <name> <time_key_1> <time_value_1> <time_key_2> <time_value_2s> ...
space function : <name> <space_key_1> <space_value_1> <space_key_2> <space_value_2s> ...
<KIND> stimulation : <name> <stim_key_1> <stim_value_1> <stim_key_2> <stim_value_2> ...

where <KIND> is volumic, boundary or pacer. Each key has default value and are optional.

  • time function control keys
    • PROFILE : select regularity of stimulation C4, C5, CINF or CONSTANT with default value C4.
    • START : time at wich stimulation starts in milliseconds, default value is 0.0.
    • END : time at wich stimulation ends in milliseconds, default value is 5000.0.
    • DURATION : in milliseconds, default value is 2.0 milliseconds.
    • PERIOD : in milliseconds, use none or -1 for non-periodic stimulations, default value is -1.
  • space function control keys
    • PROFILE : select regularity of stimulation C4, C5, CINF or CONSTANT with default value C4.
    • CENTER : center of stimulated region, three scalars with default values 0.0 0.0 0.0.
    • DIAMETER : width of stimulated region with default value 0.3.
  • volumic/boundary stimulation keys:
    • TIME : name of a time function or internal tag
    • SPACE : name of a space function or internal tag
    • AMPLITUDE : scale factor to rescale stimulation between zero and amplitude
    • COMPONENT : excited unknown of the problem (intra or extracellular potential ?)
    • ATTRIBUTE : mesh region attribute on which to apply stimulation (careful, if you set also type or origin or diameter the behaviour can not be ensured !)
  • pacer stimulation keys:
    • TIME : name of a time function or internal tag
    • ANODE : name of a space function or internal tag defining anode electrode
    • CATHODE : name of a space function or internal tag defining cathode electrode
    • INTENSITY : the current flowing through the electrodes
    • ATTRIBUTE : mesh region attribute on which to apply stimulation (careful, if you set also type or origin or diameter the behaviour can not be ensured !)

The profiles of the regular functions are the following, plotted as functions of the distance to center, normalized by the radius:

\[ \begin{tikzpicture}[x=6cm,y=5cm] \foreach \x in {0,0.5,1} { \draw (\x,0) node [below]{\x} -- +(0,5pt); } \foreach \y in {0,0.5,1} { \draw (0,\y) node [left]{\y} -- +(5pt,0); } \draw[thick] (0,0) rectangle (1.3,1.1); \node[below left] (xlabel) at (1.3,0) {$d/r$}; \draw[thick,domain=0:1,smooth,variable=\x,blue] plot ({\x},{1.+6*\x^5*(-21+70*\x-90*\x*\x+52.5*\x*\x*\x-35/3.*\x^4)}); \draw[line width=1.25,domain=0:1,smooth,variable=\x,red] plot ({\x},{(1-\x)^6*(1+6*\x+21*\x*\x+56*\x^3+126*\x^4+252*\x^5)}); \draw[line width=1pt,domain=0:0.99,smooth,variable=\x,samples=50,orange] plot ({\x},{exp(1-1/(1-\x*\x))}); \draw[line width=1.5pt,blue] (1,0) -- (1.3,0); \draw[line width=1.25pt,red] (1,0) -- (1.3,0); \draw[line width=1pt,orange] (1,0) -- (1.3,0); \draw[line width=0.75pt,green] (0,1) -- (1,1) -- (1,0) -- (1.3,0); \draw[line width=1.5pt,blue] (1.5,0.65) -- (1.6,0.65) node[right]{\texttt{C4}}; \draw[line width=1.25pt,red] (1.5,0.55) -- (1.6,0.55) node[right]{\texttt{C5}}; \draw[line width=1pt,orange] (1.5,0.45) -- (1.6,0.45) node[right]{\texttt{CINF}}; \draw[line width=0.75pt,green] (1.5,0.35) -- (1.6,0.35) node[right]{\texttt{CONSTANT}}; \end{tikzpicture} \]

The C4 and C5 functions are such that their integral on $[-1,+1]$ is equal to 1.

For example an input like

time function : t_stim PROFILE CINF START 2.5 DURATION 2.0
space function : s_stim PROFILE CINF CENTER 0.5 0.5 0.0 DIAMETER 0.05
volumic stimulation : vstim TIME t_stim SPACE s_stim AMPLITUDE 10.0 COMPONENT 1
time function : t_st START 0.5 DURATION 10.0 PERIOD 1
space function : s_st CENTER 0.1 0.2 0.5 DIAMETER 0.3
volumic stimulation : vst TIME t_stim SPACE s_stim COMPONENT 0

defines two volumic stimulation:

  • vstim : a non-periodic CINF stimulation centered at (0.5 0.5 0.0) with a diameter of 0.05. It starts at 2.5ms and lasts 2ms with an amplitude of 10 on component 1.
  • vst : a periodic C4 (default) stimulation centered at (0.1, 0.2, 0.5) with a diameter of 0.3. It starts at 0.5ms and ends at 10ms with a 1ms period. The amplitude is 1.0 (default) on component 0.

In the same way an input like

time function : t_stim PROFILE C4 START 0.5 DURATION 1.0
boundary stimulation : vstim TIME t_stim AMPLITUDE 15.0 COMPONENT 0 ATTRIBUTE 200

defines a boundary non-periodic stimulation on the first unknown on region with attribute 200 with values between 0 and 15.0 which starts at 0.5ms and lasts 1ms.




Time parameters

Time unit is the millisecond.

pde start time: <scalar>
pde end time: <scalar>

Beginning and ending time of simulation. Default start time is 0.

stop at complete activation : <yes,1,true>

For cardiac solvers, it is possible to stop the stimulation once all the points have seen an action potential, even before "pde end time". Of course, activation time data must be provided (see below).




PDE solver parameters

finite elements: <type>

Selects the type of finite elements to use. Currently, type can be P1 or P2.

pde solver: <type>

Selects the time discretization method of the PDE. type can be of

  • SBDF 1-4 (Stiff Backward Differentiation)
  • CN (Crank-Nicholon)
  • RK opt where opt can be an integer between 1 and 5 or Heun2 Heun3 Ralston3 SSP3 38 Ralston4 Nystrom5 Fehlberg5.
pde time step: <scalar>

The main time step of the simulation




Source terms

source term: <label1> [<scalar1>] [,<label2> [<scalar2>]]

FIXME labels ? not attributes ? For now only constant source terms are implemented. By default scalars are O, which is not really useful...




Linear algebra

linear system solver type: <type>
linear system absolute tolerance: <scalar>
linear system relative tolerance: <scalar>
linear system divergence tolerance: <scalar>
linear system max iterations: <integer>

Options of the matrix vector solver. Only iterative methods can be used for now, either CG (Conjugate Gradient) or GMRES (Generalized Minimum RESidual), the latter being suited for most of the problems, but often with less performance.

Tolerances are set to define the stopping criteria of the iterative method.




Output controls

output file name: <prefix>
snapshot time: <scalar>
post process time: <scalar>

All output files will be prefixed by the given parameter. It can contain a path to an output directory that will be created. A collection file that can be opened in Paraview is written, and data files are put in an outputFiles subdirectory.

Output will be written every snapshot time time units.

Post processed data will be computed every post process time time units. This parameter sets the temporal resolution of activation maps and all other time dependent outputs.




Cardiac outputs

save ionic current = <YES,TRUE,ON,1>
save stimulation current = <YES,TRUE,ON,1>

If active, additional data will be written with the values of ionic and/or stimulation currents

probe points: <x1> <y1> <z1> [,<x2> <y2> <z2>]

Places probe points where FIXME will be written. The z coordinates can be omitted for 2D problems. An additional file with name ending in _probes.txt will be written with time dependant data in a columns format.

activation time data: <scalar1> <scalar2> <scalar3> [...]

Values that are used to define depolarization times. For each unknown of the problem (TMP for monodomain, plus extracellular potential for bidomain), three numbers must be given:

  • the activation threshold of the unknown, defining the beginning of the APD window.
  • a minimum value of the unknown used to find the peak of the variable during AP. (FIXME useful ?)
  • a percentage XX used for computation of APDXX (as done for APD50). We assume the percentage is between 50% and 100%.
export MUSIC = <YES,TRUE,ON,1>

Writes a slightly different collection file and VTK legacy files instead of .vtu files. These files can be loaded into MUSIC.




Convergence tools

Additional parameters are required for convergence studies. These parameters can be found on the CEPS: Convergence Measurement of Implemented Methods related page.

Generated by doxygen 1.9.2