# **Simula SpringerBriefs on Computing 12**

Reports on Computational Physiology

**Kimberly J. McCabe** *Editor*

# **Computational Physiology** Simula Summer School 2021 − Student Reports

# **Simula SpringerBriefs on Computing**

# Volume 12

#### **Editor-in-Chief**

Aslak Tveito, Simula Research Laboratory, Fornebu, Norway

#### **Series Editors**

Are Magnus Bruaset, Simula Research Laboratory, Fornebu, Norway

Kimberly Claffy, San Diego Supercomputer Center, CAIDA, University of California, San Diego, San Diego, CA, USA

Magne Jørgensen, Software Engineering, Simula Research Laboratory, Fornebu, Norway

Olav Lysne, Simula Research Laboratory, Fornebu, Norway

Andrew McCulloch, Bioengineering 0412, University of California, San Diego, La Jolla, CA, USA

Fabian Theis, Institute of Computational Biology, Helmholtz Zentrum München, Neuherberg, Germany

Karen Willcox, Department of Aeronautics & Astronautics, Massachusetts Institute of Technology, Cambridge, MA, USA

Andreas Zeller, Saarbrücken, Germany

Springer and Simula have launched a new book series, *Simula SpringerBriefs on Computing*, which aims to provide introductions to select research in computing. The series presents both a state-of-the-art disciplinary overview and raises essential critical questions in the field. Published by SpringerOpen, all *Simula SpringerBriefs on Computing* are open access, allowing for faster sharing and wider dissemination of knowledge.

Simula Research Laboratory is a leading Norwegian research organization which specializes in computing. The book series will provide introductory volumes on the main topics within Simula's expertise, including communications technology, software engineering and scientific computing.

By publishing the *Simula SpringerBriefs on Computing*, Simula Research Laboratory acts on its mandate of emphasizing research education. Books in this series are published only by invitation from a member of the editorial board.

More information about this series at https://link.springer.com/bookseries/13548

Kimberly J. McCabe Editor

# Computational Physiology

Simula Summer School 2021 − Student Reports

*Editor*  Kimberly J. McCabe Simula Research Laboratory Oslo, Norway

ISSN 2512-1677 ISSN 2512-1685 (electronic) Simula SpringerBriefs on Computing ISBN 978-3-031-05163-0 ISBN 978-3-031-05164-7 (eBook) https://doi.org/10.1007/978-3-031-05164-7

Mathematics Subject Classification (2020): 65M60, 65M06, 92C10, 92C37

© The Editor(s) (if applicable) and The Author(s) 2022. This book is an open access publication. **Open Access** This book is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this book are included in the book's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the book's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

The use of general descriptive names, registered names, trademarks, service marks, etc. in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant protective laws and regulations and therefore free for general use.

The publisher, the authors and the editors are safe to assume that the advice and information in this book are believed to be true and accurate at the date of publication. Neither the publisher nor the authors or the editors give a warranty, expressed or implied, with respect to the material contained herein or for any errors or omissions that may have been made. The publisher remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This Springer imprint is published by the registered company Springer Nature Switzerland AG The registered company address is: Gewerbestrasse 11, 6330 Cham, Switzerland

## **Preface**

Since 2014, we have organized an annual summer school in computational physiology. The school starts in June each year and the graduate students spend two weeks in Oslo learning the principles underlying mathematical models commonly used in studying the heart and the brain. At the end of their stay in Oslo, the students are assigned a research project to work on over the summer. In August the students travel to the University of California, San Diego to present their fndings. Each year, we have been duly impressed by the students' progress and we have often seen that the results contain the rudiments of a scientifc paper.

Starting in the 2021 edition of the summer school, we have taken the course one step further and aim to conclude every project with a scientifc report that passes rigorous peer review as a publication in this new series called *Simula SpringerBriefs on Computing – reports on computational physiology*.

One advantage of this course adjustment is that we have the opportunity to introduce students to scientifc writing. To ensure the students get the best introduction in the shortest amount of time, we have commissioned a professional introduction to science writing by Nature. The students participate in a 2-day *Nature Masterclasses* workshop, led by two editors from Nature journals, in order to strengthen skills in high quality scientifc writing and publishing. The workshop is tailored to publications in the feld of computational physiology and allows students to gather real-time feedback on their reports.

We would like to emphasise that the contributions in this series are brief reports based on the intensive research projects assigned during the summer school. Each report addresses a specifc problem of importance in physiology and presents a succinct summary of the fndings (8-15 pages). We do not require that results represent new scientifc results; rather, they can reproduce or supplement earlier computational studies or experimental fndings. The physiological question under consideration should be clearly formulated, the mathematical models should be defned in a manner readable by others at the same level of expertise, and the software used should, if possible, be made publicly available. All reports in this series are subjected to peer-review by the other students and supervisors in the program.

We would like to express our gratitude for the very fruitful collaboration with Springer -Nature and in particular with Dr. Martin Peters, the Executive Editor for Mathematics, Computational Science and Engineering.

The editors of *Simula SpringerBriefs on Computing – reports on computational physiology*:

Oslo, Norway *Kimberly J McCabe* March 2022 *Rachel Thomas Andrew D McCulloch Aslak Tveito*

### **Acknowledgements**

The Simula Summer School in Computational Physiology is a team efort, with many scientists contributing their time to give lectures and advise projects for the students. The 2021 school was particularly challenging, as it was conducted entirely online and required creative problem solving to allow for efective project collaboration, student engagement, and knowledge transfer. We would like to thank the lecturers and project advisors for their expertise and willing participation in the course: Dr. Hermenegild Arevalo, Dr. Gabriel Balaban, Dr. Jonas van den Brink, Marius Causemann, Dr. Andrew Edwards, Dr. Nickolas Forsch, Dr. Ingeborg Gjerde, Dr. Glenn Lines, Dr. Molly Maleckar, Dr. Kimberly McCabe, Dr. Andrew McCulloch, Denis Reis de Assis, Dr. Joakim Sundnes, Abigail Teitgen, Dr. Aslak Tveito, Dr. Daniela Valdez-Jasso, Jonas Verhellen, and Dr. Samuel Wall. Administrative support for the school was provided by Dr. Kimberly McCabe, Dr. Rachel Thomas, Elisabeth Andersen, and Hanie Tampus.

We would also like to acknowledge the teaching assistants for the course, who provided much-needed support during virtual lectures: Åshild Telle, Eina Jørgensen, and Hanna Borgli. We received training and equipment relating to virtual learning from Simula Kodeskolen, and would specifcally like to acknowledge Matteus Häger, Håkon Kvale Stensland, Elin Backe Christophersen, and Marianne Aasen.

The Simula Summer School in Computational Physiology is supported through the Simula-UiO-UCSD Research PhD Training Programme (SUURPh), an endeavour funded by the Norwegian Ministry of Education and Research. Additional fnancial support is derived from SIMENT, an INTPART mobility grant from the Norwegian Research Council. The school also received funding from Digital Life Norway (DLN).

# **Contents**



#### Contents xi


# **Chapter 1 A Pipeline for Automated Coordinate Assignment in Anatomically Accurate Biventricular Models**

Lisa Pankewitz<sup>1</sup> , Laryssa Abdala<sup>2</sup> , Aadarsh Bussooa<sup>1</sup> , Hermenegild Arevalo<sup>1</sup>


**Abstract** There is an increased interest, in the feld of cardiac modeling, for an improved coordinate system that can consistently describe local position within a heart geometry across various distinct geometries. A newly designed coordinate system, Cobiveco, meets these requirements. However, it assumes the use of biventricular models with a fat base, ignoring important cardiac structures. Therefore, we extended the scope of this state-of-the-art biventricular coordinate system to work with various heart geometries which include basal cardiac structures that were previously unaccounted for in Cobiveco. First, we implemented a semi-automated input surface assignment for increased accessibility and reproducibility of assigned coordinates. Then, we extended the coordinate system to handle more anatomically accurate biventricular models including the valve planes, which are of great interest when modeling diseases that manifest themselves in the basal area. Furthermore, we added the functionality of mapping vector data, such as myocardial fber orientations, which are crucial for replicating the anisotropic electrical propagation in cardiac tissue.

#### **1.1 Introduction**

The representation of cardiac geometry independent of patient origin and the fawless transfer between diferent measuring modalities are important tools in clinical research [1, 2]. To accurately describe a local position within the heart, a robust coordinate system is required. Such a coordinate system enables a variety of applications, including the transfer of data between diferent heart geometries and comparing data produced using diferent measuring modalities, such as validating simulations with clinical data [2, 3].

A recently published biventricular coordinate system, Cobiveco, ofers a consistent and reliable approach for describing positions in biventricular heart models [2]. However, the current state-of-the-art coordinate system is limited to biventricular heart geometries, which are clipped at a specifc planar position, such that the resulting base appears completely fat. Clipping the base in this manner also clips the underlying ventricles. Although this clipping procedure remains part of the common mesh generation approach, it does not yield anatomically accurate cardiac meshes for the purpose of computer simulations. We argue that biophysical simulations of the heart should include the clipped base cardiac structures, that contain the valve openings, for more realistic results. The ventricles, together with the presence of valve planes, are important features of ventricular anatomy that can infuence cardiac electrophysiology and mechanics.

Features that are connected to the valves are critical anatomical structures, such as the papillary muscle and chordae. Any structural defects that change the shape of the ventricles and alter the activation in the aortic valve annulus can have an efect on electrical dyssynchrony or ventricular dilation. Therefore, the inclusion of valve planes in cardiac models is necessary. This is especially true when modeling certain disease phenotypes, where changes in anatomy, mechanics, and activation manifest themselves in areas closer to the valve planes. An important example of this is congenital heart defect (CHD), which is the most common birth defect worldwide [4, 5]. Heterogeneous morphology and physiology in CHD patients have been shown to complicate risk assessment of individual patients requiring anatomically accurate models. This is a use-case where the inclusion of valve planes in the biventricular models may lead to enormous improvement of the model quality as morphological changes as well as scar tissues in this patient group can be located close to the base.

#### **1.2 Methods**

In this work, we extend the open-source MATLAB implementation of Cobiveco for tetrahedral meshes to take into account anatomically more accurate biventricular meshes that include valve planes instead of a fat, clipped base. First, we provide a surface extraction tool that automatically creates input surfaces fles required for setting up the biventricular coordinate system. Then, we adapt the existing Cobiveco framework to allow for more anatomically correct geometries. Last, we extend the software to allow for mapping and transfer of vector data between diferent heart geometries.

#### **1.2.1 Semi-Automated Surface Extraction**

Existing tools extract surfaces from meshes and imaging data. Image-based surface extraction operates directly on raw clinical imaging to identify cardiac structures, efectively building a mesh with automated tagging of surfaces. Mesh-based surface extraction operates directly on the meshes and identifes cardiac structures based on the position and connectivity of vertices. However, these currently existing tools require much fne-tuning and cannot efortlessly extract surfaces based solely on a seed point and a threshold.

Therefore, we present a mesh-based surface extraction tool, which identifes cardiac structures using a minimal set of parameters. As a frst step, the cardiac mesh is converted into a graph, where its nodes encode vertex identifers and surface identifers. We leverage the use of the graph topology and apply a breadth-frst search (BFS) algorithm to fnd connected nodes.

Fig. 1.1: Angular change between two neighbouring triangular surfaces, given by .

The scope of the BFS algorithm is limited by two parameters, namely a seed point lying on the surface to be extracted and an angular threshold. The BFS algorithm performs several iterations, starting with the seed point. With each iteration, the angular change between two neighbouring triangles is computed, such that they comply with the stated threshold.

$$\cos \theta = \frac{n\_1 \cdot n\_2}{|n\_1||n\_2|} \tag{1.1}$$

The angular change is given in (1.1), where <sup>1</sup> and <sup>2</sup> are the normals to two triangular surfaces (Figure 1.1). These two triangles do not have to be direct neighbours to each other. The BFS algorithm identifes neighbours in the vicinity, using a predefned depth variable, such that the overall algorithm achieves a faster execution time. The pseudo code of the implementation is given in Algorithm 1.


**Algorithm 1** Algorithm to identify connected mesh nodes

The algorithm mimics edge detection, as used in image processing, where in our case the BFS spans throughout the mesh until sharp corners are encountered. Based on heuristics, we have identifed that an angular threshold between 0.1 and 0.2 rad is optimal in correctly identifying and extracting cardiac structures in the mesh. In the context of extending the features of Cobiveco, the algorithm can be used to extract the base of the valves and the valve plane on the epicardial surface. The extracted base is excluded from the graph and subsequent surface extraction freely applies the BFS algorithm without any angular threshold restrictions.

#### **1.2.2 Biventricular Coordinate System**

Cobiveco, a **co**nsistent **bive**ntricular **co**ordinate system, provides a reliable framework for the precise and intuitive description of the position in the heart. To consistently describe a location within the heart, the coordinate system is established with four coordinates. The coordinate system fulfls a set of desired properties, which have been described in Section 1.1. The system is based on a set of four coordinates, namely a transventricular coordinate (tr), a transmural coordinate (tm), a rotational coordinate (rt) and an apicobasal coordinate (ab). The transventricular coordinate is a binary coordinate which distinguishes between the left and right ventricle. The transmural coordinate measures the distance traveled within the transmural space, so in the free walls this refers to the distance from the epicardium to the endocardium. The rotational coordinate gives information about where you are in the heart with respect to anterior and posterior direction. In more detail, it refers to the distance traveled from the interventricular posterior junction over the the interventricular anterior junction over the septum back to the interventricular posterior junction. The rotational coordinate is set up symmetrical in the biventricular model. The apicobasal coordinate describes the distance traveled from the apex point to the base. Each coordinate tuple, consisting of the four coordinates, corresponds to exactly one point in the heart. Coordinates are normalized and range between 0 and 1. Within that range, all coordinates change linearly in space, indicating that the distance traveled is directly proportional to the change in the coordinate of interest. Furthermore, both ventricles follow the same parametrization. This is also refected in the shared apex defnition. To construct the coordinate system, only landmarks which are consistent throughout variations in diferent geometries, are chosen.

#### **1.2.2.1 Creation of the Coordinate System Cobiveco**

A detailed description of the steps involved in the creation of the original Cobiveco framework can be found in [2]. In short, the construction of the coordinate system includes eight steps and is summarized below.

As with Cobiveco 1.0, Cobiveco 2.0 requires a biventricular volume that includes a base containing the four heart valve annuli, including the connecting bridges. Besides the volume mesh, fve boundary surfaces as shown in Figure 1.3 are required as input, which is one additional surface compared to Cobiveco 1.0. The surfaces required are a *basal surface SBase*, a *basal epicardial surface SEpi,base*, an *epicardial, non basal surface SEpi, nonbase*, an *LV endocardial surface SLV* and an *RV endocardial surface SRV*. The utilities for the semi-automated input-fle generation are described in Section 1.2.1.

#### **Transventricular Coordinate (tv)**

The transventricular coordinate is calculated as described in the original publication [2].

#### **Extraction of Septal Surface and Curve**

The *septal surface SSept* and the *septal curve CSept* are extracted as described in [2].

#### **Transmural Coordinate (tm)**

The calculation of the transmural coordinate follows the same steps as in Cobiveco 1.0, but takes into account the two epicardial surfaces. As we split the epicardial surface into a non-base epicardial surface and a basal epicardial surface, the whole epicardial surface is defned be the union of both, as given in (1.2):

$$\text{LS}\_{\text{Epi}} = \text{S}\_{\text{Epi}\_{=} \text{non\\_base}} \cup \text{S}\_{\text{Epi}\_{=} \text{base}} \tag{1.2}$$

#### **Heart Axes and Apex Point**

The defnition of the heart axes and apex point mainly follows the steps described in the original publication [2]. As the defnition of the orthogonal heart axes largely depends on the truncated septal surface, the calculation of the truncation needed to be revised to take into consideration the increased curvature of the septal surface at the base of the anatomically accurate biventricular heart model. Therefore the septal surface is additionally truncated by 15% at the basal side, where the distance is based in the direction of **vLongAx**. This modifcation results in the fnal truncated septal surface *SSeptTrunc* being calculated by (1.3), where P*<sup>q</sup>* refers to the *q th* percentile.

$$S\_{\rm SepTrunc} = \left\{ \mathbf{x} \in S\_{\rm Sept} \mid \mathbf{x} \cdot \mathbf{v}\_{\rm AP} > \mathbf{P}\_{20} \left( \mathbf{x} \cdot \mathbf{v}\_{\rm AP} \right) \text{ and } \tag{1.3}$$

$$\mathbf{x} \cdot \mathbf{v}\_{\rm AP} < \mathbf{P}\_{90} \left( \mathbf{x} \cdot \mathbf{v}\_{\rm AP} \right) \right\} \text{ and }$$

$$\mathbf{x} \cdot \mathbf{v}\_{\rm long} > \mathbf{P}\_{15} \left( \mathbf{x} \cdot \mathbf{v}\_{\rm long} \right) \}$$

Since the *septal curve* needs to be split in two segments, the new geometry with a closed base requires a diferent solution than in Cobiveco 1.0 as well. Hence, we exclude the new, *basal epicardial surface* from the epicardial defnition to allow for a separation of the anterior and posterior part of the septal curve.

$$\mathcal{C}\_{\text{Sept}} = \left\{ \mathbf{x} \in \mathcal{S}\_{\text{Epi}, \text{base}} \mid \boldsymbol{\mu}\_{\boldsymbol{\nu}}(\mathbf{x}) = \mathbf{0}.\mathbf{5} \right\} \tag{1.4}$$

#### **Extraction of Ridge Surfaces**

As the more anatomically accurate biventricular models contain a closed surface at the base, we needed to modify the ridge defnition in Cobiveco 2.0. In Cobiveco 2.0 we aim to replicate the original ridge assignment but use the valve planes as guiding points resulting in a symmetric ridge. The ridge is used to provide a boundary condition for the rotational coordinate. Currently, the ridge is defned from the posterior interventricular junctions, where both ventricles symmetrically impose a boundary condition, via the mitral valve and tricuspid valve in the LV/RV septum respectively. The anterior ridge defnition is defned via the mitral valve and pulmonary valve, resulting in the ridge defnition as shown in Figure 1.2.

Fig. 1.2: Ridge defnition in Cobiveco 1.0 and Cobiveco 2.0. The anterior part of the ridge is colored in red, while the posterior part of the ridge is highlighted in grey.

#### 1.2 Methods 7

Therefore, the ridge is obtained by defning the solution to Laplace's equation with boundaries applied to be 0 on the epicardial basal surface and 1 on the septal surface. The *non-base epicardial surface* is now excluded as a boundary condition, as shown in (1.5).

$$
\mu\_{\text{Ridge}}\left(V\right) = 0 \quad \text{with} \quad \mu\_{\text{Ridge}}\left(S\_{\text{Epi}} \backslash S\_{\text{Sept}}\right) = 0 \quad \text{and} \quad \mu\_{\text{Ridge}}\left(S \mid\_{\text{Sept}}\right) = 1 \tag{1.5}
$$

The solution to Laplace's equation is calculated as shown here in (1.6).

$$
\Delta\mu\_{\text{Ridge}}(V) = 0 \quad \text{with} \quad \mu\_{\text{Ridge}}\left(S\_{\text{Epi,nonbase}} \backslash S\_{\text{Setpt}}\right) = 0 \quad \text{and} \quad \mu\_{\text{Ridge}}\left(S \mid\_{\text{Setpt}}\right) = 1 \tag{1.6}
$$

The second step remains as set up in the original Cobiveco 1.0. However, the resulting ridge cannot be applied as it is. Currently, a manual fltering step is involved as there remains a ridge within the RV in-between the tricuspid and the pulmonary valve, which is not used as a boundary condition when setting up the rotational coordinate.

#### **The Rotational Coordinate (r)**

The rotational coordinate is defned as described in [2].

#### **Computation of the Apicobasal Coordinate (ab)**

