See code in GitLab.
Author: Moritz Linder moritz.linder@kit.edu, Tobias Gerach tobias.gerach@kit.edu
cd ${TUTORIALS}/02_EP_tissue/03E_study_resolution
This example implements the classic benchmark problem for electrophysiology introduced by Niederer et al. in 2011. [Niederer2011]
The geometry of the tissue model is defined by a cuboid with dimensions 20.0 x 7.0 x 3.0 mm. Fiber directions are defined to be oriented in the X axis of the cuboid. All boundaries are assumed to have a zero-flux boundary condition. The stimulus current is delivered to a volume of 1.5 x 1.5 x 1.5 mm originating at the corner of the slab with coordinates \(x=y=z=0\).
The ionic model that is used in the benchmark is the ten Tusscher & Panfilov model [TTP] with the epicardial cell model variant.
This is set in the base parameter file alongside a single cell surface-to-volume ratio of \(0.14 \mu m^{-1}\) and the portion of the volume that is occupied by cells (=1).
To use a given set of initial values for the state variables, we load the state variable file singlecell.sv with the parameter -imp_region[0].im_sv_init
.
num_imp_regions = 1
imp_region[0].name = "Ventricle"
imp_region[0].im = tenTusscherPanfilov
imp_region[0].im_param = "flags=EPI"
imp_region[0].cellSurfVolRatio = 0.14
imp_region[0].volFrac = 1
The tissue is assumed to be transverse isotropic. Therefore, we set a single conductivity region with the intra- and extracellular conductivites in longitudinal and transversal directions.
num_gregions = 1
gregion[0].g_il = 0.17
gregion[0].g_it = 0.019
gregion[0].g_in = 0.019
gregion[0].g_el = 0.62
gregion[0].g_et = 0.24
gregion[0].g_en = 0.24
To set up the stimulus current, we define a box by providing the lower left corner \(p0 = (0,0,0)\) and the upper right corner \(p1 = (1500,1500,1500)\). The stimulation protocol is set to start at 0 s with a duration of 2 ms and a pulse strength of \(50000 \mu A cm^{-3}\).
num_stim = 1
stim[0].elec.geom_type = 2
stim[0].elec.p0[0] = 0
stim[0].elec.p0[1] = 0
stim[0].elec.p0[2] = 0
stim[0].elec.p1[0] = 1500
stim[0].elec.p1[1] = 1500
stim[0].elec.p1[2] = 1500
stim[0].pulse.strength = 35.71
stim[0].ptcl.duration = 2.00
stim[0].ptcl.start = 0.00
A quantitative comparison of the benchmark simulations is achieved by comparing activation times at certain points of the geometry and along a diagonal line from one corner of the cuboid to another. An activation time is recorded at the time the membrane potential passes through 0 mV for the first time. This is achieved by setting the following parameters:
num_LATs = 1
lats[0].measurand = 0
lats[0].all = 0
lats[0].mode = 0
lats[0].threshold = 0.
lats[0].method = 1
The experiment specific options are:
--tend
duration of the simulation in ms
--dx
resolution of the mesh in um to be used for benchmark
--dt
temporal resolution in us
--massLumping
toggle mass massLumping on (1) or use full mass matrix (0)
--plot
generate a plot of activation times along the block diagonal
--compare
run benchmark with --dx 100 200 and 500 and make comparison plot
The user can either run single experiments by calling python3 run.py
with the optional arguments --tend
, --dx
, --dt
, and --massLumping
or run the entire benchmark setup with the argument --compare
.
With this option you will run all simulations required for the benchmark and create a comparison plot for the activation times on the block diagonal with \(\Delta x = 0.5, 0.2, 0.1\) mm resolution.
Notice that the accuracy of the solution deteriorates with lower resolution meshes.
Now we will run the benchmark setting with the additional parameter --massLumping 1
to see how mass lumping affects the solution in the otherwise unchanged setup.
Mass lumping is a common numerical technique in FEM to speed up computation. However, notice how the accuracy of the solution deteriorates even faster than before. This was one of the beneficial outcomes of the benchmark, which showed that mass lumping was leading to the biggest inaccuracies among finite element codes especially with lower spatial resolutions. Therefore, you should always remember to consider the tradeoffs of this numerical technique.
References
Niederer S. A. et al., Verification of cardiac tissue electrophysiology simulators using an <i>N</i>-version benchmark, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 369, 4331-4351, 2011. [ePrint]
ten Tusscher, K. H. W. J. & Panﬁlov, A. V., Alternans and spiral breakup in a human ventricular tissue model, Am. J. Physiol. Heart Circ. Physiol. 291, H1088-H1100, 2006. [ePrint]
© Copyright 2020 openCARP project Supported by DFG and EuroHPC Contact Imprint and data protection