DREAM reentry simulation example

See code in GitLab.
Author: Tobias Gerach tobias.gerach@kit.edu, Stephanie Appel stephanie.appel@kit.edu

cd ${EXAMPLES}/03_Eikonal/03B_DREAM_reentry

Introduction

In the Diffusion Reaction Eikonal Alternant Model (DREAM) the possibility of simulating a reentry is given, since it uses the cyclical fast iterative method to compute activation times from the eikonal equation [Barrios2024]. In each iteration, converged nodes are allowed to be activated again if they fulfill certain conditions. These conditions involve the repolarization time of the individual cell at each node. The most important condition that is relevant for the user is that the DREAM uses conduction velocity restitution based on the diastolic interval. This makes it possible to respect the refractory state of the cell and to reflect the typical behavior of cardiac tissue during pacing, i.e. a reduced conduction velocity when pacing frequency is increased. By default, conduction velocity restitution in DREAM is calculated by the formula

\[v(d) = v_\mathrm{ref} \cdot \left( 1 - \rho \cdot \exp(-\frac{\log(\rho)}{\psi} (d + \kappa))\right)\]

if \(d > \theta\), where \(d\) is the diastolic interval and \(v\) is conduction velocity. Therefore, \(\theta\) determines a cutoff below which conduction is not possible anymore, since we are inside the refractory phase of the cell.

Problem setup

Mesh

First, we build a geometry for this experiment using mesher and the dynamic retagging feature of openCARP. The mesh is a simple 60mm by 60mm slab, which is discretized using triangles. In the centre of this mesh, we place a spherical region with elemTag = 2. Later, this will represent a non-conductive region with a conduction velocity of 0. With dynamic retagging, we place a conductive channel through the middle of this spherical region with elemTag = 3. By default, the mesh resolution is set to \(1000 \mu m\), but you can change it during runtime with the parameter --dx.

Ionic model

This example uses the Courtemanche ionic model for human atrial myocytes [Courtemanche1998]. Since we want to simulate a reentry, we modify ion channel conductances using the -imp_region[].im_param parameter according to [Loewe2016] to include remodelling during atrial fibrillation. To switch between a healthy and AF parameterization, use the option --AF or --no-AF, respectively.

Pacing protocol

We use a simple S1-S2 protocol to induce a reentry. The S1 stimulus is applied on the left edge of the 2D slab. The S2 stimulus is applied to a cluster of nodes 4mm around the coordinate (9, 30, 0)mm and the time is set by the user using the parameter --S2 XXX with XXX being the delay in ms after the S1 start time. The parameter --pulse_strength sets the stimulus current strength used for both S1 and S2. This same current also determines how strongly a node is still affected when a stimulus reaches it but does not trigger a full activation due to refractoriness. Setting it to 0 means that node is not affected at all, so it can recover and be reactivated by nearby tissue sooner than it should, which can prevent a conduction block from holding. The simulation warns you if this happens.

Usage

The following optional arguments are available (default values are indicated):

python3 ./run.py --help
...
Example specific options:
  --dx DX               Set resolution of the mesh in um (default: 1000.0)
  --cv CV               Set conduction velocity in mm/s in fiber direction (default: 400.0)
  --ar AR               Set anisotropy ratio perpendicular to fiber direction. (default: 1.0)
  --model {0,1,2,3}     Choose model to solve: Eikonal(0), RE+(1), RE-(2), DREAM(3) (default: 3)
  --S2 S2               Set start time of S2 stimulus (default: None)
  --pulse_strength PULSE_STRENGTH
                        Set the stimulus current strength in uA/cm^2 for both stimuli S1 and S2. (default: 20.0)
  --plotCVrest          Plot CV restitution properties (default: False)
  --AF, --no-AF         Determine if you want to use Courtemanche with remodelled ionic channels (--AF) or using the default
                        Courtemanche (--no-AF). (default: None)

After running the script, activation time, repolarization time, diffusion current and transmembrane voltage are output to the files act.igb, rpt.igb, idiff.igb, and vm.igb, respectively. Running the script with the parameter --visualize will open the transmembrane voltage file in meshalyzer after the simulation is done.

Tasks

1. Use the CV restitution properties provided above to determine an S2 stimulus time that induces a reentry.

2. Can you determine an S2 stimulus time for which a reentry is possible using the standard model parameters of Courtemanche? Use the option --no-AF to turn AF parameterization off.

Solutions to the tasks

1. According to the provided CV restitution, the minimum diastolic interval at which a stimulus can be propagated is 66.28ms.

First, check the repolarization time (RT) of nodes across the S2 stimulus area in the output file rpt.igb. Since S1 reaches different parts of this area at slightly different times, RT also varies across the stimulus area. For example, node ID 1836 (further left within the stimulus area) has an RT of 99.225ms, while node ID 1840 (center area) has an RT of 110.575ms.

Propagation through a given point only becomes possible again once the minimum diastolic interval of 66.28ms has passed since its RT. For node 1836, this is from 99.225 + 66.28 = 165.505ms onward; for node 1840, from 110.575 + 66.28 = 176.855ms onward.

Choose an S2 stimulus time late enough for the side of the stimulus region that the S1 wavefront reached earlier (here: the left side) to have already recovered. At the same time, it should be early enough that the side the wavefront reached later (here: the right side) has not yet recovered. This produces a unidirectional block, where propagation succeeds toward the side that recovered first, but is blocked toward the side that is still refractory.

Therefore, one possible solution is

python3 run.py --visualize --S2 175 --AF

2. We can approach this in the same manner as task 1. However, the default Courtemanche parameterization has a longer action potential duration and using --plotCVrest and --no-AF shows the minimum diastolic interval is now around 137ms.

Following the same approach, node 1836 (left) has an RT of 193.825ms and node 1840 (center) has an RT of 208.075ms, giving thresholds of 330.825ms and 345.075ms respectively. An S2 stimulus time at 345ms is one option to again produce a unidirectional block.

However, this does not lead to a stable reentry. By the time the wavefront circles back around, the tissue is not yet re-excitable, blocking propagation for the same timing argument as before. Therefore, without changing other parameters such as CV, the answer is no. We cannot induce a reentry with the given protocol using the default Courtemanche model.

The following command generates the figures

python3 run.py --visualize --no-AF --S2 345

Recent questions tagged experiments, examples, tissue, eikonal