# Accelerating stabilization of whole-heart models after changes in cycle length

Hassaan Bukhari, Carlos Sánchez, Esther Pueyo, Mark Potse

► To cite this version:

Hassaan Bukhari, Carlos Sánchez, Esther Pueyo, Mark Potse. Accelerating stabilization of whole-heart models after changes in cycle length. Computing in Cardiology 2022, Sep 2022, Tampere, Finland. hal-03936920

HAL Id: hal-03936920

https://inria.hal.science/hal-03936920

Submitted on 12 Jan 2023

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Accelerating stabilization of whole-heart models after changes in cycle length

Hassaan A. Bukhari¹,²,³,⁴, Carlos Sánchez¹,², Esther Pueyo¹,², Mark Potse³,⁴

¹ I3A, University of Zaragoza, IIS Aragón, Zaragoza, Spain

² CIBER en Bioingeniería, Biomateriales y Nanomedicina, Spain

³ Carmen team, Inria Bordeaux – Sud-Ouest, Talence, France

⁴ Univ. Bordeaux, IMB, UMR 5251, Talence, France

## Abstract

Parameter changes can cause long-term drift in membrane models. To reduce the cost of whole-heart simulations with such changes the stabilization can be performed in isolated-cell models, but it can then still take many beats to stabilize the full model. We hypothesized that differences in activation time leading to cycle length (CL) variability before the first beat contribute to this. To remove this variability we froze most state variables of the model until the sodium current activated.

Simulations were performed with CL 400, 500, 600 and 1000 ms and modified Ten Tusscher-Panfilov 2006 dynamics. Isolated endocardial, mid-myocardial, and epicardial cells were simulated for 1000 beats. Their final states were then copied to a model of the whole human ventricles, which was run for 5 beats, with and without freezing.

Stabilization of the full model took three to four beats. Freezing of the membrane state accelerated stabilization in some cell types but caused opposite drifts in others. Drifts were largest in the epicardial and mid-myocardial layers, and not in particular at their interfaces.

Freezing of membrane state may help to accelerate stabilization but in our scenarios other types of drift dominate and may be aggravated by freezing, as it inhibits electrotonic interactions.

## 1. Introduction

Changes in the parameters of cardiac membrane models can cause drift in model variables that lasts for hundreds of beats. This is true in particular for changes in the extra-cellular concentrations of potassium ([$K^+$]₀) and calcium ([$Ca^{2+}$]₀).

[$K^+$]₀ and [$Ca^{2+}$]₀ outside normal ranges increase the risk for life-threatening arrhythmia and sudden cardiac death and they have also been shown to affect depolarization and repolarization features of the electrocardiogram (ECG) [1–3]. This means that simulations of chronic kidney disease, in which these concentrations are abnormal, require long

stabilization periods. This is problematic for large-scale models, for which such long periods are costly. Therefore stabilization is often performed in models of isolated cells, and then the membrane state variables are copied to the large-scale model. However, the model can then still take a few beats to reach a stable behavior, which can still cause high computational cost.

Aiming to further reduce this cost, we tested a method to reach a stable behavior faster. We hypothesized that part of the drift is caused by the different cycle lengths (CL) experienced by the cells at the start of the first beat in the whole-heart model: the cells are all initialized identically, but they activate at different times, and thus experience different CL. When the CL is short, this could lead to significant perturbations in the behavior of the cells during the first beat and possibly during subsequent beats as well. To prevent this, we froze the state variables of the cells in the whole-heart model until their first activation. Only the activation gate of the sodium channel was not frozen, as it is required for activation.

## 2. Methods

**Computational Modeling** Cardiac electrical activity was simulated using isolated cells as well as a model of the whole human ventricles created from computed tomography data of a single patient [4]. Cellular electrophysiology was represented by the human ventricular membrane model of Ten Tusscher and Panfilov [5], as modified by Severi et al. [6]. The ventricular model was simulated with a monodomain reaction-diffusion model with 200 µm resolution and a time step of 0.01 ms. Simulations were run for CL of 400, 500, 600 and 1000 ms.

12-lead ECGs were computed with a lead-field method [7] using lead fields calculated in a 1-mm resolution finite-difference torso model. All simulations were performed on a cluster computer using the Propag-5 software [8] with provisions for the handling of lead fields [7] and recent corrections in the handling of model initialization from stored membrane states.
<figure><img src="image_2.png" /><figcaption>Figure 1. Human torso-heart model used for the ECG simulations.</figcaption></figure>

