<figure><img src="image_1.png" /></figure>

HAL
open science

# Continuity of Mircoscopic Cardiac Conduction in a Computational Cell-by-Cell Model

Joshua Steyer, Fatemeh Chegini, Mark Potse, Axel Loewe, Martin Weiser

<figure><img src="image_2.png" /></figure>

► To cite this version:

Joshua Steyer, Fatemeh Chegini, Mark Potse, Axel Loewe, Martin Weiser. Continuity of Mircoscopic Cardiac Conduction in a Computational Cell-by-Cell Model. CinC 2023 - 50th Computing in Cardiology Conference, Oct 2023, Atlanta (GA), United States. 10.22489/CinC.2023.385. hal-04407818

HAL Id: hal-04407818

https://inria.hal.science/hal-04407818v1

Submitted on 21 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Continuity of Microscopic Cardiac Conduction in a Computational Cell-by-Cell Model

Joshua Steyer<fnref n="1" />, Fatemeh Chegini<fnref n="2" />, Mark Potse³,⁴,⁵, Axel Loewe<fnref n="1" />, Martin Weiser<fnref n="2" />

## Abstract

Conduction velocity in cardiac tissue is a crucial electrophysiological parameter for arrhythmia vulnerability. Pathologically reduced conduction velocity facilitates arrhythmogenesis because such conduction velocities decrease the wavelength with which re-entry may occur. Computational studies on CV and how it changes regionally in models at spatial scales multiple times larger than actual cardiac cells exist. However, microscopic conduction within cells and between them have been studied less in simulations. In this work, we study the relation of microscopic conduction patterns and clinically observable macroscopic conduction using an extracellular-membrane-intracellular model which represents cardiac tissue with these subdomains at subcellular resolution. By considering cell arrangement and non-uniform gap junction distribution, it yields anisotropic excitation propagation. This novel kind of model can for example be used to understand how discontinuous conduction on the microscopic level affects fractionation of electrograms in healthy and fibrotic tissue. Along the membrane of a cell, we observed a continuously propagating activation wavefront. When transitioning from one cell to the neighbouring one, jumps in local activation times occurred, which led to lower global conduction velocities than locally within each cell.

## 1. Introduction

Being one of the main reasons for sudden cardiac death, which causes around 700,000 deaths per year in Europe alone, cardiac arrhythmias remain a major public health issue [1]. An important tissue property affecting arrhythmia vulnerability is the conduction velocity (CV) of cardiac excitation waves, which can be derived from local activation times (LATs), defined as the moment at which the excitation wave front reaches a certain point in space. The spatial distribution [2] and effect [3] of these two quantities

have been investigated in homogenised models in which a model element averages over several hundreds of myocytes. Such models, however, do not directly consider myocyte size and shape, an inhomogeneous arrangement of them, connectivity patterns, and microscopic conduction from cell to cell which has been found to play a role in arrhythmogenesis [4].

In this study, we investigate the distribution of LAT and CV values using the novel extracellular-membrane-intracellular (EMI) model, which explicitly represents these three domains in the mesh [5, 6]. To do so, we used a 3D myocyte configuration embedded in a bath for which excitation propagation is initiated for three different cases to study purely longitudinal and transversal wave propagation, as well as a combination of them. We then determined LATs both on the membrane and the intracellular space and derived local CV, which were both analysed in terms of difference between the local (subcellular and cell-to-cell) and global behaviour of the excitation propagation.

## 2. Methods

We used the EMI model described in [5], which represents the extracellular medium $\Omega_e$, intercalated discs $\Gamma_g$, the cell-to-extracellular membrane $\Gamma_m$ and the intracellular domain $\Omega_i = \bigcup_j \Omega_j$ consisting of individual myocyte domains, $\Omega_j$, explicitly in the mesh. Its model equations read