The apicobasal coordinate is calculated as described in the original Cobiveco article. Currently, however, we used the solution to Laplace's equation as a place holder, as the rotational coordinate still includes discontinuities which prohibit the assignment of the apicobasal coordinate as described in Cobiveco 1.0

#### **1.2.3 Mapping Vector Fields**

The original Cobiveco implementation has a scalar feld mapping functionality available. To map a scalar feld from the source mesh to a target mesh , it constructs a matrix ←− from the nodes of the source mesh to the nodes of the target mesh [2]. The user can choose between linear and nearest-neighbor interpolation.

Mapping vector felds is of interest since data, such as muscle fber felds, are crucial to advance the cardiovascular computational simulations feld. Here we enable the functionality of mapping such felds by treating each coordinate as a scalar feld. More specifcally, the vector feld is represented as a matrix of the nodes of the source mesh by three. Each of its columns represents the coordinate of the vector feld in each source node. The end result of the vector mapping process is shown in Figure 1.5.

#### **1.3 Results**

In this project, we have successfully founded the basis for extending Cobiveco to include more anatomically accurate biventricular models. The results are presented in three steps, namely a pre-processing, processing, and post-processing step. Each step reduces the manual manipulation of meshes, generalizes the biventricular coordinates, and enables vector data transfer, respectively.

First, we successfully implemented a semi-automated surface extraction method that uses a minimal set of parameters, based on a seed point and angular threshold, to identify structures of interest in cardiac meshes. The resulting surfaces, after extraction, are shown in Figure 1.3.

Fig. 1.3: Extracted surfaces after using BFS algorithm and angular threshold.

Second, we adapted the previous Cobiveco framework to work with biventricular geometries, which includes the four cardiac valve planes. The preliminary results of Cobiveco 2.0 are shown in Figure 1.4. The rotational coordinate sufers from inconsistencies at the septum, owing to the manual exclusion of the ridge boundary. Currently, the apicobasal coordinate is only represented by the solution to Laplace's equation. Last, we adapted the framework to include the mapping of vector data. The result for mapping synthetic data is shown in Figure 1.5.

#### **1.4 Conclusion**

In this project, we present an updated version of the consistent biventricular coordinates introduced by [2]. The pipeline can be applied to biventricular geometries for mapping scalar and vector data between diferent hearts.

Cobiveco 2.0 builds upon the original Cobiveco [2], by extending the coordinates for biventricular geometries that include the ventricular base. We aim to keep the

Fig. 1.4: Visual Comparison of the four coordinates created by Cobiveco 2.0 and Cobiveco 1.0. The apicobasal coordinate shown for Cobiveco 2.0 is represented by the solution to Laplace's equation.

Fig. 1.5: Cobiveco 2.0 has the functionality of mapping vector felds. First, the coordinates are built in the source and target biventricular geometries (top). Then the map ←− is used to map vector felds (bottom).

resulting coordinates with the same properties from the original ones: bijective, continuous (apart from the binary transventricular coordinates), normalized, complete, linear, with consistent parametrization, and consistent landmarks.

The new pipeline reduces manual mesh manipulations for surface extraction. The required surfaces can be efortlessly extracted using fewer parameters than other conventional methods. Cobiveco 2.0 enables the mapping of vector data in addition to scalar data, which is useful for computational modeling and data comparison. A remarkable application of this feature is myofber data mapping which is widely used in electrophysiology simulations. Moreover, by enabling data transfer between anatomical accurate biventricular geometries, it will be possible to validate computational models for diseases that manifest themselves in areas close to the valve planes.

The newly developed pipeline, with the inclusion of the valve planes in the cardiac model, is of special interest for studying the most common form of CHD, namely Tetralogy of Fallot (ToF). In ToF patients, the scar tissue is located close to the base, rendering electrophysiological simulations feasible with our pipeline [6, 7, 8, 9, 10, 11].

To the best of our knowledge, there is no consistent coordinate system available that can be readily applied to a four-chamber heart model. The current state-ofthe-art coordinates for the atria, Universal Atrial Coordinates (UAC) is not based on circular coordinates but rather in lateral-septal and posterior-anterior coordinates [12]. The two-dimensional framework can be extended to three dimensions by adding a transmural coordinate. Therefore, merging our improved Cobiveco pipeline with an updated version of the UAC could be used to create a consistent four-chamber heart coordinate system.

#### **1.4.1 Limitations**

The limitations of the current work include incomplete assignment of the rotational and apicobasal coordinates. Therefore, we aim to improve the ridge assignment to obtain a more symmetric ridge, which will result in a more symmetric rotational coordinate, as the ridge defnes the boundaries set in the rotational coordinate. To achieve this goal, we will modify the anterior part of the ridge in the LV to be defned via the aortic valve and not via the mitral valve as the defnition is set now. This will ensure a symmetric set up of the rotational coordinate in the RV and LV. The remaining parts of the proposed framework, however, were successfully tested in one patient-specifc geometry.

Besides, a statistical analysis of the errors using a cohort of geometries is necessary to ensure this is a reliable coordinate system. Moreover, the mappings were performed using artifcial data. Future work could include transferring experimental data between two geometries to ensure the results are physiologically consistent. Furthermore, the post-processing could also feature tensor data mapping.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Chapter 2 3D Simulations of Fetal and Maternal Ventricular Excitation for Investigating the Abdominal ECG**

Julie Johanne Uv<sup>1</sup> , Lena Myklebust<sup>1</sup> , Hamid Khoshfekr Rudsari2,3, Hannes Welle<sup>4</sup> , Hermenegild Arevalo<sup>1</sup>


**Abstract** Congenital heart disease (CHD) is a leading cause of infant death. To diagnose CHD, recordings from abdominal fetal electrocardiograms (fECG) can be used as a non-invasive tool. However, it is challenging to extract the fetal signal from fECG recordings partly due to the lack of data combining fECG recordings with a ground truth for the fetal signal, which can be obtained by using a scalp electrode during delivery. In this study, we present a computational model of a pregnant female torso, in which we simulate fetal and maternal ventricular excitation during sinus rhythm to derive fECGs, so as to enable isolated measurement of the fetal and maternal signal contributions. To extract the fetal contribution from a combined signal, we apply an adaptive fltering algorithm to wavelet transformed signals. Further development of the model may enable improvements in the recording and processing capabilities for fECGs, the reliable estimation of fetal heart rates, and possibly interpretation of fetal signal morphologies that could improve the overall diagnostic signifcance of abdominal fECGs.

#### **2.1 Introduction**

Congenital heart disease (CHD) is a common birth defect referring to abnormal function or structure of the heart [1]. It arises in early pregnancy stages during heart development [2]. With a prevalence of about 7–10 per 1000 live births, CHD is a leading cause of infant death [1, 3]. For children who survive with CHD, the odds of developing mental disabilities is 9 times higher than in the general population, while the odds of experiencing health problems that limit physical activity is 14 times higher [4].

CHD may be diagnosed during pregnancy, which has been found to decrease mortality risk and improve patient outcome following post natal surgery [5]. The most common fetal monitoring technique today is Doppler ultrasound, either through a single hand-held Doppler ultrasound transducer or cardiotocography (CTG). CTG and hand-held Doppler ultrasound provide estimates of the fetal heart rate (FHR), but no information on the fetal heart rate variability (FHRV) or more sophisticated analysis of the electrical cardiac activity [6]. In addition, CTG has technical limitations such as low sampling frequency, which limits its FHR estimation capability [7]. Since the frst commercial fetal heart monitor was launched in 1968 [8], continuous CTG monitoring during delivery has been linked to an increase in caesarean sections and instrumental births, which is a risk to mother and fetus [9], but did not have an impact on perinatal mortality and morbidity [10].

The waveform of the electrocardiogram (ECG) may contain important information about the cardiac electric activity beyond the FHR [11], which is routinely interpreted to assess cardiac function in adults [12]. Several studies have pointed out the need of ECG waveform analysis when monitoring pregnancies [7, 13]. In particular, ECG waveform analysis may be used in addition or as an alternative to CTG [14], to improve the detection of pathologic cardiac events or to provide information on uterine contraction [13].

Despite the above-mentioned advantages, abdominal ECGs show a low signal to noise ratio due to the small myocardial volume of the fetus and interference from e.g. the maternal heart, motion artifacts and muscle contractions. Therefore, postprocessing is required to extract the fetal signal contribution. This is challenging, because signal and artifacts are overlapping in the frequency domain [15, 16]. State of the art methods use adaptive fltering of the abdominal signal with a thoracic reference signal to overcome this problem [17, 18]. A lack of limited gold standard data, comprising simultaneous fECG recordings from abdominal electrodes, and an invasive scalp electrode further complicate the extraction [19]. Previous work has addressed this challenge by creating a fECG simulator in which each heart is represented by a moving dipole [20]. Herein, we simulated fetal and maternal cardiac activity using an image-based fnite element model to derive a realistic abdominal fECG.

#### **2.2 Methods**

#### **2.2.1 Geometrical mesh construction**

For our maternal and fetal heart we used a biventricular mesh of a female adult heart based on CT images [21]. We used myocardial fber orientations from the same study, which have been generated using a rule-based approach to reproduce experimental fndings [22, 23]. Both hearts were then augmented to ft into a female pregnant torso from the FEMONUM database [24, 25, 26]. For the maternal heart, the biventricular triangular surface mesh was translated, scaled and rotated by matching the surrounding torso to the pregnant torso which it was incorporated into. For the fetal heart, the same was done for a fetal torso originally included in the model from the FEMONUM database. Incorporation of biventricular surface meshes into the pregnant torso was done in Paraview [27], while volume mesh generation was done using gmsh [28]. The fnal model is a 3D fnite element mesh with ∼ 17 million nodes as shown in Fig. 2.1a. Elements are tetrahedral with sizes approximately 0.4 and 50 mm for the hearts and torso respectively.

(a) Visualisation of the fnite element mesh used in our simulations. Hearts and torso are based on models by [21] and the FEMONUM database [24, 25, 26] respectively.

(b) The extracellular potential was measured at 11 nodes on the torso to obtain fECG traces. The fgure shows thoracic (purple) and abdominal (green) electrodes. Numbers denote the diferent channels used for signal processing.

#### **2.2.2 Electrophysiological modelling**

Ionic current properties of the maternal heart were modelled using the Ten Tusscher model of human ventricular myocytes [29]. For the fetal heart a modifed version of the model, adapted to match fetal ventricular electrophysiology, was used [30]. Tissue propagation was modelled with the pseudo-bidomain approach [31] as implemented in CARPentry [32]. Following [22], intracellular and extracellular conductivities in the myocardium were set to = (0.27,0.081,0.045) and = (0.9828,0.3654,0.3654) S m−<sup>1</sup> in the longitudinal, transverse and normal direction respectively, while the torso was assigned an isotropic conductivity of 0.216 S m−<sup>1</sup> .

Sinus rhythm was simulated over 5 seconds by pacing both hearts with a 10 ms stimulus of 100 uA/cm<sup>2</sup> . The stimuli were delivered to all vertices on the endocardial surfaces through a transmembrane current, resulting in the activation pattern displayed in Fig. 2.2. Basic cycle lengths were set to 500 and 450 ms for the maternal and fetal heart, respectively. The fetal cycle length was chosen to refect the heart rate after 38 weeks of gestation [30].

Fig. 2.2: Activation map of the maternal heart during simulated sinus rhythm. Activation times are measured from stimuli onset and encoded by colour.

#### **2.2.3 Extracellular potential measurements**

In order to compute extracellular potential traces using CARPentry, the average potential was set to zero. For the fECG analysis, we located 11 virtual electrodes on the torso, considering both a 9-electrode and 8-electrode fECG setup [33]. The placement of the virtual electrodes is shown in Fig. 2.1b. Note the distinguished abdominal and thoracic electrodes which are meant to record the fECG and maternal ECG, respectively. The extracellular potential for the selected nodes was extracted. Subsequently, the voltage between all electrodes and a reference electrode on the lower abdomen was calculated from their potential diference to obtain ECG signals.

#### **2.2.4 Fetal ECG extraction using signal processing methods**

In this part, we study the signal processing methods for extracting the fECG from maternal abdominal recordings. An fECG is mixed with several disturbance sources in the maternal abdominal ECG. The primary source of disturbance is the maternal ECG which coincides with fECG in time, and frequency domains [16]. Hence, the extraction of fECG is challenging, and novel signal processing methods are demanded for the efcient diagnosis of fetal cardiac disorders and further treatment.

Diferent signal processing methods for fECG extraction are given in the literature. Rasti-Meymandi *et al.* proposed a deep learning algorithm named AECG-DecompNet which efciently extracts both fECG and maternal ECG [34]. AECG-DecompNet has two main network series; one network estimates the maternal ECG, and the other network eliminates the noise and interference. The two series of networks have shown superior results in QRS (a complex in the ECG signal) detection compared to utilizing one network only. Nevertheless, proper training of the frst network is important for robust fECG extraction. Independent component analysis (ICA) [35] and multi-ICA [36] also have been considered for fECG extraction, for which statistical models in the algorithms are challenging. Furthermore, adaptive fltering has also been used for extraction of fECG due to fast and simple implementation [37, 38, 17, 18]. However, fECGs extracted by adaptive fltering still have remaining maternal signal components.

We use an algorithm based on the combination of wavelet analysis and adaptive fltering proposed in [39]. We frst process abdominal and thoracic signals by wavelet transformation. Then, we process the detail coefcients of the wavelet transformed signals as inputs for adaptive fltering using the least mean square (LMS) fltering. We subsequently denoise the output of the adaptive flter, and fnally, take the inverse wavelet transform of the denoised coefcients.

The abdominal ECG is a non-stationary signal. We use stationary wavelet transform (SWT) to avoid the Gibbs phenomenon by removing down-sampling and up-sampling coefcients. The output of the SWT has the same length as the input signal. SWT has two sets of functions: the scaling and wavelet functions, denoted by () and Φ(), where stands for the th sample point of the signal. The functions are defned based on a chosen wavelet function. We decompose the signal () using the wavelet decomposition to approximation coefcient , and detail coefcient ,. The coefcients at the th scale are [39]

$$c\_{f,k} = < f(n), \Phi\_{f,k}(n)>,\tag{2.1}$$

$$d\_{j,k} = ,\tag{2.2}$$

with

$$ = \int\_{-\infty}^{\infty} 2^{-\frac{j}{2}} f(n) \Phi \star (2^j n - k) d n,\tag{2.3}$$

Fig. 2.3: Adaptive fltering for fECG extraction.

where ★ denotes the complex conjugation.

At the next step, we consider the detail coefcients of the SWT for adaptive fltering. Fig. 2.3 shows the adaptive fltering used in this section. We utilize the coefcients of the abdominal ECG as the original input () and the coefcients of thoracic ECG as the reference input (). The adaptive flter is based on LMS, which its update is given by [39]

$$
\omega(n+1) = \omega(n) + 2\mu e(n)\omega(n),\tag{2.4}
$$

where () and are the adaptive flter coefcients and step size. Also, () is the feedback in the adaptive flter where the coefcients of () are constantly adjusted until the output () is very close to the maternal ECG components of the abdominal ECG. The output signal () is fltered as

$$\mathbf{y}(n) = \sum\_{i=0}^{m-1} w\_i \mathbf{x}(n-i) = \mathbf{w}^T(n)\mathbf{x}(n),\tag{2.5}$$

where () denotes the transposed vector of (). Finally, we consider the output () for denoising at the next step. After denoising using wavelet transform, we take an inverse SWT to obtain the fECG.

#### **2.3 Results**

The conducted simulations resulted in distributions of the body surface potential as depicted in Fig. 2.4. When stimulated, both hearts produced regions on the torso of increased potential in the basal direction and decreased potential in the apical direction. The potential changes caused by the maternal ventricular excitation are about 10 times larger than for the fetal ventricles.

Fig. 2.4: Extracellular potential on the torso surface 5 ms after stimulation of the fetal (left) and maternal heart (right).

Signal processing of one clinical recording led to signals clearly representing fetal cardiac activity for some lead combinations (see Fig. 2.5). However, for other lead combinations, the maternal signal contribution is still pronounced and does not allow a clear detection of the fetal heart beats.

Fig. 2.6 illustrates one lead of the simulated ECG recording. Negative peaks, positive peaks and blunt positive hills indicate maternal depolarisation, fetal depolarisation and maternal repolarisation, respectively. The fltered signal still contains signifcant contributions from all of the three observed waveforms. Consequently, fetal heartbeats could not be reliably detected.

#### **2.4 Discussion**

The proposed model framework was able to produce abdominal ECGs containing maternal and fetal signal contributions. Furthermore, it was possible to simulate fetal and maternal cardiac activity in an isolated manner. This could in general be used as ground truth data for design and validation of signal processing methods. However, at the moment the simulated signals are too far from clinically observed recordings, rendering them unsuitable to tailor clinical solutions. In this study we modelled the

Fig. 2.5: Abdominal ECGs recorded in two leads and the resulting fltered fetal signal contribution. While lead 1 (a) enabled a clean fetal signal, lead 5 (b) still contained a lot of noise and maternal signal contribution after fltering.

Fig. 2.6: Simulated ECG, recorded on the upper right abdomen (top) and signal extracted by the adaptive flter.

torso as a uniform medium and do not account for the varying tissue conductivities which infuence the current propagation. Including additional tissue types such as bone or lungs in the model could lead to more realistic conduction and therefore more realistic fECG recordings.

Furthermore, ventricular activation was assumed to happen instantly on the whole endocardium, spreading to the epicardial surface. Several sophisticated models of ventricular activation sequences have been presented in literature [40, 41] and could be incorporated in the proposed model to enable closer resemblance of clinically observed activation patterns [42]. Future studies may also explore how time-dependent factors such as muscle contractions and respiratory movement can be incorporated.

Once a reasonable morphology of simulated signals is achieved, diferent artifcial noise sources could be incorporated to make them more realistic and test the robustness of the proposed fltering algorithm. Besides designing signal processing algorithms, realistic simulation of fECGs could also be used to assess the infuence of fetal orientation and electrode positioning. The latter could serve to establish standard electrode positions which do not yet exist.

#### **2.5 Conclusions**

We presented an electrophysiological model to simulate maternal and fetal cardiac activity during sinus rhythm. The resulting change in extracellular potential on the torso surface was used to derive the resulting voltages for several abdominal and thoracic ECG leads. Furthermore, a wavelet based adaptive fltering approach was used to extract the fetal contribution from an abdominal ECG recording for clinical as well as for simulated signals. The simulated signals contain maternal and fetal peaks but have lower complexity than clinical recordings, indicating further need to improve the proposed model.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Chapter 3 Ordinary Diferential Equation-based Modeling of Cells in Human Cartilage**

Kei Yamamoto1,2, Sophie Fischer-Holzhausen<sup>3</sup> , Maria P Fjeldstad<sup>2</sup> , Mary M Maleckar<sup>1</sup>


**Abstract** Chondrocytes produce the extracellular cartilage matrix required for smooth joint mobility. As cartilage is not vascularised, and chondrocytes are not innervated by the nervous system, chondrocytes are therefore generally considered non-excitable. However, chondrocytes do express a range of ion channels, ion pumps, and receptors involved in cell homeostasis and cartilage maintenance. Dysfunction in these ion channels and pumps has been linked to degenerative disorders such as arthritis. Because the electrophysiological properties of chondrocytes are difcult to measure experimentally, mathematical modelling can instead be used to investigate the regulation of ionic currents. Such models can provide insight into the fnely tuned parameters underlying fuctuations in membrane potential and cell behaviour in healthy and pathological conditions. Here, we introduce an open-source, intuitive, and extendable mathematical model of chondrocyte electrophysiology, implementing key proteins involved in regulating the membrane potential. Because of the inherent biological variability of cells and their physiological ranges of ionic concentrations, we describe a population of models that provides a robust computational representation of the biological data. This permits parameter variability in a manner mimicking biological variation, and we present a selection of parameter sets that suitably represent experimental data. Our mathematical model can be used to efciently investigate the ionic currents underlying chondrocyte behaviour.

