# **SIMULA SPRINGER BRIEFS ON COMPUTING 7** REPORTS ON COMPUTATIONAL PHYSIOLOGY

Aslak Tveito Kent-Andre Mardal Marie E. Rognes Editors

# Modeling Excitable Tissue The EMI Framework

# **Simula SpringerBriefs on Computing**

Reports on Computational Physiology

Volume 7

**Editor-in-Chief** Aslak Tveito, Fornebu, Norway

# **Series Editors**

Kimberly J. McCabe, Fornebu, Norway Rachel Thomas, Fornebu, Norway Andrew D. McCulloch, La Jolla, USA

In 2016, Springer and Simula launched an Open Access series called the Simula SpringerBriefs on Computing. This series aims to provide concise introductions to the research areas in which Simula specializes: scientifc computing, software engineering, communication systems, machine learning and cybersecurity. These books are written for graduate students, researchers, professionals and others who are keenly interested in the science of computing, and each volume presents a compact, state-of-the-art disciplinary overview and raises essential critical questions in the feld.

Simula's focus on computational physiology has grown considerably over the last decade, with the development of multi-scale mathematical models of excitable tissues (brain and heart) that are becoming increasingly complex and accurate. This sub-series represents a new branch of the SimulaSpringer Briefs that is specifcally focused on computational physiology. Each volume in this series will explore multiple physiological questions and the models developed to address them. Each of the questions will, in turn, be packaged into a short report format that provides a succinct summary of the fndings and, whenever possible, the software used will be made publicly available.

More information about this Subseries at http://www.springer.com/series/16669

Aslak Tveito · Kent-Andre Mardal · Marie E. Rognes Editors

# Modeling Excitable Tissue

The EMI Framework

*Editors* Aslak Tveito Simula Research Laboratory Fornebu, Norway

Marie E. Rognes Simula Research Laboratory Fornebu, Norway

Kent-Andre Mardal Department of Mathematics University of Oslo Oslo, Norway

ISSN: 2512-1677 ISSN: 2512-1685 (electronic) Simula SpringerBriefs on Computing ISSN 2730-7735 ISSN 2730-7743 (electronic) Reports on Computational Physiology ISBN 978-3-030-61156-9 ISBN 978-3-030-61157-6 (eBook) https://doi.org/10.1007/978-3-030-61157-6

Mathematics Subject Classifcation: 65M60, 65M06, 92C10, 92C37

© The Editor(s) (if applicable) and the author(s) 2021. 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 specifc 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 affliations.

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

# **Foreword to** *Reports on Computational Physiology*

## Dear Reader,

In 2016, Springer and Simula launched an Open Access series called the Simula SpringerBriefs on Computing. This series aims to provide concise introductions to the research areas in which Simula specializes: scientific computing, software engineering, communication systems, machine learning and cybersecurity. These books are written for graduate students, researchers, professionals and others who are keenly interested in the science of computing. We know that entering a new field of research or getting up-to-date on a new topic of interest can be very demanding and time consuming, for students and experienced researchers alike. Each volume presents a compact, state-of-the-art disciplinary overview and raises essential critical questions in the field, all in approximately 100 pages.

Simula's focus on computational physiology has grown considerably over the last decade. Our researchers collaborate with partners around the world in interdisciplinary teams to develop multi-scale mathematical models of excitable tissues (brain and heart). These models are becoming increasingly complex and accurate, in particular as they are compared to experimental and clinical data. Since 2014, the University of California, San Diego (UCSD) and Simula have organized an annual summer school in computational physiology, in which graduate students spend the second two weeks of June in Oslo learning the principles underlying mathematical models commonly used to study the heart and the brain. The students are then 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 findings. Each year, we have been impressed by the students' abilities to learn the huge amount of mathematics and physiology theory required for their research projects, the results of which often contain the rudiments of a scientific paper.

As a result of our expanding activities in this field, we have decided to publish a branch of the SimulaSpringer Briefs that is specifically focused on computational physiology. Each volume in this series will explore multiple physiological questions and the models developed to address them. Each of the questions will, in turn, be packaged into a short report format (6-10 pages) that provides a succinct summary of the findings and, whenever possible, the software used will be made publicly available. All reports in this series are subjected to peer-review. We would like to emphasise that we do not require that reports represent new scientific results; rather, they can reproduce or supplement earlier computational studies or experimental findings.

The main driver for this series is to enable the publication of project reports from the annual summer school in Computational Physiology, but we will also publish related reports that fit the overall purpose of the SimulaSpringer Briefs. Due to the Covid-19 pandemic, the summer school in 2020 had to be cancelled, and as such this first issue of the new series is not a collection of project reports. However, the topic of the first issue fits very well in the framework of computational physiology; it presents novel methods and software for the simulation of excitable cells.

It is a pleasure to thank our collaborators at SpringerNature for their superbly efficient handling of this manuscript. In particular, we are grateful for the sound advice and support from Dr. Martin Peters, the Executive Editor for Mathematics, Computational Science and Engineering. We would also like to thank Dr. Henrik Nicolay Finsberg for his excellent technical support in this project.

Fornebu, Norway *Dr. Kimberly J McCabe* September 2020 *Dr. Rachel Thomas Dr. Andrew D McCulloch Dr. Aslak Tveito*

# **Preface**

Partial differential equations (PDEs) have proved to be immensely useful in modelling Nature; virtually all fields of science have their own equations, and every field of engineering is based on mathematical models formulated in terms of PDEs. This is astonishing given the fact that no model, but the very simplest ones, can be studied using analytical (paper and pencil) techniques. Numerical computations have proved tremendously useful in order to understand models formulated in terms of PDEs, and it can be argued that the computer was invented for the purpose of solving such equations. The computer is extremely well suited to perform the huge amounts of tedious and highly repetitive computations that earlier had to be completed by humans. However, for a very long time, the computers typically available at research labs could solve only simple models. In the eighties, PDEs was almost always studied in 1 or 2 spatial dimensions; the 2D geometry was very simple and the model was most often linear and scalar. That level of computational complexity allowed analysis of qualitative properties of PDEs, but was insufficient for studying realistic models in Science and Engineering.

Over the past 30 years, we have witnessed a tremendous development in computing power both in terms of hardware, solution algorithms and software. This development has paved the way for realism in modeling; geometrical structures can now be represented with high degree of accuracy and complex systems of PDEs is applied to model the complex dynamics under consideration. This has led to greatly improved understanding throughout many branches of Science and had led to the development of considerably more accurate tools in Engineering.

Computational electrophysiology is a branch of Science that has benefited greatly from developments in computational analysis of PDEs. This development started 70 years ago with the celebrated paper by Hodgkin and Huxley (see (6)) and was followed 10 years later by the first cardiac action potential model developed by Noble (13). Since these groundbreaking papers, there have been intense efforts to understand how excitable cells works based on modeling and computations. This development has moved in tandem with new experimental techniques providing more accurate data necessary for parameterizing the increasingly complex models. A large number of membrane models (many hundreds) have been developed (see e.g. (11)), and these models have been used together with the monodomain or bidomain models (see e.g. (5; 21)) to study the electrochemical waves underpinning the contraction of the cardiac muscle. Similarly, the Cable equation (see e.g. (15)) has been extensively used to understand propagation of electrochemical signals by neurons.

Both the monodomain and the bidomain models of electrophysiology are derived based on homogenization of the cardiac tissue. In the resulting models, this means that the both the extracellular space, the cell membrane and the intracellular space are assumed to exist everywhere in the computational domain. Specifically, this means that the cardiac cell is not explicitly present in these models. This approach to modeling cardiac tissue enable analysis of phenomena on a relatively large length scale (mm), but is useless when it comes to study processes going on at a small length scale (μm). In 1993, the bidomain model was solved by Trayanova et al. ((20)) using 257 computational nodes, and at that time this was the best that technology would allow, and using a homogenized, large scale, model made perfect sense. More recently, however, models based on about 30 million mesh points are used allowing a characteristic mesh length of about 50μm which is about half the length of a human ventricular cell. Further refinement of the mesh used in the monodomain and bidomain models is not very useful since converged solutions are obtained at a quite coarse mesh (∼0.3mm, see e.g. (12; 3)). This means that technology now allows simulation at a shorter length scale than the classical models (monodomain and bidomain) are meant to represent; it is impossible to gain information at the μmlevel using these models. Specifically, the cell is impossible to represent explicitly in the classical models and that clearly limits their usefulness.

Interesting phenomena in electrophysiology take place close to the cell membrane. But since the cell is not explicitly present in the classical models, it is very difficult, if at all possible, to use such models to get a grip on what is going on in the vicinity of excitable cells. And since mesh resolution already has reached the μm-scale, it is clearly about time to introduce the cell as the building block in models of excitable tissue. In fact, this development started several years ago and has been pursued by many authors; see e.g. (7; 10; 14; 19; 16; 18; 17; 24; 1). Recently, we have followed up on these papers aiming at developing models, algorithms and software for cellbased representation of excitable tissue. We represent the extracellular (E) domain, the cell membrane (M) and the intracellular domain (I) explicitly, and therefore refer to it as the EMI-model.

In the first paper (see (23)) using the EMI model, we compared the results of the EMI model with the results of the classical Cable equations. In particular, we assessed the magnitude of the error introduced in the Cable equation by ignoring the ephaptic effects (i.e. assuming that the extracellular potential is constant). Then the same model was used to assess the effect on the action potential of placing a microelectrode array in the extracellular domain (see (2)). Furthermore, we have developed computational techniques for solving the EMI equations (see (22; 8)) and used it to study properties of the conduction velocity of electrochemical wave traversing the cardiac muscle during a heart-beat (see (9)). Quite recently, the EMI model has been extended to also account for ion concentrations in the entire computational domain using the electroneutral Kirchhoff-Nernst-Planck (KNP) model, and the resulting model is referred to as the KNP-EMI model (see (4)). This opens the possibility for analysing the effect of depression waves traversing cortical tissue.

The main of advantage of the EMI-models is the possibility for the modeler to represent local properties of the cell and to study dynamics in the vicinity of the cells at the μm–scale. This opens vast possibilities for deeper understanding of the dynamics of collections of excitable cells. However, a main disadvantage is that the model is more complex and more computationally demanding that the common monodomain, bidomain, and cable equations. The purpose of this edition of the Simula SpringerBriefs on Computing is to provide succinct introductions to various aspects of the EMI-models, the solution algorithms and the software used to study these models. These models are not straightforward to implement and we therefore think it is useful to provide software for anyone interested in using the models. Note that it is specifically not our intention here to provide substantial new contributions to developments of models, algorithms or software, but rather to aid readers by providing easy and readable accounts of this material.

In order to avoid misunderstanding, we would like to add that we do not think this is the end of the monodomain-, bidomain-, or the cable-model. These models have been extremely useful and in combinations with membrane models they basically represent our collective knowledge about how excitable cells work. Much work remains to be done with these equations and the models we suggest are far too computationally demanding to be a realistic alternative for full scale simulations of human organs. Also, again in order to avoid misunderstandings, homogenization is still with us in the EMI models; the EMI models takes the typical length scale form mm to μm, but the atomic scale is still 10000 smaller, so the models we use in both E, M and I all represent averages.

Fornebu, Norway *Aslak Tveito* September 2020 *Kent-Andre Mardal Marie Rognes*

# **References**


# **List of Contributors**

Aslak Tveito e-mail: aslak@simula.no Karoline Horgmo Jæger e-mail: karolihj@simula.no Åshild Telle e-mail: aashild@simula.no Samuel T. Wall e-mail: samwall@simula.no Joakim Sundnes e-mail: sundnes@simula.no Ada J. Ellingsrud e-mail: ada@simula.no Cécile Daversin-Catty e-mail: cecile@simula.no Marie E. Rognes e-mail: meg@simula.no Miroslav Kuchta e-mail: miroslav@simula.no Jakob Schreiner e-mail: jakob@simula.no Kent-André Mardal e-mail: kent\_and@simula.no Kristian G. Hustad e-mail: kghustad@simula.no Xing Cai e-mail: xingca@simula.no Henrik Nicolay Finsberg e-mail: henriknf@simula.no Simula Research Laboratory, Martin Linges 25, 1364 Fornebu, Norway Alessio Paolo Buccino e-mail: alessio.buccino@bsse.ethz.ch

Bio Engineering Laboratory, Department of Biosystems and Science Engineering, ETH Zurich, Switzerland

# **Contents**




# **Chapter 1**

# **Derivation of a Cell-Based Mathematical Model of Excitable Cells**

Karoline Horgmo Jæger<sup>1</sup> and Aslak Tveito1,<sup>2</sup>

**Abstract** Excitable cells are of vital importance in biology, and mathematical models have contributed significantly to understand their basic mechanisms. However, classical models of excitable cells are based on severe assumptions that may limit the accuracy of the simulation results. Here, we derive a more detailed approach to modeling that has recently been applied to study the electrical properties of both neurons and cardiomyocytes. The model is derived from first principles and opens up possibilities for studying detailed properties of excitable cells. We refer to the model as the EMI model because both the extracellular space (E), the cell membrane (M) and the intracellular space (I) are explicitly represented in the model, in contrast to classical spatial models of excitable cells. Later chapters of the present text will focus on numerical methods and software for solving the model. Also, in the next chapter, the model will be extended to account for ionic concentrations in the intracellular and extracellular spaces.

# **1.1 Introduction**

Mathematical modeling has a great potential for increasing our understanding of the physiological processes underlying the function of the body. For example, modeling of the electrical properties of excitable cells may provide insight into the complex electrical signaling involved in a number of important functions, like transfer of information through neurons and coordination of the pumping of the heart. A popular model of the conduction of electrical signals in neurons is the so-called cable equation (30), whereas the extracellular potential surrounding neurons is often modeled using the point-source or line-source approximations (6). The three aforementioned models

<sup>1</sup>Simula Research Laboratory, Norway

<sup>2</sup>Department of Informatics, University of Oslo, Norway

have been used extensively to gain insight into the function of neurons and the interpretation of measurements of the extracellular potential around neurons (6; 5; 12). Correspondingly, the conduction of electrical signals through the heart is traditionally modeled using the homogenized bidomain and monodomain models (20). These models are also widely used and have, for instance, provided insight into mechanisms of cardiac arrhythmias (26; 33).

However, despite the success of the above-mentioned classical models of excitable cells, the models have certain shortcomings that may make them inaccurate or impractical in some situations. For example, in the derivation of the cable equation, the extracellular potential is often assumed to be constant (30; 15). Therefore, changes in the extracellular potential generated by the neuron itself or by neighboring neurons (i.e., ephaptic effects) are ignored in the model. This could potentially lead to inaccuracies (16; 2; 34). In addition, the point-source and line-source approximations rely on the assumption that the extracellular space is infinite and homogeneous. Consequently, the models might not be well-suited to interpret extracellular potentials measured when large measurement electrodes are present in the extracellular space close to the neurons (3).

Moreover, the bidomain and monodomain models represent cardiac tissue in a homogenized manner, assuming that the intracellular space, the extracellular space and the cell membrane exist everywhere in the tissue. Because the geometry of the individual cells is not represented, it is very hard to use the models to study the effect of, e.g., the cell geometry or a non-uniform distribution of ion channels on the cell membrane. These properties are both believed to influence cardiac conduction, but their exact effects are not fully understood and call for further investigations (28; 21). In addition, it has been proposed that ephaptic coupling between cardiac cells might occur at small extracellular clefts located at the intercalated discs between cells (29). Since the geometry of the extracellular space is not represented in the homogenized models, it is difficult to use these models to study such ephaptic effects.

In order to account for the difficulties related to the classical models, alternative electrophysiological models have been developed (e.g., (28; 31; 24)). In this chapter, we consider one of these alternative models, referred to as the EMI model, because it explicitly represents the extracellular space (E), the cell membrane (M) and the intracellular space (I). This model has been used to study both neurons (1; 34; 3) and cardiac tissue (32; 31; 18). Because the model represents the extracellular space, the membrane and the intracellular space in a coupled manner, the model allows for representation of ephaptic effects between neurons (34) or cardiomyocytes (18). In addition, since the geometry of the extracellular space is explicitly represented, the model allows for representation of non-homogeneous extracellular surroundings (3). Furthermore, since the geometry of each cell is represented, the model allows for study of the effect of cell geometry and non-uniform distributions of ion channels on cardiac conduction properties (18).

In other words, the EMI model allows for a more detailed representation of excitable cells and tissues than classical models of computational electrophysiology. In fact,

Fig. 1.1: Illustration of an EMI model domain consisting of an extracellular domain, Ωe, a cell membrane, Γ, and an intracellular domain, Ωi.

the classical models mentioned above can be derived from the more detailed EMI model by introducing certain simplifying assumptions, see e.g., (1; 8; 11). In this chapter, however, we focus on deriving the EMI model from Maxwell's equations of electromagnetism.

# **1.2 Derivation of the EMI Model**

In this section, we present a derivation of the EMI model for excitable cells. This derivation is to a large extent based on the derivation found in (1; 17). We consider a domain separated into an extracellular part, Ω<sup>e</sup> and an intracellular part, Ωi, like illustrated in Figures 1.1 and 1.3. The cell membrane, denoted by Γ, is defined as the boundary between Ω<sup>i</sup> and Ωe. Here, we derive a model for the electrical potentials in both a domain with a single cell (Figure 1.1) and in a domain with two cells connected at an intercalated disc denoted by Γ1,<sup>2</sup> (Figure 1.3).

## **1.2.1 Fundamental Equations**

We base the derivation of the EMI model on two of the quasi-static approximations of Maxwell's equations, i.e.,

$$
\nabla \times \mathbf{E} = \mathbf{0},
\tag{1.1}
$$

$$
\nabla \times \mathbf{H} = \mathbf{J}.\tag{1.2}
$$

Here, **E** is the electric field (typically in μF/cm), **H** is the magnetic field (typically in μA/cm) and **J** is the density of free currents (typically in μA/cm2). In the quasi-static approximation of (1.2), it is assumed that free unbalanced charges are instantly balanced. The assumptions hold in the intracellular and extracellular spaces. However, we assume that charges may accumulate at the cell membrane and at the intercalated discs between cells. Therefore, we let (1.2) at these locations be replaced by the corresponding equation without the quasi-static approximation, i.e., by

$$
\nabla \times \mathbf{H} = \mathbf{J} + \varepsilon \frac{\partial \mathbf{E}}{\partial t},
\tag{1.3}
$$

where ε is the permittivity of the medium (typically in μF/cm). In addition, we assume that Ohm's law applies in the intracellular and extracellular domains. This means that

$$
\mathbf{J} = \sigma \mathbf{E},
\tag{1.4}
$$

where σ is the conductivity of the considered medium (typically in mS/cm). We also note that (1.1) implies that **E** is a conservative vector field and that it therefore can be defined as the gradient of a scalar field (9). More specifically, we can define

$$\mathbf{E} = -\nabla \boldsymbol{\mu},\tag{1.5}$$

where the scalar *u* is the electric potential (typically in mV).

## **1.2.2 Model for the Intracellular and Extracellular Domains**

In order to derive equations for the intracellular and extracellular domains, we take the divergence of both sides of (1.2) and apply the vector identity ∇ · (∇ × **H**) = 0, which holds for any vector **H** (9). This yields

$$\nabla \cdot \mathbf{J} = 0.$$

Inserting (1.4) and (1.5), we obtain the Laplace equation

$$
\nabla \cdot \sigma \nabla \mu = 0.\tag{1.6}
$$

More specifically, for the intracellular and extracellular domains, we have

$$
\nabla \cdot \sigma\_i \nabla u\_i = 0 \quad \text{in } \Omega\_i,\tag{1.7}
$$

$$
\nabla \cdot \sigma\_e \nabla \mu\_e = 0 \quad \text{in } \Omega\_e,\tag{1.8}
$$

where σ<sup>i</sup> and σ<sup>e</sup> are the intracellular and extracellular conductivities, respectively, and *u*<sup>i</sup> and *u*<sup>e</sup> are the electric potentials in the intracellular and extracellular domains, respectively.

## **1.2.3 Model for the Membrane**

In order to derive the EMI model equations for the membrane, we consider a volume element, *B*, intersected by the membrane. This volume element may be divided into an extracellular part, *B*e, and an intracellular part, *B*i, such that *B*<sup>e</sup> ∪ *B*<sup>i</sup> = *B* and *B*<sup>e</sup> ∩ *B*<sup>i</sup> = ∅, as illustrated in Figure 1.2A. In each of these domains, we assume that (1.3) holds. Taking the divergence of both sides of (1.3) and applying the vector identity ∇ · (∇ × **H**) = 0 results in

$$
\nabla \cdot \mathbf{J} = -\nabla \cdot \varepsilon \frac{\partial \mathbf{E}}{\partial t}.\tag{1.9}
$$

Integrating this equation over each of the volume elements *B*<sup>i</sup> and *B*e, we get

$$\begin{aligned} \int\_{B\_i} \nabla \cdot \mathbf{J} \, dV &= - \int\_{B\_i} \nabla \cdot \boldsymbol{\varepsilon} \frac{\partial \mathbf{E}}{\partial t} \, dV, \\ \int\_{B\_\varepsilon} \nabla \cdot \mathbf{J} \, dV &= - \int\_{B\_\varepsilon} \nabla \cdot \boldsymbol{\varepsilon} \frac{\partial \mathbf{E}}{\partial t} \, dV, \end{aligned}$$

and applying the divergence theorem (see e.g., (9)), we obtain

$$\int\_{\partial B\_i} \mathbf{J} \cdot \mathbf{n}\_{B\_i} \, dS = -\int\_{\partial B\_i} \varepsilon \frac{\partial \mathbf{E}}{\partial t} \cdot \mathbf{n}\_{B\_i} \, dS,\tag{1.10}$$

$$\int\_{\partial B\_{\varepsilon}} \mathbf{J} \cdot \mathbf{n}\_{B\_{\varepsilon}} \, dS = -\int\_{\partial B\_{\varepsilon}} \varepsilon \frac{\partial \mathbf{E}}{\partial t} \cdot \mathbf{n}\_{B\_{\varepsilon}} \, dS. \tag{1.11}$$

Here, **<sup>n</sup>**Bi and **<sup>n</sup>**Be are the outward pointing normal vectors of *<sup>B</sup>*<sup>i</sup> and *<sup>B</sup>*e, respectively. Furthermore, in (1.10) and (1.11), the left-hand side terms represent the free ionic current and the right-hand side terms represent the capacitive current.

#### **1.2.3.1 Ionic Current**

We start by considering the left-hand side of (1.10), representing the ionic current. Here, we note that the boundary ∂*B*<sup>i</sup> may be split into two parts, Γ<sup>B</sup> and ∂*B*<sup>i</sup> \ ΓB, where Γ<sup>B</sup> is the part of ∂*B*<sup>i</sup> coinciding with the membrane and ∂*B*<sup>i</sup> \ Γ<sup>B</sup> is the remaining part (see Figure 1.2A). We can then write

$$\int\_{\partial B\_i} \mathbf{J} \cdot \mathbf{n}\_{B\_i} \, dS = \int\_{\partial B\_i \backslash \Gamma\_B} \mathbf{J} \cdot \mathbf{n}\_{B\_i} \, dS + \int\_{\Gamma\_B} \mathbf{J} \cdot \mathbf{n}\_i \, dS,\tag{1.12}$$

where **<sup>n</sup>**<sup>i</sup> is the outward pointing normal vector of the membrane and **<sup>n</sup>**Bi is the outward pointing normal vector of the remaining part of *B*i, as illustrated in Figure 1.2A. At the membrane, the current density, **J**, consists of currents through different types of ion channels, pumps and exchangers located at the membrane. This current den-

Fig. 1.2: **A**: Illustration of a volume element, *B*, intersected by the membrane. The volume element is separated into an intracellular part, *B*i, and an extracellular part, *B*e. **B**: Illustration of a small volume element, *B*, located on the extracellular part of the membrane.

sity is typically denoted by *I*ion and given in units of μA/cm2. By convention, *I*ion is defined to be positive for a flux of positive ions out of the cell (i.e., in the direction of **n**i). This gives

$$\int\_{\Gamma\_B} \mathbf{J} \cdot \mathbf{n}\_i \, dS = \int\_{\Gamma\_B} I\_{\text{ion}} \, dS. \tag{1.13}$$

The boundary ∂*B*<sup>i</sup> \ Γ<sup>B</sup> is located in the intracellular domain. Here, we assume that the current density, **J**, is given by Ohm's law (1.4), such that

$$\int\_{\partial B\_i \backslash \Gamma\_B} \mathbf{J} \cdot \mathbf{n}\_{B\_i} \, dS = \int\_{\partial B\_i \backslash \Gamma\_B} \sigma\_i \mathbf{E} \cdot \mathbf{n}\_{B\_i} \, dS. \tag{1.14}$$

Inserting (1.13) and (1.14) into (1.12), we get

$$\int\_{\partial B\_i} \mathbf{J} \cdot \mathbf{n}\_{B\_i} \, dS = \int\_{\partial B\_i \backslash \Gamma\_B} \sigma\_i \mathbf{E} \cdot \mathbf{n}\_{B\_i} \, dS + \int\_{\Gamma\_B} I\_{\text{ion}} \, dS,\tag{1.15}$$

and similar arguments for the extracellular part of the membrane yield

$$\int\_{\partial B\_{\varepsilon}} \mathbf{J} \cdot \mathbf{n}\_{B\_{\varepsilon}} \, dS = \int\_{\partial B\_{\varepsilon} \backslash \Gamma\_{B}} \sigma\_{\varepsilon} \mathbf{E} \cdot \mathbf{n}\_{B\_{\varepsilon}} \, dS - \int\_{\Gamma\_{B}} I\_{\text{ion}} \, dS. \tag{1.16}$$

Note that the negative sign in front of the last term is due to the fact that **n**<sup>e</sup> = −**n**i.

#### **1.2.3.2 Capacitive Current**

For the right-hand side part of (1.10), representing the capacitive current, we again split the integral into two parts

$$\int\_{\partial B\_i} \boldsymbol{\varepsilon} \frac{\partial \mathbf{E}}{\partial t} \cdot \mathbf{n}\_{B\_i} \, dS = \int\_{\partial B\_i \backslash \Gamma\_B} \boldsymbol{\varepsilon} \frac{\partial \mathbf{E}}{\partial t} \cdot \mathbf{n}\_{B\_i} \, dS + \int\_{\Gamma\_B} \boldsymbol{\varepsilon}\_\Gamma \frac{\partial \mathbf{E}}{\partial t} \cdot \mathbf{n}\_i \, dS.$$

