# Exploring Chemotherapy-Induced Cardiotoxicity Combining A 3D Computational Model and Preclinical Cardiac Imaging Data

Javier Villar-Valero, Jesus Jairo Rodríguez Padilla, Buntheng Ly, Juan F Gomez, Mihaela Pop, Beatriz Trenor, Maxime Sermesant

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

## ► To cite this version:

Javier Villar-Valero, Jesus Jairo Rodríguez Padilla, Buntheng Ly, Juan F Gomez, Mihaela Pop, et al. Exploring Chemotherapy-Induced Cardiotoxicity Combining A 3D Computational Model and Preclinical Cardiac Imaging Data. STACOM 2024 - 15th workshop on Statistical Atlases and Computational Modeling of the Heart, Oct 2024, Marrakesh, Morocco. hal-04840698

HAL Id: hal-04840698

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

Submitted on 16 Dec 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Exploring Chemotherapy-Induced Cardiotoxicity Combining A 3D Computational Model and Preclinical Cardiac Imaging Data

Javier Villar-Valero¹[0009-0003-1889-0515], Jesus Jairo Rodríguez Padilla²,³[0000-0002-2013-1898], Buntheng Ly³, Juan F. Gomez⁴, Mihaela Pop²,⁵, Beatriz Trenor¹, and Maxime Sermesant²,³

¹ Universitat Politècnica de València, Valencia, Spain

² Centre Inria d'Université Côte d'Azur, Epione team, Sophia Antipolis, France

³ {abc,lncs}@uni-heidelberg.de

⁴ IHU-Liryc, Université de Bordeaux, Pessac, France

⁵ Valencian International University, Valencia, Spain

⁶ Sunnybrook Research Institute, Toronto, Canada

**Abstract.** Cardiotoxicity following anthracycline-based chemotherapy is becoming a critical clinical burden, leading in time to heart dysfunction such as arrhythmia and heart failure. The purpose of this computational study is to better understand the structure-function interaction and subtle longitudinal changes in the presence of progressive fibrosis induced by chemotherapeutic drugs such as doxorubicin, DOX. Specifically, here we investigated via simulations how the cardiac action potential wave propagates in the presence of collagenous fibrotic clusters generated by this drug. For this, we built isotropic and anisotropic 3D left ventricular (LV) models from in vivo contrast-enhanced MRI data (1.4mm isotropic resolution) acquired in n=2 pigs pre-DOX and post-DOX therapy (at weeks 5 in one pig, at weeks 9 in the other pig). For computations, we employed a modified biophysical Mitchell-Schaeffer model with a fast GPU-based Lattice-Boltzmann numerical method (which runs directly on MR images Cartesian grid). Overall, our results showed that several key parameters (fibrosis density, tissue conductivity, action potential duration, anisotropy) independently contribute to the perturbation of the electrical wave propagation. Future work will focus on simulating potentially lethal arrhythmia in hearts subjected to DOX therapy.

**Keywords:** Cardiotoxicity · Arrhythmia · Fibrosis · Computational Modeling · Simulation

## 1 Introduction

Drug-related treatments can have severe consequences and toxic effects on the heart muscle. Among such therapies, it has been recently shown that chemotherapy-induced cardiotoxicity prevails as a critical clinical burden by severely altering cardiac electrical and mechanical functions [1]. One of the most widely used
chemotherapeutic drugs in cancer patient is doxorubicin (DOX), which, although effective in killing cancerous cells, has been recently demonstrated to induce structural myocardial damage through a gradual deposition of collagen-based re-active fibrosis within the heart wall [2]. This leads in time to a progressive heart dysfunction (e.g. non-ischemic cardiomyopathy and abnormal heart rhythms), culminating in heart failure and potentially lethal arrhythmic events [3, 4]. In particular, the abnormal propagation of electrical wave in sub-chronic phases (i.e., in the weeks after completion of DOX treatment) in the presence of diffuse fibrosis, is not well understood.

Computer simulations of electrophysiology are nowadays widely spread ro-
bust tools that can help us virtually dissect the propagation of electrical impulse
through the heart [5]. In particular, 3D image-based models can predict the be-
haviour of action potential (AP) wave under various experimental conditions,
while integrating high level of structural details [6]. Thus, such models could
be adapted and exploited to study the consequences and gradual (possible irre-
versible) changes in heart function following chemotherapy.