#### **3.1 Introduction**

An important feature of vertebrate joints is the presence of articular cartilage, a type of connective tissue covering the end of the bones and greatly reducing the friction of bone against bone and facilitating movement. Cartilage is largely composed of proteoglycans and collagens, forming a meshwork of extracellular matrix proteins that can withstand large mechanical load while also remaining fexible compared with bone tissue [1]. As cartilage has poor regenerative properties, joint disorders like osteoarthritis are becoming an increasing health burden as the global population ages. Osteoarthritis is a progressive degradation of the cartilage in joints, and is estimated to afect 12.1% of adults over the age of 60 in the United States [2]. The progression of cartilage degradation, and what triggers it, remains poorly understood, and there are no treatments other than pain relief or joint replacement surgery [3]. Cartilage is largely made up of extracellular matrix. The cells responsible for the synthesis and turnover of this matrix are known as chondrocytes, which are embedded in small clusters in cartilaginous tissue. These cells are not innervated by the nervous system and are thus generally considered non-excitable. Additionally, as cartilage is not vascularised, they rely on difusion in order to exchange nutrients and waste material [4]. Their extracellular milieu is therefore often slightly hypoxic, and chondrocytes have low metabolic turnover rates and poor regenerative abilities [5]. This makes them prone to degenerative disorders including osteoarthritis, in which the cartilage generally is degraded at a faster rate than it is synthesised by chondrocytes [6].

Although chondrocytes are considered non-excitable, they do express a range of ion channels and ion pumps involved in the regulation of intracellular ionic concentrations, and, in turn, the membrane potential and downstream intracellular processes. However, the exploration of electrophysiological properties of chondrocytes has been limited due to their small size (of less than 10m), representing a technical challenge for electrophysiological experiments, implying large uncertainties in measurement results [7].

A mathematical model can provide an alternative, more accessible method of studying the membrane potential and ion dynamics for chondrocytes, and the construction of the model itself can be benefcial in understanding the system's behaviour. The mathematical model we present here can be seen as a reduced and abstracted representation of what is currently known about chondrocytes. Additionally, such models can be useful to develop hypotheses, to identify mechanistic details or knowledge gaps, and to guide experiments [8]. Also, mechanistic models such as this, built on biophysical knowledge, have an increasing interest in pharmacological research as they have proven themselves useful for qualitative understanding of the underlying physiology and pathophysiology [9]. Moreover, mathematical models can fnd application in drug development pipelines. Strauss *et al.* [10], for instance, demonstrate in a comprehensive review the value of *in silico* models as an integrated parts of risk assessment strategies to evaluate the proarrhythmic risk of certain drugs in humans. Finally, computational models are also essential tools to study inter- and intra-variability of cell physiology. Studying the sources of variability can provide better understanding of biological processes and aid in making meaningful predictions. One intuitive approach to investigate sources of variability is to create and calibrate a population of models, which is employed in the present study [11].

Maleckar *et al.* [12] introduce a frst generation mathematical model for the membrane potential of the chondrocyte. Their model is based on an intensive literature review and is additionally partly supported by experimental data. Both time-dependent and time-independent ion channel currents, as well as ion pump and ion exchanger currents, are represented in the model, with detailed description of the underlying model development and equation derivation. This frst generation of the chondrocyte model mainly focused on K + currents across the chondrocyte cell membrane. Furthermore, Maleckar *et al.* [13] presents a model extension by including functions for the Na<sup>+</sup> /K <sup>+</sup> pump, an active transport pump with an important role in setting the electrochemical gradient of Na<sup>+</sup> and K + across the cell membrane.

In this study, the model is expanded by including the ATP-sensitive K + (KATP) channel, a subgroup of inward-rectifying K + channels that have been identifed in human chondrocytes [14]. This channel has been found to play an important protective role against hypoxia in the cardiovascular and nervous system [15, 16, 17], which is notable given oxygen tension in cartilage is lower than in vascularised tissue. KATP channel function is sensitive to intracellular Mg2<sup>+</sup> concentrations, and channel activation is blocked by the binding of ATP [18]. Thus, when energy expenditure exceeds energy storage in the cell, intracellular ATP concentration decreases and the block of KATP channels is relieved [19], providing a direct link between cellular respiration and the membrane potential.

The previously published models were created using MATLAB [12, 13]. However, as a frst-generation computational model of chondrocyte electrophysiology, we believe it benefcial to present our model on an open-source platform, so that more users can access it without having to purchase access to MATLAB. Python is an open-source programming language that has been embraced by researchers across various branches of science, owing to its powerful libraries like Numpy [20] and Scipy [21] in addition to interactive environments (IPython, Jupyter) that enable developers to provide interactive manuals of the software. In this work we present a re-implementation of the mathematical model for the chondrocyte's membrane potential fst published in 2018 by Maleckar *et al.* [12]. Our Python implementation is freely available at https://github.com/mmaleck/chondrocyte. By replicating fgures from Maleckar *et al.* [13] we ensure a correct model transfer. The addition of KATP channel activity to the membrane dynamics of the model indicate that KATP channels may have a stabilising efect on overall K + currents in the presence of high Mg2<sup>+</sup> concentrations.

To learn more about the model's behaviour and parameter sensitivity, a population of models was created, allowing the investigation of ranges of parameter values and their efect on the steady state of the system. This investigation identifed an important role of the Na<sup>+</sup> /K <sup>+</sup> pump current in the modulation of the ionic balance of the model, as well as setting the resting membrane potential. Furthermore, the model can be expanded upon by the addition of ion channel and ion pump dynamics not described here. Thus, this readily available re-implementation of a frst generation chondrocyte model can provide insight into the dynamics of the chondrocyte membrane potential, and how this may be afected by variations in ion channel and ion pump function.

#### **3.2 Methods**

To investigate the chondrocyte membrane potential, we adopted a chondrocyte model and simulation protocols from previously published works [12, 13].

Fig. 3.1: Schematic of the model and its implemented ion channels and ion pumps. −, −, and −: background ionic activity for Na<sup>+</sup> , Cl<sup>−</sup> , and K + , respectively, : Na<sup>+</sup> /Ca2<sup>+</sup> exchanger, : Na<sup>+</sup> /H + exchanger, : Na<sup>+</sup> /K <sup>+</sup> pump, , : ATP-dependent Ca2<sup>+</sup> pump, −: delayed-rectifer K + channel, −2: two-pore K + channel, −: Ca2<sup>+</sup> activated K + channel, − : ATP-sensitive K + channel, 4: Transient Receptor Potential V-4 (a mechanosensitive cation channel). Adapted from [12] with permissions dictated by the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

The given model is defned by a set of ordinary diferential equations. Those systems are defned as

$$\frac{d\mathbf{x}}{dt} = f(\mathbf{x}(t), \boldsymbol{\mu}(t), k) \tag{3.1}$$

and describe the temporal evolution of states *x(t)*. In the used model, the states are the membrane potential (*V*) and the ion concentrations of Na<sup>+</sup> , K + , and Ca2<sup>+</sup> . Furthermore, the function *f* depends on currents *u(f*) for included ion pumps, channels and exchangers and model parameters *k*. Systems described in the form as given in Equation (3.1) need to be solved numerically. The equation which governs the chondrocyte transmembrane potential follows the Hodgkin–Huxley formalism where each cell component is represented as an electrical element, where cell membranes have a capacitance *C* and ion channels are treated as resistors. This gives a mathematical formalism for the membrane potential *V* which includes all transport processes that are electrogenic and reads as follows:

$$C\_m \frac{dV\_m}{dt} = -\sum\_{ions} I\_{ion} \tag{3.2}$$

This work presents a model extension by including the KATP currents to the previously published chondrocyte model. Furthermore, the population of models idea is used to investigate the model behaviour under diferent parameter settings. The role of the maximum Na<sup>+</sup> /K <sup>+</sup> pump current is studied in depth.

#### **3.2.1 Mathematical modelling of ATP-sensitive K**<sup>+</sup> **currents**

We have added an ATP-sensitive K+ channel to the model from Maleckar *et al.* [13] to investigate its role for the functionality and membrane potential of chondrocytes. The mathematical formulation of the KATP current − reads:

$$\mathbf{I}\_{\rm K-ATP} = \sigma \mathbf{g}\_0 \mathbf{p}\_0 \mathbf{f}\_{\rm ATP} \left(\mathbf{V}\_{\rm m} - \mathbf{E}\_{\rm K}\right) \tag{3.3}$$

where is the channel density, g is the unitary channel conductance, p<sup>0</sup> is the maximum open channel probability, and fATP is the fraction activated channels. E<sup>K</sup> denotes the equilibrium potential for K + in the given circumstances. In the previous implementation of − , a constant value was used for g0. However, the unitary conductance can be expressed as [22]

$$\mathbf{g}\_0 = \gamma\_0 \mathbf{f}\_\mathbf{M} \mathbf{f}\_\mathbf{N} \mathbf{f}\_\mathbf{T} \tag{3.4}$$

where <sup>0</sup> is the unitary conductance in the absence of intracellular Na<sup>+</sup> and Mg2<sup>+</sup> and depends on [K + ] :

$$\gamma\_0 = 35.375 \left( \frac{[\text{K}^+]\_0}{5.4} \right)^{0.24} \tag{3.5}$$

The term f<sup>M</sup> in Equation (3.4) represents the inward rectifcation generated by intracellular Mg2<sup>+</sup> ions and can be expressed by means of a Hill equation :

$$\mathbf{f\_M} = \frac{1}{1 + \frac{\left[\text{Mg}^{2+}\right]\_i}{\text{K}\_{\text{h}, \text{Mg}}}} \tag{3.6}$$

where Kh,Mg is the half-maximum saturation constant that depends on membrane potential and on [K + ]<sup>o</sup> :

$$\mathbf{K}\_{\rm h,Mg} = \mathbf{K}\_{\rm hMg}^{0} \left( \left[ \mathbf{K}^{+} \right]\_{0} \right) \exp\left( -\frac{2\delta\_{\rm Mg} \mathbf{F}}{\mathbf{RT}} \mathbf{V}\_{\rm m} \right) \tag{3.7}$$

$$\mathbf{K}\_{\rm h,Mg}^{0}\left(\left[K^{+}\right]\_{\rm o}\right) = \frac{0.65}{\sqrt{[K^{+}]\_{0} + \rm S}}\tag{3.8}$$

Here the value of the electrical distance Mg was set to 0.32, F is the Faraday constant, R is the gas constant, and T is the absolute temperature. In a similar manner, the term f<sup>N</sup> represents the inward rectifcation caused by intracellular Na<sup>+</sup> ions and is expressed as

$$\mathbf{f}\_N = \frac{1}{1 + \left(\frac{[\text{Na}^\*]\_i}{K\_{\text{h,Na}}}\right)^2} \tag{3.9}$$

$$K\_{\rm h,Na} = K\_{\rm h,Na}^{0} \exp\left(-\frac{\delta\_{\rm Na} F}{RT} \mathbf{V}\_{\rm m}\right) \tag{3.10}$$

where Na = 0.35 and K 0 <sup>h</sup>,Na <sup>=</sup> <sup>25</sup>.<sup>9</sup> was used.

#### **3.2.2 Population of Models**

In the context of our work, we use the term population of models for set model simulations with randomly varied parameter sets for ionic current conductances [11]. Generally, a population of modes is useful to investigate parameter sensitivity [23, 24], sources of variability [25], as well as emergent model behavior dependency for a wide range of parameters.

In this study, we searched for parameters that had a signifcant efect on the simulation results and the physiological relevance of the model. Table 3.1 gives an overview of the selected set of parameters, the values in the original implementation, the parameter regime used to create a population of models, and a brief parameter description.

Parameter ranges originate from parameter sampling from log normal distribution. The parameter distribution was not ftted to experimental data and therefore underlies the assumptions that the parameters' distributions are skewed and the parameter values are positive. The distribution underlies the reasoning that numerous works have shown that a log normal distribution is often a useful assumption to describe the random variation in biological samples [26]. The population of models presented in the results is created out of 100 simulation runs. Each member of the population has its characteristic set of parameters, and all parameter set were sample from log normal distributions with the mean being the parameter value and a variance = 0.15.


Table 3.1: Ensemble of parameters selected for creating a population of models.

#### **3.3 Results**

The following results section is divided in three parts. Firstly, we provide the validation for our new implementation that are used for tow following sections. In the second part, we present the simulations focused on the KATP currents and the efect of diferent Mg2<sup>+</sup> concentrations. Lastly, the parameter and model behaviour is studied based on the population of models approach.

#### **3.3.1 Validation**

We performed the validation against a previous publication, Maleckar *et al.* [13]. Figure 3.2 shows the results of temperature-dependent contribution of the Na<sup>+</sup> /K + pump electrogenic current to the chondrocyte resting membrane potential implemented both in MATLAB and Python as a replication of Figure 2 from Maleckar *et al.* [13]. Although our model includes other currents than shown in the fgure, we conclude that our new implementation is well validated.

#### **3.3.2 Results for the ATP-sensitive K**<sup>+</sup> **current**

In order to investigate the efect of − on K <sup>+</sup> dynamics in the chondrocyte, numerical simulations were performed as shown in Figure 3.3. Accurate measurements of intracellular Mg2<sup>+</sup> concentrations in chondrocytes are, to our knowledge, unavailable, hence we tested a range of initial Mg2<sup>+</sup> values. The Mg2<sup>+</sup> concentrations used for the simulation were (a) 0.1 mM (b) 1.0 mM (c) 10 mM, while three diferent K + concentrations, 5 mM, 30 mM, 70 mM, were used for each Mg2<sup>+</sup> concentration as indicated in the fgure. To further clarify the efect of diferent intracellular Mg2<sup>+</sup> concentration to the overall chondrocyte matrix, we introduce Figure 3.4 where Figure 3.4 (a) shows the steady state voltage dependence of − while

Fig. 3.2: Temperature-dependent contribution of the Na<sup>+</sup> /K <sup>+</sup> pump electrogenic current to the chondrocyte resting membrane potential at **A** 23 ◦C and **B** 37 ◦C. Unbroken lines indicate our new implementation in Python, while dotted lines indicate the published MATLAB implementation.

(b) illustrates the time-dependent changes of chondrocyte membrane potential under diferent intracellular Mg2<sup>+</sup> concentrations.

Fig. 3.3: Sensitivity of − against varying extracellular K + concentration based on diferent intracellular Mg2<sup>+</sup> concentrations. Intracellular Mg2<sup>+</sup> concentrations used for fgures are **A** 0.1 mM **B** 1.0 mM **C** 10 mM respectively, while three concentrations of extracellular K + concentrations (5 mM, 30 mM, 70 mM) are used for all fgures.

#### **3.3.3 Populations of Models**

A population of 100 models with randomly varied sets of parameters (Table 3.1) is illustrated in Figure 3.5. The steady state solutions for all species are reached within

Fig. 3.4: Intracellular Mg2<sup>+</sup> concentration dependent contribution of the ATPsensitive K + current to the chondrocyte resting membrane potential. In **A**, we illustrate steady-state voltage dependence and the Mg2<sup>+</sup> dependence of the − current, while **B** shows the hyperpolarization of chondrocyte membrane potential with diferent intracellular Mg2<sup>+</sup> concentration. The extracellular K + concentration used for **A** and **B** is fxed at 7 mM.

the simulation time. Interestingly, the steady states for both Na<sup>+</sup> and K + reach similar concentration levels regardless the set of parameters, whereas Ca2<sup>+</sup> concentrations at steady state solution are strongly afected by the set of parameters. It is notable that the behaviour of the membrane potential follows the shape of the Ca2<sup>+</sup> time curve, and is therefore also infuenced by the composition of the parameter sets.

For all 100 parameter sets, the Na<sup>+</sup> steady state solution is close to zero. However, the time point at which steady state is reached is parameter-dependent (Figure 3.5 **B** and **D**, Figure 3.6 **A** and **B**). To investigate the individual parameter efects, simulations with randomly drawn parameters for each individual parameter represented in Table 3.1 were performed. Larger parameter ranges were also investigated; Figure 3.5 , ranges from 1.375 to 2.1125, whereas for Figure 3.6, , values between 0.1 and 7.0 were tested. The simulation results indicate that the scaling factors of the maximum Na<sup>+</sup> /K <sup>+</sup> pump current ( ,) have a major efect on the membrane potential, the fnal concentrations for Na<sup>+</sup> and K + , as well as the duration required for the system to reach steady-state.

A wide range of , values resulted in a Na<sup>+</sup> depletion, which is unlikely to occur in a physiologically viable chondrocyte preparation. A step-wise scan through , (0.1 < , < 7.0) reveals three parameter regimes for the scaling of this current. For 0.1 < , < 0.6, here called the low parameter regime, Na<sup>+</sup> concentrations range from 10mM to 90mM (Figure 3.6(a)). If 0.6 < , < 5.0, a non-physiological state of Na<sup>+</sup> depletion is eventually reached, as intracellular sodium moves towards zero over time. The value of , also afects the timing of this depleted state (Figure 3.6 **B**): the smaller the scaling factor and thus the smaller the current, the later the systems reaches its steady state. Henceforth, this parameter range is referred to as the middle parameter regime. For , > 5.0, steady-

Fig. 3.5: A population of 100 models. Each simulation trajectory results from a randomly drawn parameter set. A simulation time of 2500 seconds ensures that all species reach their steady state. **A** shows the profle of membrane potential, **B** shows the profle of K + concentration, **C** shows the profle of Ca2<sup>+</sup> , and **D** shows the profle of Na<sup>+</sup> concentration.

state Na<sup>+</sup> concentrations become negative, which is biologically and numerically infeasible and causes the simulation to fail.

A population of models approach was again used to further probe the model dependence on parameter variation in for low and middle parameter regimes, this time varying initial conditions (Fig. 3.7). Figure 3.7 shows strong model dependence on initial conditions for evolving variables (panels **A**–**D**). However, while there are a variety of stable steady states at diferent parameter sets, where for example, intracellular sodium is not depleted (panel **B**), all these occur for the low parameter regime of , (blue traces, Fig. 3.7, all panels). Thus, the low parameter regime of INaK conductance (0.1 < , < 0.6), corresponds to, roughly, an INaK current of approximately 0.9–5.35 pA/pF and likely represents a physiologically relevant model regime incorporating a variety of stable steady states for evolving ionic concentrations and the membrane potential. To further investigate the overall dependence of the chondrocyte's resting membrane potential on K+ currents, we varied the conductance parameters for these at a variety of , values in the low parameter regime (Fig. 3.8). While there are a variety of stable steady states for the resting membrane potential of the chondrocyte in this example, even as , increases (moving right from panel **A** to **C**), all the values for each unique parameter set in the current model set the chondrocytes resting membrane potential somewhere between approximately -50 and -70 mV.

Fig. 3.6: The scaling of the maximum Na<sup>+</sup> /K <sup>+</sup> pump current afects the timing and concentration of the Na<sup>+</sup> steady state as well as the membrane potential. The simulation trajectory for membrane potentials and the Na<sup>+</sup> concentrations in the low and the middle parameter regime are displayed here. Simulations for , > 5.0 can not be presented because negative Na<sup>+</sup> concentrations causes the simulation to fail. **A** and **B** display simulation curves for 0.1 < , < 0.6. The steady state concentration of Na<sup>+</sup> ranges from 10mM up to 90mM. **C** and **D** shows simulation trajectories for 0.6 < , < 5.0. In all cases Na<sup>+</sup> depletion can be observed. The timing is parameter dependant.

Fig. 3.7: Initial conditions for evolving concentrations and membrane potential were varied, with 100 parameters / Simultaneously, 10 diferent values for , through the low to middle parameter regimes (0.1 - 4.9, 0.5 step) were used, for a total of 10 x 100 simulations. The colors reveal this , gradation, with dark blue as the lowest value (0.1), and deep red as the highest value (4.9). The same colors always represent the same , value but diferent overall parameter sets.