**Stabilization of the Whole-Heart Model** In order to get the whole-heart model into a stable state as quickly as possible for each CL, we found initial values for all the membrane model variables that correspond to a stable state under the given CL. This was done by running single-cell simulations for all different cell types (endocardial, mid-myocardial and epicardial, for LV and RV) for 1000 beats, and then the membrane state variables from the last beat were copied to the cells of the same type in the large-scale model. Then the heart model was simulated for 5 beats to see how long it takes before it reaches a stable state.

Two initialization methods were tested: freezing and non-freezing. In the freezing method, at the beginning of the whole-heart simulation, the membranes are “frozen”, meaning that none of their state variables can change, except those of the activation gate of the sodium (Na) current. When the Na channel opens, the other variables are reactivated. This serves to let all model nodes experience the same CL, despite their different activation times. In the non-freezing method, all variables are integrated normally throughout the whole-heart simulation.

Model stability was evaluated in terms of action potential duration at 90% repolarization (APD$_{90}$), and in terms of correlation between the ECG of the last beat and the earlier beats.

**Statistical Analysis** Pearson correlation coefficients ($\rho$) were computed to assess the strength and the effects of different beats in each of the 12-lead ECGs at different CL.

To assess how different layers of heart (middle of endocardial, midmyocardial, epicardial layers and the bound-

aries of the layers) modulated APD$_{90}$ at different CL with and without freezing, a sensitivity analysis was performed. For APD$_{90}$, the percentage of change ($D_{M;c;a_i}$) and its sensitivity ($S_{M;c;a_1,a_2}$) to changes in each of the ventricular layers were computed as follows [9, 10]:

$$D_{M;c;a_i} = \left( \frac{M_{c;a_i} - M_{\text{control}}}{M_{\text{control}}} \right) \cdot 100, \quad i \in \{1, 2\} \quad (1)$$

$$S_{M;c;a_1,a_2} = \frac{(D_{M;c;a_2} - D_{M;c;a_1})100}{a_2 - a_1}, \quad (2)$$

where $M_{c;a}$ is the value of APD$_{90}$ marker $M$ under different levels of CL (a) in each of the beat (c) calculated at CL $a_1$ and $a_2$. The values of $a_1$ and $a_2$ were taken as the minimum and maximum value of CL in each beat, respectively. $M_{\text{control}}$ is the value of APD$_{90}$ at a CL of 1000 ms, for the 5th beat.

### 3. Results

Figure 2 shows the differences in APD$_{90}$ with and without freezing for 1st and 5th beat at different layers of the heart, with varying CL (top panels). The bottom panels show APD$_{90}$ and its differences with and without freezing for 1st to 5th beats at each tested CL. The effect of freezing is more visible at shorter CL, particularly at the 1st beat. Freezing accelerated the stabilization in the endocardial and epicardial layers, but induced an opposite drift in the midmyocardial layer and its boundaries.

Figure 3 illustrates the correlation coefficient between the 12-lead ECGs of first three beats and the 5th beat, with (right panel) and without (left panel) initialization for 1000 beats as well as with and without freezing. Higher correlation was obtained (median correlation coefficients of 0.994, 0.999 and 1.000) with initialization than without initialization (0.838, 0.899 and 0.960 for the 1st, 2nd and 3rd beat with freezing, respectively).

Table 1 shows the APD$_{90}$ values at varying CL, computed for the 1st to 5th beats. Results show that APD$_{90}$ values are stable from 2nd or 3rd beat onwards. Differences in the APD$_{90}$ values between the 1st and 5th beat were larger at shorter CL, particularly for the middle of the midmyocardial layer and its interface with the endocardial layer. Therefore, drifts could be largest in these layers at shorter CL. Freezing effects are more prominent in endocardial than in other layers.

Table 2 shows the differences in the sensitivities of APD$_{90}$ with (F) and without (nF) freezing, in different layers at varying CL for the 1st to 5th beat. The highest absolute differences in the sensitivities were obtained at the 1st beat, particularly for the middle of the midmyocardial layer and its boundary with the endocardial layer. However, sensitivity differences were almost negligible from the 3rd beat onwards in the whole-heart model.
<figure><img src="image_3.png" /><figcaption>Figure 2. Top panels: Differences in APD$_{90}$ with (F) and without (nF) freezing for 1st and 5th beat at different CL. Bottom panels: APD$_{90}$ and its differences with (F) and without (nF) freezing for 1st to 5th beats at CL=400 ms (bottom left panels) and 1000 ms (bottom right panels).</figcaption></figure>

