Joakim Sundnes Glenn Terje Lines Xing Cai

Bj0rn Fredrik Nielsen Kent-Andre Mardal Aslak Tveito

Simula Research Laboratory P.O. Box 134 1325 Lysaker, Norway email: [email protected] [email protected] [email protected] [email protected] [email protected] [email protected]

Aslak Tveito has received financial support from the NFF - Norsk faglitter^r forfatter- og oversetterforening

Library of Congress Control Number: 2006927369

Mathematics Subject Classification: 35Q80, 65M55, 65M60, 92C50, 92C0510

ISBN-10 3-540-33432-7 Springer Berlin Heidelberg New York ISBN-13 978-3-540-33432-3 Springer Berlin Heidelberg New York

This work is subject to copyright. All rights are reserved, whether the whole or part of the material is concerned, specifically the rights of translation, reprinting, reuse of illustrations, recitation, broadcasting, reproduction on microfilm or in any other way, and storage in data banks. Duplication of this publication or parts thereof is permitted only under the provisions of the German Copyright Law of September 9, 1965, in its current version, and permission for use must always be obtained from Springer. Violations are liable for prosecution under the German Copyright Law.

Springer is a part of Springer Science+Business Media springer.com

© Springer-Verlag Berlin Heidelberg 2006 Printed in The Netherlands

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

Typesetting: by the authors and techbooks using a Springer LTjeX macro package Cover design: design & production GmbH, Heidelberg

Printed on acid-free paper SPIN: 11737407 46/techbooks 543210


The heart is a fantastic machine; during a normal lifetime it beats about 2.5 billion times and pumps 200.000 tons of blood through an enormous system of vessels extending 160.000 kilometres throughout the body.

For centuries, man has tried to understand how the heart works, but there remain many unsolved problems, problems that have captured the attention of thousands of researchers worldwide. There is, for example, a huge amount of research being devoted to the analysis of single heart cells. Other areas of research include trying to understand how it works as a complete muscle, and how blood flows through the heart. The entire process is extremely complex.

The history of bioelectricity can be traced back to the late eighteenth century and the experiments of Luigi Galvani. A century later, in 1887, Augustus Wallers managed to measure the electrical signal generated by the heart at the surface of the body [142]. His dog Jimmy earned a place in history by being the first to have his heart measured in this way; see Figure 1.1. In 1903 Willem Einthoven [34] developed the first commercial device for recording electrocardiograms (ECGs); see Figure 1.2.

This book is about computing the electrical activity in the heart. In order to do so, we will need mathematical models of how the electrical signals are generated in the heart. Furthermore, we will need mathematical models of how the signals are transported through the heart and distributed in the body. All these models are formulated in terms of differential equations. Based on these models, we will derive discrete models suitable for numerical simulations. We will base our methods on a finite element approach and will solve the equations that arise on parallel computers.

In order to understand this book, you will need to have a basic course in partial differential equations, and it is definitely an advantage to know a bit about finite element methods; but that is about it. We will explain all the biology you will need. We will provide a detailed introduction to parallel computing and how to solve linear systems of equations.

This book is about computational science; in fact, its main objective is to present a large project in computational science. We will begin by giving you the necessary physiological background. What is going on in the heart and how can we use a recording on the body surface to say anything about the condition of the heart? These questions are discussed briefly in Chapter 1.

In Chapter 2 we jump straight to the mathematical models. A wonderful feature of mathematics is its ability to phrase extremely complex phenomena in rather simple equations. Yet despite this ability, the equations we describe in Chapter 2

are rather complicated. These models are absolutely necessary in order to be able to compute the electrical activity on a computer. We guide you through the basic physics of the process, ranging from models of what happens within a single cell to a complete mathematical model of the electrical activity in the entire heart.

In Chapter 3 we discuss how to discretize the equations derived in Chapter 2. We will do this using the finite element method. The approach is straightforward, but matters become complicated due to the complexity of the mathematical models. Chapter 3 will motivate three specific problems that have to be considered carefully. First, the finite element method leads to large systems of linear equations. It is well known from other fields that these equations may be very challenging to solve, and in Chapter 4 we introduce and analyze preconditioning techniques to solve the linear system that arises from the finite element discretizations. Second, systems of ordinary differential equations that model the processes going on in single cells have to be solved. We discuss appropriate methods for this in Chapter 5. Third, the problem that we are going to solve is huge, and therefore has to be solved on parallel computers. We introduce such computers in Chapter 6 and demonstrate their usefulness for the problems at hand.