In this computational study, we sought to better understand how the action potential wave propagates in the presence of DOX-induces diffuse fibrosis. We first aim to use MR images acquired in a preclinical pig model of pathology before/after DOX treatment to build 3D anatomical and morphological models of the left ventricle (LV). We then employ those 3D LV models to perform several simulation-based parametric studies. Specifically, we are interested to learn how a gradual increase in the density of DOX-induced fibrosis, along with the anisotropy and regional heterogeneities in conductivity (diffusivity) and AP duration contribute to a perturbed electrical activity within the heart under pacing conditions. The novel contribution of this work lies in its comprehensive simulation of the interplay between fibrosis and electrical activity, providing insights into the progressive and potentially irreversible changes in heart function resulting from chemotherapy.

## 2 Methods

2.1 Preclinical experiment

The experiments involved n=2 juvenile Yorkshire swine weighing ~25kg prior
to starting the doxorubicin therapy, and were carried out with approval from
the research ethics board of the Animal Care Committee at our institution
(anonymous). Briefly, a 1mg/kg weekly dose of DOX was injected intravenously
throughout a four-week plan, mimicking clinical treatments in cancer patients.

Non-invasive CMR imaging was performed on a 3T MR750 scanner GE scan-
ner using a free-breathing 3D LGE method (late gadolinium enhancement). The
scans were done before the treatment (hereto referred as pre-DOX therapy time
point), as well as after the completion of treatment (post-DOX) at the following
time points: week 5 in pig 1, and week 9 in pig 2, respectively. The MR images
were acquired at 1.4 × 1.4 × 1.4 mm³ spatial resolution, using a conventional
3D inversion recovery fast gradient echo (IR-FGRE) pulse sequence and a bo-
lus of 0.2 mmol/kg Magnevist (a gadolinium-based contrast agent) sensitive to
collagenous fibrosis.

## 2.2 Left ventricular computational models from CMR images

In this work, we built isotropic (i.e., no fiber directions integrated into the do-
main) and anisotropic 3D left ventricular (LV) models from in vivo CMR images.
For each raw cine-MRI image, an expert initially performed manual contouring
of both the endocardial and epicardial surfaces of the LV, identifying the dias-
tolic phase as the most suitable for segmentation. Once the diastolic phase was
identified, a region of interest (ROI) from the posterior side of the remote my-
ocardium was selected. Then, fibrosis was identified by applying a threshold on
the intensity values on LGE images based on a clinically accepted 5SD (standard
deviation) cut-off [7]. For the anisotropic cases, synthetic fiber orientations were
assigned by linearly interpolating the helix angle from -70° on the epicardium
to +70° on the endocardium [8]. Figure 2.2 illustrates the entire workflow for
generating our computational LV models.

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Diagram of the pipeline: (left) generation of the 3D computational models from MR images with fibrosis and fiber directions integrated as well as exemplary AP waves parameterized in healthy and fibrotic zones; (right) simulation output showing activated tissue (red) within the 3D LV model and select AP waves at one voxel following a train of 5 stimuli.</figcaption></figure>

## 2.3 A modified Mitchell-Schaeffer model for LV electrophysiology

To simulate action potential propagation (AP) we used a monodomain formula-
tion of a modified Mitchell-Schaeffer (MS) model as in Djabella et. al. [9]. This
modification of the original MS model introduces an additional parameter, $\lambda$, that allows us to control the degree of excitability in a region of interest. The evolution equation of the system is given by:

$$ \frac{\partial u}{\partial t} = \nabla \cdot (D\nabla u) + \frac{hu(u-\lambda)(u_{max}-u)}{\tau_{in}} - \frac{u}{\tau_{out}} + J_{stim}(t), \quad (1) $$

$$ \frac{\partial h}{\partial t} = \epsilon(u, h)(h_{\infty}(u) - h), $$