Fig. 3.8: Simulations (100 parameter sets, varying I K,b:g\_K\_b\_bar, I K-DR: g\_K\_DR, I K2p: I\_K\_2pore\_scale, I K-Ca (BK): gBK, I K, ATP: gamma\_0\_scale over a distribution as described in Methods) were run with difering values for , in the low parameter regime range (0.1-0.6) until intracellular sodium reached steady state. Each color in panels **A**-**C** above refers to a unique parameter set, and panels **A**, **B**, and **C** show the chondrocyte resting membrane potential reaching steady state over time for increasing values of , (0.1, 0.3, and 0.5, respectively).

#### **3.4 Discussion and Conclusion**

As can be seen in Figure 3.3, higher Mg2<sup>+</sup> concentrations help stabilise − currents with respect to varying extracellular K + concentrations. Also, it can be seen from Figure 3.4 that the lower Mg2<sup>+</sup> concentration contributes to the hyperpolarisation of the membrane potential. This stabilising efect is interesting, as experimental data suggest a protective role of high Mg2<sup>+</sup> on joint health. Although the exact mechanisms behind this protection remain elusive, some hypotheses point to a role of Mg2<sup>+</sup> in reducing chondrocyte apoptosis and facilitating chondrocyte proliferation [27, 28, 29, 30]. One caveat of this paper, however, is that physiological Mg2<sup>+</sup> concentrations in chondrocytes remain unknown, and thus we have tested our model with a range of concentrations. Reliable measurements of intracellular Mg2<sup>+</sup> concentrations will be necessary in order to assess the physiological efects of Mg2<sup>+</sup> on K + currents, and the potentially protective role of Mg2<sup>+</sup> and the ATP-dependent K + channel in cartilage maintenance.

To get a complete understanding of the three parameter regimes for the scaling factor of the maximum Na<sup>+</sup> /K <sup>+</sup> pump, current further investigation is needed. It might be that the observed behaviour is a result of the chosen initial value for the Na<sup>+</sup> concentration. One could investigate this hypothesis by treating initial Na<sup>+</sup> concentrations as a parameter and scanning through a wide range of initial concentrations comparable to what has been demonstrated in Figure 3.6. But it could also be possible that the scaling parameter has this narrow parameter regime for the given model complexity and parameterisation. Additionally, it would be interesting to investigate the link between Ca2<sup>+</sup> concentration and the membrane potential further.

The lack of available biological data presents a challenge for the validation of models such as the one presented here. As our study shows, measurement data are needed to determine parameter values rather than parameter ranges leading to plausible results. Due to the high parameter uncertainty of the model, care must be taken in interpreting results in order to make useful predictions and to suggest hypotheses for further testing. A more thorough investigation into the range of appropriate ionic concentrations remains to be performed.

Here we present a model of chondrocyte electrophysiology, with particular attention to the role of K + currents in setting the steady-state membrane potential. Our preliminary results support a protective role of Mg2<sup>+</sup> in cartilage maintenance as recorded in clinical studies by potentially stabilising K + currents through the ATPdependent K + channel, although more research is required. Furthermore, we show how a population of models can be used to examine fuctuations in a range of parameters, and their interactions as the model reaches steady-state membrane potential. Such random variations in parameter values will also act to make the model more realistic and more robust to naturally occurring fuctuations as seen in biological data. Finally, the chondrocyte model has been re-implemented from MATLAB into Python to increase its accessibility, and is available as an open-source repository on Github, with demo scripts to aid interested parties in getting started.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Chapter 4 Conduction Velocity in Cardiac Tissue as Function of Ion Channel Conductance and Distribution**

Kristian Gregorius Hustad<sup>1</sup> , Ena Ivanovic<sup>2</sup> , Adrian Llop Recha<sup>3</sup> , Abinaya Abbi Sakthivel<sup>3</sup>


**Abstract** Ion channels on the membrane of cardiomyocytes regulate the propagation of action potentials from cell to cell and hence are essential for the proper function of the heart. Through computer simulations with the classical monodomain model for cardiac tissue and the more recent extracellular-membrane-intracellular (EMI) model where individual cells are explicitly represented, we investigated how conduction velocity (CV) in cardiac tissue depends on the strength of various ion currents as well as on the spatial distribution of the ion channels. Our simulations show a sharp decrease in CV when reducing the strength of the sodium (Na<sup>+</sup> ) currents, whereas independent reductions in the potassium (K1 and Kr) and L-type calcium currents have negligible efect on the CV. Furthermore, we fnd that an increase in number density of Na<sup>+</sup> channels towards the cell ends increases the CV, whereas a higher number density of K1 channels slightly reduces the CV. These fndings contribute to the understanding of ion channels (e.g. Na<sup>+</sup> and K + channels) in the propagation velocity of action potentials in the heart.

#### **4.1 Introduction**

A healthy heart rhythm is essential for the proper functioning of the cardiac pump, and requires the coordinated propagation of electrical impulses through the myocardium. The cardiac action potential is a change in the membrane potential governed by the ionic current fowing through ion channels, which are distributed along the cell membrane. Current fowing into the cell through activated sodium (Na<sup>+</sup> ) channels is responsible for the rapid upstroke of the action potential [1]. This is followed by a current infow through L-type calcium (CaL) channels and an outfow through diferent types of potassium channels (e.g. Kr and K1) leading to repolarisation, thus bringing the cell membrane to resting membrane potential [1]. More precisely, Kr is the rapid component of the delayed rectifer, and K1, the inward rectifying potassium current, which stabilises the resting membrane potential [1]. L-type calcium channels, on the other hand, are also responsible for the excitation-contraction coupling of the cardiac muscle [1].

Disorders of electrical conduction, such as slow conduction and conduction block, can lead to life-threatening arrhythmias, which occur frequently in the diseased heart. In this study, we investigated how both the strength of several transmembrane ionic currents and the spatial distribution of ion channels along the cell membrane infuence conduction velocity (CV) in cardiac tissue. Our aim is to compare the efect on CV with the use of two computational models: the monodomain model and the extracellular-membrane-intracellular (EMI) model.

The monodomain model is a classical approximation of the electrical propagation in myocardial tissue based on a homogenised mathematical model of the cell. The intra- and extracellular domains overlap and are considered continuous. As a consequence, the monodomain model ofers a good insight of large scale efects, but it is limited when sizes are reduced to a single cell. Alternatively, the EMI model represents the extracellular, the membrane, and the intracellular spaces at the expense of computational power. Therefore, one of the main advantages is the possibility to introduce changes in local cell properties (e.g. changes in ion channel density distribution on the cell membrane) that might contribute signifcantly to action potential propagation [2], [3].

#### **4.2 Models and methods**

#### **4.2.1 The monodomain model**

The monodomain model is a simplifcation of the bidomain model [4]. In the bidomain model, heart tissue is classifed into two groups or domains: extracellular and intracellular, defned by their respective electric potentials, and , and conductivities and .

Each point in the heart is considered to be in both domains. Therefore, both spaces overlap.

The physical description can be addressed using a generalisation of Ohm's Law. The current densities at each domain will be: = −∇ , and = −∇. Assuming that there are no other sources than the membrane, the conservation of charge applies, and thus: ∇( + ) = 0.

The current fowing from one domain to the other through the cell membrane is called transmembrane current, . Because the charge is conserved, then ∇ = −∇ = . Transmembrane current (Eq. 4.1) depends on the voltage drop between both domains, = −, the membrane capacitance , the ionic current, ion, and

#### 4.2 Models and methods 43

surface area-to-volume ratio of cardiac cell, .

$$I\_m = \beta\_m \left( C\_m \frac{\partial \nu}{\partial t} + I\_{\rm ion}(\nu) \right) \tag{4.1}$$

It can be shown that the set of equations governing the bidomain model are the ones expressed in Eq.4.2 and Eq.4.3.

$$\nabla \mathbf{G}\_I (\nabla \boldsymbol{\nu} + \nabla \boldsymbol{\mu}\_e) = \beta\_m \left( C\_m \frac{\partial \boldsymbol{\nu}}{\partial t} + I\_{\rm ion} (\boldsymbol{\nu}) \right) \tag{4.2}$$

$$
\nabla G\_i \nabla \nu + \nabla (G\_i + G\_e) \nabla \mu\_e = 0 \tag{4.3}
$$

Solving the bidomain equation is a computationally heavy process. Therefore, the monodomain model is often used instead. In the monodomain model, the anisotropy between extra- and intracellular spaces is assumed to be the same, i.e., their respective conductances are proportional = .

If we defne an efective conductivity, **ef** = <sup>1</sup>+ , then the bidomain equation can be simplifed and rearranged as shown in Eq.4.4.

$$\frac{\partial \nu}{\partial t} = \frac{1}{C\_m \beta\_m} \nabla G\_{\text{eff}} \nabla \nu - \frac{1}{C\_m} I\_{\text{ion}} \tag{4.4}$$

#### **4.2.2 The EMI model**

In the EMI model [5], the extracellular (E), cell membrane (M) and intracellular (I) domains, are represented explicitly as depicted in Figure 4.1. The intracellular spaces of both cells, Ω<sup>1</sup> and Ω<sup>2</sup> , are separated from the extracellular domain, Ω, by the cell membrane boundaries, Γ<sup>1</sup> and Γ2. Additionally, Γ1,<sup>2</sup> is the boundary separating the intracellular domains of two connected cells. Note that the EMI model is always solved in a three-dimensional space, as the extracellular space should be a single, connected domain.

The equation system describing the potentials in the EMI model is summarised in the following set of equations:

$$\begin{aligned} \nabla \cdot \sigma\_e \nabla \mu\_e &= 0 & \text{in } \Omega\_e, & \quad \mathfrak{n}\_e \cdot \sigma\_e \nabla \mu\_e &= -\mathfrak{n}\_i^k \cdot \sigma\_i \nabla \mu\_i^k \equiv I\_m^k & \text{at } \Gamma\_k, \\ \nabla \cdot \sigma\_i \nabla \mu\_i^k &= 0 & \text{in } \Omega\_i^k, & \quad \nu^k &= \frac{1}{\longrightarrow} (I\_m^k - I\_{i-}^k) & \text{at } \Gamma\_{\bar{\nu}}. \end{aligned}$$

$$\begin{aligned} \mu\_i^k &= 0 & \text{in } \Omega\_i^k, & \quad \nu^k &= \frac{1}{C\_m} (I\_m^k - I\_{\text{ion}}^k) & \text{at } \Gamma\_k, \\ \mu\_e &= 0 & \text{at } \partial \Omega\_e^D, & \quad \nu\_\omega &= \frac{1}{\omega\_\omega} (I\_m^k - I\_{\text{ion}}^k) & \quad \text{at } \Gamma\_k, \end{aligned}$$

$$\begin{aligned} \mathfrak{n}\_{\mathfrak{e}} - \mathfrak{o} &\quad \text{at } \partial \mathfrak{A}\_{\mathfrak{e}} \text{ } ,\\ \mathfrak{n}\_{\mathfrak{e}} \cdot \sigma\_{\mathfrak{e}} \nabla \mathfrak{u}\_{\mathfrak{e}} &= 0 \quad \quad \text{at } \partial \mathfrak{D}\_{\mathfrak{e}}^{N} \text{ } ,\\ \mathfrak{u}\_{i}^{k} - \mathfrak{u}\_{\mathfrak{e}} &= \nu^{k} \quad \quad \text{at } \Gamma\_{k} , \end{aligned} \qquad \begin{aligned} \mathfrak{n}\_{\mathfrak{e}} \cdot \sigma\_{\mathfrak{e}} \Delta\_{\mathfrak{e}} \circ \mu\_{i} &= -\mathfrak{n}\_{i}^{1} \cdot \sigma\_{i} \nabla u\_{i}^{1} \equiv I\_{1,2} \quad \quad \text{at } \Gamma\_{1,2},\\ \mathfrak{n}\_{i}^{1} - \mathfrak{u}\_{i}^{2} &= \mathcal{W} \quad \quad \text{at } \Gamma\_{1,2},\\ \mathfrak{n}\_{\mathfrak{e}} \cdot \sigma\_{\mathfrak{e}} \nabla \mathfrak{u}\_{\mathfrak{e}} &= \nu^{k} \quad \quad \text{at } \Gamma\_{1,2}, \end{aligned}$$