Our primary goal is to contribute to the development of machinery that can compute the location and geometry of a myocardial infarction. This is a huge task occupying the minds of many researchers around the world, but a technologically feasible solution still lies in the future. In Chapter 7 we initiate a discussion of how to utilize our ability to simulate the electrical activity in the heart and the body in order to detect an infarction. This is, in fact, an inverse problem; we measure the result of a process and we then compute its cause. In the present setting, the measurements are electrocardiograms and the cause is the infarcted area of the heart. The process is the transport of electrical signals in the heart and the distribution of these in the body.

We hope this book can serve as an introduction to this field for applied mathematicians and computational scientists and also for researchers in bioengineering who are interested in using the tools of computational science in their research. The book presents the authors' view of the field, with a strong emphasis on the mathematical models and how to solve them numerically. There are few new scientific results in the book. Most of the material presented here is based on already published material.

Fornebu, Norway Joakim Sundnes

March 2006 Glenn Terje Lines

Xing Cai Bj0rn Fredrik Nielsen Kent-Andre Mardal andAslak Tveito

Table of Contents

1 Physiological Background 1

1.1 The Electrocardiogram 2

1.1.1 Physics and Physiology 5

1.1.2 Cellular Electrical Activity 9

1.1.3 Signal Conduction 11

1.1.4 Diagnosis and the Inverse Problem 12

1.2 Computer Simulations 14

1.2.2 State of the Art 15

2 Mathematical Models 21

2.1 Modelling the Body as a Volume Conductor 21

2.1.1 A Volume Conductor Model 21

2.2 A Model for the Heart Tissue 24

2.2.1 Excitable Tissue 24

2.2.2 The Bidomain Model 25

2.2.3 The Monodomain Model 30

2.3 Coupling the Heart and the Body 32

2.4 Models for the Ionic Current 33

2.4.1 The FitzHugh-Nagumo Model 34

2.4.2 The Cell Membrane 36

2.4.3 The Nernst Equilibrium Potential 38

2.4.4 Models for Ionic Flux 40

2.4.5 Channel Gating 42

2.4.6 The Hodgkin-Huxley Model 44

2.4.7 A Model for Cardiac Cells 46

2.4.8 Models for Ventricular Cells 48

2.4.9 Second Generation Models 51

2.5 Summary of the Mathematical Model 53

3 Computational Models 57

3.1 The Finite Element Method for the Torso 57

3.1.1 A Simplified Model Problem 57

3.1.2 A Dipole Model of the Heart 63

3.1.3 Known Potential on the Heart Surface 67

3.2 The Heart Equations 70

3.2.1 Operator Splitting 71

3.2.2 Operator Splitting for the Monodomain Model 75

3.2.3 Operator Splitting for the Bidomain Model 78

3.3 Coupling the Heart and the Torso 82

3.4 Numerical Experiments 85

3.4.1 Convergence Tests 85

3.4.2 Simulation on a 2D Slice 92

3.4.3 Normal Propagation 93

3.4.4 Ischemia 94

4 Solving Linear Systems 99

4.1 Overview 99

4.2 Iterative Methods 100

4.2.1 The Richardson Iteration 101

4.2.2 The FDM Discretization Poisson Equation in 1D 103

4.2.3 The Richardson Iteration Revisited 106

4.2.4 Preconditioning 109

4.2.5 The Jacobi Method 110

4.2.6 The Relaxed Jacobi Method 112

4.2.7 The Exact and the Inexact Block Jacobi Methods 113

4.2.8 The Gauss-Seidel Method 114

4.2.9 The Relaxed Gauss-Seidel Method 115

4.2.10 The Symmetric Gauss-Seidel Method 115

4.2.11 The Exact and the Inexact Block Gauss-Seidel Method 115

4.3 The Conjugate Gradient Method 116

4.3.1 The CG Algorithm 116

4.3.2 Convergence Theory 120

4.3.3 Numerical Experiments 121

4.4 Multigrid 123

4.4.1 Idea 123

4.4.2 Theoretical Framework 125

4.4.3 Convergence Theory 127

4.4.4 Numerical Experiments 128

4.5 Domain Decomposition 128

4.6 Preconditioning Revisited 133

4.6.1 Idea 133

4.6.2 Spectral Equivalence 133

4.6.3 The Richardson Iteration Re-Revisited 134

4.6.4 Preconditioned Conjugate Gradient Method 135

4.6.5 Convergence Analysis Revisited 136