Here, ε<sup>Γ</sup> is the permittivity of the membrane. Following the quasi-static assumptions, we assume that the term ε <sup>∂</sup>**<sup>E</sup>** <sup>∂</sup><sup>t</sup> is negligible for the part of <sup>∂</sup>*B*<sup>i</sup> that does not coincide with the membrane. Furthermore, from (1.5), we get **<sup>E</sup>** · **<sup>n</sup>**<sup>i</sup> <sup>=</sup> −∇*<sup>u</sup>* · **<sup>n</sup>**<sup>i</sup> <sup>≈</sup> <sup>v</sup>/*d*, where

$$
\nu = \mu\_i - \mu\_e \tag{1.17}
$$

is the membrane potential and *d* is the thickness of the membrane. We assume that the membrane can be treated as a capacitor formed by two parallel plates separated by an insulator. In that case, the membrane capacitance per area is given by *C*<sup>m</sup> = εΓ/*d* (13). Therefore,

$$\int\_{\partial B\_i} \boldsymbol{\varepsilon} \frac{\partial \mathbf{E}}{\partial t} \cdot \mathbf{n}\_{B\_i} \, dS = \int\_{\Gamma\_B} \boldsymbol{\varepsilon}\_\Gamma \frac{\partial \mathbf{E}}{\partial t} \cdot \mathbf{n}\_i \, dS = \int\_{\Gamma\_B} \frac{\boldsymbol{\varepsilon}\_\Gamma}{d} \frac{\partial \mathbf{v}}{\partial t} \, dS = \int\_{\Gamma\_B} \mathbf{C}\_m \frac{\partial \mathbf{v}}{\partial t} \, dS. \tag{1.18}$$

Similar arguments for the extracellular side yield

$$\int\_{\partial B\_{\rm ef}} \mathcal{E} \frac{\partial \mathbf{E}}{\partial t} \cdot \mathbf{n}\_{B\_{\rm ef}} \, dS = - \int\_{\Gamma\_B} \mathbf{C}\_m \frac{\partial \mathbf{v}}{\partial t} \, dS,\tag{1.19}$$

where the change of sign again is due to the fact that **n**<sup>e</sup> = −**n**i.

#### **1.2.3.3 Collecting the Ionic and Capacitive Currents**

Collecting the ionic and capacitive currents by inserting (1.15)–(1.16) and (1.18)– (1.19) into (1.10)–(1.11), we obtain

$$\begin{aligned} \int\_{\partial B\_i \backslash \Gamma\_B} \sigma\_i \mathbf{E} \cdot \mathbf{n}\_{B\_i} \, dS + \int\_{\Gamma\_B} I\_{\text{ion}} \, dS &= - \int\_{\Gamma\_B} C\_m \frac{\partial \nu}{\partial t} \, dS, \\ \int\_{\partial B\_e \backslash \Gamma\_B} \sigma\_e \mathbf{E} \cdot \mathbf{n}\_{B\_e} \, dS - \int\_{\Gamma\_B} I\_{\text{ion}} \, dS &= \int\_{\Gamma\_B} C\_m \frac{\partial \nu}{\partial t} \, dS, \end{aligned}$$

which can be rewritten to

$$\int\_{\partial B\_i \backslash \Gamma\_B} \sigma\_i \mathbf{E} \cdot \mathbf{n}\_{B\_i} \, dS = - \int\_{\Gamma\_B} I\_m \, dS,\tag{1.20}$$

$$\int\_{\partial B\_{\mathbf{e}} \backslash \Gamma\_B} \sigma\_{\mathbf{e}} \mathbf{E} \cdot \mathbf{n}\_{B\_{\mathbf{e}}} \, dS = \int\_{\Gamma\_B} I\_m \, dS,\tag{1.21}$$

where the total membrane current density *I*<sup>m</sup> is defined as

$$I\_m = C\_m \frac{\partial \nu}{\partial t} + I\_{\text{ion}}.\tag{1.22}$$

Fig. 1.3: Illustration of an EMI model domain consisting of two cells, Ω<sup>1</sup> <sup>i</sup> and <sup>Ω</sup><sup>2</sup> i , connected at an intercalated disc, Γ1,<sup>2</sup> and surrounded by an extracellular domain, Ω<sup>e</sup>

We now wish to rewrite (1.20)–(1.21) to a differential form. We note that we can divide any volume element, *B*, intersecting the membrane into a purely intracellular, a purely extracellular, and a membrane intersecting part. We also know that (1.7)–(1.8) hold in the purely intracellular and extracellular parts. Therefore, we are interested in equations (1.20)–(1.21) as the size of *B* approaches zero. For example, we may consider a small extracellular volume element shaped as a cylinder, as illustrated in Figure 1.2B. As the height, Δ*h*B, of this cylinder approaches zero, the integral over ∂*B*<sup>e</sup> \ Γ<sup>B</sup> approaches the integral over ΓB, and we therefore get

$$\int\_{\partial B\_{\mathfrak{e}} \backslash \Gamma\_B} \sigma\_{\mathfrak{e}} \mathbf{E} \cdot \mathbf{n}\_{B\_{\mathfrak{e}}} \, dS \approx \int\_{\Gamma\_B} \sigma\_{\mathfrak{e}} \mathbf{E} \cdot \mathbf{n}\_{B\_{\mathfrak{e}}} \, dS.$$

Inserting this approximation into (1.21), we obtain

$$\int\_{\Gamma\_B} \sigma\_e \mathbf{E} \cdot \mathbf{n}\_{B\_e} \, dS = \int\_{\Gamma\_B} I\_m \, dS \qquad \Rightarrow \qquad \sigma\_e \mathbf{E} \cdot \mathbf{n}\_{B\_e} = I\_m,$$

and inserting (1.5) and **<sup>n</sup>**<sup>e</sup> <sup>=</sup> <sup>−</sup>**n**Be , we get

$$
\sigma\_e \nabla \mu\_e \cdot \mathbf{n}\_e = I\_m. \tag{1.23}
$$

Similar arguments for the intracellular part of the membrane yield

$$-\sigma\_i \nabla u\_i \cdot \mathbf{n}\_i = I\_m,\tag{1.24}$$

where the negative sign is due to the negative sign in (1.20). Finally, combining (1.23) and (1.24), we obtain

$$
\sigma\_{\text{e}} \nabla u\_{\text{e}} \cdot \mathbf{n}\_{\text{e}} = -\sigma\_{i} \nabla u\_{i} \cdot \mathbf{n}\_{i} = I\_{m}, \tag{1.25}
$$

where *I*<sup>m</sup> = *C*<sup>m</sup> <sup>∂</sup><sup>v</sup> <sup>∂</sup><sup>t</sup> <sup>+</sup> *<sup>I</sup>*ion and <sup>v</sup> <sup>=</sup> *<sup>u</sup>*<sup>i</sup> <sup>−</sup> *<sup>u</sup>*<sup>e</sup> (see (1.22) and (1.17)).

## **1.2.4 Model for the Intercalated Disc**

In some cases, we wish to model cells that are connected to each other, as illustrated in Figure 1.3. We then let the intercalated discs connecting the cells be represented as boundaries between the intracellular domains, like the membrane is a boundary between the intracellular and extracellular domains. Furthermore, we assume that the intercalated disc have capacitive properties like the membrane, and that gap junctions allow for currents between neighboring cells, in the same manner as ion channels allows for currents between the intracellular and extracellular spaces. Therefore, the derivation of equations for an intercalated disc follows the exact same lines as the derivation of the membrane equations. More precisely, for two connected cells, we define an intercalated disc potential, w, by

$$
\hbar \boldsymbol{\omega} = \boldsymbol{\mu}\_i^1 - \boldsymbol{\mu}\_i^2,\tag{1.26}
$$

where *u*<sup>1</sup> <sup>i</sup> and *<sup>u</sup>*<sup>2</sup> <sup>i</sup> are the electric potentials in <sup>Ω</sup><sup>1</sup> <sup>i</sup> and <sup>Ω</sup><sup>2</sup> <sup>i</sup> , respectively. In addition, we define a total intercalated disc current density, *I*1,2, by

$$I\_{1,2} = C\_{1,2} \frac{\partial \boldsymbol{w}}{\partial t} + I\_{\text{gap}},\tag{1.27}$$

where *I*gap is the current density through the gap junctions, with positive direction in the direction from Ω<sup>1</sup> <sup>i</sup> to <sup>Ω</sup><sup>2</sup> <sup>i</sup> , *<sup>C</sup>*1,<sup>2</sup> is the capacitance of the intercalated disc, and *C*1,<sup>2</sup> ∂w <sup>∂</sup><sup>t</sup> is the capacitive current density of the intercalated disc. Furthermore, following the same arguments as for the derivation of the membrane equations, we end up with an analogue to (1.25) of the form

$$
\sigma\_i^2 \nabla \mu\_i^2 \cdot \mathbf{n}\_i^2 = -\sigma\_i^1 \nabla \mu\_i^1 \cdot \mathbf{n}\_i^1 = I\_{1,2}, \tag{1.28}
$$

representing the total current density across the interface.

## **1.2.5 Models of the Ionic Currents**

Mathematical models of the ionic currents governing the membrane potential of excitable cells come in a large variety of versions; see (4) for several hundred examples. The simplest possible model is just a passive current of the form *I*ion = const · v, followed by a third order polynomial model. More realistic models tend to be more complex and are usually written on the form

$$I\_{\rm ion} = \sum\_{i=1}^{N} I\_i,\tag{1.29}$$

given in μA/cm2. Here, the individual currents can usually be written on the form *<sup>I</sup>*<sup>i</sup> <sup>=</sup> *<sup>I</sup>*i(v,*s*), where <sup>v</sup> denotes the membrane potential, given by *<sup>u</sup>*<sup>i</sup> <sup>−</sup> *<sup>u</sup>*e, and *<sup>s</sup>* denotes gating variables and ionic concentrations. The celebrated model of the action potential of a neuron presented by Hodgkin and Huxley (see (14)) can be written on this form, and so can the first model of a cardiac cell presented by Nobel (25). A comprehensive and readable introduction to models of the membrane ionic currents is given in the survey (27).

Correspondingly, the ionic currents through gap junctions between neighboring cells are often modeled by a simple passive model of the form *<sup>I</sup>*gap <sup>=</sup> const · <sup>w</sup>. More detailed models of voltage-dependent gap junction dynamics have also been introduced (see e.g., (10; 35)).

## **1.2.6 Summary of the Model Equations**

In summary, the EMI model for a single cell surrounded by an extracellular domain (as illustrated in Figure 1.1) is given by the equations (1.7), (1.8), (1.17), (1.22) and (1.25), that is

$$\nabla \cdot \sigma\_i \nabla u\_i = 0 \qquad\qquad\qquad\text{in } \Omega\_i,\tag{1.30}$$

$$\nabla \cdot \sigma\_{\epsilon} \nabla u\_{\epsilon} = 0 \qquad\qquad\qquad\text{in } \Omega\_{\epsilon},\tag{1.31}$$

$$
\sigma\_e \nabla \mu\_e \cdot \mathbf{n}\_e = -\sigma\_i \nabla \mu\_i \cdot \mathbf{n}\_i \equiv I\_m \quad \text{at } \Gamma,\tag{1.32}
$$

$$
\nu = \mu\_i - \mu\_e \qquad\qquad\qquad\text{at }\Gamma,\tag{1.33}
$$

$$\frac{\partial \mathbf{v}}{\partial t} = \frac{1}{C\_m} (I\_m - I\_{\rm ion}) \qquad \text{ at } \Gamma, \tag{1.34}$$

where *<sup>u</sup>*i, *<sup>u</sup>*<sup>e</sup> and <sup>v</sup> are the intracellular, extracellular and membrane potentials, respectively, typically given in mV. Moreover, σ<sup>i</sup> and σ<sup>e</sup> are the intracellular and extracellular conductivities, respectively (typically in mS/cm), *C*<sup>m</sup> is the membrane capacitance (typically in μF/cm2), and Γ denotes the cell membrane. The ionic currents through channels, pumps and exchangers at the membrane are denoted by *I*ion and typically given in μA/cm2.

If several cells are connected at intercalated discs, as illustrated for two cells in Figure 1.3, the system of equations must be extended to include equations for the currents between cells. For two cells, this extension consists of the equations

$$\left[\sigma\_{i}\nabla u\_{i}^{2}\cdot\mathbf{n}\_{i}^{2}=-\sigma\_{i}\nabla u\_{i}^{1}\cdot\mathbf{n}\_{i}^{1}\equiv I\_{1,2}\quad\text{at }\Gamma\_{1,2},\tag{1.35}$$

$$
\mu\_i^1 - \mu\_i^2 = w \tag{1.36}
$$

$$\mathcal{W}\_{\rm I} = \frac{1}{C\_{1,2}} (I\_{1,2} - I\_{\rm gap}) \qquad \text{at } \Gamma\_{1,2}, \tag{1.37}$$

where, as above, Γ1,<sup>2</sup> is the intercalated disc, **n**<sup>1</sup> <sup>i</sup> is the outward pointing normal vector of Ω<sup>1</sup> <sup>i</sup> , **<sup>n</sup>**<sup>2</sup> <sup>i</sup> is the outward pointing normal vector of <sup>Ω</sup><sup>2</sup> <sup>i</sup> , and *<sup>u</sup>*<sup>1</sup> <sup>i</sup> and *<sup>u</sup>*<sup>2</sup> <sup>i</sup> are the intracellular potentials (typically in mV) of Ω<sup>1</sup> <sup>i</sup> and <sup>Ω</sup><sup>2</sup> <sup>i</sup> , respectively. Furthermore, *C*1,<sup>2</sup> is the specific capacitance of the intercalated disc (typically in μF/cm2), and *I*gap is the current through the gap junctions (typically in μA/cm2).

# **1.3 Conclusion**

In the present chapter, we have derived the EMI model. The EMI model predicts electrical potentials in cells with an explicit geometrical representation and thus allows for more detail than homogenized models of excitable tissue. In the next chapter (7, Chapter 2), the model will be extended by taking ion concentration in the extracellular and intracellular spaces into account. Numerical solutions of the EMI models will be presented in (19, Chapter 4), (23, Chapter 5) and (22, Chapter 6). In these chapters the readers will also be pointed to open software that can be used to solve the EMI model.

**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.

# **References**


#### 1 Derivation of EMI 13


# **Chapter 2**

# **A Cell-Based Model for Ionic Electrodiffusion in Excitable Tissue**

Ada J. Ellingsrud1, Cécile Daversin-Catty1 and Marie E. Rognes<sup>1</sup>

**Abstract** This chapter presents the KNP-EMI model describing ion concentrations and electrodiffusion in excitable tissue. The KNP-EMI model extends on the EMI model by removing the assumption that ion concentrations are constant in time and space, and may as such be more appropriate in connection with modelling e.g. spreading depression, stroke and epilepsy. The KNP-EMI model defines a system of time-dependent, nonlinear, mixed dimensional partial differential equations. We here detail the derivation of the system and present a numerical example illustrating how ion concentrations evolve during neuronal activity.

# **2.1 Introduction and Motivation**

In this chapter, we present an extension of the EMI model, presented in (11, Chapter 1), describing ion concentrations and electrodiffusion in excitable tissue. The EMI model is based on the assumption that intra- and extracellular ion concentrations are constant in time and space. This is often a good approximation, as ion concentrations in healthy tissue typically quickly return to base levels after neuronal activity due to cellular mechanisms such as e.g. membrane pumps and glial cell buffering. However, there are scenarios where this assumption is inadequate.

Several cerebral pathologies are associated with increased neuronal activity (3), such as e.g. seizures and epilepsy (10; 6; 1), stroke (17), and spreading depression (22). In particular, periods of neuronal hyperactivity can lead to substantial variations in extracellular ion concentrations. These variations will in turn (i) influence membrane reversal potentials and (ii) generate diffusive currents. Changes in the reversal potentials, caused by local ionic shifts, may affect the dynamical properties of the

<sup>1</sup>Simula Research Laboratory, Norway

neurons (12; 16; 24). On the other hand, diffusive currents, driven by ion concentration gradients, can shift the extracellular potential (8; 3). Mathematical models addressing the aforementioned phenomena and pathologies should therefore also account for ion concentrations, their spatial and temporal gradients and associated dynamics.

In this chapter, we derive a system of time-dependent, nonlinear partial differential equations describing the distribution and evolution of ion concentrations in a geometrically-explicit representation of the intra- and extracellular domains using the electroneutral Kirchhoff-Nernst-Planck (KNP) model (21). We will refer to this model as the KNP-EMI model, see also e.g. (5).

# **2.2 Derivation of the Equations**

Let the computational domain Ω and subdomains Ωi, Ωe, and Γ be defined as in the previous chapter 1.1. For simplicity and clarity, we present the mathematical model for one intracellular region Ωi<sup>1</sup> = Ω<sup>i</sup> with membrane Γ below. We model a set *K* of intracellular and extracellular ion concentrations, and note that key ions in excitable tissue are potassium (K+), sodium (Na+), and chloride (Cl−). For each ion species *k* ∈ *K* and each region *r* ∈ {*i*, *e*}, we model the *ion concentrations c*k <sup>r</sup> : Ω<sup>r</sup> × (0,*T*] → R (mol/m3), and *electrical potentials u*<sup>r</sup> : Ω<sup>r</sup> × (0,*T*] → R (V), and additionally the *total transmembrane current density I*<sup>m</sup> : Γ × (0,*T*] → R (A/m2).

## **2.2.1 Equations in the Intracellular and Extracellular Volumes**

In the EMI model, the free current densities **J**i, **J**<sup>e</sup> (μA/cm2), c.f. (1.4), are assumed to satisfy Ohm's law. To include diffusive ion effects, we instead assume that the free current density is composed of flux density contributions **J**<sup>k</sup> <sup>r</sup> (mol/(m2s)) from different ions *k* as:

$$\mathbf{J}\_r = \sum\_{k \in K} F z^k \mathbf{J}\_r^k \quad \text{in } \Omega\_r,\tag{2.1}$$

where *z*<sup>k</sup> is the valence of ion species *k* and *F* (C/mol) is Faraday's constant. Furthermore, we assume that ions can move by diffusion and/or in response to the electrical field as charged particles. Hence, the ion flux densities are modelled as the sum of two terms: (i) the ion concentrations that are transported via electrical potential gradients σ<sup>k</sup> <sup>r</sup> ∇*u*<sup>r</sup> and (ii) the diffusive movement of ions due to ionic gradients *D*<sup>k</sup> <sup>r</sup> ∇ *c*<sup>k</sup> r :

$$\mathbf{J}\_r^k = -\sigma\_r^k \nabla \mu\_r - D\_r^k \nabla c\_r^k \quad \text{in } \Omega\_r,\tag{2.2}$$

where *D*<sup>k</sup> <sup>r</sup> (m2/s) and σ<sup>k</sup> <sup>r</sup> denote the effective diffusion coefficient and the conductivity for ion species *k* in region *r*, respectively. The conductivity σ<sup>k</sup> <sup>r</sup> depends on the concentration of ion species *k* and the diffusion coefficient *D*<sup>k</sup> <sup>r</sup> in the following manner:

$$
\sigma\_r^k = \sigma\_r^k(c\_r^k) = \frac{D\_r^k \boldsymbol{\mathcal{z}}^k}{\boldsymbol{\psi}} c\_r^k \quad \text{in } \Omega\_r. \tag{2.3}
$$

Here, the constant ψ = *RT F*−<sup>1</sup> combines Faraday's constant *F*, the absolute temperature *T* (K), and the gas constant *R* (J/(K mol)). Moreover, the bulk conductivity σ<sup>r</sup> can be expressed as:

$$
\sigma\_r = \sigma\_r(c\_r^k) = \frac{F}{\psi} \sum\_{k \in K} D\_r^k c\_r^k (z^k)^2 \quad \text{in } \Omega\_r. \tag{2.4}
$$

See e.g. (21) for a derivation of the conductivity (2.3) and the bulk conductivity (2.4). Comparing with (1.4) and (1.5), we note the dependency on the ion concentrations in the conductivity σ<sup>r</sup> in (2.3), and the second term accounting for ion diffusion in (2.2).

As in Chapter 1, we stipulate that:

$$\nabla \cdot \mathbf{J}\_i = 0 \quad \text{in } \Omega\_i,\tag{2.5}$$

$$\nabla \cdot \mathbf{J}\_{\epsilon} = 0 \quad \text{in } \Omega\_{\epsilon}. \tag{2.6}$$

Finally, conservation of ions for the bulk of each region Ω<sup>r</sup> gives that:

$$\frac{\partial [k]\_i}{\partial t} + \nabla \cdot \mathbf{J}\_i^k = 0 \qquad \text{in } \Omega\_i,\tag{2.7}$$

$$\frac{\partial [k]\_e}{\partial t} + \nabla \cdot \mathbf{J}\_e^k = 0 \qquad \text{in } \Omega\_e,\tag{2.8}$$

for *t* ∈ (0,*T*].

## **2.2.2 Membrane Currents**

We next turn to modelling the cell membrane currents and membrane potential across the interface <sup>Γ</sup>. As in Chapter 1, we introduce the membrane potential v as the jump in the electrical potential over the membrane:

$$
\nu = \mu\_i - \mu\_e \qquad \text{on } \Gamma. \tag{2.9}
$$

We also introduce the total membrane current as the combination of a capacitive current and ion specific currents:

2 A Cell-Based Model for Ionic Electrodiffusion in Excitable Tissue 17

$$I\_m = I\_{\rm cap} + I\_{\rm ion} = C\_m \frac{\partial \nu}{\partial t} + I\_{\rm ion},\tag{2.10}$$

where the total channel current *I*ion is the sum of the ion specific channel currents *I*k ion:

$$I\_{\rm ion} = \sum\_{k \in K} I\_{\rm ion}^k, \quad I\_{\rm ion}^k = I\_{\rm ion}^k(\nu, c\_\cdot^k, \dots). \tag{2.11}$$

The channel currents *I*<sup>k</sup> ion are subject to modelling, and will be discussed briefly in Section 2.2.2.1.

Using our concepts, we have that the *total ionic current density I*<sup>m</sup> : Γ × (0,*T*) → R (A/m2) across the interface Γ (from the intracellular to the extracellular domain) is given by:

$$-F\sum\_{k\in K} z^k \mathbf{J}\_e^k \cdot \mathbf{n}\_e = F \sum\_{k\in K} z^k \mathbf{J}\_i^k \cdot \mathbf{n}\_i \equiv I\_m. \tag{2.12}$$

It now remains to specify a set of interface conditions for the specific ion fluxes **J**k <sup>r</sup> · **n**<sup>r</sup> for *r* ∈ {*i*, *e*}.

Here, we propose a heuristic approach via ion specific capacitive current modelling, and note that an alternative approach is presented in (15). As for the total current, we assume that the capacitive current can be represented as a sum of ion specific contributions:

$$I\_{\rm cap} = \sum\_{k \in K} I\_{\rm cap}^k. \tag{2.13}$$

Without loss of generality, we let the ion specific capacitive current *I*<sup>k</sup> cap,<sup>r</sup> in region Ω<sup>r</sup> at the interface Γ be some fraction α<sup>k</sup> <sup>r</sup> of the total capacitive current *I*cap:

$$I\_{\rm cap,r}^k = \alpha\_r^k I\_{\rm cap}.\tag{2.14}$$

Specifically, we assume that:

$$\alpha\_r^k = \frac{D\_r^k (z^k)^2 [k]\_r}{\sum\_{l \in K} D\_r^l (z^l)^2 [l]\_r},\tag{2.15}$$

and note that <sup>k</sup> <sup>∈</sup><sup>K</sup> α<sup>k</sup> <sup>r</sup> = 1 for*r* ∈ {*i*, *e*}. By the above definitions, (2.10) and (2.12), we let the intracellular and extracellular ion fluxes across the membrane be given by:

$$\mathbf{J}\_i^k \cdot \mathbf{n}\_i = \frac{I\_{\rm ion}^k + \alpha\_i^k (I\_m - I\_{\rm ion})}{F \varepsilon^k}, \qquad -\mathbf{J}\_e^k \cdot \mathbf{n}\_e = \frac{I\_{\rm ion}^k + \alpha\_e^k (I\_m - I\_{\rm ion})}{F \varepsilon^k}, \tag{2.16}$$

for *k* ∈ *K*.

#### **2.2.2.1 Modelling Specific Ion Channels**

The membrane channel currents *I*<sup>k</sup> ion(v) for each ion species *<sup>k</sup>* are subject to modelling. These currents are typically expressed on the form:

$$I\_{\rm ion}^k(\nu) = g\_L^k(\nu - E^k),\tag{2.17}$$

where g<sup>k</sup> <sup>L</sup> is the conductivity, and *<sup>E</sup>*<sup>k</sup> is the ion specific reversal potential (or Nernst potential), given by:

$$E^k = \frac{RT}{\varepsilon^k F} \ln \frac{c\_e^k}{c\_i^k}.\tag{2.18}$$

This Nernst potential depends on the concentration ratio, whereas the Nernst potential in models without explicit modelling of ion concentrations is constant. Typical models include synaptic input currents, passive neuronal leak channels, or e.g. the Hodgkin-Huxley model (9). For more details on membrane current models and modelling, see e.g. (18).

## **2.2.3 Summary of KNP-EMI Equations**

The KNP-EMI model equations follow from inserting (2.1) into (2.5)–(2.6), combined with (2.7), (2.8), (2.9), (2.10), and (2.16), and read as follows.

For each ion species *k* ∈ *K* and each region *r* ∈ {*i*, *e*}, find the *ion concentrations c*k <sup>r</sup> : Ω<sup>r</sup> × (0,*T*] → R (mol/m3), the *electrical potentials u*<sup>r</sup> : Ω<sup>r</sup> × (0,*T*] → R (V), and the *total transmembrane current density <sup>I</sup>*<sup>m</sup> : <sup>Γ</sup> × (0,*T*] → <sup>R</sup> (A/m2) such that1:

$$\nabla \cdot (F \sum\_{k} z^{k} \mathbf{J}\_{r}^{k}) = 0 \tag{2.19}$$

$$\frac{\partial c\_r^k}{\partial t} + \nabla \cdot \mathbf{J}\_r^k = 0 \tag{2.20}$$

$$-F\sum\_{k} z^{k} \mathbf{J}\_{e}^{k} \cdot \mathbf{n}\_{e} = F \sum\_{k} z^{k} \mathbf{J}\_{i}^{k} \cdot \mathbf{n}\_{i} \equiv I\_{m} \tag{2.21}$$

$$\nu = \mu\_i - \mu\_e \tag{2.22}$$

$$\frac{\partial \mathbf{v}}{\partial t} = \frac{1}{C\_m} (I\_m - I\_{\text{ion}}) \tag{2.23}$$

where the ion flux density **J**<sup>k</sup> <sup>r</sup> is given by (2.2), and *I*ion is subject to modelling. A set of initial and boundary or compatibility conditions will close the system.

<sup>1</sup> Note that the additional negative signs in (2.19) and (2.21), compared with the corresponding equations in Chapter 1, result from our physically consistent definition of the ion flux density **J**<sup>k</sup> r as the negative gradient, cf. (2.2).

# **2.3 Numerical Solution of the KNP-EMI Equations**

The KNP-EMI model defines a complicated system of time-dependent, nonlinear, mixed dimensional partial differential equations. The number of unknowns depends on the number of ion species modelled. Some of the variables exist in the intracellular and extracellular domains, while others live on the lower-dimensional membrane. This setting is numerically challenging and calls for advanced techniques.

To solve the KNP-EMI model numerically, one may consider a finite difference scheme to approximate the time derivatives, a linearization of ion flux densities **J**<sup>k</sup> r and fractions α<sup>k</sup> <sup>r</sup> , a splitting scheme to handle active ion channel current models, and a finite element discretization in space. Such a solution algorithm is detailed in (5), and we refer the reader to this description for further details.

# **2.4 Comparing KNP-EMI and EMI during Neuronal Hyperactivity**

Neurons are negatively charged relative to their environment, with a resting membrane potential of about −70 mV. This resting potential is maintained by low concentrations of sodium ions (Na+) and high levels of potassium ions (K+) inside the cell (23). Action potentials (neuronal activity) are generated by the opening of sodium and potassium channels in the cell membranes. The ionic gradient will drive sodium into the cell and depolarize the cell membrane. Next, the potassium channels open causing an outflux of potassium which in turn repolarizes the cell.

As a result, there is a continuous need to pump potassium into the intracellular space and sodium out to the extracellular space to restore the electrochemical gradient across the cell membrane. One of the key mechanisms for this process is the Na/K/ATPase pump. The Na/K/ATPase pump actively transports 3 Na<sup>+</sup> ions out of the cell and 2 K<sup>+</sup> ions into the cell (7; 14; 20). Several pathologies are associated with increased neuronal activity, e.g. seizures and epilepsy (10; 6; 1), and spreading depression (22). In periods of neuronal hyperactivity, the Na/K/ATPase pumps may not be able to restore the concentrations to baseline levels. Consequently, the electrochemical gradients may be reduced, and silenced neuronal activity and cellular swelling may occur (13).

The ion concentration gradients observed during neuronal hyperactivity thus yields a suitable setting for illustrating differences between the KNP-EMI and the EMI frameworks. In particular, we compare the two frameworks both during normal neuronal activity (firing rate of 1 Hz) and during hyperactivity (firing rate of 50 Hz).

## **2.4.1 Model Parameters and Membrane Mechanisms**

We consider two idealized axons, represented by two parallel, rectangular domains, surrounded by extracellular space in three dimensions. The diameter of each axon is 2.0 · 10−<sup>7</sup> m, and they are separated by 1.0 · 10−<sup>7</sup> m of extracellular space. Parameter values are as listed in Table 2.1. We refer to the supplementary code for a complete description of the model set-up (4).

**KNP-EMI membrane mechanisms** The membrane mechanisms in the KNP-EMI model, cf. (2.11), are modelled using the standard Hodgkin-Huxley model (9) combined with a model for the Na/K/ATPase pump (12), the KCC2 cotransporter (24) and the NKCC1 cotransporter (24). The Na/K/ATPase pump current *I*ATP (A/m2) is modelled as:

$$I\_{\rm ATP} = I\_{\rm ATP}(c\_i^{\rm Na}, c\_e^{\rm K}) = \frac{\hat{I}}{(1 + \frac{m\_{\rm K}}{c\_e^{\rm K}})^2 (1 + \frac{m\_{\rm Na}}{c\_i^{\rm Na}})^3},\tag{2.24}$$

i

e

where ˆ*I* is the maximum pump strength and *m*<sup>K</sup> and *m*Na denote the pump threshold for extracellular potassium and intracellular sodium, respectively. Further, the transmembrane currents generated by the KCC2 cotransporter *I*KCC2 (A/m2) and the NKCC1 cotransporter *I*NKCC1 (A/m2) are modelled as:

$$I\_{\rm KCC2} = S\_{\rm KCC2} \ln(\frac{c\_i^{\rm K} c\_i^{\rm Cl}}{c\_e^{\rm K} c\_e^{\rm Cl}}),\tag{2.25}$$

$$I\_{\rm NKCC1} = S\_{\rm NKCC1} \frac{1}{1 + e^{16 - c\_e^{\rm K}}} (\ln(\frac{c\_i^{\rm K} c\_i^{\rm Cl}}{c\_e^{\rm K} c\_e^{\rm Cl}}) + \ln(\frac{c\_i^{\rm Na} c\_i^{\rm Cl}}{c\_e^{\rm Na} c\_e^{\rm Cl}})),\tag{2.26}$$

where *S*KCC2 and *S*NKCC1 are the maximal cotransporter strengths. Moreover, the cell is stimulated by prescribing a synaptic input *I*syn of the form:

$$I\_{\rm syn}^k = \mathcal{g}\_{\rm syn} H e^{\frac{\nu - t\_0}{\alpha}} (\nu - E^k), \tag{2.27}$$

where α (s) is the synaptic time constant, *H* is the Heaviside function for a small region on the left side of the axons, and <sup>g</sup>syn <sup>=</sup> <sup>1</sup>.<sup>25</sup> · <sup>10</sup>−<sup>3</sup> S/m2. In summary, the membrane channel currents for sodium, potassium and chloride are modelled as:

$$\begin{split} I\_{\mathrm{ion}}^{\mathrm{Na}}(\boldsymbol{\upsilon},\boldsymbol{c}\_{r}^{\mathrm{k}}) &= \mathbf{g}\_{\mathrm{leak}}^{\mathrm{Na}}(\boldsymbol{\upsilon} - \boldsymbol{E}^{\mathrm{Na}}) + \bar{\mathbf{g}}^{\mathrm{Na}}\boldsymbol{m}^{3}h(\boldsymbol{\upsilon} - \boldsymbol{E}^{\mathrm{Na}}) + 3I\_{\mathrm{ATP}} + I\_{\mathrm{NKCC1}} + I\_{\mathrm{syn}}^{\mathrm{Na}}, \\ I\_{\mathrm{ion}}^{\mathrm{K}}(\boldsymbol{\upsilon},\boldsymbol{c}\_{r}^{\mathrm{k}}) &= \mathbf{g}\_{\mathrm{leak}}^{\mathrm{K}}(\boldsymbol{\upsilon} - \boldsymbol{E}^{\mathrm{K}}) + \bar{\mathbf{g}}^{\mathrm{K}}\boldsymbol{n}^{4}(\boldsymbol{\upsilon} - \boldsymbol{E}^{\mathrm{K}}) - 2I\_{\mathrm{ATP}} + I\_{\mathrm{NKCC1}} + I\_{\mathrm{KCC2}} \\ I\_{\mathrm{ion}}^{\mathrm{Cl}}(\boldsymbol{\upsilon},\boldsymbol{c}\_{r}^{\mathrm{k}}) &= \mathbf{g}\_{\mathrm{leak}}^{\mathrm{Cl}}(\boldsymbol{\upsilon} - \boldsymbol{E}^{\mathrm{Cl}}) - 2I\_{\mathrm{NKCC1}} - I\_{\mathrm{KCC2}}, \end{split}$$

where, g<sup>k</sup> leak and ¯g<sup>k</sup> is the leak conductivity and the maximal conductivity for ion species *k*, respectively, the Nernst potential *E*<sup>k</sup> for ion species *k* is as described in Section 2.2.2.1, and the gating variables *m*, *h* and *n* are described by the standard Hodgkin-Huxley ODEs, see e.g. (23) for details.

**EMI membrane mechanisms** For the EMI model, we apply the standard Hodgkin-Huxley model and stimulate the cell by prescribing an input current of the form (2.27); thus, the membrane channels currents are modelled as:

$$\begin{split} I\_{\mathrm{ion}}(\boldsymbol{\nu}) &= \mathbf{g}\_{\mathrm{leak}}^{\mathrm{Na}}(\boldsymbol{\nu} - \boldsymbol{E}^{\mathrm{Na}}) + \mathbf{g}\_{\mathrm{leak}}^{\mathrm{K}}(\boldsymbol{\nu} - \boldsymbol{E}^{\mathrm{K}}) + \mathbf{g}\_{\mathrm{leak}}^{\mathrm{Cl}}(\boldsymbol{\nu} - \boldsymbol{E}^{\mathrm{Cl}}) \\ &+ \bar{\mathbf{g}}^{\mathrm{Na}} m^{3} h(\boldsymbol{\nu} - \boldsymbol{E}^{\mathrm{Na}}) + \bar{\mathbf{g}}^{\mathrm{K}} n^{4} (\boldsymbol{\nu} - \boldsymbol{E}^{\mathrm{K}}) + I\_{\mathrm{ATP}} + I\_{\mathrm{syn}}, \end{split}$$

where *E*K, *E*Na and *E*Cl are calculated by (2.18) with the initial values from the KNP-EMI model for the sodium and potassium concentrations. Similarly, the bulk conductivities σ<sup>i</sup> and σ<sup>e</sup> are calculated by (2.4), and the net current from the Na/K/ATPase pump *I*ATP is given by (2.24). Finally, there is no contribution from KCC2 and NKCC1, as both cotransporters mediate ion transport without any net charge movement across the membrane.


Table 2.1: The physical and model parameters used in the simulations. The values are collected from Sterratt et al. (23), Hodgkin et al. (9), Wei et al. (24), whereas the values marked with \* are computed by a steady state estimation.

The initial conditions for the intra- and extracellular ion concentrations, the membrane potential and the gating variables are listed in Table 2.2. At the exterior boundary, we apply no flux boundary conditions for each ion species.


Table 2.2: Initial conditions. The initial ion concentrations are chosen such that the Nernst potentials are equal to those in the Hodgkin-Huxley model (9). The membrane potential is computed by a steady state estimation.

## **2.4.2 Results and Discussion**

During normal activity, the KNP-EMI and the EMI models behave similarly, both for the membrane potential and the extracellular potential (Figure 2.1 A, B). The stimuli current depolarizes the membrane potential above the threshold for firing, and an action potential is initiated (Figure 2.1 A). Simultaneously, the extracellular potential decreases by ∼ 0.13 mV, before quickly returning to baseline (Figure 2.1 B).

During hyperactivity, the KNP-EMI and EMI models differ (Figure 2.1 C, D, E, F). In both models, repeated action potentials are triggered. But, for the KNP-EMI model, we observe changes in the membrane potential between hyperpolarization phases. In particular, we conclude that the KNP-EMI membrane resting potential increases with repeated firing: after 5 action potentials (at *t* = 90 ms) the membrane potential has a minimum value of −75 mV, which is an 9% increase from the first action potential. Eventually, the membrane is depolarized to the point where action potentials can long longer be fired (Figure 2.1 E).

The observed changes are caused by alterations in the ion concentration gradients. For each action potential, the extracellular Na<sup>+</sup> concentration decreases by 0.15 mM and the extracellular K<sup>+</sup> concentration increases by 0.16 mM (Figure 2.2 A, B). During normal activity (Figure 2.2 A, B), the ion concentrations will slowly be pumped back toward baseline levels, and the membrane potentials are not substantially affected by the small ion concentration changes. However, in the case of hyperactivity, the membrane mechanisms (i.e. pumps and cotransporters) are not able to keep up. Consequently, the extracellular Na<sup>+</sup> concentration will keep decreasing and the extracellular potassium will keep increasing, causing the cell to depolarize (Figure 2.2 C, D).

In the KNP-EMI model (Figure 2.2 A, B), we note that 7.92 % of the extracellular K<sup>+</sup> concentration is restored, and 7.3 % of the extracellular Na<sup>+</sup> concentration is restored after 100 ms. That is, the extracellular concentrations do not reach baseline levels within the simulation period. Other studies have reported that it takes on the order of minutes (0.5 minutes (19), 6 minutes (2)) before the concentrations return to baseline after neuronal activity.

# **2.5 Conclusions and Outlook**

In this chapter, we have presented a mathematical model, the KNP-EMI model, for ionic electrodiffusion in excitable tissue with an explicit representation of the intracellular, extracellular and membrane domains. For further reading on methodological aspects, we refer to (5; 15) and references therein. This model extends on the EMI model presented in Chapter 1 and may be more accurate in situations with rapid and persistent changes in ion concentrations. Moreover, the KNP-EMI framework allows for modelling ligand-gated ion channels (e.g. NMDA receptors).

The complexity of the KNP-EMI system yields a number of numerical challenges. The mere number of unknowns result in large systems of equations calling for efficient solution techniques. The nonlinearities in the system can easily lead to non-convergence and thus call for robust algorithms. Moreover, the coupling of full and lower dimensional domains and fields calls for well-posed numerical methods together with suitable simulation software. Further, the system couples different time scales: from neuronal action potentials taking place at the microscale to the slower diffusion process. In short, modelling ionic electrodiffusion in the EMI setting is an area with vast opportunities for further research.

**Acknowledgements** The authors would like to thank Geir Halnes for useful discussions and Min Ragan-Kelley for technical assistance.

Fig. 2.1: Comparison of potentials over time at fixed points in space predicted by the KNP-EMI and the EMI frameworks during normal activity (upper panels) and during hyperactivity (mid and lower panels). The membrane potentials for KNP-EMI and EMI during normal activity (**A**) and hyperactivity (**C**, **E**, **F**), and the extracellular potentials for KNP-EMI and EMI during normal activity (**B**) and hyperactivity (**D**).

Fig. 2.2: Time development of extracellular ion concentrations at a fixed point in space for the KNP-EMI framework during normal activity (upper panels) and hyperactivity (lower panels). The extracellular sodium (**A**) and potassium (**B**) concentrations during normal activity, and the extracellular sodium (**C**) and potassium (**D**) concentrations during hyperactivity.

**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.

# **References**


# **Chapter 3**

# **Modeling Cardiac Mechanics on a Sub-Cellular Scale**

Åshild Telle1, Samuel T. Wall1 and Joakim Sundnes<sup>1</sup>

**Abstract** We aim to extend existing models of single-cell mechanics to the EMI framework, to define spatially resolved mechanical models of cardiac myocytes embedded in a passive extracellular space. The models introduced here will be pure mechanics models employing fairly simple constitutive laws for active and passive mechanics. Future extensions of the models may include a coupling to the electrophysiology and electro-diffusion models described in the other chapters, to study the impact of spatially heterogeneous ion concentrations on the cell and tissue mechanics.

# **3.1 Introduction**

A vast range of models have been developed for the force development of cardiac and skeletal muscle, on the scale of a single cross bridge (10), myofilament (3), sarcomere (2), and the complete cell (13). The scales involved and the main functional units considered on each scale are schematically illustrated in Figure 3.1. Common to most existing models is the fact that they focus on a single spatial scale, and any coupling between scales is fairly crudely represented. As an example, the model by Rice et al. (13) is essentially a model of a single sarcomere (Fig. 3.1 D), which is normalized and then scaled to yield a realistic force output for cell- and tissue-level mechanics applications. Other models provide detailed descriptions of mechanisms and interactions on a molecular level (Fig. 3.1 F)(4; 3), and are able to capture many of the characteristic non-linearities of muscle cell mechanics. However, key aspects of mechanical activation and force-length relationships are still not fully understood, and they may be the result of interactions between individual sarcomeres and myofibril bundles. A few attempts have been made at modeling interactions at this scale, and

<sup>1</sup>Simula Research Laboratory, Norway

Fig. 3.1: The heart (A) is mainly composed of cardiac muscle cells, also called cardiomyocytes (B). Each cell (C) is composed of long tubes of sarcomeres (D), in which the thin and thick myofilaments overlap in layers (E). The interaction between these (F) causes the cardiac muscle to contract in a process called the cross-bridge cycle.

have shown potentially interesting emergent behaviours (2; 11). Furthermore, heart failure and other pathologies are linked with heterogeneous intracellular calcium concentration resulting from disruptions of the calcium regulation system. Describing the effect of such heterogeneities on the cell contraction and force development requires spatially resolved mechanics models on the sub-cellular scale.

Finite element models of contracting myocytes have been proposed (8; 14), and have been used to explore the impact of model assumptions, calcium heterogeneity, and boundary conditions. The model presented by Ruiz-Baier et al. (14) describes the individual myocyte as a hyperelastic material, and uses an active strain approach to describe the contraction. Both the passive and active mechanical properties are assumed to be homogeneous, but sub-cellular heterogeneities can easily be introduced. We here propose to extend the single myocyte model in (14) to include the extracellular domain, and to model collections of cells, based on similar ideas used for the electrophysiology model presented in (17; 18) and (7, Chapter 1).

# **3.2 Models and Methods**

The motion and deformation of the heart can be described by the classical theory of non-linear solid mechanics. The primary unknown in our computational model will be the displacement vector *u*, which for each *material point* describes the difference between its current and original position. We have *u* = *x* − *X*, where *X* is the original (reference) position of a point, and *x* is its position after the deformation. From the displacement vector we can define the *deformation gradient F* = ∂*x*/∂*X* = *I* + ∂*u*/∂*X*, which is an essential quantity describing the deformation of a solid. See for instance (6) for a detailed introduction to non-linear solid mechanics.