where $u$ is the normalized AP variable, $h$ is the gating variable which makes the gate open and close (corresponding to depolarization and repolarization phases), $J_{stim}$ is the stimulation current at pacing location. The other parameters and their values can be found in [8]. To numerically solve this PDE system 1, we employed a fast GPU-based Lattice-Boltzmann numerical method (which runs directly on images) [10]. We validated its efficacy by comparing it to the finite element method using FEiCS on a 3D benchmark [11], demonstrating significantly faster computation times for large, high-resolution meshes, achieving a speed-up of 140x in the simulations. The results of the calibrations can be found in the **Supplemental Material**.

To account for tissue anisotropy, propagation speed was set to be 2.5 times slower in the direction transverse to myocardial fibers; thus, we used $d = 3.5$ for the diffusion coefficient parallel to fiber directions and $d/2.5$ for the transverse and normal directions to fibers. This gives us an estimate of a conduction velocity of 0.8 m/s along the fiber direction in healthy tissue. Under these considerations, our diffusion coefficient matrix is defined as $D = \Lambda A \Lambda^T$ with

$$ A = \begin{pmatrix} d & 0 & 0 \\ 0 & \frac{d}{2.5} & 0 \\ 0 & 0 & \frac{d}{2.5} \end{pmatrix}, $$

and $\Lambda$ having the local basis information (fiber, transverse and normal directions).

Simulations were carried out on NVidia GTX 1080 ti (11 GB of RAM) GPU.

**Stimulation protocols and experiments description** For both isotropic and anisotropic cases, a series of 5 beats at 500 ms were performed from a selected point in the LV free wall. On average, the whole protocol took around 20 minutes per simulation. As our goal is to understand the effect of structure on electrophysiological function, the parameters concerning the shape of the APD were fixed in all our simulations: ($\tau_{open}$ = 120, $\tau_{close}$ = 150, $\tau_{in}$ = 0.3, $\tau_{out}$ = 6, $u_{max}$ = 1, $u_{gate}$ = 0.13, all in milliseconds with the exception of both $u_{max}$ and $u_{gate}$ which are dimensionless). In addition, a fast pacing protocol was performed with modifying the APD to assess whether this could lead to wavefront fragmentation. The protocol consisted of adding 2 pulses at 300ms and 350ms consecutively, with the same values for the model parameters except $\tau_{close}$, which was set to 250 milliseconds.
*Fibrosis modeled via reduced excitability* In this scenario, for the part of the domain exhibiting fibrosis we modify the excitability parameter ($\lambda = 0.2$) to decrease the propagation speed in these regions. Additionally, diffusion in both the longitudinal and transverse directions was reduced in fibrotic areas to 1/4th of the one in normal tissue [13, 14].

*Fibrosis modeled as non-conductive tissue* In this case, we remove the voxels of MR-defined fibrosis from the ventricular domain. Consequently, no computations were performed in the fibrotic region, essentially creating clusters of tissue acting as electrical insulators.

## 3 Results and discussion

We next present the simulation results obtained after varying different key parameters separately (i.e., anisotropy, fibrosis density, tissue diffusivity, and APD), in order to observe their independent influence on the perturbation of wavefront propagation and activation heterogeneity.

Figure 2 illustrates clear differences in the pattern of activation and AP wave propagation immediately following stimulation and thereafter, in the LV models from both pigs when the fiber directions were integrated (i.e. anisotropic cases) or not (isotropic cases), into the 3D computational domains. For instance, for pig 1, as the speed of wave travels 2.5 faster in the longitudinal direction (parallel to the muscle fibers), the initial activation (50ms) is ellipsoidal in shape and, as a result, the wavefront reached the base at a later time point (200ms) compared to the isotropic scenario where it reached the base at 150ms. Such differences are important for realistically reproducing the AP wave spread through the heart, underscoring the need to use anisotropic model. Therefore, in the following sections, we restrain from using isotropic models.

Figure 3 presents the simulation results obtained for both pre-DOX and post-DOX at weeks 5 (pig 1) as well as pre-DOX and post-DOX at weeks 9 (pig 2). Most importantly, by week 9, fibrosis deposition had progressed substantially, which resulted in a noticeable wave perturbation due to larger obstacles. The most significant differences were observed at 9 weeks post-DOX (pig 2), where excitability ($\lambda = 0.2$) and diffusivity (1/4 of the healthy value) were reduced.