<table><caption>Table 1. APD$_{90}$ at different CL and different layers (middle of endocardial, mEndo, midmyocardial, mMid, epicardial, mEpi, the boundary of endocardial, bEndo, and the boundary of endocardial with midmyocardial, endo-mid, layers) after freezing, computed from the 1st to the 5th beat.</caption><thead><tr><th rowspan="3">beat</th><th colspan="4">mEndo</th><th colspan="4">mMid</th><th colspan="4">mEpi</th><th colspan="4">bEndo</th><th colspan="4">endo-mid</th></tr><tr><th colspan="4">cycle length (ms)</th><th colspan="4">cycle length (ms)</th><th colspan="4">cycle length (ms)</th><th colspan="4">cycle length (ms)</th><th colspan="4">cycle length (ms)</th></tr><tr><th>400</th><th>500</th><th>600</th><th>1000</th><th>400</th><th>500</th><th>600</th><th>1000</th><th>400</th><th>500</th><th>600</th><th>1000</th><th>400</th><th>500</th><th>600</th><th>1000</th><th>400</th><th>500</th><th>600</th><th>1000</th></tr></thead><tbody><tr><td>1</td><td>253</td><td>268</td><td>279</td><td>308</td><td>268</td><td>284</td><td>296</td><td>329</td><td>239</td><td>252</td><td>262</td><td>289</td><td>254</td><td>270</td><td>282</td><td>309</td><td>266</td><td>283</td><td>295</td><td>328</td></tr><tr><td>2</td><td>253</td><td>268</td><td>279</td><td>308</td><td>275</td><td>288</td><td>299</td><td>332</td><td>242</td><td>254</td><td>264</td><td>290</td><td>254</td><td>270</td><td>281</td><td>310</td><td>272</td><td>286</td><td>298</td><td>330</td></tr><tr><td>3</td><td>252</td><td>267</td><td>279</td><td>308</td><td>274</td><td>288</td><td>300</td><td>333</td><td>242</td><td>255</td><td>264</td><td>291</td><td>254</td><td>270</td><td>281</td><td>309</td><td>272</td><td>287</td><td>298</td><td>331</td></tr><tr><td>4</td><td>252</td><td>267</td><td>279</td><td>308</td><td>274</td><td>289</td><td>300</td><td>334</td><td>242</td><td>255</td><td>264</td><td>291</td><td>254</td><td>270</td><td>281</td><td>309</td><td>272</td><td>287</td><td>299</td><td>331</td></tr><tr><td>5</td><td>252</td><td>267</td><td>279</td><td>308</td><td>274</td><td>289</td><td>301</td><td>334</td><td>242</td><td>255</td><td>264</td><td>291</td><td>254</td><td>270</td><td>281</td><td>309</td><td>272</td><td>287</td><td>299</td><td>332</td></tr></tbody></table>

## 4. Discussion

We investigated a “membrane freezing” method to accelerate the stabilization of whole-heart models after initialization with membrane states that were pre-stabilized with isolated-cell simulations. We found that freezing of the membrane state accelerated stabilization in some cell types but caused opposite drifts in others and did not significantly accelerate stabilization.

The idea behind the freezing method was that differ-

ences in apparent CL at the start of the first beat would be the most important cause of short-term drift. Although this certainly seems to be one of the causes, our results show that other factors are at least as important. One possible explanation is that when cells of different types are first coupled to each other, electrotonic interactions will change the behaviors of the cells and thus invalidate their initialization state, computed for a slightly different functional cell type. Thus, we would expect that the largest drifts occur at the interfaces between layers of different cell types.
<figure><img src="image_4.png" /><figcaption>Figure 3. Median and 25th/75th percentiles of the correlation coefficient ($\rho$) between the 12-lead ECGs of the first three beats and the 5th beat, without (left panel) and with (right panel) initialization for 1000 beats in the whole-heart. Blue and red boxplots show the effects with and without freezing, respectively. Each purple dot represents an individual ECG lead.</figcaption></figure>