$$ 
\left\{ 
\begin{array}{ll}
\nabla \cdot (\sigma_i \nabla \phi_i) = 0 & \text{in } \Omega_i, \\
\nabla \cdot (\sigma_e \nabla \phi_e) = 0 & \text{in } \Omega_e, \\
V_m = \phi_i - \phi_e & \text{on } \Gamma_m, \\
w_g = \phi_j - \phi_k \quad \text{where } j \neq k & \text{on } \Gamma_g, \\
C_m \dot{V}_m + I_{\text{ion}} = -n_i^T \sigma_i \nabla \phi_i = n_e \sigma_e \nabla \phi_e & \text{on } \Gamma_m, \\
C_m \dot{w}_g + I_g = -n_j^T \sigma_i \nabla \phi_j = n_k \sigma_i \nabla \phi_k & \text{on } \Gamma_g,
\end{array}
\right. 
\quad (1) 
$$

where $\phi$ denotes the potential in the corresponding domain, while $V_m$ describes the membrane voltage between $\Omega_i$ and $\Omega_e$, and $w_g$ between two adjacent cell domains,

<footnote n="1">¹ Karlsruhe Institute of Technology, Institute of Biomedical Engineering, Karlsruhe, Germany</footnote>

<footnote n="2">² Zuse Institute Berlin, Berlin, Germany, ³ LIRYC, Pessac, France, ⁴ Inria Centre Bordeaux, Talence, France, ⁵ Bordeaux Institute of Mathematics, Talence, France</footnote>
$\Omega_j$ and $\Omega_k$. $n$ denotes the unit outer normal of the respective domain and the conductivities of the intra- and extracellular domains are $\sigma_i = 0.3$ S/m and $\sigma_e = 2.0$ S/m, and the membrane capacitance is $C_m = 10^{-4}$ F/m$^2$. To model membrane dynamics, we used the Aliev-Panfilov model [7] with its ion current given by

$$I_{\text{ion}}(V_m, w_g) = g_a V_m (V_m - a)(V_m - 1) + V_m w_g, \quad (2)$$

and the gating dynamics being governed by

$$\frac{\partial w_g}{\partial t} = -\frac{1}{4}\varepsilon(V_m, w_g)(w_g + g_s V_m(V_m - a - 1)), \quad (3)$$

where $\varepsilon(V_m, w_g) = \varepsilon_0 + \mu_1 w_g / (V_m + \mu_2)$, $g_a = g_s = 8.0$, $a = 0.1$, $\varepsilon_0 = 0.01$, $\mu_1 = 0.07$ and $\mu_2 = 0.3$. Being purely phenomenological, the normalised state variables of the dynamics $V_m$ and $w_g$ are restricted to [0,1] a.u. Similar to [8], we used a linear (Ohmic) relationship between the current through the gap junctions and the respective voltage between the two cells:

$$I_g = \frac{w_g}{R_g}, \quad (4)$$

where $R_g = 4.5 \cdot 10^{-4} \Omega m^2$. This system was solved using the Kaskade 7 finite element toolbox [9].

The 3D mesh [10] is shown in Fig. 1. Its dimensions are 400 µm × 240 µm × 240 µm and it consists of 507648 tetrahedral mesh elements. The intracellular domain, $\Omega_i$ comprises $n_x \times n_y \times n_z = 7 \times 3 \times 7 = 147$ myocytes.

<figure><img src="image_4.png" /><figcaption>Figure 1: Myocyte arrangement used for the simulations embedded in $\Omega_e$ (semi-transparent grey box) and an exemplary single myocyte (top) with edges indicated in white. Note that colours are used to discriminate between individual myocytes, not varying electrophysiological properties.</figcaption></figure>

The myocytes' largest extension (i.e. preferential orientation) is along the x-axis, which is thus the direction of fastest propagation and may be referred to as longitudinal direction. For transversal direction, we consider the y-axis. The intercalated discs with the gap junctions are located at $x = 150, 250$ µm along the longitudinal axis, such that we have two planes ($y, z$) with gap junction properties given by Eq. 4.

In total, we ran three different simulations with different initial conditions leading to different excitation dynamics: a) initial depolarisation of all leftmost myocytes along the x-axis (causing an almost purely longitudinal propagation and referred to as the **left** case), b) depolarisation of all lowermost cells (purely transversal propagation, **bottom** case) and c) depolarisation of one cell in the corner of the myocyte configuration (longitudinal and transversal propagation, **single** case).

To quantify the spatial propagation of the excitation along the cell membrane (quantified via $V_m(x, t)$) and the intracellular potential ($\phi_i(x, t)$), where $x \in \mathbb{R}^3$ denotes the position, we chose $V_m^{th} = 0.5$ a.u and $\phi_i^{th} = 0.45$ a.u.s thresholds to determine the local activation times (LATs). Using a projection of these LAT values onto the primary axis of propagation, i.e. either $x$ or $y$, we then defined longitudinal (‖) and transversal (⊥) CVs as:

$$CV_{\parallel}(x) = \left| \frac{\partial \text{LAT}(x)}{\partial x} \right|^{-1}, \quad (5)$$

$$CV_{\perp}(x) = \left| \frac{\partial \text{LAT}(x)}{\partial y} \right|^{-1}. \quad (6)$$

### 3. Results

As already mentioned, LATs are obtained with the thresholds imposed on the time evolution of $V_m$ and $\phi_i$. An exemplary snapshot of $\phi_i(x, t)$ is given in Fig. 2 for a purely longitudinal propagation. By determining the nearest neighbour nodes, we compared the LAT values obtained via $V_m$ and $\phi_i$.