Figure 4 shows the temporal evolution of the electrical activation of the ventricular volume, with the time steps represented in milliseconds after the last stimulus of a 5-pulse train. In the first row, the propagation of the action potential in the pre-Doxo case is illustrated. We observed that in this case, the propagation proceeded normally, with a continuous wavefront, and the LV is fully activated in the absence of fibrosis (pre-DOX).

The next two rows correspond to the simulation of fibrosis treated as less excitable tissue ($\lambda = 0.2$) and as a non-conducting tissue (d=0), respectively. In the simulation of fibrosis as less excitable tissue, the wavefront is observed to fragment after the stimulus (t=50 ms). The propagation throughout the ventricle is more discontinuous than in the pre-DOX case, and, when the ventricle is
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** AP wave propagation at different time points in the LV MRI-based models built from pig 1 (at week 5) and pig 2 (at week 9), with fiber directions (anisotropy) integrated and without fiber directions (isotropic case). Red indicates depolarized tissue and blue corresponds to the repolarized phase.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** The effect of fibrosis density on AP wave propagation for both pigs: (left) 3D anatomical LV models with fibrosis zone integrated; (right) simulated cardiac activation maps pre-/post-DOX presented in short-axis corresponding to the cutting planes.</figcaption></figure>
fully activated, the fibrotic zones remain inactive. In addition, the fibrotic zone repolarizes earlier than in pre-DOX, showing more heterogeneities in the APDs. To study the differences between cases, error maps of the differences have been computed and quantified, as detailed in the **Supplemental Material**.

For the simulation case where fibrosis was treated as nonconducting tissue, the wavefront also fragments after the stimulus, but in this case the fibrosis acts as a functional block without introducing further heterogeneity in the wavefront propagation.

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Comparison of simulations obtained for pig 2: pre-DOX (top row) vs post-DOX at w9 in two cases: when the fibrotic tissue is treated with a modification of $\lambda$ (middle row); and when fibrotic region treated as non-conductive tissue (bottom row).</figcaption></figure>

Figure 5 shows the simulation results for pig 2 after the introduction of changes in the APD. In the second row (APD increased in the fibrotic zone), we observed that the propagation wavefront fragments into two after the stimulus ($t=100\text{ms}$). This is not seen in the first row, which corresponds to the pre-DOX simulation without fibrotic zone.

The fragmentation of the wavefront into two waves may be indicative of tissue arrhythmicity, as the introduction of these structural changes together with an early stimulus may result in a unidirectional block that triggers a reentrant arrhythmia.

# 4 Conclusion and future work

Doxorubicin administration results in a gradual deposition of reactive fibrosis. Such type of fibrosis was identified in LV models built from LGE imaging data
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Simulation of fast pacing. Case without fibrosis and w9 with increased APD in the fibrosis zone (epicardial and short-axis views).</figcaption></figure>

acquired in swine after DOX therapy. Our computational study showed that several key parameters such as: fibrosis density, tissue conductivity, action potential duration, and anisotropy, independently contribute to a local perturbation of electrical wave propagation in the presence of DOX-induced fibrosis.

The use of the modified Lattice-Boltzmann method (LBM), based on GPU, has allowed simulations with high computational speed. This technique offers remarkable flexibility in the calibration of the parameters, which facilitates their adjustment to real cases. This robust method will enable us to efficiently perform protocols for studying arrhythmia inducibility, which in future lines will allow an exhaustive study of doxorubicin-induced cardiotoxicity. Our results highlight the ability of LBM to investigate structure-function interactions in hearts with chemotherapy-induced fibrosis, providing a valuable tool for the study of potentially lethal arrhythmias in hearts undergoing DOX therapy.

**Acknowledgments.** This work was supported by the SimCardioTest project, which has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No. 101016496; by the French government through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002; by a Canadian CIHR project grant (PJT) 153212; by Grant PRE2020-091849 [MCIN/AEI/10.13039/501100011033] and "ESF Investing in your future"; Grant PID2019-104356RB-C41 [MCIN/AEI/10.13039/501100011033]; PID2022-136273OA-C33 and PID2022-140553OB-C41 [MICIU/AEI/10.13039/501100011033] and by ERDF/EU]; and by the Barcelona Supercomputing Center [IM-2021-1-0001 and IM-2021-3-0001]. The authors are grateful to the OPAL infrastructure from Université Côte d'Azur for providing resources and support.
**Disclosure of Interests.** The authors have no competing interests to declare that are relevant to the content of this article.