$$\begin{aligned} u\_i^k - u\_\varepsilon &= \nu^k & \text{at } \Gamma\_k, \\ s\_i^k &= F^k & \text{at } \Gamma\_k, \\ &\end{aligned} \qquad \begin{aligned} \cdots & \quad \cdots & \quad \cdots & \quad \cdots & \cdots \\ &\text{w = \frac{1}{\Gamma\_{1,2}} (I\_{1,2} - I\_{\mathfrak{g}}^k) & \quad \text{at } \Gamma\_{1,2}, \end{aligned}$$

Fig. 4.1: Two-dimensional schematic of the diferent domains for two connected cells in the EMI model. Adapted from [5] with permissions dictated by the Creative Commons CC BY license (https://creativecommons.org/licenses/by/4.0/).

As in the case of the monodomain model, potentials of intracellular, extracellular and transmembrane domains are denoted by , and , respectively. In addition to the ion current, ion, and the transmembrane current, , defned at the cell membrane, Γ, a gap junction current, g, and a transmembrane current, 1,2, are defned at the interface between two cells, Γ1,2, in order to model the gap junction dynamics.

The gap junction between neighbouring cells, Γ1,2, is modeled as a passive membrane with a constant resistance, . Its electric dynamics are described using the currents and 1,2, and the potential drop at the junction, .

Furthermore, and 1,<sup>2</sup> are the transmembrane and the gap junction capacitances. and are the conductivities of both intracellular and extracellular domain, whereas and are the outward pointing normal vectors of the inner and outer cell domains, respectively.

The boundary of the extracellular domain, Ω, is further divided into two parts: one corresponding to the Dirichlet boundary conditions, Ω , and other one corresponding to the Neumann boundary conditions, Ω . Additionally, the index can take the values 1 or 2 depending on which cell is described. Two cells were used to introduce the EMI model, however the model can easily be scaled up to consider more cells in the system.

Finally, represents a collection of additional state variables introduced in the membrane model, whereas (, ) represents the ordinary diferential equations describing the dynamics of the additional state variables.

Since there is no analytical solution to the EMI model, a numerical solution is

#### **4.3 Results**

We use code based on [6] to solve the monodomain and EMI models using fnite diference method discretisation, and study the infuence on CV. The simulations are run over a domain size of 2000µm × 40µm for the monodomain model and 1956µm×40µm×30µm for the EMI model, both with a spatial resolution of Δ = Δ = Δ = 2µm. For the EMI model the cells are arranged in a line such that all cells are connected in the x direction, similar to Figure 4.1. Each cell is comprised of fve disjoint subdomains, as depicted in Figure 4.3, and their extents are listed in Table 4.2. Regarding the time domain, the system evolves during 5 ms in steps of 0.01 ms.

The base model representing the cell membrane is described in [7]. Furthermore, the values of the most relevant parameters used in both monodomain and EMI model simulations are compiled in Table 4.1 and Table 4.2, respectively.


Table 4.1: Relevant parameter values used in the monodomain simulations.


Table 4.2: Relevant parameter values used in the EMI simulations.

Our study aims to investigate the CV dependence with ion channel properties from two diferent perspectives: when ion channel conductances change, and when ion channel distributions along the cell membrane is modifed.

First, for exploring the relation between CV and ion channel conductance, we focused particularly on the Na<sup>+</sup> , K1, Kr and CaL channels. The nominal values for the Na, K1 and Kr channel conductances are 12.6, 0.37 and 0.025 mSµF −1 , respectively, while the nominal value for CaL channel conductance is 0.12 nLµF <sup>−</sup><sup>1</sup> ms−<sup>1</sup> as specifed in [7]. Every conductance was varied from 20% to 150% of their respective nominal value by sweeping an adjustment factor from 0.2 to 1.5 in steps of 0.1.

Figure 4.2 shows the resultant CV dependence with each channel conductance in both models, the monodomain model and the EMI model.

Fig. 4.2: CV dependence with channel conductance of Na, CaL, K1 and Kr channels simulated using (a) the monodomain and (b) EMI models.

In order to address the second part of our study, we change the uniform distribution of ion channels into a non-uniform distribution along the cell membranes. Changes in local properties of cells can only be implemented using the EMI model by allowing movement of ion channels towards the cell ends, i.e., the region that is closest to the Γ1,<sup>2</sup> domain (see Figure 4.1). To run the simulations, we considered two types of channels, Na<sup>+</sup> and K1, and we explored their four corner distribution states, i.e, when both Na<sup>+</sup> and K1 are uniformly distributed, when both type of channels are completely shifted towards the cell end (see Figure 4.3), and a combination of these two. The resulting CV for each case is compiled in Table 4.3.


Table 4.3: CV (cm/s) with uniform and non-uniform channel distribution.

Fig. 4.3: In our simulations of a non-uniform distribution of ion channels, the ion channel density was increased in the brown areas near the cell ends (Ω<sup>W</sup> and ΩE) and decreased elsewhere (ΩN, Ω<sup>S</sup> and ΩO).

#### **4.4 Discussion**

#### **4.4.1 Infuence of ion channel conductance on CV**

Figure 4.2 shows that the CV increases monotonously as a function of the sodium channel conductance, Na, and this dependency is observed for both models.

To assess the compatibility between the monodomain and EMI model, the calculated CV points for Na<sup>+</sup> in Figure 4.2 are ftted into a function of the form (Na) = · Na, where and are constants, using a non-linear least-squares method. The curve adjustment is shown in Figure 4.4. For the constant , we obtain = 0.294 and = 0.3 for the EMI model and the monodomain model, respectively. Thus, consistent results were obtained with both models.

Furthermore, Figure 4.2 shows that CV remains almost constant when sweeping K + and CaL channels conductances. Therefore, varying the strength of these ion channels did not lead to signifcant changes in CV.

#### **4.4.2 Infuence of ion channel distribution**

From Table 4.3, when both K1 and Na<sup>+</sup> channels are uniformly distributed, the CV reaches 55.1 cm/s. However, when all Na<sup>+</sup> channels are placed at the cell ends (the coupling junction area between neighbouring cells) while keeping K1 channels uniformly distributed, then the CV experienced an increment of around

Fig. 4.4: Fitting the CV to the function (gNa) = · Na.

9% with respect to the previous situation. Conversely, when K1 channels are moved towards the cell ends and Na<sup>+</sup> channels are kept uniformly distributed along the cell membrane, CV decreased around 2%.

These results suggest that Na<sup>+</sup> channel distribution contributes signifcantly to the overall CV, whereas the efect of K1 distribution is relatively small. While Na<sup>+</sup> channels increased the CV as their accumulation to the cell end increased, the opposite efect was observed with a non-uniform distribution of K1 channels. Additionally, we deduce that both contributions are asymmetrical. The displacement of Na<sup>+</sup> channels causes a major impact on CV compared to the displacement of K1 channels. When both channels are non-uniformly distributed, the CV increases about 6.5% (see Table 4.3). This CV is the result of the displacement of Na<sup>+</sup> channels, which largely increases CV, mitigated by the displacement of K1 channels, which slightly decreases CV.

Immunohistochemical studies revealed that about 50% of the Na<sup>+</sup> channels are located in the membranes of the intercalated discs [8]. In the diseased heart with reduced gap junctional coupling, action potential propagation can be maintained through a mechanism known as ephaptic coupling. A prerequisite for ephaptic effects to occur is a high density of Na<sup>+</sup> channels at the intercalated disc, where the intermembrane distance between two adjacent cells is small (< 30 nm, [9]).

#### **4.5 Conclusions**

In this study, we investigated the infuence of ion channels, their conductance and their physical distribution along the cell membrane, on CV. To that end, two diferent models were used: the monodomain model and the EMI model. While the former ofers a good insight of large scale efects by reducing the cell model complexity, the latter allows the implementation of changes in local cell properties, at the expense of increased computational efort.

Regarding the CV dependence on ion channel conductance, the study focused particularly on Na<sup>+</sup> , K + (K1 and Kr) and CaL channels. Both models showed that Na<sup>+</sup> channel conductance strongly infuences CV, whereas the efect of the other ion channel conductances on CV is negligible. Moreover, the changes in CV as a result of modifying Na<sup>+</sup> and K1 channels distribution on the cell membranes, were explored with the use of the EMI model. The simulation results suggest that the infuence of these channels to the CV is opposed and asymmetrical. The infuence is considered opposed because the CV increases, when Na<sup>+</sup> channels are moved towards the cell ends, but decreases in the case of K1 channels being located at the cell ends. Furthermore, the efect on CV is asymmetrical because the movement of Na<sup>+</sup> channels along the cell membrane causes a substantial modifcation of the CV, of around 9%, compared to the transfer of K1 channels, which accounts for a variation of 2%.

**Acknowledgements** The illustrations in Figure 4.1 and Figure 4.3 were created by Karoline Horgmo Jæger and reused with permission. The authors are grateful to Jæger for proofreading the manuscript and for providing code solving the monodomain and EMI models, and to Aslak Tveito for discussions about our fndings in this study.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Chapter 5 Computational Prediction of Cardiac Electropharmacology - How Much Does the Model Matter?**

John Dawson<sup>1</sup> , Anna Gams<sup>2</sup> , Ivan Rajen<sup>3</sup> , Andrew M Soltisz<sup>4</sup> , Andrew G Edwards<sup>5</sup>


**Abstract** Animal data describing drug interactions in cardiac tissue are abundant, however, nuanced inter-species diferences hamper the use of these data to predict drug responses in humans. There are many computational models of cardiomyocyte electrophysiology that facilitate this translation, yet it is unclear whether fundamental diferences in their mathematical formalisms signifcantly impact their predictive power. A common solution to this problem is to perform inter-species translations within a collection of models with internally consistent formalisms, termed a "lineage", but there has been little efort to translate outputs across lineages. Here, we translate model outputs between lineages from Simula and Washington University for models of ventricular cardiomyocyte electrophysiology of humans, canines, and guinea pigs. For each lineage-species combination, we generated a population of 1000 models by varying common parameters, namely ion conductances, according to a Guassian log-normal distribution with a mean at the parameter's species-specifc default value and standard deviation of 30%. We used partial least squares regression to translate the infuences of one model to another using perturbations to calculated descriptors of resulting electrophysiological behavior derived from these parameter variations. Finally, we evaluated translation fdelity by performing a sensitivity analysis between input parameters and output descriptors, as similar sensitivities between models of a common species indicates similar biological mechanisms underlying model behavior. Successful translation between models, especially those from diferent lineages, will increase confdence in their predictive power.

#### **5.1 Introduction**

Preclinical drug development relies on animal-derived data for determining drug efcacy and toxicity. This is especially the case for interactions within cardiac physiology, yet there is a crisis in translating these data to human outcomes and clinical utility or liability [1]. This crisis is consequential: one third of all discontinued drugs are withdrawn from the market on account of poor safety, with cardiovascular toxicity being the leading cause of both post-approval and preclinical withdrawal [2]. In other words, cardiotoxic reactions as predicted in animal models are preventing a signifcant portion of drugs from advancing to human trials, and drugs that manage to pass approval are having unanticipated cardiotoxicity in humans. To both understand the diferences between non-human preclinical species as experimental models, and to better translate those animal-based screening results to human outcomes, many groups have developed computational models that recapitulate the cardiomyocyte electrophysiology of each species [3, 4, 5, 6]. These models are an efective and inexpensive way to predict mechanisms of action for drug activity and toxicities and have signifcant utility for human medicine.

However, there are sources of variability that hamper translating discoveries made in one computational model to another. In the case of inter-species translation, the most obvious are the signifcant and characteristic phenotypic diferences between species that refect the meaningful and real diferences in their cardiac electrophysiology. Additionally, the mathematical formulations used to describe ion channel open probabilities, gating mechanisms, or other model components chosen by the model's developers, may difer between models of the same species. Laboratories often develop a computational model for one species that uses their previous formulations to parameterize models for another species. This has created collections of models sharing common mathematical formulations, but for which parameters are varied in order to capture key physiologic diferences among diferent species. These distinct collections are often termed "lineages", and it is unclear whether their difering formulations impact cross-lineage predictions. Finally, all existing models can only represent an idealized phenotype of their respective species, further complicating accurate model-based species-translations. That is, they capture the functional properties of the average cardiomyocyte of that species. Inherent biological heterogeneity and individual variability refected in experimental measurements are not represented within the single set of parameter values used in a computational model. This valuable biological information is unused and needs to be accounted for when developing and efectively translating results between these models. This is particularly true for cardiac safety screening, where rare but lethal events are a critical outcome. Generating populations of models whose parameters vary in a way that represents both intra- and inter-individual variability is one potential method for overcoming this issue.

Prior work has attempted to translate simulated drug efects between computational models of cardiac electrophysiology, with varied success, but such work has typically only been conducted within a single model lineage or using an idealized model representing average behavior. For example, Tveito *et al.* demonstrated methods for translating drug efects from animal experiments to computational models of ventricular myocytes for a panel of pro-arrhythmic drugs [7]. Population-based modeling has recently been developed to overcome this drawback by generating sets of models from which statistical analyses can be harnessed. For example, by simulating multiple iterations of induced-pluripotent stem cell cardiomyocyte models, Gong & Sobie [8] built populations of cell models and successfully translated their interactions with drugs onto human cardiomyocyte models through application of a Partial Least Squares Regression Analysis [9]. This same approach can be generalized to models of the same species but arising from diferent lineages, or translating across both species and lineages. The only requirement is that all models have a set of common (and corresponding) parameters that can be perturbed identically to simulate the same population-level variation.

Here we have applied these approaches to systematically understand the impact of mixing models from diferent lineages on translation performance. In particular, we sought to understand whether diferences across model lineages imposed similar challenges to translation as the intrinsic electrophysiological diferences across species. We found this to not be the case, and report that, at least for two distinct model lineages, species-translation works comparatively well across lineages as within. Importantly, this general fnding did not hold for translation of action potential duration (APD) measures, and this may be due to diferent sensitivities of the model lineages to variations in the maximal transport rate of the <sup>+</sup> -Ca2<sup>+</sup> exchanger and <sup>+</sup> - + -ATPase.

#### **5.2 Methods**

#### **5.2.1 Models of Cardiac Electrophysiology**

We modeled the action potential of a single ventricular cardiomyocyte for human, canine, and guinea pig species using the Washington University [3, 4, 5] (WU) and Simula [6] lineages. These models describe currents through voltage-gated ion channels in the form = ( −), where is the channel conductance, the membrane potential, the channel's equilibrium potential, and the channel's open probability which can be a function of either membrane potential or time-dependent gating variables. Figure 5.1 illustrates the currents, fuxes, and compartments available to each model and Table 5.1 lists which of these each model represents. All Simula models share the same formalisms and only difer in exact parameter values, while the WU models, though sequentially infuenced by one another, were individually derived and feature notable diferences in the ionic currents and cellular compartments they represent. Major similarities and diferences include:

• All models split the sarcoplasmic reticulum (SR) into junctional and network compartments and the intracellular space into bulk cytosol and dyadic subspace. The two exceptions are that Simula models feature a third cytosolic compartment,

Fig. 5.1: Membrane currents (I), 2<sup>+</sup> fuxes (J), and intracellular compartments for the (A) Simula model lineage, from [6], and (B) WU model lineage, from [5]. All WU models include the 2<sup>+</sup> bufers calmodulin (CMDN) and troponin (TRPN) while only the WU human and canine models include calsequestrin (CSQN), anionic SR and sarcolemmal 2<sup>+</sup> binding sites (BSR, BSL), and CaMKII which infuences the gating of ion channels marked with a black ring in panel B. Both images in this fgure are reproduced with permissions dictated by the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

the subsarcolemma (SSL), and the WU Guinea Pig model features only one cytosolic compartment.


Baseline parameter and initial condition values were taken from each model's original publication. We generated populations of 1000 diferent parameter confgurations for each model by multiplying the baseline value of a selection of common parameters (1, , , , , , , , , ) by a scaling factor sampled from the Guassian log-normal distribution with mean of 1 and 30% standard deviation; the same scaling factors were used for each model.

Fig. 5.2: (A) Membrane potential and (B) CaT traces depicting model features of resting membrane potential (RMP), maximum upstroke velocity (dV/dt), voltage amplitude (V), action potential duration to 30, 50, and 90% repolarization (APD30, 50, 90), diastolic [2<sup>+</sup> ] (CaTd), maximum CaT velocity (dCaT/dt), time and value of maximum CaT (CaT, t), time to 50% CaT decay (CaT50), and CaT time constant ( CaT).

#### **5.2.2 Feature Extraction**

Descriptors or features of each model's resulting electrophysiology were measured during steady state activity by applying a stimulus of -80mV for 20ms to the model 500 times at a frequency of 1Hz. Features were measured as their average over the fnal 10 stimulations, and included maximum upstroke velocity, resting membrane potential (RMP), voltage amplitude, 2<sup>+</sup> transient (CaT) amplitude, maximum velocity, time constant, time to peak, and time to 50% decay, diastolic [2<sup>+</sup> ], and APD at 30, 50, and 90% repolarization (Figure 5.2). Model confgurations were omitted from subsequent analysis if a convergent solution to their system of equations could not be reached or if their steady state behavior featured a voltage amplitude less than 5mV, resting membrane potential greater than -20mV, APD90 standard deviation greater than 10%, or change in CaT amplitude of more than 2% for total SR [2<sup>+</sup> ], cytosolic [2<sup>+</sup> ], or cytosolic [<sup>+</sup> ].


Table 5.1: Currents represented in each model.

#### **5.2.3 Sensitivity Analysis and Translation**

Sensitivity analysis is a quantifcation of the correlation between input model parameters and resulting electrophysiological features and thus determines the functional dependence of model features on specifc ionic conductances. Translation correlates features between pairs of models and thus translates the electrophysiological response of diferent models to identical parameter perturbations. Projection to latent structures, also referred to as partial least squares (PLS), is a form of matrix decomposition used for sensitivity analysis and translations of the models. The nonlinear iterative partial least squares algorithm [9] was used to calculate PLS regression coefcients, referred to as the B matrix from the matrix decomposition Y = XB. The decomposition for sensitivity analysis in matrix format is = , and for translations is = . For sensitivity analysis, the set of confgurations for a model and the corresponding set of features are used as the and matrices resulting in a matrix for each of the six models in the two lineages. For translations, there are ffteen matrices, one for each pair of models.

The matrix represents the linear relationships between and , meaning we cannot determine features from parameters or translate without loss of information quantifed by the R2 values computed on the diference between the actual data values and the values predicted from multiplying a data set with its applicable matrix. In practice, the matrix can be utilized on parameter confguration vectors outside of the population to predict and translate linearly approximated feature vectors without solving the lumped conductance ordinary diferential cardiac electrophysiological model, simplifying the analysis of other datasets.

#### **5.3 Results**

#### **5.3.1 Model Translation**

Translation performance, as measured by the average <sup>2</sup> of all feature translations for each model combination, was generally high with all average <sup>2</sup> values greater than 0.75 and most close to 1 (Figure 5.3A-B); model lineage, species, or translation direction did not appear to signifcantly afect this metric. When averaging 2 s based on whether the translation was across or with lineages, inter-lineage translations only modestly underperformed compared to intra-lineage translations with respective 2 values of 0.859 ± 0.038 and 0.793 ± 0.073 (Figure 5.3C). Focusing on a subset of the translations for all combinations of WU, Simula, human, and canine models, this general trend no longer holds (Figure 5.3D-G). Instead, all <sup>2</sup> values are well above 0.8 except for measures of APD (<0.6) which was translated well for only the Simula canine-human translation (Figure 5.3G).

#### **5.3.2 Translation Discrepancies**

A likely source of this discrepancy may come from irreconcilable diferences in AP morphology between human and canine models (Figure 5.4A). Regardless of lineage, canine model APs tend to feature the classic notch-dome shape due to transient repolarization from inactivation of depolarizing inward sodium currents and activation of the transient outward potassium and sodium-calcium exchanger currents [10]. The WU canine model had the largest notch, causing some confgurations to reach 30% of repolarization immediately following peak voltage, thus resulting in biphasic distributions of APD (Figure 5.4B-D). With this behavior entirely absent from human models, linear translation between the two species was signifcantly impeded (Figure 5.4E-G).

Sensitivity analyses of this subset of models revealed that the sign of APD regression coefcients for and were the same for all models (Figure 5.5A-F, J-L) except WU canine (Figure 5.5G-I). This diference indicates that the WU canine model will have an entirely opposite change to APD given the same perturbation to these conductances. Such a fundamental diference in model behavior may also provide an explanation for poor inter-species translation.

Fig. 5.3: Translation performance as measured by their resulting <sup>2</sup> values. (A-B) Mean ± standard deviation <sup>2</sup> values over all feature translations for each combination of human (h), canine (c), and guinea pig (gp) models. (C) Mean ± standard deviation <sup>2</sup> values for translations within and across lineages for all species. (D-G) <sup>2</sup> values for each feature of WU and Simula canine-human translation.

#### **5.4 Discussion**

Computational tools are increasingly seen as a key bridge for making the inferential jumps between electrophysiologic data collected among the major cardiac cell types used in drug screening (i.e. rodent, rabbit, canine, and human stem cell-derived cardiac myocytes). Because models of these cell types contain the basic properties of each cell type and diferences between them, both linear and non-linear methods have been developed to translate among them, and importantly, to adult human cardiac cells. In this way, one holy grail of computational cardiac pharmacology is to leverage these models to reliably predict clinical outcomes of drugs based on data collected in other species or human cell lines.

Due to the inherent biological variability, it is important to use a technique that captures this variability. To this end, we developed a systematic approach of modifying the values of parameters common to all investigated models to obtain a population of 1000 confgurations for each electrophysiological model. Resulting features of model behavior extracted from each variation of the model, such as APD and CaT velocity, were used to translate between models where the behavior of one model was predicted based on the measured behavior of another. Applying

Fig. 5.4: (A) Steady state action potential traces generated using baseline parameter values for human (grey) and canine (black) models from both WU (solid) and Simula (dashed) lineages. (B-D) Histograms depicting the distributions of APD30 (29.1±48.7ms), 50 (191.0±29.8), and 90 (236.0±30.8) from the WU canine model measured as the average of the last 10 stimulations following 490s of 1Hz pacing. (E-G) Predicted APD30, 50, and 90 values for the WU human model using the regression coefcients from the WU canine – WU human translation analysis. Translation 2 values for APD30, 50, and 90 are respectively 0.5268, 0.5014, and 0.5523. Black dashed lines are the identity line.

the NIPALS method for translation provided adequate predictive power as indicated by high <sup>2</sup> values on average. Unexpectedly, there was not a large dependence of translation performance on model lineage, suggesting that their specifc formulations might not be as large a source of inherent model behavior as initially thought. It appears that interspecies diferences play a more signifcant role in afecting translation performance than the model's formulation.

As an illustrative example, we focused on canine to human translation between and within both lineages. All features were predicted well with the exception of ADP, except in the case of Simula canine to human translation which performed uniquely well. Sensitivity analyses revealed that the APD of WU canine models reacted oppositely to perturbations of conductances and compared to all other models, indicating that blockers of the <sup>+</sup> -2<sup>+</sup> or <sup>+</sup> - + -ATPase exchanger

Fig. 5.5: Normalized sensitivity coefcients, grouped by APD30, 50 and 90 (columns), for the Simula canine (A-C), Simula human (D-F), WU canine (G-I), and WU human (J-L) models.

would have an opposite efect on the WU Canine model compared to the other models we investigated. Ultimately, the poor performance of APD translation implies that translation of electrophysiological behavior between species would inaccurately predict pathologies of APD such as QT interval prolongation and cardiotoxicity afecting the repolarization phase.

#### **5.5 Conclusion**

Here, we have proposed a ventricular cardiomyocyte model translation and parameter sensitivity analysis. We demonstrated that lineage dependence is not as strong as initially hypothesized, thus interspecies translation can be accurately performed between models of diferent origin. Most model features were robustly translated within and between model lineages, however, there were shortcomings in APD translation between human and canine models. Sensitivity analysis identifed and as potential candidates for the translation inconsistency where the WU canine model had an opposite APD reaction given the same perturbation to conductances of those transporters. This anomaly requires further investigation to ascertain whether it is due to that model's specifc formulation or simply a general characteristic of the species.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Chapter 6 A Computational Study of Flow Instabilities in Aneurysms**

Nanna Berre<sup>1</sup> , Gabriela Castro<sup>2</sup> , Henrik Kjeldsberg<sup>3</sup> , Rami Masri<sup>4</sup> , Ingeborg Gjerde<sup>3</sup>


**Abstract** The majority of hemorrhagic strokes are caused by cerebral aneurysms, harboured in a large portion of the human population. Currently, it is unclear what constitutes a dangerous aneurysm, and the risk of rupture of such aneurysms is challenging to quantify. Previous studies have shown that fow dynamics play an important role in the development of aneurysms. However, there is varying consensus on blood fow patterns and stability. In this work, we formulate a Reynolds–Orr method to quantify the stability of blood fow in arteries. Applying this method to blood fow in four diferent arterial models, we fnd that the blood fow therein is unstable under physiological conditions. We show the most unstable eigenmodes for each of these models and discuss how they potentially could help explain the initiation and growth of aneurysms.

**Keywords:** fow instability, aneurysm, perturbation, Reynolds–Orr Method

#### **6.1 Introduction**

Stroke is a leading cause of death, accounting for approximately 11% of all deaths worldwide in 2015 [1], and is currently the main cause of disabilities [2]. As reported in a large cohort study in Japan, 85% of hemorrhagic strokes are caused by cerebral aneurysms [3]. A cerebral aneurysm is an outward bulging of an artery, thought to form due to a localized weakness in the arterial wall. As an estimate, 5-8% of the population harbour an aneurysm [4, 5], with an annual risk of rupture at 1-2% [6]. The rupture of such aneurysms result in non-traumatic subarachnoid hemorrhage which often leads to death or disability [4].

Unruptured aneurysms are typically detected from neurovascular imaging obtained when diagnosing some unrelated issue. Due to the increasing availability of such imaging, clinicians are also detecting an increasing number of unruptured cerebral aneurysms [4]. This observation prompts difcult clinical decisions on whether one should monitor the aneurysm without treatment or proceed with surgical intervention [4]. This raises the question of what quantitative measures reliably determine rupture risk [7].

The dynamics of blood fow are thought to play a signifcant role in the initiation and propagation of lesions along an artery [8]. In particular, computational fuid dynamics (CFD) has been used to calculate the wall shear stresses (WSS) of the blood fow in an artery. Several CFD studies have successfully used abnormal WSS as a marker in retroactively classifying aneurysms according to their risk of rupture [9]. However, there has been observed signifcant diferences between CFD results among diferent scientifc groups [10]. In particular, the prediction of fow patterns and fow instability has been found to vary, which further impacts the computation of e.g. the WSS.

In simulations of blood fow in aneurysms, the fow itself is typically assumed laminar. However, experimental evidence has shown transitional and turbulent fow to occur in blood vessels such as the aorta [11] and carotid artery [12], and in saccular aneurysms in dogs [13] and humans [14, 15]. The presence of turbulent-like fow could signifcantly change the magnitude of the WSS.

In this work, we investigate fow instabilities in patient-specifc geometries. More precisely, we use the Reynolds–Orr method to quantify the most unstable perturbation one can make to patient-specifc aneurysm geometries. The Reynolds–Orr method, as revisited by Scott in [16], provides an exact relation between a basefow and the kinetic energy of a perturbation. The analysis itself is nonlinear, with no linear approximations, and yields a well posed linear symmetric eigenvalue problem that can be solved via standard methods. Combining the open-source softwares FEniCS [17] and SLEPc [18], we fnd that the basefow and corresponding eigenproblem can be implemented in less than 150 lines of code. The eigenvector solutions represent the most unstable perturbations one can make to the basefow, while the eigenvalue indicates the growth rate of the resulting instability.

The main contributions of this report are as follows:


The rest of the paper is organized as follows: Section 6.2.1 introduces the Navier– Stokes equations with Dirichlet and traction boundary conditions. The variational formulation for this problem is also presented. The derivation of the kinetic energy relation is presented in Section 6.2.2 where the eigenvalue problem is recalled, followed by its discretization in Section 6.2.3. In Section 6.2.4, we present our methodology for solving the eigenvalue problem, followed by the results in Section 8.3. Finally, we present a discussion and potential future work in Section 8.4.

#### **6.2 Methods**

#### **6.2.1 Basefow equations**

Fig. 6.1: Two arterial models, harbouring a terminal aneurysm (left) and a saccular aneurysm (right). The aneurysm is indicated in red.

Let Ω ⊂ R <sup>3</sup> denote an open, bounded fow domain with boundary Ω. Fig. 6.1 shows the two fow domains considered. The domain Ω is assumed to have an infow boundary Γ in ⊂ Ω and an outfow boundary Γ out ⊂ Ω. We denote by Γ the union of the infow and outfow boundaries: Γ = Γin ∪Γ out .

In Ω we then consider the steady incompressible Navier–Stokes equations for the fuid velocity and pressure feld :

$$-\nu \Delta \mathfrak{u} + \mathfrak{u} \cdot \nabla \mathfrak{u} + \nabla p = f \quad \text{in} \quad \Omega,\tag{6.1a}$$

$$
\nabla \cdot \mathbf{u} = 0 \quad \text{in} \quad \Omega,\tag{6.1b}
$$

$$
\mu = \mathbf{0} \quad \text{on} \quad \partial \Omega / \Gamma. \tag{6.1c}
$$

In the equations above, is the kinematic viscosity, is a given external force, and (6.1c) is a no-slip boundary condition to be applied on a subset of the boundary. We augment (6.1a) - (6.1c) by the following traction boundary conditions on Γ, namely

$$(\nu \nabla \mathfrak{u} - pI)\mathfrak{n} = p^i \mathfrak{n}, \quad \text{on} \quad \Gamma^i \text{ for } i \in \{\text{in}, \text{out}\}, \tag{6.2}$$

where in and out are given data for the infow and outfow boundaries, respectively.

Let , ∈ 1 (Ω) 3 . For the variational formulation, we frst defne

$$a(\mathfrak{u}, \mathfrak{v}) = \nu \int\_{\mathfrak{Q}} \nabla \mathfrak{u} : \nabla \mathfrak{v},\tag{6.3}$$

$$c\left(\mathfrak{u},\mathfrak{u},\mathfrak{v}\right) = \int\_{\mathfrak{Q}} \left(\mathfrak{u} \cdot \nabla \mathfrak{u}\right) \cdot \mathfrak{v}.\tag{6.4}$$

For ∈ 1 (Ω) 3 and ∈ 2 (Ω), we defne

$$b(\mathfrak{u}, p) = \int\_{\mathfrak{U}} \nabla \cdot \mathfrak{u} \, p. \tag{6.5}$$

We also defne the following function space:

$$V\_{\Gamma} = \{ \mathbf{v} \in H^1(\mathfrak{Q})^d : \mathbf{v}|\_{\partial \mathfrak{Q} \backslash \Gamma} = \mathbf{0} \}.$$

Then, the variational problem reads: fnd (, ) ∈ <sup>Γ</sup> × 2 (Ω) such that for all (, ) ∈ <sup>Γ</sup> × 2 (Ω) the following holds.

$$a(\mathfrak{u}, \mathfrak{v}) + c(\mathfrak{u}, \mathfrak{u}, \mathfrak{v}) - b(\mathfrak{v}, p) = \int\_{\Omega} f \cdot \mathfrak{v} + \sum\_{i \in \{\text{in,out}\}} \int\_{\Gamma^i} p^i \mathfrak{n} \cdot \mathfrak{v},\tag{6.6a}$$

$$b(\mathfrak{u}, q) = 0.\tag{6.6b}$$

#### **6.2.2 Flow perturbations and instability**

In this section, we will derive an eigenvalue problem that describes diferent perturbations one can make to this basefow. This derivation is non-trivial with pressure boundary conditions. For this reason, we from now on consider Dirichlet boundary conditions for the fux for the entire boundary. In particular, let (, ) be the solution of the following time-dependent Navier–Stokes equations:

$$
\partial\_t \mathfrak{u} - \nu \Delta \mathfrak{u} + \mathfrak{u} \cdot \nabla \mathfrak{u} + \nabla p = f \quad \text{in} \quad \mathfrak{Q} \times (0, T], \tag{6.7a}
$$

$$
\nabla \cdot \mathbf{u} = 0 \qquad \text{in} \quad \Omega \times (0, T], \tag{6.7b}
$$

$$
\mathfrak{u} = \mathfrak{g} \qquad \text{on} \quad \partial \mathfrak{Q} \times (0, T], \tag{6.7c}
$$

$$
\mathfrak{u} = \mathfrak{u}^0 \quad \text{on} \quad \mathfrak{Q} \times \{0\}. \tag{6.7d}
$$

In the above, 0 is the solution of the basefow equations (6.1a)-(6.1c). We derive the kinetic energy relation for the velocity solving (6.7a)-(6.7d). The derivation presented here closely follows the derivation in [20]. Then, we state the eigenvalue problem which results from the energy relation. Let now (, ) solve the same equations (6.7a)-(6.7d) with the same boundary condition , but diferent initial data 0 , such that <sup>0</sup> ≠ 0 . Defne (, ) as the diference of the two solutions:

$$\mathbf{w} = \mathbf{u} - \mathbf{w}, \quad o = p - q. \tag{6.8}$$

Thus, (, ) satisfy:

$$
\partial\_t \nu - \nu \Delta \nu + (\mathbf{u} \cdot \nabla \mathbf{u} - \mathbf{w} \cdot \nabla \mathbf{w}) + \nabla o = \mathbf{0} \qquad \qquad \text{in} \quad \Omega \times (0, T], \tag{6.9a}
$$

$$
\nabla \cdot \mathbf{v} = 0 \qquad\qquad\text{in} \quad \mathfrak{Q} \times (0, T], \qquad\qquad(6.9b)
$$

$$\nu = \mathbf{0} \qquad\qquad\text{on}\quad\partial\Omega \times (0, T],\qquad(6.9c)$$

$$\mathbf{w}^{0} = \mathbf{u}^{0} - \mathbf{w}^{0} \quad \text{on} \quad \mathfrak{Q} \times \{0\}. \tag{6.9d}$$

Following [20], we decompose the nonlinear terms as follows

$$\begin{split} \boldsymbol{u} \cdot \nabla \boldsymbol{u} - \boldsymbol{w} \cdot \nabla \boldsymbol{w} &= \boldsymbol{u} \cdot \nabla \boldsymbol{u} - \boldsymbol{u} \cdot \nabla \boldsymbol{w} + \boldsymbol{u} \cdot \nabla \boldsymbol{w} - \boldsymbol{w} \cdot \nabla \boldsymbol{w} \\ &= \boldsymbol{u} \cdot \nabla \boldsymbol{v} + \boldsymbol{v} \cdot \nabla \boldsymbol{w} = \boldsymbol{u} \cdot \nabla \boldsymbol{v} + \boldsymbol{v} \cdot \nabla \boldsymbol{u} - \boldsymbol{v} \cdot \nabla \boldsymbol{v} . \end{split} \tag{6.10}$$

Multiplying (6.9a) by , using (6.10), and integrating over Ω yields

$$a(\partial\_t \mathbf{v}, \mathbf{v}) + a(\mathbf{v}, \mathbf{v}) + c(\mathbf{u}, \mathbf{v}, \mathbf{v}) + c(\mathbf{v}, \mathbf{u}, \mathbf{v}) - c(\mathbf{v}, \mathbf{v}, \mathbf{v}) + b(\mathbf{v}, o) = 0. \tag{6.11}$$

Defne the following function space:

$$\mathcal{V} \coloneqq \{ \boldsymbol{\nu} \in H^1(\Omega)^3 : \boldsymbol{\nabla} \cdot \boldsymbol{\nu} = 0 \quad \text{and} \quad \boldsymbol{\nu} = \mathbf{0} \quad \text{on} \quad \partial \Omega \}. \tag{6.12}$$

It can be observed from Lemma 20.1 in [16] that for any ∈ we have

$$c(\omega, \nu, \nu) = 0 \quad \forall \nu \in \mathcal{V}.$$

In addition, since ∇ · = 0, we have

$$b(\nu, o) = \int\_{\Omega} \nabla \cdot \nu o = 0.$$

Thus, we obtain the following

$$\begin{split} \frac{1}{2} \frac{\partial}{\partial t} \int\_{\Omega} |\mathbf{v}|^{2} &= -\nu \int\_{\Omega} |\nabla \mathbf{v}|^{2} - \int\_{\Omega} (\mathbf{v} \cdot \nabla \mathbf{u}) \cdot \mathbf{v} \\ &= -\nu \int\_{\Omega} |\nabla \mathbf{v}|^{2} - \frac{1}{2} \int\_{\Omega} \mathbf{v}^{l} (\nabla \mathbf{u} + \nabla \mathbf{u}^{l}) \mathbf{v}. \end{split} \tag{6.13}$$

From (6.13), we can make the following observation. The fow is energy unstable at time = 0 if there exists a <sup>0</sup> ∈ such that

$$-\frac{1}{2}\int\_{\Omega} \nu\_0^l (\nabla \boldsymbol{\mu}\_0 + \nabla \boldsymbol{\mu}\_0^l) \nu\_0 - \nu \int\_{\Omega} |\nabla \boldsymbol{\nu}\_0|^2 > 0. \tag{6.14}$$

To simplify, we defne

$$A\_{\nu} = \frac{B\_{\mathfrak{u}}(\nu, \nu)}{a(\nu, \nu)} \quad \text{with} \quad B\_{\mathfrak{u}}(\nu, \nu) = \frac{1}{2} \int\_{\Omega} \nu^{t} (\nabla \mathfrak{u} + \nabla \mathfrak{u}^{t}) \mathfrak{w} \dots$$

In the above defnition, the form (, ) is given in (6.3) and we recall that for ∈ :

$$a(\nu, \nu) = \nu \int\_{\Omega} \nabla \nu : \nabla \nu = \nu ||\nabla \nu||^2.$$

Using the above defnitions and expressions, we obtain an equivalent instability condition to (6.14). Namely, the fow is unstable at time = 0, if there is a <sup>0</sup> ∈ such that

<sup>0</sup> < −1.

In addition, we denote

$$
\lambda = \inf\_{0 \neq \boldsymbol{\nu} \in \boldsymbol{V}} \lambda\_{\boldsymbol{\nu}}.\tag{6.15}
$$

It is clear that if ≥ −1 then the fow is stable. On the other hand, if < −1, then the fow is unstable. In particular, from (6.13), there exists a such that

$$\left(\frac{1}{\|\nabla\nu\|^2}\right)\frac{\partial}{\partial t}\int\_{\Omega}|\nu|^2 = -2\nu(1+\lambda\_{\nu}) > 0.$$

We also recall Poincare's inequality which states that

$$\|\|\mathbf{v}\|\|^2 \le C\_P \|\|\nabla \mathbf{v}\|\|^2 \quad \forall \mathbf{v} \in \mathbf{V}.$$

Hence, we have the following relation:

$$\frac{\partial}{\partial t} \log \left( \int\_{\Omega} |\mathbf{v}|^2 \right) \ge -\frac{2\nu}{C\_P} (1 + \lambda\_\nu).$$

From the above, one can deduce that the most negative leads to the most unstable mode and this value is what we seek our computations.

It is shown in [20] that the solution of (6.15) solves the following eigenvalue problem. Find (, ) ∈ (R,) such that

$$B\_{\mathfrak{u}}(\nu, \mathfrak{w}) = \lambda a(\nu, \mathfrak{w}) \quad \forall \mathfrak{w} \in V. \tag{6.16}$$

For more details on the derivation of this eigenvalue problem, we refer to sections 5.2 and 5.3 in [20].

#### **6.2.3 Discretization**

We introduce the following discrete polynomial spaces to approximate solutions to the basefow equations (6.1a)-(6.1c) and to the eigenvalue problem (6.16). Let ℎ denote the space of <sup>0</sup> piecewise polynomials of degree on a regular mesh of the domain Ω. Defne the vector valued polynomial space

$$\mathcal{V}\_{\Gamma,h}^k = \{ \mathbf{v} \in (\mathcal{V}\_h^k)^3 \, : \, \mathbf{v}|\_{\partial \Omega \backslash \Gamma} = \mathbf{0} \} \dots$$

We approximate the solutions to the variational form resulting from the basefow equations: Find (ℎ, ℎ) ∈ ( 2 <sup>Γ</sup>,ℎ, 1 ℎ ) such that for all (ℎ, ℎ) ∈ ( 2 <sup>Γ</sup>,ℎ, 1 ℎ ):

$$a(\mathfrak{u}\_h, \mathfrak{v}\_h) + c(\mathfrak{u}\_h, \mathfrak{u}\_h, \mathfrak{v}\_h) - b(\mathfrak{v}\_h, p\_h) = \int\_{\Omega} f \cdot \mathfrak{v}\_h + \sum\_{i \in \{\text{in,out}\}} \int\_{\Gamma^i} p^i \mathfrak{n} \cdot \mathfrak{v}\_h,\tag{6.17a}$$

$$b(\mathfrak{u}\_h, q\_h) = 0.\tag{6.17b}$$

To approximate the solutions to the eigenvalue problem (6.16), we defne

$$\mathcal{V}\_h^k = \{ \boldsymbol{\nu} \in (\mathcal{V}\_h^k)^3 : \boldsymbol{\nu}|\_{\partial \Omega} = \mathbf{0} \} \dots$$

We fnd (, ) ∈ (R, 2 ℎ ) such that

$$B\_{\mathfrak{U}\_{h}}(\boldsymbol{\nu}\_{h}, \boldsymbol{\mathcal{w}}\_{h}) = \lambda a(\boldsymbol{\nu}\_{h}, \boldsymbol{\mathcal{w}}\_{h}) \quad \forall \boldsymbol{\mathcal{w}}\_{h} \in \mathcal{V}\_{h}.\tag{6.18}$$

#### **6.2.4 Computational Methodology**

As discussed in the introduction, the aim of this work is to numerically investigate fow instabilities in aneurysms. We investigate this by studying the fow through four diferent arterial models.

Fig. 6.2: Streamlines and vector feld in the model with (left) and without (right) a saccular aneurysm, scaled and colored by the velocity magnitude.

To compute the most unstable fow modes for each of these models, we have followed the following procedure:


Fig. 6.3: Streamlines and vector feld in the model with (left) and without (right) a terminal aneurysm, scaled and colored by the velocity magnitude.

The models have been collected from the Aneurisk database [22], and have been manually manipulated using Meshmixer [23]. To solve the basefow equations (6.1a)- (6.2) and the eigenvalue problem (6.16), we used *FlowInstabilities* [19], an opensource CFD solver developed to investigate fow instabilities.

For the basefow, we considered a laminar, steady state regime where the pressure drop was scaled so that the blood fow velocity is similar to the fow peak at systole. During the peak of the ventricular systole, the maximum distension of the artery wall occurs, and results in a greater diameter with less complacency. This lends validity to our assumption of rigid artery walls. The blood is considered to behave as a Newtonian fuid [24], with a kinematic viscosity of = 3.5 · 10−<sup>6</sup> m<sup>2</sup> /s and density of = 1060 kg/m<sup>3</sup> .

#### **6.3 Results**

Figures 6.2 and 6.3 show the numerically computed basefows. In addition to the vector feld , we also show the streamlines, i.e., the lines parallel to the velocity vector. Figure 6.2 shows the basefow computed for the domains with and without a saccular aneurysm. The pressure drop across the domain was set to 0.001 mmHg, which resulted in a basefow with peak magnitude of 50 cm/s. This corresponds well to values cited in the literature [25]. A small circulation is observed inside the aneurysm. This may indicate low values of WSS, which heightens the risk of disruption [26].

Figure 6.3 shows the basefow computed for the domains with and without a terminal aneurysm. The pressure drop across the domain was set to 0.0001 mmHg, which resulted in a basefow with peak magnitude of 10 cm/s. This is slightly lower than values cited in the literature. For increased pressure drops, however, the Newton solver did not converge. Examining the streamlines and the vector feld of Fig. 6.3a and Fig. 6.3b, one sees a large recirculation in the area where the two diferent inlets meet. It is possible that for a higher pressure drop, no steady state solution to the Navier-Stokes equation exists for this branched geometry. For the model harbouring an aneurysm (Fig. 6.3a), the streamlines show a small fow circulation inside of it.

The most unstable eigenmodes for these basefows are presented in Fig. 6.4 and Fig. 6.5. The fgures also give the corresponding eigenvalue; an increase in absolute eigenvalue implies an increase in the kinetic energy of the instability. On the top row we present the most unstable mode, and on the bottom row the second most unstable mode, with its respective eigenvalue.

Fig. 6.4: Most (top) and second (bottom) most unstable eigenmodes with corresponding eigenvalue, , for the model with a terminal aneurysm (left) and without (right).

First considering the terminal aneurysm case in Fig. 6.4, we see that for both geometries (Fig. 6.4a and Fig. 6.4b) the most unstable eigenmode has an eigenvector directed to the point of the impingement of the jet located at the bifurcation. The corresponding eigenvalues are similar in magnitude. This is not unexpected as the basefows for (6.4a) and (6.4b) are quite similar. Considering the impact of these perturbations on the basefow, the most unstable perturbation points towards the aneurysm. Recalling that an aneurysm is caused by a weakening of the blood vessel wall, this may indicate that the fow perturbations and resulting instabilities play a role in aneurysm rupture. The eigenvectors of the second mode of instability for both cases (Fig. 6.4c and Fig. 6.4d) are similarly directed to the same location, near the entrance of the aneurysm.

Fig. 6.5: Most (top) and second (bottom) most unstable eigenmodes with corresponding eigenvalue, , for the model with a saccular aneurysm (left) and without (right).

For the saccular geometry (Fig. 6.5), we fnd that the most unstable mode, for both the models - with (Fig. 6.5a) and without aneurysm (Fig. 6.5b) - are located downstream, i.e. near the outlet. This was similarly found for the most unstable fow modes for fow past a cylinder [27]. The most unstable mode is similar for the domains with and without an aneurysm. The second most unstable mode, conversely, shows a diferent behaviour for the domain with aneurysm. Considering Fig. 6.5c we observe a fow mode that is instead similar to the most unstable fow mode for the artery with the terminal aneurysm. Again, the eigenvector points toward the inside of the aneurysm. It is interesting to note that the second most unstable mode for the healthy geometry occurs in the constriction of the artery, see Fig. 6.5d. This could indicate that fow instabilities are the reason why the aneurysm occurs in this region.

#### **6.4 Discussion**

In this computational study, we used the theory of kinetic-energy instability [16] to analyze the most unstable fow perturbations one can make to blood fow in four diferent arterial models. The results indicate that blood fow in aneurysms can be unstable under physiological conditions. Moreover, the magnitude of the eigenvalue is seen to increase in the domains containing an aneurysm. This indicates that the kinetic energy of the instability increases if the domain contains an aneurysm. It is known, however, that the magnitude of the eigenvalue increases with the size of the domain, as discussed by [27]. Thus it is difcult to say if this efect is due to domain shape changes or simply occurs as the domain with aneurysm is larger.

For future work, it would be interesting to add the most unstable perturbations we computed to the basefow and see the efect of this perturbation over time. In particular, it would be interesting to examine what efect the resulting perturbations has on the computation of the WSS.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Chapter 7 Investigating the Multiscale Impact of Deoxyadenosine Triphosphate (dATP) on Pulmonary Arterial Hypertension (PAH) Induced Heart Failure**

Kristen Garcia<sup>1</sup> , Marcus Hock<sup>1</sup> , Vikrant Jaltare<sup>1</sup> , Can Uysalel<sup>2</sup> , Kimberly J McCabe<sup>3</sup> , Abigail Teitgen<sup>1</sup> , Daniela Valdez-Jasso<sup>1</sup>


**Abstract** 2-deoxy-ATP (dATP) is a myosin activator known to improve cardiac contractile force [1]. *In vitro* studies have shown that dATP alters the calcium transient profle in addition to the kinetics of the cross-bridge cycle [2]. Furthermore, *in vivo* studies of transgenic mice with increased production of dATP show elevated left ventricular systolic function [3]. Pulmonary arterial hypertension (PAH) is a rare disease of the pulmonary vasculature in which pressure overload in the right ventricle results in reduced contractile function and right heart failure [4]. We hypothesize that dATP may have a therapeutic efect on PAH-induced heart failure, by improving contractile function and restoring cardiac output and ejection fraction. However, because the efects of dATP cannot easily be assessed experimentally, we propose using a computational multiscale modeling approach to predict cardiac function. By altering parameters in an existing multiscale biventricular cardiac model [5], we were able to reproduce end-systolic and end-diastolic pressures and volumes that refect the PAH condition, as well as healthy hearts. dATP was simulated by adjusting parameters in the model at the molecular and cellular levels based on experimental data [1], allowing us to predict the efects of dATP on PAH at the organ level. Our results show that the molecular efects of dATP can increase cardiac output and restore ejection fraction in PAH conditions, though at the cost of elevated mean arterial pressure, and may provide a new approach to treating this disease. Our multiscale modeling approach paves the way for further studies mapping out cardiovascular mechanics. As novel therapeutics continue to be discovered, their application and mechanism can be further explored through these computational models.

#### **7.1 Introduction**

Pulmonary arterial hypertension (PAH) is a rare, yet aggressive disease with a 3-year survival rate of 54.9% for humans [6]. This life-threatening disease is characterized by pulmonary vascular remodeling, is represented by the narrowing of blood vessels in the lungs, and is indicated by a mean pulmonary arterial pressure (mPAP) above 20 mmHg [7]. PAH results in increased blood pressure which causes pressure overload in the right ventricle (RV) and leads to chamber remodeling [8]. The function and morphology of the RV has been found to have a very high correlation with the outcome of PAH. The RV initially adapts to the pressure overload by increasing muscle contractility and wall thickness in order to maintain cardiac output (CO), however, over time will maladapt and right heart failure (RHF) will occur leading to premature death. RHF is associated with increased flling pressures and dilation, along with reduced CO, ejection fraction (EF), and contractility [9]. Therapies such as vasodilator medications can alleviate the symptoms and improve quality of life for a short period of time [10], but do not work long term and eventually the RV will fail, leading to death.

2-deoxy-adenosine triphosphate, also known as dATP, is being studied as a potential treatment for heart failure [3]. dATP is a naturally occurring nucleotide that has been previously found to enhance crossbridge binding and cycling kinetics to improve contractility in muscle, when used to replace adenosine triphosphate (ATP) as the energy source [11]. Molecular modeling approaches have demonstrated that dATP induces allosteric changes in the myosin S1 fragment that lead to more favorable electrostatics and accelerated cycling kinetics [12, 13]. Furthermore, previous *in vivo* studies have demonstrated using a gene therapy approach to overexpress ribonucleotide reductase (RNR) [14] to restore function in animal models of heart failure, but focusing on the function of the left ventricle.

In this paper, we use a multi-scale modeling approach to determine if dATP can be used as a potential therapeutic for PAH. We combine cellular level models that include cross-bridge cycling and electrophysiology with an organ level lumped circulatory model and rat experimental data to predict cardiovascular function of the RV. We will discuss our *in silico* investigation of using dATP as a potential therapeutic for PAH with the use of multiscale models.

#### **7.2 Methods**

#### **7.2.1 Cell Level Changes**

Tewari *et al.* [15] have developed a cell level model of cardiomyocyte function described using a system of ordinary diferential equations. While computationally relatively efcient, this model incorporates the essential components such as mitochondrial energetics, thin flament activation, and detailed description of crossbridge kinetics. Furthermore, we have coupled this model to an electrophysiology model from Morotti *et al.* to produce an excitable model [16], while additionally incorporating a more detailed model of the SERCA pump [17]. As such, this model provides a means to replicate experiments via detailed simulations. In this work, we change specifc model parameters to account for physical changes caused by the presence of dATP, as compared to ATP.

#### **7.2.1.1 The SERCA Pump and Calcium transients**

The Sarco/Endoplasmic Reticulum ATP-ase (SERCA) is an active pump that transports calcium-ions from the cytosol to the sarcoplasmic reticulum (SR). SERCA removes calcium-ions from the cytosol and contributes to maintaining the resting calcium-ion concentration in the 50 - 100 nM range [18]. In normal conditions, the SERCA pump consumes energy from the hydrolysis of 1 ATP molecule to transport 2 calcium-ions across the SR membrane [19, 17]. In addition, the SERCA also transports 2 H<sup>+</sup> out of the SR for every 2 Ca2<sup>+</sup> fowing inside the SR [17, 20, 21].

Korte *et al.* [2] have shown that intracellular calcium transients undergo changes in the presence of dATP. The most notable feature of this dATP-induced efect is the shortening of the calcium decay time – defned as the time taken for the calcium-driven fuorescence indicator to reach half its peak value on stimulating the cardiomyocyte with a specifc frequency – without having a signifcant efect on its peak amplitude. Korte *et al.* suggest that altered function of SERCA due to energetic efects of dATP on the pump rate could potentially explain this phenomenon but the mechanism still remains debated.

In the present study, we focused on the H<sup>+</sup> /Ca2<sup>+</sup> countertransport in SERCA to explain the increased Ca2<sup>+</sup> pumping preferentially in the decay phase of the calciumtransient. Previous studies have reported changes in the H<sup>+</sup> /Ca2<sup>+</sup> countertransport due to conformational changes in the SERCA brought about by pH [22] and the presence of dATP [23]. To study the changes to H<sup>+</sup> /Ca2<sup>+</sup> countertransport, we varied the parameter in the SERCA model from Tran et al. This parameter gives the stoichiometric coefcient of H<sup>+</sup> binding for state transitions <sup>4</sup> → <sup>5</sup> and <sup>8</sup> → <sup>9</sup> (see Figure 7.1). Based on the results from Tran *et al.*, we varied this parameter by ±10% to ft the experimental calcium transient from Korte *et al.*

#### **7.2.1.2 Cross-bridge cycling kinetics**

Tewari *et al.* [15] use a cross-bridge model composed of six states (super relaxed, non-permissive, permissive, weakly bound, strongly bound pre-powerstroke, and post-powerstroke). By altering the rate-constant parameters that describe transitions between states, the kinetics and force generation of the cross-bridge cycle can be modifed to replicate the behavior observed in dATP conditions. In this work, we use experimental measurements from Force-pCa experiments and slack re-stretch ( ) of isolated rat trabeculae in the presence of ATP or dATP [1] to re-parameterize

Fig. 7.1: SERCA pump reaction schematic. Figure adapted from Tran *et al.*, 2009 [17] with permission from Elsevier and Copyright Clearance Center (License number 5232480317559).

our cross-bridge model, in addition to using Brownian Dynamics (BD) simulation results from [13].

Fig. 7.2: (a) Force development measurements from experimental data and model simulations. ATP conditions shown in black, and dATP conditions in blue. Forces normalized to ATP experimental measurements at pCa 4. (b) Steady state normalized force development at diferent calcium concentrations in the presence of 5 mM ATP (black) and dATP (blue) from the model and experimental data.

Based on BD simulation results, the association rate of the myosin head to actin was increased from 250 −<sup>1</sup> to 567 −<sup>1</sup> . Similarly, the dissociation rate was adjusted proportionately to remain thermodynamically constrained. In order to reduce overftting, we manually optimized the fewest number of rate constants possible that could still characterize the change in behavior as seen by dATP for the force-pCa and slack re-stretch experiments. While this model does not reproduce exact measurements from the experimental data, relative changes between conditions provide an efective means to compare behavior. Figure 7.2 shows that at diferent calcium concentrations (pCa 4 and 5.5) dATP increases the rate of force redevelopment as compared to ATP. The exact rate measurements vary between the experimental data and model results, but the relative change in the rate constant is nearly proportional at both calcium concentrations (Table 7.1).

The model was also able to predict changes to the steady state force production in force-pCa experiments when in the presence of dATP compared to ATP. Most notably, the model predicted pCa<sup>50</sup> and maximal force for dATP was signifcantly higher, matching the experimental results (Table 7.2). In total, only fve new parameters were manually optimized to match dATP model behavior with experimental, and an additional three reverse rates adjusted to maintain proper cycle thermodynamics (Table 7.3).

Table 7.1: Slack Re-stretch Measurements.




\* p < 0.05 diference between ATP and dATP condition.

#### **7.2.2 Organ Level Model**

Previously developed models from Tewari *et al. (2016)*, Bazil *et al. (2016)*, Gao *et al. (2019)*, and a whole-organ cardiac mechanics model created by Lumens *et al. (2009)* are embedded in a lumped circulatory model that is used to simulate whole-body


Table 7.3: Cross-bridge Model Parameter Changes.

<sup>∗</sup>Note:\* denotes thermodynamically constrained.

cardiovascular function. TriSeg model created by Lumens *et al. (2009)* functions as a biventricular model which simulates left- and right- ventricular mechanics. Data from each rat is matched to the model which is run for 120 heart beats. Outputs of this model include values of mean arterial pressure and ejection fraction, along with plots for predicted right ventricular PV loops for the individual animal.

Experimental research around PAH is being conducted on Sprague Dawley rats in order to investigate right ventricular function throughout the progression of this disease. PAH is induced in rats following a Sugen-Hypoxia protocol, in which Sugen5416 (a vascular endothelial growth factor receptor (VEGFR) inhibitor) is injected into the rats following a 3-week period of hypoxia (10% oxygen). The rats are then removed from hypoxia and placed back into normoxia. The rats undergo open-chest terminal surgeries at diferent weeks throughout the disease progression. For the purpose of this project, we were focused on week 21 healthy and PAH induced male rats. During these open-chest surgeries, a pressure volume (PV) catheter is inserted into the apex of the RV to create *in vivo* PV loops. Several PV loops are recorded in steady-state conditions, with the average of these loops taken for the purpose of our analysis. The same process can be repeated while the PV catheter is in the left ventricle, however RV function and morphology is an important indicator for the severity of PAH making this ventricle the focus of this work.

Using the organ level model to simulate the pressure and fows in the whole-body [5], the goal was to replicate the experimental PV loops found during open-chest surgeries. The known values from the surgery day for blood pressure, heart rate, body weight, LV weight, RV weight, stroke volume, CO, LV diastolic volume, EF, and volume were inputted into the model and can be found in the supplementary information Table 0.8. We then found a number of adjustable parameters in the TriSeg (Heart) model that could be adjusted to replicate the experimental data PV loops. For both the healthy and PAH case, the stifness of series element (), LV/RV/Septum wall volumes ( , ,), and the LV/RV/Septal midwall reference surface areas ( , , ) were adjusted to defne the mass and geometry of the heart for each individual animal. Table 0.4 shows the relative changes of each of these parameters for the purpose of this work. The stifness in the PAH model is 3-fold greater than that of the healthy case, representing a change that occurs in the RV during this disease. There are no original values for the , , parameters as these vary rat to rat.

We also adjusted inputs for each rat, such as the on rate constant for super relaxed state (), model parameter for force dependent super relaxed transition ( ),


Table 7.4: Adjusted Parameters for TriSeg (Heart) Model.

<sup>∗</sup>LVW is the weight of the left ventricle and RVW is the weight of the right ventricle.

ATP hydrolysis rate ( ), and the resistances across PAH and healthy rats (, ). These adjusted values, along with the inputs for each rat can be found in the supplementary information. Although, we were not able to replicate the loops exactly, we were able to get representative loops from the model that can be used to show relative changes from healthy to diseased PV loops. We were able to generate representative PV loops using the model that demonstrate relative changes between healthy and diseased conditions.

#### **7.3 Results**

We studied the impact of altered H<sup>+</sup> /Ca2<sup>+</sup> countertransport on intracellular calciumtransient by varying the parameter from the Tran *et al.* [17] SERCA model. We found that an 8% increase in this parameter (from 2 → 2.16) was able to reproduce the shortening of decay-time observed in Korte *et al.* [2] without afecting the amplitude of the calcium-transient (Figure 7.3a). This change is consistent with the fndings from Tran *et al.* which predicts ±10% [17] change to during conformational changes to the SERCA pump. A caveat, however, is that we were able to observe a decrease in the decay-time of the calcium-transient by 33.4% as compared to about that of ≈ 50% in Korte *et al.* [2] (see Figure 7.3b). This discrepancy could likely be due to diference in concentrations of metabolites like MgATP and MgADP between the experiment and the model and needs further investigation. Taken together, these results indicate that altered H<sup>+</sup> /Ca2<sup>+</sup> countertransport could potentially explain the results from Korte *et al.* [2] of the shortening of intracellular calcium transient decay-time without afecting the amplitude.

We combined the optimized parameters from the coupled SERCA and electrohysiology with the altered sarcomere XB model to predict functional changes in cardiac muscle twitch. Because the XB cycling parameterization was carried out in 100% dATP conditions, while calcium handling measurements were carried out in approximately 2% dATP concentration, the optimized parameters are used in two diferent simulated conditions as seen in Figure 7.4. When the XB cycle parameters are used to simulate the 2% dATP case, they are scaled 2% linearly interpolated to the new ftted parameters in Table 7.3. The parameter remains fxed at 2.16 in both cases. Comparable with experimental twitch observations, our results show that there is increased force production in the presence of dATP, as a result of altered XB kinetics. Furthermore, the twitch has an accelerated relaxation in the presence of dATP compared to ATP, due to adjustments in SERCA calcium handling parameters.

Fig. 7.3: (a) Calcium transient for 1 Hz stimulation for diferent values of . (b) Comparison of calcium transient from model prediction with Korte *et al., 2011* data.

Fig. 7.4: Simulated cardiomyocyte twitch in the presence of ATP (black) and dATP (blue). The dotted line represents 2% simulated dATP and solid line represents 100% simulated dATP.

Mean arterial pressure (MAP) otherwise known as mPAP, is the determining metric for whether PAH has developed (MAP > 20 mmHg) or not (MAP < 20 mmHg). Other important factors when determining the function of the heart include EF, SV, and CO. Compared to the healthy case, in heart failure the MAP increases,

Fig. 7.5: Experimental rat right ventricle pressure volume loops for (a) Healthy and (b) PAH cases.

Fig. 7.6: (a) Model healthy rat right ventricle pressure volume loop. (b) Model PAH rat right ventricle pressure volume loops, without dATP treatment (yellow) and with dATP treatment (blue).

Table 7.5: Experimental versus Model Outputs for a Healthy Rat.



Table 7.6: Experimental versus model outputs for a PAH Rat.

accompanied by a decrease in EF, SV, and CO. Although we were unable to replicate the experimental PV loops perfectly using the model, we were able to demonstrate the increase in MAP and decrease in EF, SV, and CO as shown in Table 0.5 and Table 0.6. Parameter tuning around calculating the EF requires further investigation, as the values calculated were 1.5-fold greater for the model in the healthy case and almost 3-fold greater in the PAH case.

Organ scale simulations in the presence of dATP showed variable therapeutic potential to treat PAH depending on dosage (Figures 7.5b and 7.6). When PV loops were simulated in the presence of 2% dATP, there was a reduction of MAP, likely due to the accelerated diastolic relaxation, however, the cardiac output in addition to ejection fraction were also reduced (Table 7.7). When the PAH condition was simulated with 100% dATP, cardiac output was restored and ejection fraction improved signifcantly, yet this came at the cost of further elevated MAP. The elevated MAP indicates dATP may not be an ideal therapy to treat PAH. Future work will focus on identifying a critical dosage of dATP to take advantage of the improved diastolic relaxation and increase cardiac output.



#### **7.4 Discussion and Conclusion**

In this paper, we have demonstrated that the cellular level models provide a framework to efectively simulate the behavior of dATP. Considering the interplay between our model and experimental data, there is a quantitative agreement of the theoretical results with the experimental observations, suggesting that our model captures the essential physics involved in the mechanics of PAH induced heart failure. However, there is a variation between experimental data and mathematical model. From these results, we can conclude the model replicates the experimental data with better precision for the PAH rat case compared with the healthy rat, although further parameter tuning is needed to advance the capabilities of the whole-organ model. Furthermore, the cell level dATP parameterization may require more detailed scaling approaches to accurately refect how contraction changes as a function of dATP concentration. The next opportunity, in our view, lies in exploring this behaviour using the addition of noise terms by transforming our ODE model to Stochastic Diferential Equation (SDE) model. A stochastic model will provide a method to study the beat to beat variability, and better characterize model stability.

From a modeling standpoint, series element resistance and cross-bridge kinetics play important roles in the mechanics of dATP on PAH induced heart failure. On the other hand, our model is not sensitive to variation in stifness parameters and compliance parameters which are proximal aortic compliance, systemic arterial compliance, systemic venous compliance, and pulmonary arterial compliance. Because PAH is characterized by stifening of the right ventricle, this is a signifcant limitation, and a more detailed model may be necessary to more accurately represent PAH. We showed that while this model can individually characterize the efects of PAH and dATP, it appears that dATP may have variable efcacy depending on the PAH conditions. As such, dATP may not be an ideal candidate to treat PAH, specifcally given that dATP increases force and pressure in a high pressure disease. Additionally animal specifc simulations are necessary for a robust prediction of dATP on function. This multiscale modeling approach shows that dATP is a promising therapeutic to treat PAH, however further detailed modeling and likely *in vivo* studies are necessary next steps.

#### **7.5 Acknowledgements**

The authors would like to acknowledge Prof. Andrew McCulloch, Dr. Kimberly McCabe, and Abby Teitgen for providing their critical comments and feedback for the project.

#### **7.6 Supplementary Information**


Table 7.8: Experimental Inputs for Rats on 1..

Table 7.9: \_ \_\_ Inputs.


#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Chapter 8 Identifying Ionic Channel Block in a Virtual Cardiomyocyte Population Using Machine Learning Classifers**

Bjørn-Jostein Singstad<sup>1</sup> , Bendik Steinsvåg Dalen<sup>2</sup> , Sandhya Sihra<sup>3</sup> , Nickolas Forsch<sup>4</sup> , Samuel Wall<sup>4</sup>


#### **Abstract**

Immature cardiomyocytes, such as those obtained by stem cell diferentiation, have been shown to be useful alternatives to mature cardiomyocytes, which are limited in availability and difcult to obtain, for evaluating the behaviour of drugs for treating arrhythmia. *In silico* models of induced pluripotent stem cell-derived cardiomyocytes (iPSC-CMs) can be used to simulate the behaviour of the transmembrane potential and cytosolic calcium under drug-treated conditions. Simulating the change in action potentials due to various ionic current blocks enables the approximation of drug behaviour. We used eight machine learning classifcation models to predict partial block of seven possible ion currents (, , , 1, , and ) in a simulated dataset containing nearly 4600 action potentials represented as a paired measure of transmembrane potential and cytosolic calcium. Each action potential was generated under 1 pacing. The Convolutional Neural Network outperformed the other models with an average accuracy of predicting partial ionic current block of 93% in noise-free data and 72% accuracy with 3% added random noise. Our results show that and current block were classifed with high accuracy with and without noise. The classifcation of , <sup>1</sup> and current block showed high accuracy at 0% noise, but showed a signifcant decrease in accuracy when noise was added. Finally, the accuracy of and classifcation were relatively lower than the other current blocks at 0% noise and also showed a signifcant drop in accuracy when noise was added. In conclusion, these machine learning methods may present a pathway for estimating drug response in adult phenotype cardiac systems, but the data must be sufciently fltered to remove noise before being used with classifer algorithms.

#### **8.1 Introduction**

Drugs that act on the cardiovascular system may cause severe arrythmias, therefore animal or *in vitro* models are widely used for testing to validate the drug efect. The use of human induced pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs) and microphysiological systems represents a new era in assessing drug-induced cardiotoxicity and screening of drug efects in microtissues [1, 2]. Even though hiPSC-CMs are one of the most developed and well-characterized model systems derived from induced pluripotent stem cells (iPSCs), there are still some limitations. Intraindividual variability, the lack of reproducibility of results across multiple laboratories, and a long maturation process are some of the biggest challenges in this feld [2]. A recent study found the maturation process to take up to 8 months [3]. Another signifcant limitation is the small sample size that is being used in current studies [2], which limits the feasibility of using machine learning models. In this paper we show how simulated action potentials (AP) from *in silico* models can be used to train supervised machine learning models to classify partial ionic current blocks. Despite using simulated data sampled from a virtual population of hiPSC-CMs, this study demonstrates the feasibility of using machine learning models to classify ionic current block from action potential recordings.

#### **8.2 Methods**

This section describes the data used in this study and some pre-processing applied to the data. The diferent model architectures tested and the processes of model selection and hyperparameter tuning are described. Finally, the diferent scoring metrics used to analyze classifer performance are presented and explainable AI is discussed.

#### **8.2.1 Data**

The data set used in this study contained 4582 simulated action potentials from hiPSC-CMs under both normal conditions (control case, not drug-treated) and drugtreated conditions. Partial block of ionic currents was simulated in the cell model by reducing the maximal conductances of the corresponding channels. Each signal under drug-treated conditions was subtracted from the paired signal under normal conditions such that the resulting signal represented the change in the action potential due to drug treatment.

The seven current blocks that were simulated in this study were , , , 1, , and . Figure 8.1 presents an overview of an action potential with diferent degrees of ion current block.

Fig. 8.1: Action potentials measured in terms of transmembrane voltage and cytosolic calcium with diferent percentages of ion current block. The f ion current block (shown in row three and four in the third column) is not present in the dataset used in this study.

Each simulated action potential was generated under 1 pacing and was represented as two signals: transmembrane voltage () and cytosolic calcium (Ca2<sup>+</sup> ). The two signals were represented as time-series with a length of 200, with 5ms between each time point.

The ground truth labels in the data set assumed that the ion channel was either blocked (1) or not blocked (0). Considering all seven ion currents, this gives a theoretical total of 128 unique combinations. All of the possible combinations were represented in the data set except for the case in which all currents are completely open, resulting in an actual total of 127 combinations.

$$2^{\mathcal{T}} - 1 = 127$$

The combinations were not equally represented in the data set; the least represented combination of ion current blocks had 27 examples, while the most represented combination had 40 examples. The mean and standard deviation were 36.1 and 2.1 respectively, which makes the data set slightly imbalanced and can translate to imbalanced classifer training.

#### **8.2.2 Preprocessing**

#### **8.2.2.1 Noise**

To provide more realistic signals as input to the classifcation models, random Gaussian noise was added to the simulated and Ca2<sup>+</sup> signals. Noise (within a certain percentage of the standard deviation of the mean across all signals in the data set) was incorporated into the raw signal by adding a vector of random noise of equal length as the signal to the signal itself. Listing 8.1 shows a code snippet of the noise generation function. 1

Listing 8.1: The Python function that was used to add noise to the ion current signals.

**d e f** a d d \_ n oi s e (X, p e r c e n t = 5 . 0 ) : s t d = np . nanmean (X, a x i s = 0 ) . s t d ( ) n o i s e = np . random . n o rmal ( 0 , st d ,X. s h a p e )∗ p e r c e n t / 1 0 0 X \_ n oi se = X + n o i s e **r e t u r n** X \_ n oi se

Noise levels of 0, 1, and 3 percent were added to all data and tested with all classifer models in this study.

#### **8.2.2.2 Normalizing**

After the addition of noise, the signals were normalized using min-max normalization. Equation 8.1 shows how , which is the -th element in X, can be normalized between 0 and 1 using the largest value in X, , and the smallest value, .

$$
\tilde{\chi}\_i = \frac{\mathbf{x}\_i - \mathbf{x}\_{min}}{\mathbf{x}\_{max} - \mathbf{x}\_{min}} \tag{8.1}
$$

<sup>1</sup> The rest of the code developed in this study is available here: https://github.com/ SSCP2021-group-9/hiPSC-CMs-ionic-current-block-detecion

#### **8.2.2.3 Subtract drug signals from control signals**

After the addition of noise and the normalization of the signals between 0 and 1, the signal from each ion's blocked case was subtracted from the signal from the control case. The same was done for the Ca2<sup>+</sup> signal.

#### **8.2.2.4 and Ca2**<sup>+</sup> **concatenation**

Finally the and the Ca2<sup>+</sup> signals were concatenated, meaning both and Ca2<sup>+</sup> for one simulated AP were appended into a single array. Figure 8.2 shows an example of the concatenation for all current blocked signals with and without noise.

(a) Signal derived from transmembrane voltage with Kr block and 0% noise.

(b) Signal derived from transmembrane voltage with Kr block and 3% noise.

(c) Signal derived from cytosolic calcium with Kr block and 0% noise.

(d) Signal derived from cytosolic calcium with Kr block and 3% noise.

Fig. 8.2: Signals derived from the transmembrane voltage and cytosolic calcium signals after the preprocessing steps for the channel. Figures (a) and (c) show the signals without noise, while (b) and (d) show the same signals with noise.

#### **8.2.3 Multi-label classifcation methods**

This section explains the three methods used to transform single-label classifcation methods into multi-label methods. These methods were implemented using Scikitmultilearn [4].

#### **8.2.3.1 Binary relevance**

The binary relevance method is considered to be the simplest strategy for problem transformation. It converts a multi-label problem into several independent binary classifcation problems [5].

#### **8.2.3.2 Classifer chains**

Similar to the binary relevance method, the classifer chain trains one model per class [6]. However, a notable diference between the methods is that all previous predictions are taken in addition to the input features to predict the next model. The -th classifer in the chain is trained on the original input features in addition to the −1 classifcation results.

#### **8.2.3.3 Label Powerset**

A label powerset problem transformation converts a multi-label classifcation into a multi-class problem with one multi-class classifer trained on all unique label combinations found in the training data [4]

#### **8.2.4 Model architectures**

#### **8.2.4.1 Gaussian Naive Bayes**

Naive Bayes (NB) is a series of classifcation methods that relies on Bayes theorem to make predictions, and Gaussian NB (GNB) is a variant of NB in which the data is assumed to follow a normal distribution.

If the input signal is **x**, the likelihood of the signal belonging to a class is

$$p\left(C\_k \mid \mathbf{x}\right) = \frac{p\left(C\_k\right)p\left(\mathbf{x} \mid C\_k\right)}{p\left(\mathbf{x}\right)} \propto p\left(C\_k\right) \prod\_{i=1}^n p\left(\mathbf{x}\_i \mid C\_k\right). \tag{8.2}$$

Since it would be the same for each class, (**x**) can be ignored. The proportion of training samples belonging to the class is (). Further, the mean value and the standard deviation for a class can be used in the formula for Gaussian probability distribution to calculate ( | ):

$$p\left(\mathbf{x}\_{i}\mid\mathcal{C}\_{k}\right) = \frac{1}{\sigma\_{k}\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{\mathbf{x}\_{i}-\mu\_{k}}{\sigma\_{k}}\right)^{2}}\tag{8.3}$$