A characteristic feature of the heart and other muscles is that they contract and deform even in the absence of external loads. The overall deformation and mechanical state of the heart depends both on this active contraction and on the passive mechanical properties of the tissue. There are two main approaches for modeling the coupling of active and passive mechanics in cardiac tissue, often referred to as *active strain* and *active stress*. Both approaches are based on modeling the active and passive contributions separately, then combining them into a complete coupled model.

In the active strain approach, the active-passive coupling is incorporated through a multiplicative decomposition of the deformation gradient *F* into active and passive components, *F* = *F*p*F*a. Here, *F*<sup>a</sup> represents an active deformation governed by the cell state, and *F*<sup>p</sup> is a passive elastic deformation which ensures compatibility with loads and kinematic boundary conditions. The active stress approach is based on an additive split of the stress tensor into its active and passive components. In terms of the first Piola-Kirchhoff stress tensor *P*, the stress is written as *P* = *P*<sup>p</sup> + *P*a, where *P*<sup>a</sup> is a function of the cellular activation state and *P*<sup>p</sup> is a standard elastic stress derived from a strain energy function.

Both of these approaches have their strengths and weaknesses. In general, the active strain approach is considered to be more suitable for deriving mathematically wellbehaved constitutive laws, while the active stress concept is more easily coupled to biophysically detailed models of cell contraction.

## **3.2.1 Fundamental Equations**

In this study we will primarily use the active stress approach, but for completeness we also present the equations arising from the active strain approach. This model can be derived as a direct extension of the single myocyte model in (14), using a similar approach as in (17; 18) to consider both the intra- and extracellular domains:

Fig. 3.2: Illustration of the intra- and extracellular domains for a single cell and its surroundings.