4.6.6 Variable Coefficients 137

4.6.7 Numerical Experiments 138

4.7 The Monodomain Model 140

4.7.1 Multigrid 141

4.7.2 Numerical Experiments 141

4.7.3 Domain Decomposition 142

4.8 The Bidomain Model 142

5 Solving Systems of ODEs 149

5.1 Simple ODE Solvers 149

5.1.1 The Euler Methods 149

5.1.2 Stability Analysis for the Euler Methods 151

5.2 Higher-Order Methods 154

5.2.1 Multistep Methods 154

5.2.2 Runge-Kutta Methods 156

5.3 Solving Nonlinear Equations 160

5.3.1 Newton's Method 161

5.3.2 Newton's Method for Higher-Order Solvers 162

5.4 Automatic Time Step Control 164

5.5 The Cell Model Equations 168

5.5.1 Explicit versus Implicit Methods 168

5.5.2 Simulation Results 171

6 Large-Scale Electrocardiac Simulations 175

6.1 The Electrocardiac Simulator 175

6.1.1 The Numerical Strategy 176

6.1.2 Software Components of the Electrocardiac Simulator 176

6.2 Requirements for Large-Scale Simulations 177

6.2.1 The Memory Requirement 178

6.2.2 Realistic Estimates for Memory and Time Usage 179

6.3 Introduction to Parallel Computing 181

6.3.1 Hardware and Programming Models 181

6.3.2 Division of Work and the Resulting Overhead 182

6.3.3 Speedup and Parallel Efficiency 184

6.4 Two Simple Examples 186

6.4.1 Adding Two Vectors 186

6.4.2 Inner Product 188

6.5 Domain-Based Parallelization 191

6.5.1 Division of Data for FDM 192

6.5.2 Division of Data for FEM 194

6.5.3 Summary of Principles 196

6.6 Explicit FDM in Parallel 197

6.6.1 Discretization by Finite Differences 197

6.6.2 A Sequential Program 198

6.6.3 Parallelization 200

6.7 Parallel Conjugate Gradient Iterations 201

6.7.1 Conjugate Gradient Revisited 202

6.7.2 Parallel CG and FDM 203

6.7.3 Parallel CG and FEM 203

6.8 Domain Decomposition as Parallel Preconditioners 207

6.8.1 Additive Schwarz Preconditioner 208

6.8.2 Parallel DD for the Monodomain Equation 210

6.8.3 Parallel DD for the Bidomain Equations 210

6.8.4 Extension to the Torso 211

6.9 Parallelizing Electrocardiac Simulations 212

6.9.1 The Overall Simulation Process 212

6.9.2 Partitioning the Domains 213

6.9.3 Straightforward Parallelization Tasks 214

6.9.4 Solving the Block Linear System in Parallel 215

6.10 Simulation on a Realistic Geometry 215

6.11 Summary 216

7 Inverse Problems 219

7.1 A Simple Example 221

7.1.1 Problem Formulation 222

7.1.2 Fourier Analysis 224

7.1.3 Ill-Posedness 225

7.1.4 Discretization and an Output Least Squares Formulation of the Problem 229

7.1.5 Regularization Techniques 231

7.2 The Classical Inverse Problem of Electrocardiography 238

7.2.1 Mathematical Formulation 239

7.2.2 A Linear Problem 241

7.2.3 Discretization 243

7.2.4 The Time-Dependent Problem 246

7.2.5 Numerical Experiments 249

7.3 Computing the Location and Orientation of a Dipole 253

7.3.1 Preliminaries 255

7.3.2 The Inverse Problem 256

7.3.3 Parameterizations Leading to Linear Problems 257

7.3.4 A Numerical Example 261

7.3.5 Parameterizations Leading to Nonlinear Problems 265

7.3.6 A Numerical Experiment 270

7.4 Computing the Size and Location of a Myocardial Infarction 271

7.4.1 Modelling Infarctions, the Direct (Forward) Problem 272

7.4.2 Modelling Infarctions, the Inverse Problem 274

7.4.3 Differentiation of the Objective Function 276

7.4.4 Numerical Experiments 281

A Color Figures 287

B Rate Functions and Ionic Currents 297

B.1 The Hodgkin-Huxley Model 297

B.2 The Noble Model for Purkinje Cells 297

B.3 The Beeler-Reuter Model 298

C Coefficients for the Implicit RK Solvers 300

Bibliography 301

Index 308

Chapter 1

0 0

Post a comment