<figure><img src="image_5.png" /><figcaption>Figure 2: Exemplary snapshot of a purely longitudinal propagation of $\phi_i(t)$ initialised by depolarising all leftmost myocytes.</figcaption></figure>
In all of the three simulations, the Pearson correlation was always $\rho > 0.999$, which is why we here look at both LAT datasets combined. The projections of these LAT values onto these axes are given in Fig. 3. For all cases, we fitted a linear function to a) the data within each “column” of cells (perpendicular to the propagation direction, green lines) and b) the complete data set (red line). The upper left panel shows the projection of LAT values of the central and right myocyte columns onto the x-axis (i.e. the axis of propagation). Both values grow linearly in propagation direction distance with a jump at the intercalated disc ($x = 250 \mu m$). The upper right panel illustrates the distribution of LAT values along the transversal axis obtained from the “bottom” case. Note that, in contrast to the fully longitudinal propagation, the separate myocyte planes overlap, since the cells are intertwined both along the y and z direction (Fig. 1). Finally, for the depolarisation of one cell in the corner, we can project the LATs onto both the longitudinal and transversal direction, given in the bottom row of Fig. 3, respectively. While the LATs in transversal direction show linear increase in propagation direction, both within the myocyte planes and globally, the longitudinal data for this initialisation does not show a clear trend and is broadly distributed.

<figure><img src="image_6.png" /><figcaption>Figure 3: Spatial distribution of LAT values on the membranes as projections onto the axis of propagation. Green lines indicate the linear regression of LAT values within one cell layer (see Fig. 2) and the red line the linear regression of all LATs along this axis. Dark grey background indicates intercalated disc regions, which are wider in the right column due to intertwining of myocytes along the y-axis. Bottom row: Diagonal excitation starting from a single myocyte in the bottom left corner of the mesh.</figcaption></figure>

We now quantify the time delay of the activation propagation when transitioning from one myocyte to the next via

the gap junctions in the intercalated disc. Both the mean and standard deviation of the LAT delays, $\Delta$LAT, are given in Tab. 1.

<table><caption>Table 1: Mean and standard deviation of LAT delays on the membrane during transition from one cell layer to the neighbouring one.</caption><thead><tr><td>propagation</td><td>ΔLAT (ms)</td></tr></thead><tbody><tr><td>fully long.</td><td>0.8</td></tr><tr><td>fully trans.</td><td>1.6 ± 0.9</td></tr><tr><td>long. (single)</td><td>-0.1 ± 0.3</td></tr><tr><td>trans. (single)</td><td>1.5 ± 0.7</td></tr></tbody></table>

The local and global CVs defined as the inverse slopes of the single myocyte column regressions and the inverse of the global slope, according to Eq. 6, are given in Tab. 2. Overall, the global CVs are smaller than the locally determined ones.

<table><caption>Table 2: Mean and standard deviation of the single myocyte layer LAT slopes and the global one along propagation direction $x_i$.</caption><thead><tr><td>propagation</td><td>loc. CV̄ (m/s)</td><td>glob. CV (m/s)</td></tr></thead><tbody><tr><td>fully long.</td><td>0.0839 ± 0.0006</td><td>0.06</td></tr><tr><td>fully trans.</td><td>0.06 ± 0.03</td><td>0.01</td></tr><tr><td>long. (single)</td><td>1.7 ± 2.5</td><td>1.6</td></tr><tr><td>trans. (single)</td><td>0.08 ± 0.05</td><td>0.01</td></tr></tbody></table>

## 4. Discussion

We have analysed local CV values in the EMI model [5, 6], which explicitly represents intra- and extracellular media, as well as the membrane separately in the mesh. As opposed to what would be possible with homogenised models, such as the frequently used monodomain or bido main models [11], this allowed studying the the local behaviour of the potential propagation. We have shown that LATs undergo jumps when the wave transitions from one myocyte to the next. CVs, calculated by the inverse slope of the linearly fitted LAT values reveal that the potential and excitation propagation within a cell is faster than globally, which in turn can be explained by the time delay in the order of a few milliseconds caused at each cell-to-cell transition. This is in accordance with what has been theorised e.g. in [12] and shown in 2D simulations using the EMI model [13]. The absolute value of global and local CV is a result of the choice of model parameters and gap junction distribution [14]. The setup used here motivated by
previous work led to CV values one to two orders of magnitude smaller than typically observed in healthy human tissue [2, 15]. There are several limitations to this analysis. First, the method to calculate CVs used here is not suitable for the longitudinal propagation in the single cell initialisation case due to multimodal LATs, as can be seen in Fig. 3, and will have to be done node-wise. Furthermore, we assumed the intracellular medium to be homogeneous without considering potential conduction anisotropy due, to e.g., the nucleus or other organelles in the cell. The findings presented here provide a first insight into electrophysiological excitation propagation behaviour (CV) on a microscale quantitatively. If the model parameters are adjusted such that physiological CV values are obtained, it may form a basis for future work analyzing for example the genesis of intracardiac electrograms in such microstructural models [16] or extending the setup to pathological scenarios allowing for mechanistical studies of the role of microscopic cell arrangement and distribution of, e.g., fibrosis.