$$\begin{array}{llll} a & \mathbf{\nabla \cdot P\_{i} = 0, & b & : & P\_{i} = \frac{\partial \Psi\_{i}}{\partial F\_{i}}, & c & : & F\_{i} = F\_{i}^{P} F\_{i}^{a}, & \text{in } \Omega\_{i}, \\ d & \mathbf{: } \nabla \cdot P\_{e} = 0, & e & : P\_{e} = \frac{\partial \Psi\_{e}}{\partial F\_{e}}, & & & \text{in } \Omega\_{e}, \\ f & \mathbf{: } u\_{i} = u\_{e}, & g & : n\_{i} \cdot P\_{i} = n\_{e} \cdot P\_{e}, & & & \text{on } \Gamma, \\ h & \mathbf{: } n\_{e} \cdot P\_{e} = 0, & & & \text{on } \partial \Omega\_{e, \Gamma}, \\ i & \mathbf{: } u = 0, & & & \text{on } \partial \Omega\_{e, \Gamma}. \end{array} \tag{3.1}$$

Here, Ω<sup>i</sup> and Ω<sup>e</sup> are the intra- and extracellular domains, respectively, Γ is the interface between the domains, with the normal vector *n*<sup>i</sup> pointing out of the intracellular domain and *n*<sup>e</sup> out of the extracellular domain. Finally, ∂Ωe,<sup>T</sup> and ∂Ωe,<sup>D</sup> are the parts of the outer boundary ∂Ω<sup>e</sup> subject to traction- and displacement boundary conditions, respectively. See Figure 3.2 for a sketch of a typical computational domain, including a single cell and its immediate surroundings. Following (14), we here apply the active strain approach to incorporate active contraction of the myocyte, where the intracellular deformation gradient *F*<sup>i</sup> is decomposed as described above. The passive part is assumed to be hyper-elastic and derived from a strain energy function, see for instance (6) for details. A common choice for the active part is *F*a <sup>i</sup> <sup>=</sup> *dia*g((<sup>1</sup> <sup>−</sup> <sup>γ</sup>),(<sup>1</sup> <sup>−</sup> <sup>γ</sup>) <sup>−</sup>1/2,(1 − γ) <sup>−</sup>1/2), where γ describes the fiber contraction and is a function of the cell activation state. For a more detailed introduction and discussion of active strain models, we refer to (1).

The active stress model is the most widely used approach for modeling coupled active and passive mechanics on tissue level, and this is the approach we will employ in the subsequent numerical experiments. In the present context the active stress model involves a decomposition of the intracellular first Piola-Kirchhoff stress *P*<sup>i</sup> into a passive elastic part *P*<sup>p</sup> <sup>i</sup> and an active part *<sup>P</sup>*<sup>a</sup> <sup>i</sup> . The passive stress is derived from a strain energy function in the usual way, while the active stress is a function of the cell activation state. For the simplified model considered here we write the active stress as a function of time and the local fiber stretch λ, but the approach can easily be extended to include detailed biophysical models of the contractile mechanisms. The full active stress model may be written as

32 Telle et al.

$$\begin{array}{llll} a: \,\nabla \cdot P\_{i} = 0, & b: \, P\_{i} = \frac{\partial \Psi\_{i}}{\partial F} + P\_{i}^{a}(t, \lambda), & \text{in } \Omega\_{i}, \\ c: \,\nabla \cdot P\_{e} = 0, & d: \, P\_{e} = \frac{\partial \Psi\_{e}}{\partial F\_{e}}, & \text{in } \Omega\_{\mathrm{e}}, \\ e: \,\boldsymbol{u}\_{i} = \boldsymbol{u}\_{e}, & f: \,\boldsymbol{n}\_{i} \cdot P\_{i} = \boldsymbol{n}\_{e} \cdot P\_{e}, & \text{on } \Gamma, \\ \text{g}: \,\boldsymbol{n}\_{e} \cdot P\_{e} = 0 & \text{on } \partial \Omega\_{\mathrm{e},T}, \\ h: \,\boldsymbol{u} = 0 & \text{on } \partial \Omega\_{\mathrm{e},D}. \end{array} \tag{3.2}$$

Both approaches treat the extracellular domain in the same way, as a passive hyperelastic material governed by a strain energy function Ψe. As given by (3.1) f-g and (3.2) g-h we assume continuity of stresses *P*i, *P*<sup>e</sup> and displacements *u*i,*u*<sup>e</sup> across the cell membrane Γ, implying that the membrane itself has no stiffness. The outer boundary Ω<sup>e</sup> is assumed to be stress free, with Dirichlet conditions applied to parts of the boundary to avoid rigid body motion. Models for the active stress *P*<sup>a</sup> <sup>i</sup> come in many forms, including simple phenomenological models as well as detailed biophysical models of cell electro-mechanics (12; 13). For the present study we apply a simple model where the active stress is derived from a (pseudo-) strain energy in the same way as the passive stress:

$$P\_i^a = \frac{\partial \Psi\_i^a}{\partial F}.\tag{3.3}$$

Here, Ψ<sup>a</sup> <sup>i</sup> is given by

$$
\Psi\_i^a = \frac{T\_{active}(t)}{2} \lambda^2,
$$

where λ = ||*Fe*1|| is the stretch in the so-called *fiber direction* (i.e. the main orientation of the muscle cells), defined by the unit vector *e*1, and *T*active(*t*) is a prescribed function defining the active contractile force as a function of time.

## **3.2.2 Specific Model Choices**

In this section we describe specific choices of the constitutive laws describing active and passive material properties in the models above, to arrive at a complete model that can be solved for the deformations and stresses. As noted above, we will in the following only consider the active stress model, given by (3.2). For the strain energy defining the passive stress-strain relationships we have applied a model from (19), which belongs to the family of models first presented by Guccione et al. (5). The same form of strain energy is used in the intra- and extracellular domains, but we allow the material parameters to be different. Both domains are modeled as nearly incompressible, with volume changes during deformations controlled by a penalty term. We have

$$\Psi\_i = C\_i(e^{Q\_i} - 1) + \kappa (J \ln J - J + 1) \qquad \qquad x \in \Omega\_i,\tag{3.4}$$

$$\Psi\_e = C\_e(e^{Q\_e} - 1) + \kappa (J \ln J - J + 1) \qquad \qquad x \in \Omega\_e,\tag{3.5}$$

where *Q*i,*Q*<sup>e</sup> are functions depending on components of the Green-Lagrange strain tensor *E* = <sup>1</sup> <sup>2</sup> (*F*<sup>T</sup> *<sup>F</sup>* <sup>−</sup> *<sup>I</sup>*):

$$\begin{aligned} \mathcal{Q}\_j = b\_{f,j} E\_{11}^2 + b\_{t,j} (E\_{22}^2 + E\_{33}^2 + E\_{23}^2 + E\_{32}^2) \\ + b\_{fs,j} (E\_{12}^2 + E\_{21}^2 + E\_{13}^2 + E\_{31}^2). \end{aligned} \tag{3.6}$$

Furthermore *C*j, *b*<sup>f</sup> ,j, *b*t,j, and *b*f s,j, for *j* = *i*, *e* are material parameters characterizing the material's stiffness to the various strain modes, κ is a penalty parameter that controls the volume changes, and *J* = det *F*. For a fully incompressible deformation we have *J* = 1, and in our nearly incompressible model we tune the parameter κ to keep *J* ≈ 1.

In its most general form, the materials described by (3.4)-(3.6) are are transversely isotropic, which is a special case of orhtotropic materials. While an orthotropic material has different mechanical properties in three different directions, a transversely isotropic material is isotropic in planes normal to a characteristic direction. Passive cardiac tissue is known to behave as an orthortopic material (9), with the three directions dictated by the orientation and organization of the myocytes. However, a transversely isotropic material is shown to be a good approximation, with material isotropy in planes normal to the fiber direction, the main orientation of the muscle cells. The details of the intra- and extracellular material behavior in our micro-structural model are less well-studied, and the degree of anisotropy has not been characterized. From the microstructure of the contractile apparatus occupying most of the intracellular space (see Figure 3.1) it is natural to assume anisotropic behavior, but the exact degree of aniostropy is not known. As a starting point, we set the intracellular material parameters to

$$b\_{f,i} = 8, \qquad b\_{l,e} = 2, \qquad b\_{fs,e} = 4. \tag{3.7}$$

For the extracellular space we assume isotropic material behaviour, setting

$$b\_{f,e} = b\_{t,e} = b\_{fs,e} = 1.\tag{3.8}$$

The bulk compressibility was set to κ = 1000 kPa in both domains, while we explored different values of the scaling parameters *C*<sup>i</sup> and *C*e, to be specified below.

For the active stress model defined in (3.3) we have used a pre-computed transient tension *T*active(*t*) as shown in Figure 3.3. The curve was computed using the model of Rice et al. (13) with default parameters, which outputs a normalized force. This value was then scaled such that the peak value reaches 2 kPa, giving a reasonable contractile stress for our application.

Fig. 3.3: Transient tension *T*active(*t*) over time (left), first computed in (13), then scaled to give values on a reasonable scale. In the intracellular domain the active tension is homogeneously set to this value; in the extracellular domain there is no such tension, implemented as being set to zero for all time steps.

## **3.2.3 Numerical Methods**

The problem defined by (3.2) is solved with the displacement *u* as the primary unknown. To solve the system with the finite element method, it is convenient to formulate it as a single PDE defined over the entire domain Ω = Ω<sup>i</sup> ∪ Ωe. Such a formulation is not possible for the strong form of the PDEs, so we first need to derive the weak form of the equations. Starting with (3.2)a, we define a suitable vector function space *V*(Ωi) defined over the intracellular domain, multiply the equation with a test function v <sup>∈</sup> *<sup>V</sup>*(Ωi) and integrate by parts, to arrive at a weak formulation

$$\int\_{\Omega\_{\bar{t}}} P\_{\bar{t}} \cdot \nabla \mathbf{v} \, d\mathbf{x} - \int\_{\Gamma} (n\_{\bar{t}} \cdot P\_{\bar{t}}) \mathbf{v} = \mathbf{0}.$$

This equation is to be satisfied for all <sup>v</sup> <sup>∈</sup> *<sup>V</sup>*(Ωi). Performing the same steps for the extracellular domain, and using the boundary condition (3.2)g on the outer boundary, we get

$$\int\_{\Omega\_e} P\_e \cdot \nabla \nu d\mathbf{x} - \int\_{\Gamma} (n\_e \cdot P\_e) \nu = 0.1$$

This equation should hold for all test functions <sup>v</sup> <sup>∈</sup> *<sup>V</sup>*(Ωe), where *<sup>V</sup>*(Ωe) is a suitable space of functions defined over the domain Ωe. Using similar arguments as in (15), we can define a function space *V*(Ω) as the set of functions defined over Ω that belong to both *V*(Ωi) and *V*(Ωe) and are continuous over Γ. With this definition, we may add the two weak forms above to obtain

$$\int\_{\Omega\_i} P\_i \cdot \nabla \boldsymbol{\nu} \, d\mathbf{x} - \int\_{\Gamma} (\boldsymbol{n}\_i \cdot \boldsymbol{P}\_i) \boldsymbol{\nu} + \int\_{\Omega\_e} P\_e \cdot \nabla \boldsymbol{\nu} \, d\mathbf{x} - \int\_{\Gamma} (\boldsymbol{n}\_e \cdot \boldsymbol{P}\_e) \boldsymbol{\nu} = 0.1$$

Which is to be satisfied for all <sup>v</sup> <sup>∈</sup> *<sup>V</sup>*(Ω). Since *<sup>n</sup>*<sup>e</sup> <sup>=</sup> <sup>−</sup>*n*i, the surface integrals over Γ cancel because of (3.2)f. We can also use (3.2)e to define a single displacement field over Ω, and we are left with the following weak form: Find *u* ∈ *V*(Ω) such that

Fig. 3.4: **A**: Volume element of one single cell; lines indicate cross section area. **B**: Cross section along longitudial direction of the cell. **C**: Volume element, 5 x 5 cells; lines indicate cross section area. **D**: Cross section along longitudial direction of the cells.

$$\int\_{\Omega} P \cdot \nabla \mathbf{v} \, d\mathbf{x} = 0,\tag{3.9}$$

is satisfied for all v <sup>∈</sup> *<sup>V</sup>*(Ω). with *<sup>P</sup>* defined by (3.2)b and (3.2)d in the respective domains.

# **3.3 Results**

We here present a number of numerical experiments to illustrate the general behavior of the models defined above. The code is implemented using FEniCS, and an archieved version of the code is available, see (16).

For the simulations we used two different meshes; one representing a single cell and one representing a sheet of five by five cells, see Figure 3.4. Both meshes include subdomains defining the intra- and extracellular domains, separated by the cell membrane. To avoid rigid body motion, we keep a few points in the middle fixed. The rest of the boundary is kept unloaded to allow free contraction of the cells.

For each experiment we calculated the Green-Lagrange strain tensor *E* and the Cauchy stress tensor σ, given by

$$\sigma = \begin{cases} |F|^{-1} P\_i F^T & x \in \Omega\_i \\ |F|^{-1} P\_e F^T & \text{otherwise.} \end{cases}$$

On matrix form we can can write these out as

$$E = \begin{bmatrix} E\_{11} \ E\_{12} \ E\_{13} \\ E\_{21} \ E\_{22} \ E\_{23} \\ E\_{31} \ E\_{32} \ E\_{33} \end{bmatrix} \qquad \sigma = \begin{bmatrix} \sigma\_{11} \ \sigma\_{12} \ \sigma\_{13} \\ \sigma\_{21} \ \sigma\_{22} \ \sigma\_{23} \\ \sigma\_{31} \ \sigma\_{32} \ \sigma\_{33} \end{bmatrix}$$

and for each of these we present plots for the first and the middle components, (*E*11, *E*22, σ11, σ22), which characterize strain and stress in the fiber and cross-fiber directions.

Fig. 3.5: Tracking points, for which we evaluate functions of interest across various experiments. The points are uniformly distributed on a line from one corner to the middle, in the xy-direction, corresponding to the cross-section shown in Figure 3.4. Two of them are both located in the extracellular subdomain, and one should expect them to show different patterns than the three located in the intracellular subdomain.

We first considered a single cell, and simulated contraction over a single cardiac cycle with homogeneous active force applied throughout the cell. For this simulation we chose parameter values *C*<sup>e</sup> = *C*<sup>i</sup> = 0.5. The results are presented in Figure 3.6, where we observe that the deformation follows the expected pattern of a contraction in the longitudinal direction of the cell. Furthermore, in spite of the homogeneous applied active stress we see slight spatial variations in the deformation state, resulting from the discontinuity of active force across the cell membrane.

Similar patterns are observed in the simulation of the sheet of 25 cells, shown in Figure 3.7. In this experiment the same active stress transient through the intracellular domain of all the cells, with the same material parameters. We still observe spatial variations in the deformation pattern – each cells is affected by mechanical deformation around them.

We then considered two cases where we kept all parameters but one fixed, exploring the choices of material stiffness parameters *C*<sup>e</sup> and *C*i. The results are presnted in

Fig. 3.6: First and middle components of the Cauchy stress tensor σ and Green-Lagrange strain *E*, for a single cell. The plots to the left shows values plotted over time, for the first 500 ms (out of 1000), following tracking points as shown in Figure 3.5. The plots to the right shows values over a cross-section as shown in Figure 3.4, as the active tension reaches it's peak value. The grey rectangle indicates initial configuration.

Figures 3.8 and 3.9. These simulations were again performed on a mesh representing a single cell, with active force applied as described above. For the first experiment we kept *C*<sup>e</sup> fixed at 0.5, changing *C*i; that is, we let the material stiffness in the extracellular domain remain the same while increasing the stress/strain scaling parameter in the intracellular domain. As *C*<sup>i</sup> increases the material becomes stiffer, and for the same active stress applied, one should expect less contraction. This can indeed be observed; both components of the Cauchy stress tensor (in magnitude) and the strain tensor decreases everywhere, and for the last three parameter choices there is almost no difference in deformation. On the other hand, we still apply an active stress in the intracellular domain, and we observe that the strain close to the membrane doesn't change much even if it changes everywhere else.

For the next experiment we changed to keeping *C*<sup>i</sup> = 0.5 constant, while increasing *C*e. We can observe higher Cauchy stress for the first component, and lower Cauchy stress for the second component, with increasing values of *C*e. The strain decreases for both components. This is exactly as expected – in one end of the spectrum, having

Fig. 3.7: First and middle components of the Cauchy stress tensor σ, and and Green-Lagrange strain *E*, for 5 x 5 cells. Values are plotted over the cross-section as shown in 3.4, as the active tension reaches it's peak value. The grey rectangle indicates initial configuration.

*C*<sup>e</sup> = 0.5, one would expect the extracellular subdomain to not affect the intracellular domain as it's rather flexible. For a given tension in the intracellular domain, it will just move along quite easily, while the overall behaviour in the whole domain is governed by the contraction inside the cell. As *C*<sup>e</sup> increases, the material is modeled as stiffer and hence constrain the movement more. For very high values the material is so stiff that it hardly moves, efficiently keeping the membrane close to fixed.

# **3.4 Discussion**

We have presented a general framework for modeling cardiac mechanics on a subcellular scale, by extending a model of the type defined in (14) to the extracellular domain. A series of preliminary numerical experiments demonstrate that the model behaves as expected, with the discontinuity across the cell membrane giving rise to spatially varying deformation fields even though both the active stress and other model parameters are spatially homogeneous over the intracellular domain.

The main purpose of this work was to present the model framework and to illustrate the general behaviour of the model, while more detailed investigations and model extensions are left for future studies. A complete list of model limitations and potential extensions would be too extensive to present here, but it is worth commenting on a few of the most obvious ones. First, the model derivation above included a number of simplifying assumptions on the mechanical properties of the cell membrane and

Fig. 3.8: First and middle components of the Cauchy stress tensor σ and Green-Lagrange strain *E*, for a single cell, as we vary the parameter *C*i, which defines the stiffness of the material in the intracellular domain. Panel A shows spatial variation over a cross-section of the cell (see Figure 3.4), at peak. Panel B shows how the value, at peak, changes in given tracking points (see Figure 3.5).

Fig. 3.9: First and middle components of the Cauchy stress tensor σ and Green-Lagrange strain *E*, for a single cell, as we vary the parameter *C*e, which defines the stiffness of the material in the extracellular domain. Panel A shows spatial variation over a cross-section of the cell (see Figure 3.4), at peak. Panel B shows how the value, at peak, changes in given tracking points (see Figure 3.5).

the two domains. The continuity of stress across the cell membrane implies that the membrane itself has no stiffness, which is obviously incorrect, but it may be a reasonable assumption for many applications. The impact of different membrane mechanical properties should be explored further in a future study. Similarly, both the intra- and extracellular domains are assumed to be hyperelastic materials, which is probably a fairly crude approximation of the actual behaviour. In reality both of these domains are complex compositions of fluids and various embedded proteins structures, and the material behavior is most likely quite complex. Visco-elastic material models could potentially be a more accurate description than the hyper-elastic models applied here, but the required level of detail and material model complexity remains to be determined. Finally, we have here assumed that both domains are initially in a stress-free resting state, while experiments have shown that the extracellular matrix shrinks considerably when the myocytes are removed. Thus indicates that the resting state is actually an equilibrium state with non-zero stress in both domains, and accurately capturing the overall mechanical behaviour may require including this pre-stress in the model.

In general, the level of detail and complexity of the model formulation will be dictated by the application. Some applications may require further development of the model along the lines suggested above, while for studies of a more qualitative nature the simplest version would be sufficient. One obvious application of the developed model framework, where a fairly simple model would probably give interesting results, is to study the impact of heterogeneities in calcium concentration and mechanical properties on the contractile properties of cells and tissue.

**Acknowledgements** The simulations were performed on resources provided by UNINETT Sigma2 - the National Infrastructure for High Performance Computing and Data Storage in Norway. The meshes used for the simulations were generated using software provided by Miroslav Kuchta.

**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.

# **References**


# **Chapter 4**

# **Operator Splitting and Finite Difference Schemes for Solving the EMI Model**

Karoline Horgmo Jæger1, Kristian Gregorius Hustad1,2, Xing Cai1,<sup>2</sup> and Aslak Tveito1,<sup>2</sup>

**Abstract** We want to be able to perform accurate simulations of a large number of cardiac cells based on mathematical models where each individual cell is represented in the model. This implies that the computational mesh has to have a typical resolution of a few μm leading to huge computational challenges. In this paper we use a certain operator splitting of the coupled equations and show that this leads to systems that can be solved in parallel. This opens up for the possibility of simulating large numbers of coupled cardiac cells.

# **4.1 Introduction**

In recent publications (31; 30; 13) we have shown that a cell-based model is useful for accurately representing the electrophysiology of excitable cells. Traditionally, excitable tissue is simulated based on homogenized models where the cells are not explicitly resolved, see e.g., (26; 7). In the cell-based model, we explicitly represent both the extracellular space (E), the cell membrane (M) and the intracellular space (I), and it is therefore referred to as the *EMI model*. Similar approaches to modeling excitable tissue have been used by several authors; see e.g., (2; 18; 25; 22; 24; 23; 11; 16; 34).

The EMI model is solved, numerically, using an operator splitting scheme which results in two steps; a non-linear system of ordinary differential equations (ODEs) to be solved in each computational node (i.e, degree of freedom) placed on the cell membrane, and a linear system of algebraic equations coupling the discrete Laplace equations of E and I with continuity requirements of the current over M. The spatial

<sup>1</sup>Simula Research Laboratory, Norway

<sup>2</sup>Department of Informatics, University of Oslo, Norway

resolution used in the discretization of the model is usually between 1 μm and 4 μm, thus only 1 mm3 of tissue leads to more than 107 computational nodes. For an adult human cardiac cell, with a resolution of 2 μm, the number of computational nodes per cell (including the associated extracellular space) is about 6000 (see (30), Table 7). Thus, for a limited number of cells, the linear system coupling all the discrete Laplace equations is manageable. In fact, the system was solved using Matlab for up to 16,384 cells, with about 9.8×10<sup>7</sup> computational nodes, see (30).

However, not only the sheer size of the linear system is a challenge, also the properties of the linear system are unusual. In scientific computing, one of the most well-studied problems is solution of linear systems arising from the discretization of elliptic boundary value problems; see e.g., (5; 21; 8). Unfortunately, the EMI system does not naturally fall into the category of elliptic boundary value problems that can be solved using well-developed numerical machinery. It is therefore of importance to develop a splitting strategy for the EMI model that leads to sub-problems of the elliptic type. In (14), we showed that such a splitting can indeed be achieved. Here, we will review this convenient way of splitting the EMI model and show how to solve the system numerically using a finite difference method. Moreover, we will use the numerical scheme to assess the conduction properties in a small collection of cells where a sub-group of the cells are ischemic. Furthermore, we will present a parallel implementation of the splitting strategy, based on using open-source numerical libraries. This code is considerably faster than the existing Matlab code, and well suited for shared-memory parallel computers.

# **4.2 The EMI Model**

We model the electrical properties of collections of cardiac cells using the EMI model introduced in (24; 1; 2; 31; 30). In Figure 4.1 we show the computational domains in the case of two coupled cells. Here, Ω<sup>1</sup> <sup>i</sup> and <sup>Ω</sup><sup>2</sup> <sup>i</sup> denote the intracellular domains, and Ω<sup>e</sup> denotes the extracellular space. The cell membranes are denoted by Γ<sup>1</sup> and Γ2, respectively. The intercalated disc at the intersection between Ω<sup>1</sup> <sup>i</sup> and Ω2 <sup>i</sup> , allowing for currents between the cells, is denoted by Γ1,2. With this notation at hand, the EMI model takes the following form:

∇ · σi∇*u*<sup>k</sup> <sup>i</sup> <sup>=</sup> 0 in <sup>Ω</sup><sup>k</sup> <sup>i</sup> , *<sup>n</sup>*<sup>e</sup> · <sup>σ</sup>e∇*u*<sup>e</sup> <sup>=</sup> <sup>−</sup>*n*<sup>k</sup> <sup>i</sup> · <sup>σ</sup>i∇*u*<sup>k</sup> <sup>i</sup> <sup>≡</sup> *<sup>I</sup>* k <sup>m</sup> at Γk, ∇ · <sup>σ</sup>e∇*u*<sup>e</sup> <sup>=</sup> 0 in <sup>Ω</sup>e, <sup>v</sup><sup>k</sup> <sup>t</sup> = <sup>1</sup> Cm (*I* k <sup>m</sup> − *I* k ion) at Γk, *u*<sup>e</sup> = 0 at ∂Ω<sup>D</sup> <sup>e</sup> , *u*<sup>k</sup> <sup>i</sup> <sup>−</sup> *<sup>u</sup>*k˜ <sup>i</sup> <sup>=</sup> <sup>w</sup><sup>k</sup> at <sup>Γ</sup>k,k˜, *n*<sup>e</sup> · σe∇*u*<sup>e</sup> = 0 at ∂Ω<sup>N</sup> <sup>e</sup> , *<sup>n</sup>*k˜ <sup>i</sup> · <sup>σ</sup>i∇*u*k˜ <sup>i</sup> <sup>=</sup> <sup>−</sup>*n*<sup>k</sup> <sup>i</sup> · <sup>σ</sup>i∇*u*<sup>k</sup> <sup>i</sup> <sup>≡</sup> *<sup>I</sup>*k,k˜ at <sup>Γ</sup>k,k˜, *u*k <sup>i</sup> <sup>−</sup> *<sup>u</sup>*<sup>e</sup> <sup>=</sup> <sup>v</sup><sup>k</sup> at <sup>Γ</sup>k, <sup>w</sup><sup>k</sup> <sup>t</sup> = <sup>1</sup> <sup>C</sup>g (*I*k,k˜ <sup>−</sup> *<sup>I</sup>* k gap) at Γk,k˜, *s*k <sup>t</sup> = *F*<sup>k</sup> at Γ<sup>k</sup> .

Fig. 4.1: **A**: Two-dimensional version of the EMI model domain in the case of two connected cells. Here, the cells Ω<sup>1</sup> <sup>i</sup> and <sup>Ω</sup><sup>2</sup> <sup>i</sup> , with cell membranes denoted by Γ<sup>1</sup> and Γ2, respectively, are connected to each other by the intercalated disc, Γ1,2, and surrounded by an extracellular space, denoted by Ωe. **B**: Two-dimensional illustration of the geometry used for a single cell. The intracellular domain of each cell is composed of five subdomains ΩO, ΩW, ΩE, ΩS, and ΩN. The sizes of the subdomains are specified in Table 4.1.

The model is stated for cell number *k*, and ˜ *k* denotes one of the six neighboring cells (in 3D: north, west, south, east, above, below). In the model, *u*e, *u*<sup>k</sup> <sup>i</sup> , and <sup>v</sup><sup>k</sup> <sup>=</sup> *<sup>u</sup>*<sup>k</sup> <sup>i</sup> <sup>−</sup>*u*<sup>e</sup> denote the extracellular, intracellular, and transmembrane potentials, respectively. Also, <sup>w</sup><sup>k</sup> is the potential difference across the intercalated disc1, <sup>Γ</sup>k,k˜, and <sup>σ</sup><sup>i</sup> and <sup>σ</sup><sup>e</sup> denote intracellular and extracellular conductivities, whereas *C*<sup>m</sup> and *C*<sup>g</sup> represent the specific capacitance of the membrane and the intercalated disc, respectively. Furthermore, *n*e, *n*<sup>k</sup> <sup>i</sup> , and *<sup>n</sup>*k˜ <sup>i</sup> represent the outward pointing unit normal vectors of Ωe, Ω<sup>k</sup> <sup>i</sup> and <sup>Ω</sup>k˜ <sup>i</sup> , respectively. A homogeneous Dirichlet boundary condition is applied at the outer extracellular boundary in the *x*-direction (∂Ω<sup>D</sup> <sup>e</sup> ), and a homogeneous Neumann boundary condition is applied at the boundary in the y- and *<sup>z</sup>*-directions (∂Ω<sup>N</sup> <sup>e</sup> ). The parameters used in the computations below are summarized in Table 4.1. The properties of the cell membrane and the gap junctions are represented by *F*, *I*ion and *I*gap. In the computations reported below, we use the Grandi et al. model(9), to model the dynamics of the membrane (*F* and *I*ion), and for the gap junctions we use the simple passive model *I*<sup>k</sup> gap <sup>=</sup> <sup>w</sup><sup>k</sup> /*R*g.

## **4.2.1 Operator Splitting Applied to the EMI Model**

As mentioned above, a key step in solving the EMI model is to split the equations into parts that can be solved using standard tools. In (14), we derived a splitting scheme that leads to two key numerical challenges: Non-linear systems of ODEs to be solved

<sup>1</sup>Note that w<sup>k</sup> is defined specifically for each cell.

#### 4 Operator Splitting and FD Schemes 47


Table 4.1: Parameter values used in the simulations, based on (13). For parameters of the Grandi model, see (9).


at each computational node located at the cell membranes, and a series of elliptic equations; see Algorithm 1. All the differential equations involved in Algorithm 1 are of classical type and can be solved using well-established numerical methods. In our present implementation, we apply a straightforward finite difference scheme (see e.g., (29) for an elementary introduction to finite differences) for the elliptic equations and the Forward Euler method with a substepping time step Δ*t*ODE for solving the non-linear ODEs modeling the membrane dynamics (see (30)). However, it is worth observing that elliptic equations can as well be solved using the finite element method, or a finite volume method, thus allowing for more flexible and adaptive meshes.

# **4.3 Simulating the Effect of a Region of Ischemic Cells**

In order to demonstrate an application of Algorithm 1 above, we consider a collection of cells where a fraction of the cells are ischemic. This is known to perturb the electrical conduction and may lead to arrhythmias; see e.g., (33; 28; 19; 27). This problem has been carefully studied using homogenized models (mostly the monodomain model), but here we will show that the ischemic regions also have local effects when only very few cells are considered. In Figure 4.2, we consider a collection of cells organized in a two-dimensional mesh of 22×12 cells. The cells are modeled using the Grandi model with parameters as stated in (9). Within the domain, 8×6 of the center cells are ischemic in the sense that the extracellular potassium concentration surrounding these cells is increased from 5.4 mM to 10 mM. For the ischemic cells, we use the steady-state values of the state variables for the increased extracellular potassium concentration as initial conditions, and for the remaining cells, we use the steady-state values of the default Grandi model. In addition, we run the simulation for 5 ms before stimulation.

In the simulation results we observe that the ischemic region slows down conductions and thus perturbs the wave in the intracellular potential moving from left to right. This is consistent with the result obtained in (6) (cf. Figure 5), where the monodomain model was used. Such perturbations are known to be arrhythmogenic and have been observed several times in numerical experiments; see e.g., (33; 19; 15; 3; 6). Here, we observe that such perturbations can be initiated locally when only a few cells are subject to surroundings with elevated potassium concentration.

Fig. 4.2: Extracellular potential (left) and intracellular potential (right) at four different points in time in an EMI model simulation with an ischemic region in the center of the domain, marked by the purple rectangle. The parameter values used in the simulation are given in Table 4.1.

# **4.4 A Scalable Implementation of the Splitting Scheme**

In expectation of future simulations of excitable issues that may involve a huge number of cells, we see the need of a scalable implementation of the new splitting scheme, so that it can run efficiently on parallel computers. One specific criterion is that the computation time should grow linearly with the number of cells involved. Additionally, the design goals of this new code should also include independence of proprietary software (such as Matlab) and plug-and-play of the different numerical components. This section will present a preliminary version of such a scalable implementation.

## **4.4.1 The Linear System for the Intracellular Potential**

The main benefit of the new the splitting scheme is that the intracellular Laplace equations (one per cell) are decoupled from the extracellular Laplace equation, as stated in Algorithm 1. If we assume a constant intracellular conductivity σ<sup>i</sup> and that each cell is of the same shape and size, as shown in Figure 4.1, the matrices arising from a standard finite difference discretization of the intracellular Laplace equations for the individual cells will be mostly identical. There are only a small number of unique intracellular matrices, depending on whether there is a neighboring cell connected to each of the intercalated discs. It is thus unnecessary to compute an intracellular matrix for each cell. Instead, the cells that have the same neighbor connectivity situation can share the same intracellular matrix. This not only reduces the memory usage of an implementation, but also improves data reuse in the caches of a computer. Moreover, since the number of computational nodes per intracellular domain is relatively small (each intracellular domain has about 5300 degrees of freedom for the simulations used in this chapter), it is very efficient to use a direct solver each time an intracellular Laplace problem needs to be solved. Specifically, the LU factorization of each unique intracellular matrix *A*<sup>I</sup> can be pre-calculated, which renders the solution of *A*I*u*¯<sup>k</sup> <sup>i</sup> <sup>=</sup> *<sup>b</sup>*<sup>k</sup> <sup>i</sup> per cell to be merely invoking the forward-backward substitution procedure. Parallelism of the computation mainly arises from the fact that the intracellular Laplace equations can be solved independently of each other, while limited parallelism also exists within each forward-backward substitution.

## **4.4.2 The Linear System for the Extracellular Potential**

For the overall extracellular Laplace problem, which can be huge depending on the spatial resolution and the total number of cells, an iterative solver is more appropriate. Take for instance the case of 128×128 cells. The corresponding discrete extracellular Laplace equation has 107, 202, 214 degrees of freedom. Independent of the spatial resolution and the number of cells involved, the extracellular matrix *A*<sup>E</sup> arising from a standard finite difference discretization is symmetric and positive-definite (some care is needed to discretize the boundary conditions on the membranes). The resulting linear system *A*<sup>E</sup> *u*¯<sup>e</sup> = *b*<sup>e</sup> is thus a perfect candidate for the conjugate gradients (CG) method with an algebraic multigrid (AMG) preconditioner. Under optimal conditions, an AMG preconditioner requires a constant number of iterations to reach convergence independent of the linear system size, although the number of grid levels inside the AMG preconditioner may increase with the system size. Parallelism readily exists in iterative solvers, with several software libraries providing parallel implementations of CG and AMG.

## **4.4.3 The Non-Linear ODE System for the Membrane Potential**

For solving the non-linear ODE system per computational node on the membranes, a straightforward and often very efficient numerical strategy is the Forward Euler method with a substepping time step Δ*t*ODE. Since the non-linear ODE system on each membrane node is independent of the others, the ODE computation possesses the most ample parallelism.

## **4.4.4 The Implementation**

The Python programming language has been chosen for the implementation, mostly because of its flexibility for interfacing with numerical software libraries written in performance-friendly languages such as C and C++. We have used the ctypes module from the standard Python library for this purpose. The choice of Python also simplified a partial translation from the existing MATLAB code developed in (14).

We have chosen the SuperLU library (17) for performing the LU factorization of the intracellular matrices and the subsequent forward-backward substitution, via the bindings that are provided by SciPy (32). For the extracellular Laplace equation, we have used the ViennaCL library (20) for its implementation of CG and AMG. The CG iterations are by default configured to terminate when a tolerance of 10−<sup>5</sup> is reached. The AMG preconditioner has been configured to use the maximum independent set (MIS), see (4), as the coarsening algorithm and smoothed aggregation as the interpolation algorithm. For the ODE part, the Gotran automated code generator (10) has been used to translate the Grandi cell model into C code, callable from the Python side.

## **4.4.5 Parallelization**

The Python implementation currently relies on the adopted numerical libraries (SuperLU and ViennaCL) for an implicit parallelization of the PDE computation through multi-threading. This form of parallelization suits for shared-memory parallel computers, such as laptops or servers that use multicore CPUs. Multi-threading of the ODE computation is also enabled by inserting OpenMP compiler directives into the C code that is generated automatically by Gotran. The advantage of this implicit parallelization is that the user does not have to care about parallelization-specific coding. The downside is that all the computations have to run on a shared-memory system. It is possible to achieve the more general parallelization that targets distributed-memory parallel computers, which will be a task for future work.

## **4.4.6 Performance Results**

The simulations in this section were run on a dual-socket server with two 32-core AMD EPYC 7601 CPUs, each with 8-channel memory operating at 2666 MT/s. The number of OpenMP threads was set by default to the number of logical cores, equaling 128. Moreover, the environment variable OMP\_PROC\_BIND=TRUE was set to prevent the threads from migrating between the cores (which typically leads to unnecessary performance loss).

Table 4.2 shows the average solution time per time step for the 10 first time steps, where all the parameters are as prescribed in Table 4.1. The number of cardiac cells is doubled in the *<sup>x</sup>* and y directions for each row, and we observe that the time per cell remains fairly constant, indicating that the time to solution is a linear function of the number of cardiac cells simulated.


Table 4.2: Average solution time per time step for the E, M and I domains.

# **4.5 Software**

The Matlab code used to compute the solutions shown in Figure 4.2 and the Python code discussed in Section 4.4 can be found at https://github.com/KGHustad/ emi-book-2020-splitting-code. An archived version (12) is also available.

# **4.6 Conclusion**

In this chapter we have presented a numerical scheme for solving the EMI equations using operator splitting. The scheme allows for parallel solution of individual cells combined with a global solution of the equation modeling the extracellular potential. The latter is well suited for using optimal linear solvers such as AMG. The overall code scales linearly with the number cells and thus allows for simulation of a large number of cells. It remains to be seen how well this will work for very large numbers of cells; this is subject for further work.

**Acknowledgements** The research presented in this chapter has benefited from the Experimental Infrastructure for Exploration of Exascale Computing (eX3), which is financially supported by the Research Council of Norway under contract 270053.

**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.

# **References**


sights from a multi-scale mechanistic model of the human heart. PLoS computational biology 10(11):e1003891


# **Chapter 5**

# **Solving the EMI Equations using Finite Element Methods**

Miroslav Kuchta1, Kent-André Mardal1,<sup>2</sup> and Marie E. Rognes<sup>1</sup>

**Abstract** This chapter discusses 2 × 2 symmetric variational formulations and associated finite element methods for the EMI equations. We demonstrate that the presented methods converge at expected rates, and compare the approaches in terms of approximation of the transmembrane potential. Overall, the choice of which formulation to employ for solving EMI models becomes largely a matter of desired accuracy and available computational resources.

# **5.1 Introduction**

In this chapter, we present different weak formulations and corresponding finite element methods for solving the EMI equations as presented in (7, Chapter 1) over a physiological cell Ω<sup>i</sup> and its membrane Γ surrounded by an extracellular space Ω<sup>e</sup> and a time interval (0,*T*] for some time *T* > 0. This coupled, time-dependent, and typically nonlinear system of equations can be targeted numerically by an operator splitting scheme, see e.g (8, Chapter 4). Such an approach, combined with for instance an implicit Euler discretization in time, gives the following stationary and linear, but still coupled system of equations to be solved at each time-step: find the potentials *u*<sup>e</sup> = *u*e(*x*), *u*<sup>i</sup> = *u*i(*x*) (and current *I*<sup>m</sup> = *I*m(*x*)) such that

<sup>1</sup>Simula Research Laboratory, Norway

<sup>2</sup>Department of Mathematics, University of Oslo, Oslo, Norway

$$-\nabla \cdot (\sigma\_e \nabla u\_e) = 0 \qquad\qquad\qquad\text{in } \Omega\_e,\tag{5.1a}$$

$$-\nabla \cdot (\sigma\_i \nabla u\_i) = 0 \qquad\qquad\qquad\text{in } \Omega\_{\bar{t}}\qquad\qquad(5.1b)$$

$$
\sigma\_e \nabla u\_e \cdot \mathbf{n}\_e = -\sigma\_i \nabla u\_i \cdot \mathbf{n}\_i \equiv I\_m \qquad \qquad \text{on } \Gamma,\tag{5.1c}
$$

$$
\mu\_i - \mu\_e = \nu \tag{5.1d}
$$

$$\mathbf{v} - C\_m^{-1} \Delta t I\_m = f \tag{5.1e}$$

where Δ*t* > 0 denotes a time step size, and **n**<sup>e</sup> (resp. **n**i) denotes the outward pointing normal on Γ when viewed from Ω<sup>e</sup> (resp. Ωi). In our (implicit Euler) time discretization context, the known right-hand side *f* of (5.1e) combines the previous transmembrane potential solution, v0, and the evaluation of the ionic current, *<sup>I</sup>*ion, into *<sup>f</sup>* <sup>≡</sup> <sup>v</sup><sup>0</sup> <sup>−</sup> *<sup>C</sup>*−<sup>1</sup> <sup>m</sup> Δ*t I*ion.

We assume that the potential is grounded on part of the external boundary Γ<sup>D</sup> <sup>e</sup> and that the remaining external boundary Γ<sup>N</sup> <sup>e</sup> is insulated. These assumptions give the boundary conditions:

$$
\mu\_e = 0 \qquad\qquad\qquad\text{on }\Gamma^D\_{e^-},\tag{5.2a}
$$

$$
\sigma\_e \nabla \mu\_e \cdot \mathfrak{n}\_e = 0 \qquad\qquad\qquad\text{on } \Gamma^N\_{e^-}.\tag{5.2b}
$$

This geometrical setting is illustrated in Figure 5.1.

Fig. 5.1: (Left) Illustration of the geometric setting for the single cell EMI problem. (Right) Sample meshes for our benchmark problem (5.17). The boundary facets of the intracellular mesh Ωi,<sup>h</sup> form the membrane mesh Γh.

*Remark 5.1* We remark that a single cell model is here considered for simplicity. Indeed the formulations to be studied below can be similarly derived for the intercalated model (collections of connected cells). Formulations for a number of disconnected cells are then practically identical to the case considered here.

*Remark 5.2* If (5.2) is considered without any Dirichlet boundary data, i.e. |Γ<sup>D</sup> <sup>e</sup> | = 0, then only the transmembrane potential is fixed and the intracellular and extracellular potentials are determined only up to a single, common constant.

The EMI equations (5.1) set a rich scene for numerical exploration and can be solved in a multitude of ways. In this chapter, we will derive 2×2 different weak formulations (each defining a finite element method) of this system. The two first formulations (in Section 5.2) compute the intracellular and extracellular potentials as the main unknowns. These are referred to as *primal formulations*. The latter two formulations (in Section 5.3) additionally introduce the current densities **J**<sup>i</sup> = −σi∇*u*<sup>i</sup> and **J**<sup>e</sup> = −σe∇*u*<sup>e</sup> as independent unknowns. These are referred to as *mixed formulations*. We compare finite element discretizations of the primal and mixed formulations with respect to the approximation of the transmembrane potential v in Section 5.4. This choice is motivated by the observation that v is closely coupled to the membrane dynamics as discussed in Chapter 1.

## **5.1.1 Preliminaries: Function Spaces and Norms**

The EMI equations (5.1) define a multi-dimensional<sup>1</sup> PDE system coupling unknown fields defined over cellular domains and fields defined over the cell membrane, which can be viewed as a lower-dimensional manifold. Identifying the right function spaces for the different unknown fields is key to defining well-posed weak formulations of these equations. We here present suitable Sobolev spaces for this setting; the reader is referred to e.g. (3; 5) for more material and careful formalizations.

Let Ω be a bounded, polygonal domain in R<sup>d</sup> for *d* = 2, 3. We denote the space of square-integrable functions over Ω by *L*2(Ω), and let *H*1(Ω) be the Sobolev space of functions in *L*2(Ω) with weak derivatives in *L*2(Ω). The space **H**(div,Ω) contains vector-valued functions **v** : Ω → R<sup>d</sup> such that **v** ∈ **L**2(Ω) and ∇ · **v** ∈ *L*2(Ω). In general, when clear from the context, the domain will be omitted from the notation.

The *<sup>L</sup>*2-inner product and norm for *<sup>u</sup>*, v <sup>∈</sup> *<sup>L</sup>*2(Ω) is written as

$$(\mu, \nu)\_{0, \Omega} = \int\_{\Omega} \mu \nu \, \mathrm{d}x, \quad \|\nu\|\_{0, \Omega}^2 = \int\_{\Omega} \nu^2 \, \mathrm{d}x.$$

Similarly, we define the *<sup>H</sup>*1-norm as v <sup>2</sup> <sup>1</sup>,<sup>Ω</sup> <sup>=</sup> <sup>v</sup> <sup>2</sup> <sup>0</sup>,<sup>Ω</sup> <sup>+</sup> <sup>∇</sup><sup>v</sup> <sup>2</sup> <sup>0</sup>,<sup>Ω</sup> for <sup>v</sup> <sup>∈</sup> *<sup>H</sup>*1(Ω), and the **H**(div)-norm as **v** <sup>2</sup> div,<sup>Ω</sup> <sup>=</sup> **<sup>v</sup>** <sup>2</sup> <sup>0</sup>,<sup>Ω</sup> <sup>+</sup> ∇ · **<sup>v</sup>** <sup>2</sup> <sup>0</sup>,Ω.

<sup>1</sup> PDEs coupling fields over domains of different topological dimensions are often referred to as mixed-dimensional PDEs. To avoid the confusion-inducing term mixed-dimensional mixed in the subsequent sections, we instead use the term multi-dimensional in this chapter.

For Γ ⊆ ∂Ω, we define the constrained spaces *H*<sup>1</sup> <sup>Γ</sup>(Ω) <sup>=</sup> v <sup>∈</sup> *<sup>H</sup>*1(Ω) | v <sup>=</sup> 0 on <sup>Γ</sup> , and **H**Γ(div,Ω) = {**v** ∈ **H**(div,Ω) | **v** · **n** = 0 on Γ} where **n** is the (outward pointing) normal vector of Γ. Finally the spaces *H*1/2(Γ) and *H*−1/2(Γ) are the trace spaces of *H*<sup>1</sup> and **H**(div) respectively (6, Ch. 1., 2.). Here, the spaces will be considered with the norm defined in terms of fractional powers of the Helmholtz operator, see e.g. (4), i.e.

$$\|\|\mu\|\|\_{s}^{2} = (\mu, (-\Delta + I)^{s}\mu)\_{0,\Gamma}, \; \mu \in C^{\infty}(\Gamma).$$

We remark that in the following experiments the fractional norm is evaluated using the eigenvalue decomposition of −Δ + *I* as detailed in (11).

# **5.2 Primal Formulations**

We present two primal formulations of the stationary EMI system (5.1) with the boundary conditions given by (5.2): one *single-dimensional* formulation and one *multi-dimensional*formulation. The difference in the intra- and extracellular potential across the cell membrane Γ sets up a potential jump, the transmembrane potential v, c.f. (5.1d). Due to this jump, one cannot define a global, differentiable potential *<sup>u</sup>* <sup>∈</sup> *<sup>H</sup>*1(Ω<sup>i</sup> <sup>∪</sup> <sup>Ω</sup>e) such that *<sup>u</sup>*|Ωi <sup>=</sup> *<sup>u</sup>*<sup>i</sup> and *<sup>u</sup>*|Ωe <sup>=</sup> *<sup>u</sup>*e. Instead, we seek *<sup>u</sup>*<sup>i</sup> <sup>∈</sup> *<sup>H</sup>*1(Ωi) and *u*<sup>e</sup> ∈ *H*1(Ωe) separately. In the single-dimensional formulation, these are the only unknown fields, while in the multi-dimensional formulation, we keep *I*<sup>m</sup> as an additional unknown.

## **5.2.1 Single-Dimensional Primal Formulation**

Define the function spaces

$$V\_i = H^1(\Omega\_i), \quad V\_e = H^1\_{\Gamma\_e^D\_e}(\Omega\_e). \tag{5.3}$$

To derive a weak formulation of (5.1), multiply (5.1a) by a test function <sup>v</sup><sup>e</sup> <sup>∈</sup> *<sup>V</sup>*e, (5.1b) by another test function <sup>v</sup><sup>i</sup> <sup>∈</sup> *<sup>V</sup>*i, and integrate the divergence by parts. This yields the variational formulation: find *u*<sup>e</sup> ∈ *V*e, *u*<sup>i</sup> ∈ *V*<sup>i</sup> satisfying

$$
\int\_{\Omega\_{\rm ef}} \sigma\_{\rm e} \nabla \mu\_{\rm e} \cdot \nabla \upsilon\_{\rm e} \, \mathrm{d}x - \int\_{\Gamma} \sigma\_{\rm e} \nabla \mu\_{\rm e} \cdot \mathbf{n}\_{\rm e} \upsilon\_{\rm e} \, \mathrm{d}s = 0,\tag{5.4a}
$$

e

$$\int\_{\Omega\_{\bar{l}}} \sigma\_{\bar{l}} \nabla u\_{\bar{l}} \cdot \nabla v\_{\bar{l}} \, \mathrm{d}x + \int\_{\Gamma} \left( -\sigma\_{\bar{l}} \nabla u\_{\bar{l}} \cdot \mathbf{n}\_{\bar{l}} \right) v\_{\bar{l}} \, \mathrm{d}s = 0. \tag{5.4b}$$

for all <sup>v</sup><sup>e</sup> <sup>∈</sup> *<sup>V</sup>*e, <sup>v</sup><sup>i</sup> <sup>∈</sup> *<sup>V</sup>*i. In the bracketed term of (5.4b), we recognize the membrane current density *I*<sup>m</sup> as defined by (5.1c), and similarly, the interface contribution in the corresponding extracellular equation (5.4a) hides −*I*m. Combining (5.1e) and (5.1d), we find that

$$I\_m = C\_m(\Delta t)^{-1} \left( (\mu\_i - \mu\_e) - f \right). \tag{5.5}$$

After substituting (5.5) into (5.4), the single-dimensional primal weak form of (5.1) reads: find *u*<sup>i</sup> ∈ *V*<sup>i</sup> and *u*<sup>e</sup> ∈ *V*<sup>e</sup> such that

$$\begin{aligned} \int\_{\Omega\_{\varepsilon}} \sigma\_{e} \nabla u\_{e} \cdot \nabla v\_{e} \, \mathrm{d}x + \int\_{\Gamma} C\_{m}(\Delta t)^{-1} u\_{e} v\_{e} \, \mathrm{d}s - \int\_{\Gamma} C\_{m}(\Delta t)^{-1} u\_{i} v\_{e} \, \mathrm{d}s &= \\ & - C\_{m}(\Delta t)^{-1} \int\_{\Gamma} f v\_{e} \, \mathrm{d}s, \\ \int\_{\Omega\_{i}} \sigma\_{i} \nabla u\_{i} \cdot \nabla v\_{i} \, \mathrm{d}x + \int\_{\Gamma} C\_{m}(\Delta t)^{-1} u\_{i} v\_{i} \, \mathrm{d}s - \int\_{\Gamma} C\_{m}(\Delta t)^{-1} u\_{e} v\_{i} \, \mathrm{d}s &= \\ & C\_{m}(\Delta t)^{-1} \int\_{\Gamma} f v\_{i} \, \mathrm{d}s, \end{aligned} (5.6)$$

for all <sup>v</sup><sup>e</sup> <sup>∈</sup> *<sup>V</sup>*<sup>e</sup> and <sup>v</sup><sup>i</sup> <sup>∈</sup> *<sup>V</sup>*i.

We remark that (5.6) can be viewed as a coupling of two Poisson problems with a Robin boundary condition on Γ. The well-posedness of the problem is then discussed in (10, Chapter 6). Finally, the transmembrane potential can here be computed from its definition (5.1d) as a difference of the computed potentials.

## **5.2.2 Multi-Dimensional Primal Formulation**

An alternative formulation can be derived by keeping *I*<sup>m</sup> as a separate unknown field. Since Γ is of a different (lower) dimension than Ωi,Ωe; and as *I*<sup>m</sup> : Γ → R while *u*<sup>i</sup> : Ω<sup>i</sup> → R, *u*<sup>e</sup> : Ω<sup>e</sup> → R, we will refer to this as a multi-dimensional primal formulation. Observe that (5.4) now yields two equations for three unknowns *u*<sup>i</sup> ∈ *V*i, *u*<sup>e</sup> ∈ *V*e, and *I*<sup>m</sup> ∈ *Q*:

$$\begin{aligned} \int\_{\Omega\_{\boldsymbol{e}}} \boldsymbol{\sigma}\_{\boldsymbol{e}} \nabla \boldsymbol{u}\_{\boldsymbol{e}} \cdot \nabla \boldsymbol{v}\_{\boldsymbol{e}} \, \mathrm{d}\boldsymbol{x} - \int\_{\Gamma} I\_{m} \boldsymbol{v}\_{\boldsymbol{e}} \, \mathrm{d}\boldsymbol{s} &= \boldsymbol{0}, \quad \forall \, \boldsymbol{v}\_{\boldsymbol{e}} \in V\_{\boldsymbol{e}}, \\ \int\_{\Omega\_{\boldsymbol{i}}} \boldsymbol{\sigma}\_{\boldsymbol{i}} \nabla \boldsymbol{u}\_{\boldsymbol{i}} \cdot \nabla \boldsymbol{v}\_{\boldsymbol{i}} \, \mathrm{d}\boldsymbol{x} + \int\_{\Gamma} I\_{m} \boldsymbol{v}\_{\boldsymbol{i}} \, \mathrm{d}\boldsymbol{s} &= \boldsymbol{0}, \quad \forall \, \boldsymbol{v}\_{\boldsymbol{i}} \in V\_{\boldsymbol{i}}. \end{aligned}$$

However, the missing equation can be obtained from (5.5). Let

$$\mathcal{Q} = H^{1/2}(\Gamma), \quad \mathcal{Q}^\* = H^{-1/2}(\Gamma). \tag{5.7}$$

We remind the reader that if Γ is a co-dimensional 1 subset of Ω then trace operations from Ω to Γ, *Tu* = *u*|Γ, *u* ∈ *C*(Ω) and *T*nτ = τ|<sup>Γ</sup> · **n**, τ ∈ **C**(Ω), formally have the following mapping properties *T* : *H*1(Ω) → *H*1/2(Γ) and *T*<sup>n</sup> : **H**(div,Ω) →

*H*−1/2(Γ). Hence, let *j*<sup>m</sup> be a a test function from *Q*∗. We shall then enforce that (5.5) holds in the weak sense:

$$\int\_{\Gamma} (\boldsymbol{\mu}\_{i} - \boldsymbol{\mu}\_{e}) j\_{m} \, \mathrm{d}s - \int\_{\Gamma} \Delta t \boldsymbol{C}\_{m}^{-1} I\_{m} j\_{m} \, \mathrm{d}s = \int\_{\Gamma} f j\_{m} \, \mathrm{d}s, \quad \forall \, j\_{m} \in \mathcal{Q}^{\*}.$$

In turn, the multi-dimensional primal formulation of (5.1) reads: find *u*<sup>i</sup> ∈ *V*i, *u*<sup>e</sup> ∈ *V*e, *I*<sup>m</sup> ∈ *Q*<sup>∗</sup> such that

$$\begin{aligned} \int\_{\Omega\_e} \sigma\_e \nabla u\_e \cdot \nabla v\_e \, \mathrm{d}x - \int\_{\Gamma} v\_e I\_m \, \mathrm{d}s &= 0, \\ \int\_{\Omega\_i} \sigma\_i \nabla u\_i \cdot \nabla v\_i \, \mathrm{d}x + \int\_{\Gamma} v\_i I\_m \, \mathrm{d}s &= 0, \\ \int\_{\Gamma} -u\_e j\_m \, \mathrm{d}s + \int\_{\Gamma} \mu\_i j\_m \, \mathrm{d}s - \int\_{\Gamma} \Delta t C\_m^{-1} I\_m j\_m \, \mathrm{d}s &= \int\_{\Gamma} f j\_m \, \mathrm{d}s, \end{aligned} \tag{5.8}$$

for all <sup>v</sup><sup>i</sup> <sup>∈</sup> *<sup>V</sup>*i, <sup>v</sup><sup>e</sup> <sup>∈</sup> *<sup>V</sup>*<sup>e</sup> and *<sup>j</sup>*<sup>m</sup> <sup>∈</sup> *<sup>Q</sup>*∗. We remark that (5.8) is closely related to the Babuška problem for enforcing boundary conditions by Lagrange multipliers (1) and the mortar finite element method, see e.g. (13). With regards to evaluation of the transmembrane potential, we note that v can be post-computed in several ways: from (5.1d) (as for the single-dimensional primal formulation (5.6)) or from *I*<sup>m</sup> and (5.1e).

# **5.3 Mixed Formulations**

We now turn to consider *mixed formulations* of the EMI system (5.1). Let us (re)introduce the current densities

$$\mathbf{J}\_{i} = -\sigma\_{i}\nabla\mu\_{i}, \quad \mathbf{J}\_{e} = -\sigma\_{e}\nabla\mu\_{e} \tag{5.9}$$

and the global field **<sup>J</sup>** on <sup>Ω</sup> <sup>=</sup> <sup>Ω</sup><sup>i</sup> <sup>∪</sup>Ω<sup>e</sup> such that **<sup>J</sup>**|Ωi <sup>=</sup> **<sup>J</sup>**<sup>i</sup> and **<sup>J</sup>**|Ωe <sup>=</sup> **<sup>J</sup>**e. In general, we use the convention that for a scalar or vector field *u* defined on Ω, the restriction on Ω<sup>i</sup> and Ω<sup>e</sup> is denoted by *u*<sup>i</sup> and *u*e, respectively.

With these definitions, (5.1a)–(5.1c) become: find the current densities **J**i, **J**<sup>e</sup> (or **J**) and the potentials *u*i,*u*<sup>e</sup> (or *u*) satisfying

$$
\sigma\_e^{-1} J\_e + \nabla u\_e = 0 \qquad\qquad\qquad\text{on } \Omega\_e,\tag{5.10a}
$$

$$
\sigma\_i^{-1} \mathbf{J}\_i + \nabla u\_i = \mathbf{0} \tag{5.10b}
$$

−∇ · **J** = 0 in Ω, (5.10c)

$$J\_e \cdot \mathfrak{n}\_e + J\_i \cdot \mathfrak{n}\_i = 0 \qquad\qquad\qquad\text{on }\Gamma.\tag{5.10d}$$

We refer to (5.10) together with (5.1d)–(5.1e) as the mixed EMI system with boundary conditions given by (5.2). Weak formulations of the mixed form can enjoy improved conservation properties and stability properties (3). In particular, approximations of **J** may be computed such that they are exactly divergence free, cf. (5.10c).

The continuity condition (5.10d) ensures that the normal component of **J** is continuous on Γ. We remark that **v** ∈ **H**(div,Ω) implies continuity of **v** · **n** on Γ. Moreover, we observe that (5.10c) involves only divergence of the field **J**. It is therefore sufficient to seek **<sup>J</sup>** in *<sup>S</sup>* <sup>=</sup> **<sup>H</sup>**ΓN e (div,Ω). Note that in contrast to the primal formulation, here the Neumann boundary condition (5.2b) is enforced as an essential condition; that is, it is included in the construction of the function space *S*.

## **5.3.1 Single-Dimensional Mixed Formulation**

Let

$$S = \left\{ \mathbf{J} \in H\_{\Gamma\_e^N}(\text{div}, \Omega); \mathbf{J} \cdot \mathfrak{n} \in L^2(\Gamma) \right\}, \quad V = L^2(\Omega). \tag{5.11}$$

To derive a weak form of the mixed EMI system, consider a test function τ ∈ *S*. Taking the dot product of (5.10a), (5.10b) with τi, τe, integrating and applying integration by parts then yields

$$\begin{split} \int\_{\Omega\_{\varepsilon}} \boldsymbol{\sigma}\_{e}^{-1} \mathbf{J}\_{e} \cdot \boldsymbol{\tau}\_{e} \, \mathrm{d}\mathbf{x} - \int\_{\Omega\_{\varepsilon}} \boldsymbol{u}\_{e} \nabla \cdot \boldsymbol{\tau}\_{e} \, \mathrm{d}\mathbf{x} + \int\_{\Gamma} \boldsymbol{u}\_{e} \boldsymbol{\tau}\_{e} \cdot \mathbf{n}\_{e} \, \mathrm{d}\mathbf{s} &= - \int\_{\Gamma\_{e}^{D}} \boldsymbol{u}\_{e} \boldsymbol{\tau}\_{e} \cdot \mathbf{n}\_{e} \, \mathrm{d}\mathbf{s}, \\ \int\_{\Omega\_{i}} \boldsymbol{\sigma}\_{i}^{-1} \mathbf{J}\_{i} \cdot \boldsymbol{\tau}\_{i} \, \mathrm{d}\mathbf{x} - \int\_{\Omega\_{i}} \boldsymbol{u}\_{i} \nabla \cdot \boldsymbol{\tau}\_{i} \, \mathrm{d}\mathbf{x} + \int\_{\Gamma} \boldsymbol{u}\_{i} \boldsymbol{\tau}\_{i} \cdot \mathbf{n}\_{i} \, \mathrm{d}\mathbf{s} &= \mathbf{0}. \end{split}$$

Observe that by continuity of the normal component of the test function (τ<sup>i</sup> ·**n** = τ<sup>e</sup> ·**n** on Γ), and the identity **n**<sup>e</sup> = −**n**i, the integrals on Γ can be added, resulting in ∫ <sup>Γ</sup>(*u*<sup>i</sup> <sup>−</sup> *<sup>u</sup>*e)<sup>τ</sup> · **<sup>n</sup>**i. Moreover, using (5.5), the membrane term can be rewritten as ∫ Γ *C*−1 <sup>m</sup> Δ*t***J** · **n**<sup>i</sup> + *f* τ · **n**i. In turn, we arrive at the variational problem: find **J** ∈ *S*, *u* ∈ *V* such that

$$\begin{split} \int\_{\Omega} \boldsymbol{\sigma}^{-1} \mathbf{J} \cdot \boldsymbol{\tau} \, \mathrm{d}\mathbf{x} + \int\_{\Gamma} \mathrm{C}\_{m}^{-1} \Delta t \mathbf{J} \cdot \mathbf{n}\_{i} \boldsymbol{\tau} \cdot \mathbf{n}\_{i} \, \mathrm{d}\mathbf{s} - \int\_{\Omega} \boldsymbol{u} \nabla \cdot \boldsymbol{\tau} \, \mathrm{d}\mathbf{x} &= - \int\_{\Gamma} f \boldsymbol{\tau} \cdot \mathbf{n}\_{i} \, \mathrm{d}\mathbf{s}, \\ - \int\_{\Omega} q \nabla \cdot \mathbf{J} \, \mathrm{d}\mathbf{x} &= 0, \end{split} \tag{5.12}$$

for all <sup>τ</sup> <sup>∈</sup> *<sup>S</sup>* and *<sup>q</sup>* <sup>∈</sup> *<sup>V</sup>*, with <sup>σ</sup> defined naturally as <sup>σ</sup>|Ωi <sup>=</sup> <sup>σ</sup><sup>i</sup> and likewise for <sup>Ω</sup>e. Note that due to the extra trace regularity of the trial/test space *S* all the terms in (5.12), and in particular the interface term <sup>∫</sup> <sup>Γ</sup> *<sup>C</sup>*−<sup>1</sup> <sup>m</sup> Δ*t***J** · **n**iτ · **n**<sup>i</sup> d*s*, are well defined. Without the extra regularity, i.e. if *<sup>S</sup>* <sup>=</sup> **<sup>H</sup>**ΓN e (div,Ω), this would not be the case.

We remark that (5.12) is a Γ-perturbed mixed formulation of the Poisson problem (see e.g. (3; 9) for more details on mixed formulations of the Poisson problem). Considering the task of approximating the transmembrane potential, we observe that v can be computed in two ways, as for the multi-dimensional primal formulation. Indeed, in addition to the identity <sup>v</sup> <sup>=</sup> *<sup>u</sup>*<sup>i</sup> <sup>−</sup> *<sup>u</sup>*e, cf. (5.1d), equation (5.1e) can be used since *I*<sup>m</sup> = **J** · **n**<sup>i</sup> is readily available.

## **5.3.2 Multi-Dimensional Mixed Formulation**

As for the primal formulations, the multi-dimensional mixed formulation is obtained by keeping the interface term as an explicit unknown field. Let

$$S = H\_{\Gamma\_{\mathfrak{s}}^N}(\text{div}, \Omega), \quad V = L^2(\Omega), \quad W = H^{1/2}(\Gamma). \tag{5.13}$$

To complete the formulation, the equation to be enforced weakly by test functions w <sup>∈</sup> *<sup>W</sup>* is the membrane dynamics condition (5.1e) written in the form

$$J \cdot \mathbf{n}\_i - C\_m(\Delta t)^{-1} \nu = -C\_m(\Delta t)^{-1} f \qquad \text{on } \Gamma.$$

The final multi-dimensional mixed weak formulation then reads: Find the current densities **<sup>J</sup>** <sup>∈</sup> *<sup>S</sup>*, potentials *<sup>u</sup>* <sup>∈</sup> *<sup>V</sup>*, and transmembrane potential v <sup>∈</sup> *<sup>W</sup>* such that

$$\begin{aligned} \int\_{\Omega} \boldsymbol{\sigma}^{-1} \mathbf{J} \cdot \boldsymbol{\tau} \, \mathrm{d} \mathbf{x} - \int\_{\Omega} \boldsymbol{u} \nabla \cdot \boldsymbol{\tau} \, \mathrm{d} \mathbf{x} + \int\_{\Gamma} \boldsymbol{v} \boldsymbol{\tau} \cdot \mathbf{n}\_{i} \, \mathrm{d} \mathbf{s} &= \mathbf{0}, \\ - \int\_{\Omega} q \nabla \cdot \mathbf{J} \, \mathrm{d} \mathbf{x} &= \mathbf{0}, \\ \int\_{\Gamma} \boldsymbol{w} \mathbf{J} \cdot \mathbf{n}\_{i} \, \mathrm{d} \mathbf{s} - \int\_{\Gamma} \boldsymbol{C}\_{m} (\boldsymbol{\Delta t})^{-1} \boldsymbol{v} \boldsymbol{w} \, \mathrm{d} \mathbf{s} &= -\mathbf{C}\_{m} (\boldsymbol{\Delta t})^{-1} \int\_{\Gamma} \boldsymbol{f} \boldsymbol{w} \, \mathrm{d} \mathbf{s}, \end{aligned} \tag{5.14}$$

for all <sup>τ</sup> <sup>∈</sup> *<sup>S</sup>*, *<sup>q</sup>* <sup>∈</sup> *<sup>V</sup>* and w <sup>∈</sup> *<sup>W</sup>*. Note that (5.14) is defined on the standard **<sup>H</sup>**(div) space, cf. (5.12), as the formulation no longer contains the troublesome interface term ∫ <sup>Γ</sup> *<sup>C</sup>*−<sup>1</sup> <sup>m</sup> <sup>Δ</sup>*t***<sup>J</sup>** · **<sup>n</sup>**i<sup>τ</sup> · **<sup>n</sup>**<sup>i</sup> <sup>d</sup>*s*. With regards to the approximation of <sup>v</sup> in formulation (5.14), observe that no post-processing is required to obtain this quantity. This is contrast to the previous three formulations. We remark that (5.14) is closely connected to the Babuška problem for the mixed Poisson equation (2).

# **5.4 Finite Element Spaces and Methods**

To solve the primal and mixed, single- and multi-dimensional weak formulations numerically, we approximate the continuous function spaces by discrete finite element spaces. Each choice of formulation and finite element space defines a finite element method for solving the EMI system.

Now, let Ω<sup>h</sup> be a mesh of the domain Ω = Ω<sup>i</sup> ∪ Ω<sup>e</sup> with characteristic mesh size *h*, which conforms to Γ in the sense that no element of Ω<sup>h</sup> has its interior intersected by Γ. The meshes Ωe,<sup>h</sup> and Ωi,<sup>h</sup> of the extracellular and intracellular domains are formed as non-overlapping subsets of the cells of Ωh. As a consequence, the mesh Γ<sup>h</sup> of Γ is formed by the facets of elements of Ωh, cf. Figure 5.1. We remark that the single-dimensional primal formulation allows for independent discretizations of Ωi, Ω<sup>e</sup> as well as Γ.

The choice of the finite element spaces plays a crucial role for the stability of the different discrete formulations. In particular, for the saddle-point systems, the spaces must be compatible in the sense of Babuška-Brezzi and satisfy discrete inf-sup conditions, see e.g. (3). For the primal formulations (5.4) and (5.8), we seek discrete unknowns and test functions in

$$V\_{i,h} = P\_1(\Omega\_{i,h}) \subset V\_i, \quad V\_{e,h} = P\_1(\Omega\_{e,h}) \subset V\_e, \quad \underline{Q}\_h = P\_1(\Gamma\_h) \subset \underline{Q}, \tag{5.15}$$

where *P*<sup>1</sup> denotes the space of continuous piecewise linears (defined relative to the relevant mesh). With these spaces, we expect linear convergence with the mesh size *h* for all variables in *H*1-norms and quadratic convergence in the *L*2-norm.

For the mixed formulations (5.10) and (5.12), we seek discrete unknowns and test functions in

$$S\_h = RT\_0(\Omega\_h) \subset S, \quad V\_h = P\_0(\Omega\_h) \subset V, \quad W\_h = P\_0(\Gamma\_h). \tag{5.16}$$

Here *RT*<sup>0</sup> denotes the lowest order Raviart-Thomas finite element spaces and *P*<sup>0</sup> denotes the space of piecewise constants defined relative to the relevant mesh. These spaces satisfy the relevant stability conditions, and we expect linear convergence of all unknown fields in their respective natural norms.

# **5.5 Numerical Results**

## **5.5.1 Comparison of Convergence between Formulations**

In order to compare the properties of the different formulations, and in particular their numerical stability and accuracy, we consider a manufactured solution test case with a smooth analytical solution. We define Ω = [0, 1] <sup>2</sup> and Ω<sup>i</sup> = [0.25, 0.75] <sup>2</sup> with |Γ<sup>N</sup> <sup>e</sup> <sup>|</sup> <sup>=</sup> 0. For simplicity, let <sup>σ</sup><sup>i</sup> <sup>=</sup> 1, <sup>σ</sup><sup>e</sup> <sup>=</sup> 2, *<sup>C</sup>*<sup>m</sup> <sup>=</sup> 1, <sup>Δ</sup>*<sup>t</sup>* <sup>∈</sup> 1, 10−<sup>4</sup> and consider the exact solution

$$\begin{aligned} u\_e &= \sin\left(\pi(\mathbf{x} + \mathbf{y})\right), \\ u\_i &= \frac{u\_e}{\sigma\_i} + \cos\left(\pi(\mathbf{x} - \frac{1}{4})(\mathbf{x} - \frac{3}{4})\right) \cos\left(\pi(\mathbf{y} - \frac{1}{4})(\mathbf{y} - \frac{3}{4})\right), \end{aligned} \tag{5.17}$$

which corresponds to (Δ*t* dependent) right hand sides *f* = *u*<sup>i</sup> − *u*<sup>e</sup> − Δ*t I*m. Note, that with (5.17) both v - 0 and *I*<sup>m</sup> - 0. We discretize the domain by a uniform mesh by dividing the unit square into *n* × *n* squares and dividing each square by the (left) diagonal into isosceles triangles of size *h*, cf. Figure 5.1. To compare the dimensionality of the different formulations, Table 5.1 lists the dimensions of the four different finite element pairings over these meshes.


Table 5.1: Dimensions of the different finite element spaces for uniform refinements of the unit square. The first row corresponds to a mesh of Ω with *n* = 16, i.e. having 2 × 16 × 16 cells.

Fig. 5.2: Convergence properties of (left) primal formulations (5.6)–(5.8) and (right) mixed formulations (5.12)–(5.14). The EMI system (5.1) is solved with the exact solution given by (5.17) and Δ*t* = 1. Filled symbols correspond to single-dimensional formulations. The number associated with each line indicates the convergence rate obtained from a least squares fit of the corresponding data.

On a series of meshes and for a set of different timesteps, we compute <sup>2</sup> the different approximations of all four finite element methods. We then evaluate the approximation error by evaluating the difference between (higher order interpolants of) the exact solution and the approximations in appropriate norms. Figure 5.2 shows the errors of the different formulations (with Δ*t* = 1). In the primal formulations, *u* refers to the global potential, i.e. *<sup>u</sup>*<sup>i</sup> <sup>=</sup> *<sup>u</sup>*|Ωi , *<sup>u</sup>*<sup>e</sup> <sup>=</sup> *<sup>u</sup>*|Ωe . As the formulations seek the approximations *u*<sup>i</sup> ∈ *H*1(Ωi), *u*<sup>e</sup> ∈ *H*1(Ωe) the error is considered in the natural (broken) norm *u* <sup>1</sup> = ( *u*<sup>i</sup> <sup>2</sup> <sup>1</sup> <sup>+</sup> *<sup>u</sup>*<sup>e</sup> <sup>2</sup> 1 ) <sup>1</sup>/2. We observe that all the quantities converge linearly in their respective natural norms, as expected. In particular, the errors in the current density *<sup>I</sup>*<sup>m</sup> in (5.8) and the transmembrane potential <sup>v</sup> in (5.14) are reported in the fractional norms *H*−1/<sup>2</sup> and *H*1/2, respectively. The former is computed by first interpolating the error into the space of continuous piecewise cubic polynomials on <sup>Γ</sup><sup>h</sup> while for <sup>v</sup> <sup>−</sup> <sup>v</sup><sup>h</sup> the *<sup>P</sup>*<sup>1</sup> element is used for error interpolation. Without this higher-order approximation of the error, i.e. if the error is computed in the same space as the discrete solution, we observe quadratic convergence.

Finally, we note that the primal formulations yield identical approximations of *u* cf. Figure 5.2 (left). Similarly, the mixed formulations give identical approximations of (*u*, **J**) cf. Figure 5.2 (right). Considering for comparison the error in the potential in the *L*2-norm, it can be seen that the primal formulations are more accurate than the mixed formulations. The same experiments for Δ*t* = 10−<sup>4</sup> give similar approximation results. However, it is not true that these conclusions hold in the limit of Δ*t* approaching 0, see e.g. Chapter 6.

## **5.5.2 Post-Processing the Transmembrane Potential**

With the exception of the multi-dimensional mixed formulation (5.14), the transmembrane potential v in the remaining EMI formulations is computed by postprocessing. In (5.6) the approximation <sup>v</sup><sup>h</sup> can be obtained by interpolating the difference *u*i,<sup>h</sup> − *u*e,<sup>h</sup> onto e.g. the space of continuous piecewise linear functions over Γh. This procedure can, of course, be used in the other formulations as well. However (5.8) and (5.12) also offer an alternative approach. In the multi-dimensional primal formulation, the discrete membrane current density, *I*m,<sup>h</sup> is computed in the space *<sup>P</sup>*1(Γh) of continuous piecewise linear functions on <sup>Γ</sup>h. In turn, <sup>v</sup><sup>h</sup> can be computed (in the same space) as a projection of *C*−<sup>1</sup> <sup>m</sup> Δ*t I*m,<sup>h</sup> + *f* . The same formula can be applied in the single-dimensional mixed formulation since the current density can be evaluated as **<sup>J</sup>**<sup>h</sup> · **<sup>n</sup>**i. We recall that in (5.12) the natural space for <sup>v</sup><sup>h</sup> is the space of (discontinuous) piecewise constant functions on Γ<sup>h</sup> however.

Convergence of the transmembrane potential obtained by the different formulations and post-processing strategies is shown in Figure 5.3 for the same test case as previ-

<sup>2</sup> The code used to produce results in this chapter is available at https://github.com/MiroK/ emi-book-fem and archived at (12).

ously. We observe that the primal formulations yield quadratic convergence and that the single-dimensional primal (5.6), and multi-dimensional primal formulation (5.8) are practically identical. The discrete potentials obtained by solving the mixed formulations converge linearly with the projection method in the single-dimensional mixed formulation yielding the most accurate <sup>v</sup>h. In particular, the approximation is better than that of the multi-dimensional mixed formulation for this test case. Computing the potential in the single-dimensional mixed formulation (5.12) by interpolating *u*i,<sup>h</sup> − *u*e,<sup>h</sup> leads to poorer approximation than for the multi-dimensional mixed formulation. By comparing the results for two different time steps, we observe that the rates do not change considerably if Δ*t* is modified.

Fig. 5.3: Approximation of the transmembrane potential by the different EMI formulations. (Left) Δ*t* = 1, (right) Δ*t* = 10<sup>−</sup>4. Postprocessing by projection (using the current density) is indicated by *<sup>I</sup>*m. In multi-dimensional mixed formulation <sup>v</sup><sup>h</sup> is obtained by solving (5.14). Interpolation of *u*i,<sup>h</sup> − *u*e,<sup>h</sup> is used in other formulations. The final number indicates the convergence rate.

# **5.6 Conclusions and Outlook**

All four finite element formulations provide a converging approximation to the stationary problem (5.1). This system is a key building block in any operator splitting algorithm for the time-dependent EMI equations (1.30). The formulations provide solutions which differ by accuracy as well as computational cost, cf. Figures 5.2–5.3 and Table 5.1. The formulations also differ in robustness of their approximation properties with respect to Δ*t*. This issue is however beyond the scope of this chapter, and the interested reader is referred to the discussion in Chapter 6. In terms of coupling with the membrane dynamics the single/multi-dimensional primal and single-dimensional mixed formulations require post-processing. However, all approaches discussed in Section 5.5.2 are easy to implement. Therefore, the choice of which formulation to employ in solving the EMI model is largely a matter of desired accuracy and available computational resources.

**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.

# **References**


# **Chapter 6 Iterative Solvers for EMI Models**

Miroslav Kuchta<sup>1</sup> and Kent-André Mardal1,<sup>2</sup>

**Abstract** This chapter deals with iterative solution algorithms for the four EMI formulations derived in (17, Chapter 5). Order optimal monolithic solvers robust with respect to material parameters, the number of degrees of freedom of discretization as well as the time-stepping parameter are presented and compared in terms of computational cost. Domain decomposition solver for the single-dimensional primal formulation is discussed.

# **6.1 Introduction**

Spatial discretization of EMI models describing a few cells with a complex/realistic geometry or a large collection of cells leads to linear systems with considerable number of unknowns. In our largest simulations, we may have linear systems involving billions of unknowns at millions of time steps. It is the purpose of this chapter to discuss how such systems can be solved efficiently with available algorithms.

Let us denote the system size as *N* and let *h* be a typical grid/mesh size. The complexity of a solution algorithm can then be analyzed in terms of how the computational time grows with *N*. For instance, a naive Gaussian elimination would perhaps scale as O(*N*3) and for linear system involving billions of unknowns such an approach is therefore not feasible during a life-time. Here, we shall aim for algorithms that are *order optimal*, i.e. their solution time scales linearly, O(*N*), with respect to the number of unknowns. Clearly, O(*N*) is order optimal in the sense that this is also the complexity of writing the results to file. In general, direct solvers, like for instance LU, do not scale linearly in *N*. Here, we will use Krylov solvers like the Conjugate

<sup>1</sup>Simula Research Laboratory, Norway

<sup>2</sup>Department of Mathematics, University of Oslo, Oslo, Norway

Gradient (CG) or the Minimal Residual (MinRes) methods. Both methods are known to provide efficient computations when combined with proper preconditioning techniques. Preconditioners based on multigrid and/or domain decomposition methods have been shown to be order-optimal for a variety of problems, but the theory for the EMI problem is currently limited. Efficient solvers for the EMI-problem are discussed in the following.

We say that a solver for the *transient* EMI model (1.30) is order optimal if the solution time grows linearly in Δ*t* <sup>−</sup>1, where Δ*t* is the time step. Considering the *stationary* problem (5.1), which is to be solved at every step of the temporal loop, it is clear that order optimality of the transient solver requires that the solver of the linear system does not degenerate for small (or large) Δ*t*, i.e. that it is robust with respect to the time stepping.

In this chapter we discuss solution algorithms for the linear systems due to the finite element discretization of (5.1) which are robust in *h* and well as Δ*t*. Two types of approaches will be considered. Monolithic approaches, where all the unknowns are solved for at once, are the subject of Section 6.2. Section 6.3 then concerns a domain decomposition approach where one iterates between the intra/extra-cellular subproblems. The solvers will be compared in terms of robustness and cost, however, only serial performance will be addressed. We remark that parallel scalable solvers suitable for the mixed formulations of the EMI models (5.12) (5.14) are the balancing domain decomposition methods, see e.g. (24; 25). For the elliptic single-dimensional primal formulation (5.6), the FETI domain decomposition methods, e.g. (13; 18) could be used. Moreover, to simplify the exposition the focus shall be on a single cell model. We remark that all the algorithms presented further can be generalized in a rather straightforward manner to models with multiple *disconnected* cells. However, construction of robust monolithic solvers for *multi-dimensional* formulations for collections of *connected* cells is out of the scope of this manuscript.

# **6.2 Monolithic Solvers**

Let A<sup>h</sup> *x*<sup>h</sup> = *L*<sup>h</sup> be a linear system due to discretization of the continuous problem A*x* = *L* in *W* , where *W* is a dual space to some Hilbert space *W* and A : *W* → *W* . Note that in case of (5.1) the continuous operator A depends on material parameters σi, σe, *C*<sup>m</sup> as well as the time step size Δ*t*. Here the focus is placed on the latter dependence.

The monolithic solvers considered here are preconditioned Krylov methods where the preconditioner B : *W* → *W* shall be constructed such that the number of iterations required for solving problems BhA<sup>h</sup> *x*<sup>h</sup> = B<sup>h</sup> *L*<sup>h</sup> is bounded in *h* and Δ*t*. A constructive framework for establishing B is operator preconditioning (19), where the structure of the preconditioner reflects mapping properties of A as an isomorphism between suitably chosen Hilbert space and its dual space. In particular, if A : *W* → *W* is an isomorphism then a Riesz map with respect to the inner product of *W* is a suitable preconditioner. As explained in (19), multilevel or domain decomposition algorithms are equivalent to Riesz maps for standard elliptic problems. Furthermore, these schemes also work in a fractional setting (1); a property which will be exploited in this chapter. We shall illustrate the framework briefly below using the singledimensional formulations (5.6), (5.12) as an example.

In the following we follow the notation of Chapter 5. In particular, · 0,<sup>Ω</sup> denotes the *L*<sup>2</sup> norm of a scalar or vector field over Ω, while · 1,Ω, · div,<sup>Ω</sup> are the *H*1(Ω) and **H**(div,Ω) norms respectively. For a Hilbert space *W* other than *L*2, *H*1, **H**(div) the norm is denoted as · <sup>W</sup> while (·, ·)<sup>W</sup> denotes the inner product on *W*. Moreover, we let (·, ·)<sup>Ω</sup> be the *L*<sup>2</sup> inner product. If the domain is clear the subscript will be omitted. By (·, ·) we shall also denote a duality pairing between a Hilbert space and its dual. The dual of an operator *A* with respect to the *L*<sup>2</sup> inner product is then denoted as *A* . Finally, let us introduce scaled, sum and intersection spaces, which will be required for well-posedness of some of the formulations. If *W*,*Q* are Hilbert spaces and *a* an arbitrary positive real number, the scaled space *aW*, the sum space *W* + *Q* and the intersection space *W* ∩ *Q* are Hilbert spaces with norms *x* aW = *a*· <sup>W</sup> ,

$$\|\|\mathbf{x}\|\|\_{W+\mathcal{Q}} = \inf\_{\substack{\mathbf{w}+\mathbf{q}=\mathbf{x} \\ \mathbf{w}\in W, \mathbf{q}\in \mathcal{Q}}} \sqrt{\|\|\mathbf{w}\|\|\_{W}^{2} + \|\|q\|\|\_{\mathcal{Q}}^{2}} \quad \text{and} \quad \|\mathbf{x}\|\|\_{W\cap\mathcal{Q}} = \sqrt{\|\|\mathbf{x}\|\|\_{W}^{2} + \|\|\mathbf{x}\|\|\_{\mathcal{Q}}^{2}}.$$

## **6.2.1 Single-Dimensional Primal Solvers**

In order to simplify the analysis let σi, σ<sup>e</sup> and *C*<sup>m</sup> be positive constants and let us consider a slightly modified (cf. the underlined term) single-dimensional primal formulation (5.6): Find *u*<sup>i</sup> ∈ *V*<sup>i</sup> = *H*1(Ωi), *u*<sup>e</sup> ∈ *V*<sup>e</sup> = *H*<sup>1</sup> ΓD e (Ωe) such that for all <sup>v</sup><sup>e</sup> <sup>∈</sup> *<sup>V</sup>*e, <sup>v</sup><sup>i</sup> <sup>∈</sup> *<sup>V</sup>*<sup>i</sup>

$$\begin{split} (\sigma\_{\varepsilon}\nabla\mu\_{\varepsilon}, \nabla\upsilon\_{\varepsilon}) + \frac{\mathcal{C}\_{m}}{\Delta t}(T\_{\varepsilon}\mu\_{\varepsilon}, T\_{\varepsilon}\upsilon\_{\varepsilon}) - \frac{\mathcal{C}\_{m}}{\Delta t}(T\_{i}\mu\_{i}, T\_{\varepsilon}\upsilon\_{\varepsilon}) &= \frac{\mathcal{C}\_{m}}{\Delta t}(f, \upsilon\_{\varepsilon}), \\ -\frac{\mathcal{C}\_{m}}{\Delta t}(T\_{\varepsilon}\mu\_{\varepsilon}, T\_{i}\upsilon\_{i}) + (\sigma\_{i}\nabla\mu\_{i}, \nabla\upsilon\_{i}) + \underbrace{(\mu\_{i}, \upsilon\_{i})}\_{\Delta t} + \frac{\mathcal{C}\_{m}}{\Delta t}(T\_{i}\mu\_{i}, T\_{i}\upsilon\_{i}) &= -\frac{\mathcal{C}\_{m}}{\Delta t}(f, \upsilon\_{i}). \end{split} (6.1)$$

Here, *T*e, *T*<sup>i</sup> are the trace operators *T*i*u*<sup>i</sup> = *u*<sup>i</sup> |<sup>Γ</sup> and *T*e*u*<sup>e</sup> = *u*<sup>e</sup> |Γ. Letting *W* = *V*<sup>e</sup> ×*V*i, *u* = (*u*e,*u*i) the problem (6.1) can be stated as A*u* = *L* in *W* with the operator A and functional *L* defined as

$$\mathcal{H} = \begin{pmatrix} -\nabla \cdot (\boldsymbol{\sigma}\_e \nabla) + \frac{\boldsymbol{C}\_m}{\Delta t} \boldsymbol{T}\_e' T\_e & -\frac{\boldsymbol{C}\_m}{\Delta t} \boldsymbol{T}\_e' T\_i \\\ -\frac{\boldsymbol{C}\_m}{\Delta t} \boldsymbol{T}\_i' T\_e & -\nabla \cdot (\boldsymbol{\sigma}\_i \nabla) + \boldsymbol{I} + \frac{\boldsymbol{C}\_m}{\Delta t} \boldsymbol{T}\_i' T\_i \end{pmatrix} \text{ and } \boldsymbol{L}(\boldsymbol{\nu}) = \frac{\boldsymbol{C}\_m}{\Delta t} \begin{pmatrix} (\boldsymbol{f}, \boldsymbol{\nu}\_e) \\ -(\boldsymbol{f}, \boldsymbol{\nu}\_i) \end{pmatrix}. \tag{6.2}$$

Note that the underlined lower order term in (6.1) corresponds to *I* in A.

To precondition (6.1) using the operator preconditioning framework we proceed to show that A : *W* → *W* is an isomorphism. This statement follows from the Lax-Milgram theorem, see e.g. (6, Ch 2.7). Let · <sup>W</sup> be *some* norm of the space *W* (we shall see shortly that different norms can be considered leading to different preconditioners). Hence coercivity of A, i.e.

$$\text{There exists } \alpha\_\* > 0 \text{ such that } \alpha\_\* \|u\|\_W^2 \le (\mathcal{R}u, u), \quad \forall u \in W \tag{6.3}$$

and boundedness of A, i.e.

$$\text{There exists } \alpha^\* > 0 \text{ such that } (\mathcal{H}u, \nu) \le \alpha^\* \|u\|\_W \|\nu\|\_W, \quad \forall u, \nu \in W \qquad (6.4)$$

e

need to be shown. While the details of the proof are beyond the scope of the current text, we remark that using the space *W* = *H*<sup>1</sup> ΓD (Ωe) × *H*1(Ωi), we obtained the bound

$$\rho(\mathcal{H}u, u) \le \max\left(\sigma\_e + 2C\_e \frac{C\_m}{\Delta t}, \sigma\_i, 1, 2C\_i \frac{C\_m}{\Delta t}\right) \left\|u\right\|\_{W}^2\tag{6.5}$$

by trace, Cauchy-Schwarz, and Young inequalities. Hence, the operator A is bounded. Note, however, that the constant α<sup>∗</sup> depends on Δ*t* and, in particular, it blows up for small time steps. Then, the lower bound is

$$\begin{split} \mathbb{P}(\mathcal{H}u, u) &= \sigma\_{\varepsilon} \|\nabla u\_{\varepsilon}\|\_{0}^{2} + \sigma\_{i} \|\nabla u\_{i}\|\_{0}^{2} + \|u\_{i}\|\_{0}^{2} + \frac{C\_{m}}{\Delta t} \|T\_{\varepsilon}u\_{\varepsilon} - T\_{i}u\_{i}\|\_{0}^{2} \\ &= \min\left(1, \sigma\_{i}, \sigma\_{\varepsilon}\right) \|u\|\_{W}^{2} \end{split} \tag{6.6}$$

and the coercivity thus holds with constant α<sup>∗</sup> = min (1, σi, σe) which is independent of the time step. The condition number of the preconditioned system, using a preconditioner deduced from *W*, is then <sup>α</sup><sup>∗</sup> <sup>α</sup><sup>∗</sup> and from this theoretical consideration we may then expect the number of iterations to grow linearly as Δ*t* decreases.

The exact preconditioner based on *W* is block diagonal opearator

$$\mathcal{B} = \begin{pmatrix} -\nabla \cdot (\sigma\_i \nabla) & 0 \\ 0 & -\nabla \cdot (\sigma\_e \nabla) \ +I \end{pmatrix}^{-1} \tag{6.7}$$

Observe that the preconditioner is efficient in the sense that its evaluation mounts to solving two smaller subproblems posed on Ω<sup>i</sup> and Ω<sup>e</sup> respectively. However, from our analysis we expect the performance of B to deteriorate for small time steps. Indeed, this behavior is confirmed by results<sup>1</sup> in Table 6.1. Therein it can be seen

<sup>1</sup> All the experiments in this chapter are conducted using the manufactured problem (5.17) considered in the convergence study in Chapter 5. In particular, we use uniform meshes obtained by dividing the unit square (Ωe <sup>∪</sup> <sup>Ω</sup>i) first into <sup>n</sup><sup>2</sup> h squares each of which is then subdivided to 2 isosceles triangles with diameter h. The finite element discretization is as described in Chapter 5. We set <sup>σ</sup>i <sup>=</sup> 1, <sup>σ</sup>e <sup>=</sup> <sup>2</sup>.2 and <sup>C</sup>m <sup>=</sup> 1.

The discrete linear systems <sup>A</sup>h <sup>x</sup>h <sup>=</sup> <sup>L</sup>h are solved iteratively using the CG or MinRes methods. Implementation of the methods is provided by (2). The solvers are always started from a random initial vector. As a convergence criterion the relative preconditioned residual norm is required to

that the iterations stabilize with *h* for Δ*t* > 10−<sup>4</sup> while for smaller time steps finer meshes are needed to attain the bounds.

In an attempt to find a more robust preconditioner let us observe that estimates (6.5), (6.6) show that *u* → (A*u*,*u*) <sup>1</sup>/<sup>2</sup> defines a different norm on *W*. Let now *u* <sup>W</sup> = (A*u*,*u*) <sup>1</sup>/2. Then (A*u*,*u*) = *u* <sup>2</sup> <sup>W</sup> and (A*u*, <sup>v</sup>) ≤ (A*u*,*u*) <sup>1</sup>/2(Av, v) <sup>1</sup>/<sup>2</sup> so that the conditions (6.3) and (6.4) hold with α<sup>∗</sup> = 1, α<sup>∗</sup> = 1. In another words, a Riesz map preconditioner with respect to the inner product (*u*, <sup>v</sup>)<sup>W</sup> <sup>=</sup> (A*u*, <sup>v</sup>) is independent of Δ*t*. As the Riesz map in this case is in fact A−<sup>1</sup> its exact evaluation (by LU) is not feasible. However, for the purpose of preconditioning, it suffices to replace the mapping by a spectrally equivalent operator. Then, provided that the equivalence is robust in Δ*t* the approximate operator will lead to iterations bounded in time step and the discretization parameter. We remark that since the preconditioner takes the entire A into account, including the off-diagonal terms cf. block diagonal operator (6.7), we shall refer to it as monolithic.


Table 6.1: Primal single-dimensional formulation. Number of CG iterations using (left) the diagonal operator (6.7) and (right) the Riesz map with respect to the A induced norm. Neither preconditioner is robust for Δ*t* 1.

Since (5.6) leads to symmetric positive-definite matrices, cf. the coercivity condition (6.3), we shall here use algebraic multigrid to construct the approximation of A<sup>−</sup>1. More precisely, the action of the operator is realized by single V-cycle of BoomerAMG (23). In Table 6.1 it can be seen that this choice leads to *h* bounded CG iterations for Δ*t* > 10<sup>−</sup>8. However, there is a clear sensitivity of the bound to Δ*t*. Thus BoomerAMG approximations of A−<sup>1</sup> are not Δ*t*- robust. In fact, estimates of the condition number for the preconditioned problems with *n*<sup>h</sup> = 2<sup>8</sup> are 1.2, 1.5, 2.0, 9.2 and 53 for Δ*t* = 1, 10−2, 10−4, 10−<sup>6</sup> and 10−<sup>8</sup> respectively.

be less than 10−<sup>12</sup> in magnitude. Unless specified otherwise the preconditioners <sup>B</sup>h are evaluated exactly by LU. Finally, condition number estimates of the preconditioned linear systems are obtained by using iterative Krylov-Schur solver from (10) applied to the generalized eigenvalue problem <sup>A</sup>h <sup>x</sup>h <sup>=</sup> <sup>λ</sup>h <sup>B</sup>−<sup>1</sup> h <sup>x</sup>h. The reported condition number is then max|λ<sup>h</sup> |/min|λ<sup>h</sup> <sup>|</sup>.

The source code used for the experiments can be found on https://github.com/MiroK/ emi-book-solvers and is arxived at (15).

## **6.2.2 Single-Dimensional Mixed Solvers**

We next consider preconditioners for the single-dimensional mixed formulation (5.12): Find **E** ∈ *S*, *u* ∈ *Q* such that

$$\begin{aligned} \left( (\boldsymbol{\sigma}^{-1} \boldsymbol{J}, \boldsymbol{\tau}) + \frac{\boldsymbol{\Lambda}t}{\boldsymbol{C}\_{m}} (T\_{\mathfrak{n}} \boldsymbol{J}, T\_{\mathfrak{n}} \boldsymbol{\tau}) - (u, \nabla \cdot \boldsymbol{\tau}) = -(f, T\_{\mathfrak{n}} \boldsymbol{\tau}), \quad \forall \boldsymbol{\tau} \in S, \\ (q, \nabla \cdot \boldsymbol{J}) = 0, \quad \forall q \in Q, \end{aligned}$$

where *<sup>S</sup>* <sup>=</sup> **<sup>H</sup>**ΓN e (div,Ω), *<sup>Q</sup>* <sup>=</sup> *<sup>L</sup>*2(Ω) and *<sup>T</sup>***<sup>n</sup>** is the normal trace operator *<sup>T</sup>***n**<sup>τ</sup> <sup>=</sup> <sup>τ</sup> · **<sup>n</sup>**i. Letting *W* = *S*×*Q*, *x* = (**J**,*u*) the single-dimensional mixed formulation is equivalent to the problem A*x* = *L* in *W* with

$$\mathcal{A} = \begin{pmatrix} \sigma^{-1}I + \frac{\Delta t}{C\_m}T\_\mathbf{n}^\prime T\_\mathbf{n} & -\nabla \\ \nabla \cdot & 0 \end{pmatrix} \quad \text{and} \quad L(\mathbf{x}) = -(f, T\_\mathbf{n}J). \tag{6.8}$$

Note that in (6.8) the membrane term (*T***<sup>n</sup> J**,*T***n**τ) is weighted by Δ*t*/*C*m, cf. (6.2). Thus, unlike in the single-dimensional primal formulation, the term does not dominate for small time steps.

In order to apply the operator preconditioning framework to (6.8) the operator A needs to be shown to be an isomorphism. To this end we consider A as an abstract operator over *W* = *S* × *Q* with the form

$$\mathcal{A} = \begin{pmatrix} A \ B' \\ B \ 0 \end{pmatrix} \quad \text{where} \quad \begin{array}{l} A:S \to V', \\ B:S \to Q', \end{array} \tag{6.9}$$

and apply the Brezzi theory (7). This leads us to the potential preconditioners for the multi-dimensional mixed problem

$$\mathcal{B}\_1 = \begin{pmatrix} \sigma^{-1}I - \nabla \nabla \cdot \mathbf{0} \\ \mathbf{0} & I \end{pmatrix}^{-1} \quad \text{and} \quad \mathcal{B}\_2 = \begin{pmatrix} \sigma^{-1}I - \nabla \nabla \cdot \frac{\Delta t}{C\_m} T\_\mathbf{n}' T\_\mathbf{n} & \mathbf{0} \\ \mathbf{0} & I \end{pmatrix}^{-1}, \quad (6.10)$$

which are the Riesz mappings with respect to the inner products of the spaces

$$W\_1 = S \times \mathcal{Q} \quad \text{and} \quad W\_2 = S \cap \sqrt{\frac{\Delta t}{C\_m}} N \times \mathcal{Q},\tag{6.11}$$

where *N* =  **<sup>v</sup>** <sup>∈</sup> *<sup>S</sup>*; **<sup>v</sup>** · **<sup>n</sup>** <sup>0</sup>,<sup>Γ</sup> <sup>&</sup>lt; <sup>∞</sup> . Note that in *W*<sup>2</sup> we enforce additional regularity on the normal trace since *T***<sup>n</sup>** maps *S* to *H*−1/2(Γ), e.g. (9, Ch. 2). Moreover the traces are considered in a weighted space.

A consequence of the mapping properties of the normal trace operator is the fact that the term (*T***<sup>n</sup> J**,*T***n**τ) cannot be bounded using the **H**(div) norm and in turn the Brezzi conditions are not satisfied on *W*1. We remark that if such a bound were possible the boundedness constant would depend on Δ*t*, cf. (6.5). As the Brezzi conditions do not hold on *W*<sup>1</sup> we expect preconditioner B<sup>1</sup> to perform poorly. Indeed, Table 6.2 shows that the condition numbers of B1,hA<sup>h</sup> are not bounded in *h*. The condition numbers for fixed *h* can be seen to grow with the time step. In Table 6.3 the ill-posedness of the problem in *W*<sup>1</sup> manifests as unstable iterations or failure of MinRes to converge. We note that for small Δ*t* the iterations appear bounded in *h*.


Table 6.2: Conditioning of single-dimensional mixed formulation. (Left) Posing the problem in *H*(div) × *L*<sup>2</sup> violates Brezzi conditions. (Right) Preconditioner based on *W*<sup>2</sup> in (6.11) is inf-sup stable with the inf-sup constant depending on Δ*t* for Δ*t*/*C*<sup>m</sup> > 1.

Posing (6.8) in *W*<sup>2</sup> it can be shown that the Brezzi conditions hold with the constants independent of the time step.

We shall not prove the validity of the Brezzi theory here. Instead, Table 6.2 offers numerical evidence that the discrete condition is satisfied. Therein condition numbers of the B<sup>2</sup> preconditioned problems are reported and boundedness in *h* can be observed. Moreover, it can be seen that the inf-sup constant depends on Δ*t*, in particular, the bound goes to 0 as Δ*t* grows. However, on the subspace  (**J**,*u*) ∈ *<sup>W</sup>*2; (*u*, <sup>1</sup>)0,Ωi <sup>=</sup> <sup>0</sup> the inf-sup condition holds independent of Δ*t*. The single run-away mode appears to have no effect on the MinRes iterations, cf. Table 6.3, see also (20).


Table 6.3: MinRes iterations with preconditioned single-dimensional mixed formulation using preconditioners (6.10). Failure to converge within 500 iterations is denoted by –. Convergence of B<sup>2</sup> seems unaffected by the single Δ*t* unbounded inf-sup mode.

Considering order optimality of the B<sup>2</sup> based solver we recall that the iterations in Table 6.3 were run with the exact preconditioner. In particular, the *S* ∩ *N* Riesz map was computed by LU. In practical computations such realization is not feasible and scalable approximation is needed. For standard **H**(div) inner product, i.e. the one used in B1, the Riesz map can be efficiently realized by multigrid algorithms (14). However, in the authors' experience this approach does not work equally well with the *S* ∩ *N* inner product.

We have seen in Sections 6.2.1 and 6.2.2 that construction of order optimal monolithic solvers for the single-dimensional EMI formulations presents a challenge. In primal formulationΔ*t* robustness of the monolithic approach was problematic. For the mixed formulation, order optimality required a specialized solver for the Riesz mapping over the subspace of **H**(div). These two problems are addressed by solvers for multidimensional formulations.

## **6.2.3 Multi-Dimensional Solvers**

In this section we consider the construction of preconditioners for the operators

$$\mathcal{R}\_{p} = \begin{pmatrix} -\nabla \cdot (\boldsymbol{\sigma}\_{e} \boldsymbol{\nabla}) & 0 & -T\_{e}' \\ 0 & -\nabla \cdot (\boldsymbol{\sigma}\_{i} \boldsymbol{\nabla}) & T\_{i}' \\ -T\_{e} & T\_{i} & -\frac{\Delta t}{C\_{m}}I \end{pmatrix} \quad \text{and} \quad \mathcal{R}\_{m} = \begin{pmatrix} \boldsymbol{\sigma}^{-1}I - \boldsymbol{\nabla} & T\_{\mathbf{n}}' \\ \boldsymbol{\nabla} \cdot & 0 & 0 \\ T\_{\mathbf{n}} & 0 & -\frac{C\_{m}}{\Delta t}I \end{pmatrix}, \quad (6.12)$$

which induce respectively the multi-dimensional primal weak formulation (5.8) and the multi-dimensional mixed weak formulation (5.14). In order to discuss wellposedness let

$$W\_p = H\_{\Gamma\_e^D}^1(\Omega\_e) \times H^1(\Omega\_i) \times H^{-1/2}(\Gamma) \cap \sqrt{\frac{\Delta t}{C\_m}} L^2(\Gamma) \tag{6.13}$$

and

$$W\_m = H\_{\Gamma\_e^N}(\text{div}, \Omega) \times L^2(\Omega) \times H^{1/2}(\Gamma) \cap \sqrt{\frac{C\_m}{\Delta t}} L^2(\Gamma). \tag{6.14}$$

We remark that the fractional space *H*−1/2, in which the membrane current density *<sup>I</sup>*<sup>m</sup> in (5.6) is sought, and *<sup>H</sup>*1/2, which is the space of the transmembrane potential <sup>v</sup> in (5.14), are dictated by the mapping properties of the trace operators *T*e, *T*<sup>i</sup> and *T***n**. For example, as *T*<sup>i</sup> : *H*1(Ωi) → *H*1/2(Γ), *I*<sup>m</sup> ∈ (*H*1/2) so that the term (*T*i*u*i, *I*m) is bounded. Note also that only the spaces involving Γ are now Δ*t*-dependent, cf. e.g. *W*<sup>2</sup> in (6.11).

Due to the Δ*t*-weighted (penalty) terms in A<sup>p</sup> and A<sup>m</sup> the operators cannot be established as isomorphisms on *W*<sup>p</sup> and *W*<sup>m</sup> by straightforward application of the Brezzi theory. Instead, A<sup>p</sup> for Δ*t*/*C*<sup>m</sup> ≤ 1 and A<sup>m</sup> for Δ*t*/*C*<sup>m</sup> ≥ 1 can be analyzed using the framework for saddle point systems with penalty, see (4, Ch. 3.4). A crucial assumption then is that the system without the penalty term satisfies the Brezzi conditions. While here the well-posedness shall be demonstrated only by numerical experiments, let us point out an important difference between the operators. To this end let *c* > 0 be a constant and consider a piecewise constant potential *u* such that *u*<sup>i</sup> = *c*, *u*<sup>e</sup> = 0 and set **J** = **0**. Then for any τ ∈ **H**(div) by divergence theorem

$$c(\nabla \cdot \tau, \mu)\_{\Omega} + (T\_{\mathfrak{n}}\tau, \upsilon) = (\nabla \cdot \tau, \mu)\_{\Omega\_i} + (T\_{\mathfrak{n}}\tau, \upsilon) = c(T\_{\mathfrak{n}}\tau, 1) + (T\_{\mathfrak{n}}\tau, \upsilon)$$

and in turn, if v <sup>=</sup> <sup>−</sup>*<sup>c</sup>* we have that

$$
\begin{pmatrix}
\sigma^{-1}I - \nabla \ T\_{\mathbf{n}}' \\
\nabla \cdot & 0 & 0 \\
T\_{\mathbf{n}} & 0 & 0
\end{pmatrix}
\begin{pmatrix}
\mathbf{J} \\
\boldsymbol{u}, \\
\boldsymbol{\nu}
\end{pmatrix} = 
\begin{pmatrix}
c(T\_{\mathbf{n}}\boldsymbol{\tau}, \boldsymbol{1}) + (T\_{\mathbf{n}}\boldsymbol{\tau}, \boldsymbol{\nu}) \\
0 \\
0
\end{pmatrix} = 
\begin{pmatrix}
\mathbf{0} \\
0 \\
0
\end{pmatrix}.
$$

Thus the operator A<sup>m</sup> without the penalty term is singular with a one dimensional kernel spanned by vector (**0**,*u*,−*c*) and the Brezzi conditions do not hold. We shall see implications of this property on the convergence of the iterative solvers shortly. We remark that A<sup>p</sup> without the penalty term is non-singular.

Assuming that A<sup>p</sup> : *W*<sup>p</sup> → *W* <sup>p</sup> is an isomorphism we consider as a preconditioner for the single-dimensional primal formulation the operator

$$\mathcal{B}\_{p} = \begin{pmatrix} -\nabla \cdot (\sigma\_{e} \nabla) & 0 & 0 \\ 0 & -\nabla \cdot (\sigma\_{i} \nabla) + I & 0 \\ 0 & 0 & (-\Delta + I)^{-1/2} + \frac{\Delta t}{C\_{m}}I \end{pmatrix}^{-1} \tag{6.15}$$

Observe that the preconditioner consists of (approximate) solvers for three decoupled subproblems posed on Ωe, Ω<sup>i</sup> and Γ. Moreover, the problems on extra and intracellular domains are standard elliptic operators for which efficient (black-box) multigrid techniques exist, e.g. (23). The problem on the interface is then less standard as it concerns a fractional Helmholtz operator. However, efficient multilevel solvers have been established e.g. in (5; 1). We remark that if the discrete spaces for *<sup>I</sup>*<sup>m</sup> (or <sup>v</sup>) contain only a few thousands of degrees of freedom an eigenvalue realization of the fractional preconditioner is feasible cf. (16). As the interfacial spaces here are small, cf. Table (5.1), we use further the spectral approach.


Table 6.4: Fractional preconditioner (6.15) for multi-dimensional primal formulation. (Left) Condition numbers are bounded in *h*. Growth in Δ*t* for Δ*t* > 1 is caused by a single mode. (Right) MinRes iterations counts are bounded in *h* and time step.

#### 6 EMI Preconditioners 79

Using (6.15) Table 6.4 reports the condition numbers of the preconditioned problems Bp,hAp,h. It can be seen that the results are bounded in *h*, thus providing a numerical evidence for A<sup>p</sup> being an isomorphism on *W*<sup>p</sup> in (6.13). Moreover, it can be seen that the conditioning of the problem deteriorates with Δ*t* large. However, in this case it is only a single mode *u*<sup>i</sup> = 0, *u*<sup>e</sup> = const, *I*<sup>m</sup> = 0 which causes the blowup. Considering the MinRes iteration counts, the presence of this mode seems to have almost no impact on boundedness in *h* or Δ*t*.

We finally come back to the multi-dimensional mixed formulation. Based on the space *W*<sup>m</sup> in (6.14) let a preconditioner for the multi-dimensional mixed formulation be

$$\mathcal{B}\_m = \begin{pmatrix} \sigma^{-1}I - \nabla \nabla \cdot \mathbf{0} & \mathbf{0} \\ \mathbf{0} & I & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \left( -\Delta + I \right)^{1/2} + \frac{C\_m}{\Delta t} I \end{pmatrix}^{-1} \tag{6.16}$$

We remark that the preconditioner uses a standard **H**(div) inner product, cf. B<sup>2</sup> in (6.10), which can be efficiently realized by multigrid methods, e.g. (14). Note also that the fractionality of the Laplacian is 1/2, cf. −1/2 of the multi-dimensional primal preconditioner (6.15). In addition to the previously mentioned multilevel methods, problems (−Δ*u* + *u*) <sup>s</sup> *x* = *b* for 0 < *s* < 1 can be efficiently solved by a number of approaches, see (3) and references therein.

Using (6.16) Table 6.5 reports the condition numbers of the preconditioned multidimensional mixed formulation. The conditioning can be seen to be stable in *h*, while in agreement with the limit singularity property, there is a growth with Δ*t*. Given that only a single mode is responsible for the lack of Δ*t*-boundedness and recalling results of Table 6.4 or 6.3 we might expect that also here the MinRes solver will not be affected. However, Table 6.6 shows that this is not the case. In fact, as Δ*t* grows the iterations become unstable in *h*. Figure 6.1 then shows a typical convergence behavior of the solver. We see that the relative preconditioned residual norm is quickly reduced to about 10−<sup>9</sup> in approximately 30 iterations (regardless of the mesh resolution). Afterwards the convergence stalls. We conclude that for robust preconditioning the nullspace of A<sup>m</sup> must be addressed.


Table 6.5: Conditioning of the multi-dimensional mixed formulation. (Left) Problem considered on *W*<sup>m</sup> in (6.14) yields operator A<sup>m</sup> with a singular limit as Δ*t* → ∞. (Right) Formulation including an additional constraint on the transmembrane potential <sup>∫</sup> <sup>Γ</sup> <sup>v</sup> <sup>d</sup>*<sup>S</sup>* <sup>=</sup> 0 yields *<sup>h</sup>* and <sup>Δ</sup>*<sup>t</sup>* robust condition numbers.

Motivated by the observation that A<sup>m</sup> becomes singular in the limit Δ*t* = ∞ with a mode **<sup>J</sup>** <sup>=</sup> **<sup>0</sup>**, *<sup>u</sup>*<sup>e</sup> <sup>=</sup> 0, *<sup>u</sup>*<sup>i</sup> <sup>=</sup> *<sup>c</sup>*, <sup>v</sup> <sup>=</sup> <sup>−</sup>*<sup>c</sup>* in the kernel, let *<sup>W</sup>*<sup>0</sup> <sup>e</sup> = *W*<sup>e</sup> × R and let us define a new solution operator for the multi-dimensional mixed formulation A<sup>0</sup> <sup>m</sup> : *W*<sup>0</sup> <sup>e</sup> → (*W*<sup>0</sup> e ) and its preconditioner B<sup>0</sup> <sup>m</sup> : (*W*<sup>0</sup> e ) → *W*<sup>0</sup> <sup>e</sup> as

$$
\mathcal{R}\_m^0 = \begin{pmatrix}
\sigma^{-1}I - \nabla & T\_\mathbf{n}' & 0 \\
\nabla \cdot & 0 & 0 & 0 \\
T\_\mathbf{n} & 0 & -\frac{C\_m}{\Delta t}I \ I \\
0 & 0 & I & 0
\end{pmatrix},
\mathcal{B}\_m^0 = \begin{pmatrix}
\sigma^{-1}I - \nabla \nabla \cdot 0 & 0 & 0 \\
0 & I & 0 & 0 \\
0 & 0 \left(-\Delta + I\right)^{1/2} + \frac{C\_m}{\Delta t}I & 0 \\
0 & 0 & 0 & \mu I
\end{pmatrix}^{-1},
\tag{6.17}
$$

where μ = min (1,*C*m/Δ*t*). Note that A<sup>0</sup> <sup>m</sup> includes an additional unknown, a single scalar, which enforces (1, <sup>v</sup>)<sup>Γ</sup> <sup>=</sup> 0. The new constraint thus eliminates constant transmembrane potential and in turn the new operator is non-singular2. Considering results reported in Tables 6.5, 6.6 it can be seen that the new preconditioned formulation leads to condition numbers and iteration counts, which are bounded in both the mesh size and the time step. Note that the extra unknown has only a small impact on the number of iterations compared to the B<sup>m</sup> preconditioner.


Table 6.6: MinRes iterations for the multi-dimensional mixed formulation. No convergence in 500 iterations is indicated by –. (Left) For large Δ*t* the iterations are unstable since operator A<sup>m</sup> becomes singular in the limit as Δ*t* → ∞. (Right) Constraining transmembrane potential the operator A<sup>0</sup> <sup>m</sup> does not have the limit singular property.

$$(1, I\_m + I\_{\ell m})\_{\Gamma} = (1, I\_m)\_{\Gamma} = (1, -\nabla \cdot (\sigma \tau\_{\ell} \nabla \mu\_{\ell}))\_{\Omega\_{\ell}} = 0.1$$

Thus we obtain a conservation relation <sup>C</sup>m <sup>d</sup> <sup>d</sup>t (1, <sup>v</sup>)<sup>Γ</sup> <sup>=</sup> 0.

<sup>2</sup> A physical motivation for the constraint can be found in considering the EMI interface equation <sup>C</sup>m <sup>∂</sup><sup>v</sup> <sup>∂</sup>t <sup>=</sup> <sup>I</sup><sup>m</sup> <sup>+</sup> <sup>I</sup>ion on <sup>Γ</sup>. By integrating the left-hand side we obtain <sup>C</sup><sup>m</sup> <sup>d</sup> <sup>d</sup>t (1, <sup>v</sup>)Γ. Recall that <sup>I</sup>m <sup>=</sup> <sup>σ</sup>i <sup>∇</sup>ui · **<sup>n</sup>**i on <sup>Γ</sup> and −∇ · (σi <sup>∇</sup>ui) <sup>=</sup> 0 in <sup>Ω</sup>i. Then, setting <sup>I</sup>ion <sup>=</sup> 0 and integrating the right-hand side we have

Fig. 6.1: (Left) Convergence of the MinRes method for (6.12)-based multidimensional mixed formulation with Δ*t* = 10−<sup>8</sup> and problem from Table 6.6. For all mesh resolutions the norm is reduced in cca. 30 iterations. Afterwards the reduction stalls. (Right) Convergence of the Robin/Neumann domain decomposition algorithm. For small Δ*t* the algorithm becomes sensitive to mesh size and can even diverge.

# **6.3 Domain Decomposition Solvers**

Having seen that monolithic solvers for the EMI equations can be sensitive to spatial and temporal resolution we next briefly discuss robustness of the non-overlapping domain decomposition (DD) approach. The focus here shall only be on the singledimensional primal formulation and the Robin/Neumann DD algorithm in the form presented in (12, Chapter 4). In particular, our implementation shall not include any coarse space, cf. (21), or a preconditioner, see e.g. (8) for interpretation of DD as Steklov-Poincaré operators in fractional Sobolev spaces.

For the sake of self-containedness we review the variational formulation of the Robin/Neumann algorithm. Let *V*<sup>e</sup> = *H*<sup>1</sup> ΓD e (Ωe), *<sup>V</sup>*<sup>i</sup> <sup>=</sup> *<sup>H</sup>*1(Ωi) and <sup>v</sup>0, *<sup>u</sup>*<sup>0</sup> <sup>e</sup> be the given initial transmembrane and extracellular potentials. A single iteration of the DD algorithm then produces a new approximation v<sup>1</sup> in three steps. (i) We find *u*<sup>i</sup> ∈ *V*<sup>i</sup> such that

$$(\sigma\_i \nabla \mu\_i, \nabla \upsilon\_i) + \frac{C\_{\text{m}}}{\Delta t}(T\_i \mu\_i, T\_i \upsilon\_i) = -\frac{C\_{\text{m}}}{\Delta t}(f, \upsilon\_i) + \frac{C\_{\text{m}}}{\Delta t}(T\_e \mu\_e^0, T\_i \upsilon\_i), \quad \forall \upsilon\_i \in V\_i.$$

(ii) Using the computed intracellular potential the new extracellular potential *u*<sup>e</sup> ∈ *V*<sup>e</sup> is found such that

$$(\sigma\_e \nabla \mu\_e, \nabla \upsilon\_e) = (-\sigma\_i \nabla \mu\_i \cdot n\_i, T\_e \upsilon\_e), \quad \forall \upsilon\_e \in V\_e.$$

(iii) Finally <sup>v</sup><sup>1</sup> <sup>=</sup> *<sup>T</sup>*i*u*i−*T*e*u*e. The iterations continue by assigning v<sup>0</sup> <sup>=</sup> v<sup>1</sup> and *<sup>u</sup>*<sup>0</sup> <sup>e</sup> = *u*<sup>e</sup> until convergence. which in the following is determined as <sup>v</sup>k−vk−<sup>1</sup> <sup>0</sup>/ vk−<sup>1</sup> <sup>0</sup> <sup>&</sup>lt;, where  = 10<sup>−</sup>5. We remark that the tolerance was chosen such that in the refinement studies, where convergence of DD was observed, the approximation error in the *H*<sup>1</sup> × *H*<sup>1</sup> norm decreased with the expected order as O(*h*) (we recall that *P*<sup>1</sup> − *P*<sup>1</sup> elements were used).

To test robustness of the DD solver we consider the setup of the single cell experiment used with the monolithic approaches in Section 6.2. Here the Robin problem (i) and the Neumann problem (ii) shall be solved by LU to eliminate effects of inexact subdomain solvers. We remark that due to the conforming triangulation, see Figure 5.1, and *P*<sup>1</sup> discretization of both *V*<sup>i</sup> and *V*<sup>e</sup> the discrete transmembrane potential is computed simply by interpolation. Figure 6.1 plots evolution of the potential difference v<sup>k</sup> <sup>h</sup> <sup>−</sup>vk−<sup>1</sup> <sup>h</sup> <sup>0</sup>/ vk−<sup>1</sup> <sup>h</sup> <sup>0</sup> with DD iterations for several mesh resolutions and time step values Δ*t* ≤ 1. It can be seen that for Δ*t* = 1 the algorithm converges in about 5 iterations irrespective of the spatial discretization. For smaller timesteps convergence is delayed on finer meshes and for Δ*t* = 10−<sup>4</sup> the iterations diverge.

Considering the results in Figure 6.1 we conclude that in the form presented here domain decomposition is not a robust algorithm for the single-dimensional primal formulation of the EMI equations. However, due to its simplicity (relative to e.g. the multi-dimensional formulations) and speed (see Section 6.4), DD might be the method of choice in practical cardiac modeling. It is therefore natural to ask whether the divergence conditions are likely to arise in real applications. To address this question we perform a scaling analysis based on the membrane dynamical condition *C*<sup>m</sup> <sup>∂</sup><sup>v</sup> <sup>∂</sup><sup>t</sup> <sup>∼</sup> <sup>σ</sup>i∇*u*<sup>i</sup> · **<sup>n</sup>**i. Letting *<sup>L</sup>* be a characteristic length scale and <sup>∇</sup> <sup>=</sup> *<sup>L</sup>*<sup>∇</sup> the equation becomes (*LC*m/σi) <sup>∂</sup><sup>v</sup> <sup>∂</sup><sup>t</sup> <sup>=</sup> <sup>∇</sup> *u*<sup>i</sup> · **n**<sup>i</sup> and the pre-factor *T* = *LC*m/σ<sup>i</sup> can be seen to have the unit of seconds. Following (22) let us insert *C*<sup>m</sup> = 1μF/cm2, σ<sup>i</sup> = 10mS/cm, *L* = 100μm, where the length scale is determined by the size of a typical cell used by the authors'. We remark that this a sensible choice given the setup of our experiments. As a result *T* = 10<sup>−</sup>6s and thus Δ*t* = 1 in the experiments reported here corresponds to a time step of 10−<sup>6</sup> seconds. Note that this is the finest time step considered in (11). Furthermore, therein the spatial resolution is 2μm while here with *n*<sup>h</sup> = 2<sup>8</sup> and *L* = 100μm the mesh size is approximately 0.4μm. In summary, the conditions for divergence of the DD algorithm discussed here can be encountered not far away from the parameter regime in (11).

# **6.4 Solver Comparison**

To complete our discussion of EMI solvers we finally address the speed of the different algorithms. Recall that until this point results for all solvers, but the monolithic single-dimensional primal one, were obtained using LU in the construction. Such an implementation, however, is not scalable. Here we show that the proposed algorithms are order optimal if LU is replaced by suitable multilevel methods.

Referring to the legend of Figure 6.2 we shall compare 8 different methods. Singledimensional primal formulation (sp) shall be solved either with the BoomerAMG(A) preconditioner (mono), or diagonal preconditioner (6.7). In the latter case all the blocks of the diagonal preconditioner are approximated by single algebraic multigrid V-cycle. The single-dimensional mixed (sm) formulations shall use the B<sup>1</sup> preconditioner (6.10) with the leading block approximated by **H**(div) multigrid HypreAMS of (14). We recall that the solver in general is not independent of the discretization, cf. Table 6.3, however, HypreAMS does not work well for the robust preconditioner B2. Finally, the multi-dimensional primal (mp) and multi-dimensional mixed (mm) formulations shall use (6.15) and (6.16) with the Ωi, Ω<sup>e</sup> and Ω subproblems of the preconditioners approximated by multigrid (BoomerAMG in case of (6.15) and HypreAMS for (6.16)). The fractional operators will be computed exactly. In addition, two implementations of Robin/Neumann domain decomposition algorithm are considered with the subproblems solved either exactly by LU or by 4 V-cycles of BoomerAMG. Finally, a reference solution time is provided by timings of exact solution of the single-dimensional formulation (sp-LU).

We compare the algorithms using the single cell setup of the previous experiments with Δ*t* = 10−<sup>3</sup> where this value was chosen with the intention to not favor any of the methods. Using tolerance  = 10−<sup>3</sup> for DD and 10−<sup>12</sup> for CG/MinRes, Figure 6.2 reports solution times<sup>3</sup> of the algorithms for 23 <sup>≤</sup> *<sup>n</sup>*<sup>h</sup> <sup>≤</sup> <sup>2</sup>10. It can be seen that with the exception of LU-based solvers all the methods are indeed order optimal. The monolithic approach (sp-mono) and the multi-dimensional primal (mp) solver are then the fastest solvers while solvers for the mixed formulations (sm-B1, mm) are the slowest. For *n*<sup>h</sup> = 210 the solution times of the solvers were cca. 26s, 140s, 640s, 570s respectively. We remark that the mixed formulations have approximately 5 times more unknowns compared to the single-dimensional ones. Thus the cost per degree of freedom is comparable between the two approaches. Considering the scalable DD implementation (DD-AMG) the timing for *n*<sup>h</sup> = 210 is cca. 160s. The method thus has a similar cost to multi-dimensional primal formulation.

<sup>3</sup> The timings include setup times of the preconditioners for monolithic methods. For DD the subproblems were assembled and their solvers constructed only once.

Fig. 6.2: Solution times (including preconditioner setup) of the solvers for the EMI system (5.1) in terms of number of unknowns *N* or mesh size *h*. Single cell model with Δ*t* = 10−<sup>3</sup> is used. All but multi-dimensional primal solver with LU (sp-LU) and LU-based domain decomposition solver (DD-LU) are order optimal. The former scales as *N*3/2. Legend is shared between the figures.

**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.

# **References**


# **Chapter 7**

# **Improving Neural Simulations with the EMI Model**

Alessio Paolo Buccino1,2, Miroslav Kuchta3, Jakob Schreiner3, and Kent-André Mardal3,<sup>4</sup>

**Abstract** Mathematical modeling of neurons is an essential tool to investigate neuronal activity alongside with experimental approaches. However, the conventional modeling framework to simulate neuronal dynamics and extracellular potentials makes several assumptions that might need to be revisited for some applications. In this chapter we apply the EMI model to investigate the ephaptic effect and the effect of the extracellular probes on the measured potential. Finally, we introduce reduced EMI models, which provide a more computationally efficient framework for simulating neurons with complex morphologies.

# **7.1 Introduction**

In recent years, huge efforts and resources have been spent in computational modeling of neuronal activity. For example, the Blue Brain Project (18; 16)(https: //bbp.epfl.ch/nmc-portal/welcome) has constructed and distributed several hundreds of biophysically detailed cell models (*multi-compartment models*) from rat sensory cortex. A similar effort is being conducted by the Allen Institute of Brain Science, whose cell-type database (8) (https://celltypes.brain-map.org/) includes hundreds of cell models both from mice and even from humans. As the experimental data used become more comprehensive and available, these models are expected to become elaborated and more accurate in reproducing neuronal dynamics. However, the modeling framework which is commonly used to simulate these

<sup>1</sup>Bio Engineering Laboratory, Department of Biosystems and Science Engineering, ETH Zurich, Switzerland

<sup>2</sup>Center for Integrative Neuroplasticity (CINPLA), Faculty of Mathematics and Natural Sciences, University of Oslo, Oslo, Norway

<sup>3</sup>Simula Research Laboratory, Norway

<sup>4</sup>Department of Mathematics, University of Oslo, Oslo, Norway

multi-compartment models makes several key assumptions that might be violated for certain applications.

The most widely used approach to simulate neuronal dynamics of neurons is the *cable equation*. The solution of this equation enables one to compute transmembrane currents for each of the model compartment. In order to simulate extracellular potentials, we can use *volume conduction theory* and sum the individual contributions of the currents to the electric potential at any point in space (6). Whereas the use of this modeling framework has been the gold standard to simulate neuronal activity for decades (19; 6), there are some important assumptions that need to be discussed:


In applications where any of the assumptions listed above may be violated the EMI model (10, (1.30)) provides a suitable modeling framework. In particular, the geometry of the neuron (and the extracellular space) is accurately represented. Here we will show how the model is convenient to study both the ephaptic coupling of neurons (Section 7.3) and the effects of extracellular probes on the recorded electric potentials (Section 7.4). However, the detailed representation of the geometries makes EMI much more computationally intense than the standard modeling framework (21; 3). This limits the complexity of the simulation mainly to simple neuronal morphologies, such as ball-and-stick neurons (3). In order to target more realistic morphologies Section 7.5 discusses the reduced EMI model where the (three-dimensional) intracellular space is lumped to a curve that is the centerline of the neuron.

# **7.2 EMI Simulations of Neurons using the** neuronmi **Python Package**

Before discussing the applications let us first briefly introduce a Python package called neuronmi1, which has been used for the simulation results reported further. neuronmi provides a high-level application programming interface (API) to enable users to easily set up and run EMI simulations of neurons.

The workflow of the neuronmi package consists of two parts. First, a mesh needs to be generated. This is done with the generate\_mesh function, that uses gmsh (7) as backend. With this function the user can choose different kinds of neurons to place in the mesh and optionally place a probe device in the extracellular space. Mesh resolution and sizes of the bounding box can also be adjusted, as well as parameters of the neurons and the probe. In the following code example, we create a mesh with a ball–and–stick neuron (bas) and a microwire probe. By default, the center of the soma is at (0,0,0) μm, the dendrite extends in the positive z direction and the axon in the negative z direction.

```
import neuronmi
mesh_folder = neuronmi.generate_mesh(neurons='bas',
                                      probe='microwire')
```
Once a mesh is generated, the EMI simulation can be invoked with the simulate\_emi function, which implements the finite element method for the multi-dimensional mixed formulation (14, (5.14)) following the discretization proposed in (20). Through a set of parameters, the user can stimulate the neuron with a synaptic input, a step current, or a pulse current. Alternatively, the probe contacts can be used to stimulate the neuron extracellularly. By default, the neuron receives a synaptic input on its dendrite. The user can probe electric potentials *u* at any point in the mesh, while transmembrane currents *<sup>i</sup>* and membrane potentials v are available at facets on the neuron surface. The full solutions are also saved as pvd or xdmf files. The simulation is run as follows:

```
u, i, v = neuronmi.simulate_emi(mesh_folder,
                                u_probe_locations=points_v,
                                i_probe_locations=points_i,
```
v\_probe\_locations=points\_v)

<sup>1</sup>https://github.com/MiroK/nEuronMI

Several parameters can be set to customize the mesh generation and the simulation. For further details, we refer to the package documentation (https://neuronmi. readthedocs.io/en/latest/).

# **7.3 Investigating the Ephaptic Effect between Neurons**

Neurons are surrounded by the electrically conductive extracellular space. Groups of neurons create fluctuations in the local extracellular electrical field. These fluctuations in turn influence the intracellular electrical field through the *ephaptic* effect. Ephaptic coupling cannot influence neurons at rest, however, it can affect the spike timings of a neuron receiving suprathreshold stimulus.

We will illustrate how the EMI model can be used to compute the ephaptic coupling between two ball–and–stick neurons embedded in an extracellular space. The simulation is based on the neuronmi package detailed above. One of the neurons is stimulated with a synaptic input which elicits an action potential. The intracellular potential in the other neuron is sampled. We ran several experiments increasing the distance between the neurons.

Fig. 7.1: Two ball–and–stick neurons embedded in an extracellular space (left) and the increase in the intracellular potential due to ephaptic coupling (right).

The intracellular potential is sampled in the centre of the ephaptically stimulated neuron (right) to measure the strength of the synaptic coupling. The deflection is 4.7 μV when the neurons are 5 μm apart and decreases to 3.2 μV when they are 40 μm apart with the soma of the stimulated neuron adjacent to the axon of the other neuron.

While in these simulations we only showed the effect of a single spike on an adjacent neuron, the occurrence of synchronous activity in populations of neurons can cause larger degrees of ephaptic coupling. The neuronmi package enables users to instantiate several neurons in the mesh and to define their inputs (which can be synchronous), hence allowing in principle to investigate ephaptic effects at the population level.

# **7.4 Investigating the Effect of Measuring Devices on Extracellular Potentials**

While the presence of recording devices is usually ignored in the computation of extracellular potentials, recent findings suggest that newly developed silicon-based devices, or Multi-Electrode Arrays (MEAs), have a strong effect on the measured signals (3).

Using neuronmi, which is built on the EMI model, one can easily incorporate the neural devices in the mesh and investigate their effects on the recorded signals. To demonstrate this, we built meshes with a simple ball-and-stick neuron and different types of neural probes in its vicinity:


We ran EMI simulations using the meshes with and without the probe in the extracellular space (Figure 7.2A-B-C show the meshes with the probe), and we compared the obtained extracellular action potentials - EAPs (Figure 7.2D-E-F, blue without probe - orange with probe). The probes were placed in the extracellular space at a distance of 30 μm from the center of the neuron.

Fig. 7.2: Effect of different probes on the recorded potentials. (A-B-C) Meshes including a neuron and a microwire (A), a neuron and a Neuronexus probe, a neuron and a Neuropixels probe (C). (D-E-F) Extracellular action potentials computed at the electrodes' location without (blue) and with the probe (orange) in the extracellular space. Large MEAs seem to strongly affect the recorded signals (E-F).

Microwire probes do not affect the recorded potentials, with an EAP peak of −21.63 μV without the probe and of −20.53 μV with the probe (Figure 7.2D). However, when recording with silicon MEAs, the extracellular potentials are strongly affected. For the Neuronexus probe (Figure 7.2D), the EAP peak without the probe is −30.47 μV, while with the probe it is −56.09 μV (peak ratio=1.84). For the Neuropixels probe (Figure 7.2E), the EAP peak without the probe is −32.73 μV, while with the probe it is −63.63 μV (peak ratio=1.94). The probe effect is probably due to the insulating properties of the silicon probes, which act as *electrical walls* for the generated currents. For further details and analysis we refer to (3).

# **7.5 Reduced EMI Model**

In the examples considered thus far the problem geometry was simple allowing for computations on a serial desktop computer. In order to apply the EMI framework to realistic neurons two challenges need to be addressed: representation of neurons in the form of a finite element mesh and efficient solvers for large linear systems due to the EMI equations. However, even with the order optimal algorithms discussed in (12, Chapter 6) and efficient mesh generators for neuron surface geometries, see e.g. (17), the computational cost of the 3*D*-3*D* EMI models remains large (compared to the conventional approaches). As a computationally feasible alternative we shall next discuss the 3*D*-1*D* models.

Topological order reduction is a modeling technique used e.g. in reservoir simulations (4) or studies of tissue perfusion (5), which exploits geometrical properties of the system in order to derive its reduced model. Viewing a dendrite (branch) as generalized cylinder with length *L* and radius *R* we observe that *R L* and that in a typical domain of interest the neuron's volume is negligible compared to its surroundings. This property motivates a reduced representation of the neuron in terms of a (one dimensional) curve, the centerline, along with a function *R*, which provides radius of the crossection at each point of the line. An illustration of the concept can be seen in Figure 7.3. Thus, referring to the notation of Chapter 5, Ω<sup>i</sup> is reduced to a line while Ω<sup>e</sup> newly occupies the entire domain, i.e. Ω = Ωe. In turn, the reduced EMI model presents a coupling between three dimensional extracellular space and the one-dimensional intracellular space. We remark that the membrane is one-dimensional as well.

In order to apply order reduction to the EMI model we consider the singledimensional primal formulation (5.6). Note that therein, the coupling on the membrane Γ requires that both *u*<sup>e</sup> and *u*<sup>i</sup> are restricted from Ω<sup>e</sup> and Ω<sup>i</sup> respectively by the dedicated trace operator. In a reduced model Ω<sup>i</sup> = Γ, Γ = Λ and thus *u*<sup>i</sup> needs not to be restricted. On the other hand, restriction of *u*<sup>e</sup> to *curve* Λ can no longer be realized as a trace since such an operation is not well-defined for *H*<sup>1</sup> functions, see e.g. (5). To define a value of the extracellular potential on Λ let us introduce an averaging operator

$$\left| \Pi u(\mathbf{x}) = \vec{u}(\mathbf{x}) = \left| C\_R(\mathbf{x}) \right|^{-1} \int\_{C\_R(\mathbf{x})} u(\mathbf{y}) \, \mathrm{d}S, \quad \mathbf{x} \in \Gamma, \mu \in H^1(\Omega).$$

Here *<sup>C</sup>*R(*x*) <sup>=</sup> {<sup>y</sup> <sup>∈</sup> <sup>Γ</sup>; (<sup>y</sup> <sup>−</sup> *<sup>x</sup>*) ⊥ *<sup>n</sup>*(*x*)} with *<sup>n</sup>*(*x*) being the unit tangent vector of <sup>Λ</sup> at *x*, cf. Figure 7.3. Thus, Π*u*<sup>e</sup> is computed by sampling *u*<sup>e</sup> on the original (two-

Fig. 7.3: Reduced representation of neurons. (Left) Neuron Ω<sup>i</sup> is collapsed to centerline Λ. Extracellular potential defined on Ω is evaluated at Λ by averaging *u*<sup>e</sup> on the cylindrical surface Γ. (Right) Surface mesh and centerline representation of RatS1-6-39 neuron generated by AnaMorph (17). The dendritic part satisfies *R L* property, while the reduced model assumptions do not hold on the soma.

dimensional) neuron surface. However, in practical computations we assume that Γ has a circular cross section so that |*C*R(*x*)| = 2π*R*(*x*).

Using Π the reduced single-scale primal formulation of the EMI model reads: Find *u*<sup>e</sup> ∈ *H*<sup>1</sup> ΓD (Ω), *<sup>u</sup>*<sup>i</sup> <sup>∈</sup> *<sup>H</sup>*1(Λ) such that for all <sup>v</sup><sup>e</sup> <sup>∈</sup> *<sup>H</sup>*<sup>1</sup> ΓD (Ω), <sup>v</sup><sup>i</sup> <sup>∈</sup> *<sup>H</sup>*1(Λ)

$$\begin{split} \int\_{\Omega} \sigma\_{e} \nabla u\_{e} \cdot \nabla v\_{e} \, \mathrm{d}x + \int\_{\Lambda} 2\pi R \frac{C\_{m}}{\Delta t} \bar{u}\_{e} \bar{v}\_{e} \, \mathrm{d}s - 2\pi R \int\_{\Lambda} \frac{C\_{m}}{\Delta t} \bar{v}\_{e} u\_{i} \, \mathrm{d}s \\ &= \int\_{\Lambda} 2\pi R f \bar{v}\_{e} \, \mathrm{d}s \\ - \int\_{\Lambda} 2\pi R \frac{C\_{m}}{\Delta t} v\_{i} \bar{u}\_{e} \, \mathrm{d}s + \int\_{\Lambda} \pi R^{2} \sigma\_{i} \nabla u\_{i} \cdot \nabla v\_{i} \, \mathrm{d}s + \int\_{\Lambda} 2\pi R \frac{C\_{m}}{\Delta t} u\_{i} v\_{i} \, \mathrm{d}s \\ &= - \int\_{\Lambda} 2\pi R f \bar{v}\_{i} \, \mathrm{d}s. \end{split} (7.1)$$

e

Here, the factors 2π*R* arise in reducing the integration domain from Γ to Λ and similarly for π*R*<sup>2</sup> and Ωi. Thus, defining the reduced specific membrane capacitance *C*<sup>m</sup> = 2π*RC*<sup>m</sup> and the reduced intracellular conductivity σ<sup>i</sup> = π*R*2σ<sup>i</sup> the operator form of (7.1) can be seen to be (6.2) with the new restriction operators *T*<sup>i</sup> = *I* and *T*<sup>e</sup> = Π. Note also that in (7.1) the function *f* , which characterizes the membrane dynamics, is defined on Γ.

For the proof of well-posedness of (7.1) as well as a detailed discussion of modeling assumptions, which allow for model reduction from (5.6) the reader is referred to (4). Moreover the reduced multi-dimensional primal formulation (5.8) is analyzed in (13).

e

To assess the reduced model (7.1) the surface mesh of a rat neocortex neuron RatS1-6- 39 from the NeuroMorpho database (2) has been generated using AnaMorph (17). The neuron has been embedded into a (box-shaped) domain Ω such that |Ω|/|Ω<sup>i</sup> | ∼ 100 with the resulting geometry meshed by gmsh. The full 3*D*-3*D* single-dimensional primal formation has been used to compute the response to a 1 ms synaptic stimulus of 50 nA. Referring to Figure 7.4 the lower branch close to node number 4 has been stimulated. Using the centerline representation of the neuron the response has also been computed with (7.1). We remark that *P*1-*P*<sup>1</sup> elements were used with both formulations and that spaces were setup on conforming meshes. In particular, with (7.1) the discretization of Λ consisted of the edges of the cells of Ω. However, such a geometrically conforming discretization is not required in the reduced model. In fact, the meshes of Λ and Ω can be independent, see (4). The reduced model then resulted in 4587 unknowns, which is to be contrasted with 18248 unknowns due to (5.6). In turn, the simulation time using the reduced model is about 110 seconds while the full EMI model required cca. 340 seconds to complete.

The two models are compared in Figure 7.4 which shows values of the computed intracellular potentials at different points along the centerline. In general, there is qualitative agreement between the model predictions. However, the reduced model can be seen to tend to underestimate both the minima and the maxima, while the excitation occurs faster compared to the full model. More precisely, the peak potentials due to the full model at points 2-5 were {27.64, 26.09, 22.87, 12.64} mV with occurrences after {2.08, 2.11, 2.15, 2.57} ms. For the reduced model the maxima {22.31, 21.65, 19.39, 16.02} mV were recorded at {1.56, 1.58, 1.51, 1.86} ms. In addition to the intracellular potentials, the extracellular potentials were compared by sampling 6.52μm away from the soma center (node 1 in Figure 7.4). We remark that the soma radius was 5.71μm. It can be seen that with −2.99μV < *u*<sup>e</sup> < 1.07μV for (7.1) and −1.99μV < *u*<sup>e</sup> < 0.67μV for (5.6) the reduced model overestimates the extrema. As with the extracellular potentials there is a temporal shift in the response; the negative peak is observed at 1.26ms, respectively 1.87ms.

While results of the comparison in Figure 7.4 shall be viewed as preliminary we argue that they illustrate sufficiently the potential of reduced EMI models. In particular, the reduced model is able to capture qualitatively the properties of the full EMI simulations. However, clear differences, especially in the temporal shifts of the peaks, have been observed. In the future we aim to investigate if suitable scaling of the stimulus and/or the parameters of the membrane ODEs can reduce the prediction error. In addition, the modeling error of the reduced model shall be evaluated similar to (15). More specifically, the soma, being approximated as a sphere in the 3*D*-3*D* model, cannot be represented as a slender cylinder (unlike the dendrites and axons). Thus the model reduction assumptions are not met on the soma. While localized in space, the effect of this error on temporal predictions should be analyzed. In turn, improved reduced models, which take into account the spherical nature of the soma, e.g. in construction of averaging operators, might be needed.

Fig. 7.4: Comparison of full 3*D*-3*D* and reduced 3*D*-1*D* EMI models. RatS1-6-39 neuron is stimulated close to node 4 in the dendritic part (plotted in blue). Intracellular potentials (nodes 2-5) are measured on the centerline with node 2 being the soma (in red, dashed line indicates the radius) center. Extracellular potentials are compared in node 1 next to soma. Axis of each response plot is anchored at the center next to the measurement point. The full and reduced EMI models provide qualitatively similar predictions.

# **7.6 Conclusions**

In this chapter, we have showcased some applications in which the EMI model can be a viable alternative to standard modeling frameworks in order to investigate aspects of neuronal activity and recorded signals. We introduced an open-source software package named neuronmi to easily assemble meshes including simplified neurons and probes in the extracellular space. Finally, we performed preliminary investigations into accuracy and solution cost of a reduced 3*D*-1*D* EMI model suitable for simulating complex neuronal morphologies.

**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.

# **References**


# **Index**

#### **A**

active force 28 active strain 26 active stress 26, 27 approximation 59

#### **C**

cable equation 77 capacitive current 6 cardiomyocyte 25, 27 Cauchy stress tensor 32 conjugate gradient 62 convergence 56 cross-fiber 32

#### **D**

deformation 32 deformation gradient 26, 27 displacement 28 domain decomposition 72

## **E**

ephaptic 77 ephaptic effect 79

#### **F**

fiber direction 28, 32 finite element method 48

#### **G**

gap junctions 10

Green-Lagrange strain tensor 29

#### **H**

Hodgkin-Huxley 16

#### **I**

intercalated disc 9 ionic currents 9 ischemia 42

#### **L**

Lagrange multiplier 53

#### **M**

Maxwell's equations 3 minimal residual 62 mixed formulation 50 mortar 53

#### **N**

Na/K/ATPase pump 17 Nernst potential 16

## **O**

Ohm's law 4 operator splitting 40

#### **P**

Piola-Kirchhoff stress tensor 26 preconditioner 62 primal formulation 51

The Author(s) 2021 99 A. Tveito et al. (eds.), *Modeling Excitable Tissue*, Simula SpringerBriefs on Computing 7, https://doi.org/10.1007/978-3-030-61157-6

## **R**

reduced model 76

#### **S**

splitting algorithm 41 strain energy function 26–28 stress 26, 28

**T**

transmembrane potential 49

## **W**

weak form 30