The input is assumed to belong to the class with the highest likelihood [7].

This study used 3 versions of the GNB classifer: one version using the binary relevance method, another using the classifer chain method, and the last version using the label powerset method. The classifers were implemented using sci-kit learn and scikit-multilearn [8, 4]. During training of these 3 GNB models the hyperparameters shown in Table 8.1 were tuned.


Table 8.1: Naive Bayes hyperparameters.

#### **8.2.4.2 Support Vector Classifer**

A Support Vector Machine (SVM) is a powerful method in machine learning. It is capable of doing linear and nonlinear classifcation, regression, and even outlier detection. A SVM used for classifcation purposes is sometimes called a Support Vector Classifer (SVC). SVCs are particularly well suited for classifcation of complex but small or medium-sized data sets. The mathematics behind the SVM/SVC relies on the defnition of hyperplanes and the defnition of a margin which separates classes of variables [9].

In this study the SVC was implemented in the label powerset algorithm using sci-kit learn. Table 8.2 shows the hyperparameters tuned in this study and the values used in the search space.


Table 8.2: Support Vector Machine hyperparameters.

#### **8.2.4.3 XGBoost**

XGBoost is a scalable machine learning system for tree boosting [10]. This algorithm has shown promising results in many classifcation tasks. In this paper the XGBoost classifer is used inside the previously described label powerset algorithm. The model was implemented using a Python package developed by Chen, T. & Guestrin, C. 2016 [10].