References

1. Von Hoff, DD; Rozencweig, M; Piccart, M. The cardiotoxicity of anticancer agents. Semin. Oncol. 1982, 9, 23–33.

2. Tanaka R, Umemura M, et al; Reactive fibrosis precedes doxorubicin-induced heart failure through sterile inflammation. ESC Heart Fail. 2020 Apr;7(2):588-603.

3. Christiansen, S.; Autschbach, R. Doxorubicin in experimental and clinical heart failure. European. J. Cardiothorac. Surg. 2006, 30, 611–616.

4. Herrmann, J. Adverse cardiac effects of cancer therapies: Cardiotoxicity and arrhythmia. Nat. Rev. Cardiol. 2020, 17, 474–502.

5. M Sermesant 1, H Delingette, N Ayache; An electromechanical model of the heart for image analysis and simulation; IEEE Trans Med Imaging, 2006 May; 25(5):612-25.

6. Balaban G, Costa CM, et al; 3D Electrophysiological Modeling of Interstitial Fibrosis Networks and Their Role in Ventricular Arrhythmias in Non-Ischemic Cardiomyopathy. IEEE Trans Biomed Eng. 2020 Nov; 67(11), 3125-3133.

7. Karim R, Bhagirath P, Claus P, James Housden R, Chen Z, Karimaghaloo Z, Sohn HM, Lara Rodríguez L, Vera S, Albà X, Hennemuth A, Peitgen HO, Arbel T, González Ballester MA, Frangi AF, Götte M, Razavi R, Schaeffter T, Rhode K. Evaluation of state-of-the-art segmentation algorithms for left ventricle infarct from late Gadolinium enhancement MR images. Med Image Anal. 2016 May;30:95-107.

8. Pop, M., Sermesant, M., Liu, G., Relan, J., Mansi, T., Soong, A., ... & Wright, G. A. (2012). Construction of 3D MR image-based computer models of pathologic hearts, augmented with histology and optical fluorescence imaging to characterize action potential propagation. Medical image analysis, 16(2), 505-523.

9. Djabella, K., Landau, M., & Sorine, M. (2007, December). A two-variable model of cardiac action potential with controlled pacemaker activity and ionic current interpretation. In 2007 46th IEEE Conference on Decision and Control (pp. 5186-5191). IEEE.

10. Cedilnik, N., Pop, M., Duchateau, J., Sacher, F., Jaïs, P., Cochet, H., & Sermesant, M. (2023). Efficient Patient-Specific Simulations of Ventricular Tachycardia Based on Computed Tomography-Defined Wall Thickness Heterogeneity. Clinical Electrophysiology, 9(12), 2507-2519.

11. Niederer SA, Kerfoot E, Benson AP, Bernabeu MO, Bernus O, Bradley C, Cherry EM, Clayton R, Fenton FH, Garny A, Heidenreich E, Land S, Maleckar M, Pathmanathan P, Plank G, Rodríguez JF, Roy I, Sachse FB, Seemann G, Skavhaug O, Smith NP. Verification of cardiac tissue electrophysiology simulators using an N-version benchmark. Philos Trans A Math Phys Eng Sci. 2011 Nov 13;369(1954):4331-51.

12. Lopez-Perez, A., Sebastian, R., Izquierdo, M., Ruiz, R., Bishop, M., Ferrero, J. M. (2019). Personalized Cardiac Computational Models: From Clinical Data to Simulation of Infarct-Related Ventricular Tachycardia. Frontiers in Physiology, 10.

13. Gomez JF, Cardona K, Martinez L, Saiz J, Trenor B. Electrophysiological and structural remodeling in heart failure modulate arrhythmogenesis. 2D simulation study. PLoS One. 2014 Jul 23;9(7):e103273.

14. Mendonca Costa, C., Plank, G., Rinaldi, C. A., Niederer, S. A., & Bishop, M. J. (2018). Modeling the Electrophysiological Properties of the Infarct Border Zone. Frontiers in Physiology, 9.