<table><caption>Table 2. Results of the differences in the sensitivity analysis, $S_{m;p;a_1,a_2}$, with ($F$) and without ($nF$) freezing for different values of CL, at different layers of heart in each beat from minimum CL ($a_1$) to maximum CL ($a_2$).</caption><thead><tr><th rowspan="2">beat</th><th colspan="5">$|S_{nF_{APD_{90}}}-S_{F_{APD_{90}}}|$</th></tr><tr><th>mEndo</th><th>mMid</th><th>mEpi</th><th>bEndo</th><th>endo-Mid</th></tr></thead><tbody><tr><td>1</td><td>0.38</td><td>0.50</td><td>0.34</td><td>0.39</td><td>0.60</td></tr><tr><td>2</td><td>0.11</td><td>0.10</td><td>0.06</td><td>0.01</td><td>0.05</td></tr><tr><td>3</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.06</td><td>0.05</td></tr><tr><td>4</td><td>0.05</td><td>0.00</td><td>0.00</td><td>0.04</td><td>0.05</td></tr><tr><td>5</td><td>0.05</td><td>0.00</td><td>0.00</td><td>0.04</td><td>0.00</td></tr></tbody></table>

Interestingly, this turned out not to be the case either: the largest drifts occurred in the middle of the mid-myocardial layer. Nevertheless, it is possible that the freezing method delayed the stabilization by inhibiting the effects of electrotonic interactions until the first activation.

A shortcoming of our study is that we did not systematically investigate the relation between drift and activation time, the freezing method being expected to be most effective for late-activated areas which, without freezing, experience the largest deviation from the CL for which they were initialized.

We conclude that to understand and possibly prevent short-term drift in heterogeneous heart models we must better understand the effects of electrotonic coupling between cells with different intrinsic types.

## Acknowledgements

This work was supported by projects ERC-StG 638284 (ERC), PID2019-105674RB-I00 (Ministerio de Ciencia e Innovación) and Marie Skłodowska-Curie grant 764738 (European Commission) and by European Social Fund (EU) and Aragón Government through BSICoS group T39_20R and project LMP94_21. This work was granted access to the HPC resources of IDRIS under the allocation 2021-A0110307379 made by GENCI. M. Potse was supported by the French National Research Agency, grant reference ANR-10-IAHU04-LIRYC.

## References

[1] Weiss JN, et al. Electrophysiology of hypokalemia and hyperkalemia. Circ Arrhythm Electrophysiol Mar. 2017;10:1–10.

[2] Noordam R, et al. Effects of calcium, magnesium, and potassium concentrations on ventricular repolarization in unselected individuals. J Am Coll Cardiol June 2019; 73(24):3118–3131.

[3] Gardner JD, et al. ECG diagnosis: The effect of ionized serum calcium levels on electrocardiogram. Perm J 2014; 18(1):e119–e120.

[4] Kania M, et al. Prediction of the exit site of ventricular tachycardia based on different ECG lead systems. In 2017 Computing in Cardiology (CinC). September 2017; 1–4.

[5] Ten Tusscher KHWJ, et al. Alternans and spiral breakup in a human ventricular tissue model. Am J Physiol Heart Circ Physiol Sept. 2006;291(3):H1088–H1100.

[6] Severi S, et al. From in vivo plasma composition to in vitro cardiac electrophysiology and in silico virtual heart: the extracellular calcium enigma. Philos Trans A Math Phys Eng Sci Jun. 2009;367(1896):2203–23.

[7] Potse M. Scalable and accurate ECG simulation for reaction-diffusion models of the human heart. Front Physiol 2018;9.

[8] Krause D, et al. Hybrid parallelization of a large-scale heart model. In Keller R, Kramer D, Weiss JP (eds.), Facing the Multicore-Challenge II, volume 7174 of Lecture Notes in Computer Science. Berlin: Springer, 2012; 120–132.

[9] Romero L, et al. Impact of ionic current variability on human ventricular cellular electrophysiology. Am J Physiol Heart Circ Physiol Oct. 2009;297(4):1436–1445.

[10] Bukhari HA, et al. Characterization of T wave amplitude, duration and morphology changes during hemodialysis: Relationship with serum electrolyte levels and heart rate. IEEE Trans Biomed Eng 2021;68(8):2467–2478.