Table 8.3 contains the tuned hyperparameters and the values used in the search space.


Table 8.3: XGBoost hyperparameters.

#### **8.2.4.4 Feed Forward Neural Network**

Feed forward neural networks (FFNN), also called fully connected NN and multilayer perceptrons, are the simplest form of neural networks. They consist of many nodes ordered into layers. The frst layer is the input layer, which would contain the input features used in making a prediction. The last layer is the output layer, where each node represents the prediction for each of the labels. There are several layers between the input and output layers called hidden layers. Each node is connected to all the nodes in the previous layer, meaning that the value of a node depends on all the nodes in the previous layer. The strength of the connection depends on the internal weights and biases of the network, which are static for all predictions. The FFNN is trained using a process called backpropagation, which updates the internal variables using the gradient of the loss for a prediction [11].

An example of the structure of a FFNN can be seen in Figure 8.3.

The neural network was implemented in Python using Sci-kit learn [8]. Table 8.4 shows the hyperparameters that were tuned during the model development.

#### **8.2.4.5 Convolutional Neural Network**

The convolutional neural network used in this study was inspired by the the encoder model used in Fawaz HI et al 2019 [11]. The architecture of the encoder model is shown in Figure 8.4. The CNN model had 7 output neurons, equaling the number of

Fig. 8.3: An example of the structure of a feed forward neural network. Each node in one layer is connected to all the nodes in the next layer. Here we have input features, 2 hidden layers with 6 nodes each, and 4 output nodes.