## Acknowledgments

This work was supported by the European High-Performance Computing Joint Undertaking EuroHPC under grant agreement No 955495 (MICROCARD) co-funded by the Horizon 2020 programme of the European Union (EU), the German Federal Ministry of Education and Research and the French National Research Agency ANR.

## References

[1] Chugh SS. Sudden cardiac death in 2017: Spotlight on prediction and prevention. International Journal of Cardiology. (2017);237:2-5.

[2] Verma B, Oesterlein T, Loewe A, Luik A, Schmitt C, Dössel O. Regional conduction velocity calculation from clinical multichannel electrograms in human atria. vol. 92. Elsevier; 2018. p. 188-96.

[3] Azzolin L, Eichenlaub M, Nagel C, Nairn D, Sanchez J, Unger L, et al. Personalized ablation vs. conventional ablation strategies to terminate atrial fibrillation and prevent recurrence. Europace. 2023 2;25(1):211-22.

[4] de Bakker JMT, Coronel R, Tasseron S, Wilde AAM, Opthof T, Janse MJ, et al. Ventricular tachydia in the infarcted, Langendorff-perfused human heart: Role of the arrangement of surviving cardiac fibers. Journal of the American College of Cardiology. 1990;15(7):1594-607.

[5] Tveito A, Mardal KA, Rognes ME. Modeling Excitable Tissue: The EMI Framework. Simula Springer Briefs on Computing. Springer International Publishing; 2020.

[6] Jæger KH, Edwards AG, Giles WR, Tveito A. From Millimeters to Micrometers: Re-introducing Myocytes in Models of Cardiac Electrophysiology. Frontiers in Physiology. 2021;12:763584.

[7] Aliev RR, Panfilov AV. A simple two-variable model of cardiac excitation. Chaos, Solitons & Fractals. 1996;7(3):293-301.

[8] Chegini F, Steinke T, Weiser M. Efficient adaptivity for simulating cardiac electrophysiology with spectral deferred correction methods. In: Submitted; 2022.

[9] Weiser M, Chegini F. Higher-order time integration using spectral deferred correction method (SDC) in a cell by cell discretization of cardiac excitation. Martin Weiser, Fatemeh Chegini; 2022. Available from: https://doi.org/10.35097/716.

[10] Potse M, Cirrottola L, Froehly A. A practical algorithm to build geometric models of cardiac muscle structure. In: 8th European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS). Oslo, Norway; 2022.

[11] Plank G, Loewe A, Neic A, Augustin C, Huang YL, Gsell MAF, et al. The openCARP simulation environment for cardiac electrophysiology. Comput Methods Programs Biomed. 2021 6;208:106223.

[12] Keener J, Sneyd J. Mathematical Physiology: I: Cellular Physiology. Interdisciplinary Applied Mathematics. Springer New York; (2010).

[13] Jæger K, Tveito A. Efficient, cell-based simulations of cardiac electrophysiology; The Kirchhoff Network Model (KNM). NPJ systems biology and applications. 2023 06:9;25.

[14] Rosilho de Souza G, Pezzuto S, Krause R. Effect of Gap Junction Distribution, Size, and Shape on the Conduction Velocity in a Cell-by-Cell Model for Electrophysiology. In: Bernard O, Clarysse P, Duchateau N, Ohayon J, Viallon M, editors. Functional Imaging and Modeling of the Heart. vol. 13958 of Lecture Notes in Computer Science. Cham: Springer; 2023. p. 117-26.

[15] Lubrecht JM, Grandits T, Gharaviri A, Schotten U, Pock T, Plank G, et al. Automatic reconstruction of the left atrium activation from sparse intracardiac contact recordings by inverse estimate of fibre structure and anisotropic conduction in a patient-specific model. EP Europace. 2021 Mar;23:i63-70.

[16] Steyer J, Diaz LPM, Unger LA, Loewe A. Simulated Excitation Patterns in the Atria and Their Corresponding Electrograms. In: Bernard O, Clarysse P, Duchateau N, Ohayon J, Viallon M, editors. Functional Imaging and Modeling of the Heart. Cham: Springer Nature Switzerland; 2023. p. 204-12.

Address for correspondence:

Joshua Steyer
Institute of Biomedical Engineering
Karlsruhe Institute of Technology (KIT)
Fritz-Haber-Weg 1, 76131 Karlsruhe, Germany
publications@ibt.kit.edu