Table 8.4: Neural Network hyperparameters.

classes in the dataset. A Sigmoid activation function was used in the fnal layer and binary cross-entropy was used as the loss function (Equation 8.4). The model was implemented in Python using Keras and Tensorfow [12, 13]. The hyperparameters tuned during training are shown in Table 8.5.

$$\mathcal{L}(\mathbf{y}, \tilde{\mathbf{y}}) = - \left( \mathbf{y} \log \tilde{\mathbf{y}} + (1 - \mathbf{y}) \log \left( 1 - \tilde{\mathbf{y}} \right) \right) \tag{8.4}$$


Table 8.5: Convolutional Neural Network hyperparameters.

Fig. 8.4: The architecture of an encoder model. Each block illustrates a layer or a mathematical operation. The fgure is obtained from B.Singstad et al. 2021 [14] with permission from the author.

#### **8.2.4.6 Recurrent Neural Network**

We developed a recurrent neural network (Figure 8.5) using a bidirectional long short-term memory (LSTM) block, a normal LSTM block, a dense layer with seven sigmoid outputs (one for each class), and a binary cross-entropy function as the loss function (Equation 8.4). The model was implemented in Python using Keras and Tensorfow [12, 13]. The hyperparameters that were tuned during training are shown in Table 8.6.

Fig. 8.5: The architecture of the recurrent neural network. Each block illustrates a layer or a mathematical operation.


Table 8.6: Recurrent Neural Network hyperparameters.

#### **8.2.5 Model selection and hyperparameter tuning**

The best models and hyperparameters were selected using nested cross-validation. A 10-fold stratifed cross-validation was used as the outer loop and a regular 3 fold cross-validation was used as the inner loop. The number of inner loop crossvalidations was equal to the number of possible combinations of hyperparameter settings in the search space. After being selected in the inner loop, the optimal hyperparameters were applied to the whole training data set in the outer loop.

#### **8.2.6 Scoring and metrics**

Accuracy, recall and precision were used to assess the performance of the eight models in this study.

#### **8.2.6.1 Accuracy**

The accuracy for each single channel (Equation 8.5) and the average accuracy across all channels (Equation 8.6) are both reported for each model. The accuracy score compares the predicted label (ˆ) with the true label (). In Equations 8.5 and 8.6, is the number of samples in the data set.

$$
accuracy(\mathbf{y}, \hat{\mathbf{y}}) = \frac{1}{n\_s} \sum\_{i=0}^{n\_s - 1} 1 \cdot (\hat{\mathbf{y}}\_i = \mathbf{y}\_i) \tag{8.5}
$$

In Equation 8.6, is the number of classes from which the mean is calculated.

$$
bar{v}\_{c}(\mathbf{y}, \mathbf{j}) = \frac{1}{n\_c} \sum\_{j=0}^{n\_c - 1} \frac{1}{n\_s} \sum\_{i=0}^{n\_s - 1} 1 \cdot (\mathbf{j}\_{ci} = \mathbf{y}\_{ci}) \tag{8.6}$$

#### **8.2.6.2 Recall and precision**

Precision and recall are two scoring metrics that give a number between 0 and 1. The precision score calculates the number of positive predictions that actually belong to the positive class, as seen in Equation 8.7. Recall, on the other hand, gives a score that quantifes the ratio of positive predictions made to all positive labels in the data set, as seen in Equation 8.8.

$$\text{Precision} = \frac{\text{TP}}{\text{TP} + \text{FP}} \tag{8.7}$$

$$\text{Recall} = \frac{\text{TP}}{\text{TP} + \text{FN}} \tag{8.8}$$

#### **8.2.7 Explainable AI**

Advanced AI models, such as CNN models, have been traditionally considered "black boxes" because until recently it has been impossible to explain their predictions [15]. However, a new discipline called explainable AI is making the decision basis more transparent. Explainable AI has the potential to reveal new features that could be of clinical importance. For example, explainable AI has been used to reveal new features on the electrocardiogram (ECG) that can potentially be used to detect particular gene mutations (Figure 8.6) [16]. We hypothesized that similar techniques would be helpful to highlight patterns in the data set used in this study, so local model agnostic explanations (LIME) were used to explain the predictions of the most advanced models.

#### **8.2.7.1 LIME (Local Interpretable Model-Agnostic Explanations)**

An example of a popular Python framework that can be used to implement local model agnostic explanations is LIME. By observing the classifcation result of a model, LIME makes small changes to the input data and trains a linear and explainable surrogate model. This process reveals the relative importance of the diferent properties in the input data [18].

In this study the linear surrogate model was trained to explain the CNN model. The surrogate model was trained on the same data set as the CNN model and validated on the same model as the CNN model.

#### **8.3 Results**

Figure 8.7 shows the results of the outer loop in the nested cross-validation, which was used to assess the accuracy of the models. The CNN model achieved the highest accuracy across all channels at three levels of noise with an average accuracy of 93% at 0% noise, 81% at 1% noise and 72% at 3% noise.

Fig. 8.6: Explainable AI techniques used on ECGs to explain features such as PR interval (a), QT-interval (b), QRS segment (c), J-point elevation (d), T-wave amplitude (e), R-wave amplitude (f) and heart rate (g). This fgure is taken from [17] with permission from the authors and the publisher as determined by the Creative Commons CC BY license (https://creativecommons.org/licenses/).

#### **8.4 Discussion**

This study shows that signals derived from simulated action potentials can be used as inputs to train a supervised machine learning model and to classify partial block of specifc ionic currents. The model that performed best in this study was a Convolutional Neural Network classifer, which achieved an average accuracy across all ion currents of 93% at 0% noise, 81% accuracy at 1% noise and 72% accuracy at 3% noise.

The fndings shown in Figure 8.7 state that the CNN model had an accuracy close to 100% in classifying partial block of , 1, , and when no noise was added. The 1, and currents showed a signifcant drop in accuracy when noise was added. Conversely and showed little or no decrease in accuracy when noise was added. It seems that the classifcation of and is

Fig. 8.7: Classifcation accuracy for each individual ion channel represented as boxplots. The values shown represent the 10-fold cross-validated accuracy.

less vulnerable to noise than the 1, and currents. It is possible that, as shown in Figure 8.1, and current block have a relatively large impact on the morphology of the action potential, while blocking 1, and has less impact on the action potential. This observation may support the hypothesis that classifcation of ion current block is increasingly vulnerable to noise with decreasing impact of the current on action potential morphology.

The mean cross validated accuracy of and current block classifcation were 84% and 72% even at 0% noise. This was somewhat surprising because the changes in the action potential when and currents are blocked seem to be equally large or larger than action potentials with and block in Figure 8.1. Experiments in this study show that the class activation for is high for , and blocks. The block also shows high activation in signals labeled with , and .

(a) block interpreted as block with a probability of 0.90.

(b) block interpreted as block with a probability of 0.98.

Fig. 8.8: Activation map of the CNN classifer model for a time series input corresponding to sample cases with partial block of the current only. The green line is the input to the classifer, representing the change in voltage and intracellular calcium due to partial block of . Red regions along the time series indicate a positive activation for block (label is 1), while blue regions indicate negative activation for block (label is 0).

The high accuracy of the best classifer model, despite a large degree of variability within the virtual hiPSC-CM population, demonstrates the robustness of this machine learning approach. Figure 8.8 shows two sample cases of partial block where the CNN classifer model correctly predicted the class with a probability of over 96%. Although the intracellular calcium signals have a diferent morphology, especially in the rate of repolarization, the classifer identifed the correct channel block for both cases with a high degree of confdence. This is likely due to the similarities in the input time series, − and 2<sup>+</sup> - 2<sup>+</sup> (green line, Figure 8.8). Furthermore, the classifer was able to make a correct prediction based on an input with small values relative to the original voltage and intracellular calcium signals; the control and drug-treated traces were nearly indistinguishable to the human eye. This suggests the classifer may be robust for cases with both small and large degrees of channel block.

Figure 8.8 also shows the activation map of the CNN classifer model for both cases of partial block, with red regions of the time series corresponding to positive activation for block and blue regions corresponding to negative activation for block. The darkest, highest frequency of red regions are found in the late phase of repolarization for both cases, suggesting that this region is the most important for determining class.

Figure 8.8a shows that the no-block classifcation is most pronounced for the voltage signal, while the block classifcation is most pronounced for the calcium signal. In Figure 8.8b the opposite is shown; the no-block classifcation is most pronounced for the calcium signal, while the block classifcation is most pronounced for the voltage signal. These conficting results could suggest that both signals are relatively equal in importance when making a prediction.

The activation map provides meaningful insight into the input of the classifer and can help guide strategies for improving data collection and post-processing, ultimately improving prediction of the specifc drug efect in question.

#### **8.5 Conclusion**

The purpose of this study was to investigate whether machine learning models are capable of detecting partial block of specifc ionic currents from the change in morphology of an action potential. The results of this investigation show that a CNN model achieved the highest average accuracy in classifying ionic current block. The study also shows how the classifer performance is afected due to various levels of noise in the data. Although this study is based on simulated action potentials, the fndings show the importance of fltering noise from measured signals, which may have implications in future studies that use action potential recordings of *in vitro* hiPSC-CMs. The study also showed that explainable AI methods can provide meaningful insight into the input of the classifer and can help guide strategies for improving both data collection and post-processing.

#### **References**


in Phospholamban Gene Mutation Carriers. *Circulation: Arrhythmia and Electrophysiology*, 14(2), February 2021. Number: 2 Place: United States.


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.