Computing the electrical activity in the human heart

1

• Computing the electrical activity in the human
  heart
• Aslak Tveito
• Glenn T. Lines, Joakim Sundnes, Bjørn Fredrik
  Nielsen, Per Grøttum, Xing Cai, and Kent Andre
  Mardal
• Simula Research Laboratory/
• The University of Oslo

2

• Motivation/Background
• The mathematical model
• Solution strategy
• Solving linear systems
• CPU-models
• Results
• Future work

3
The heart
For each heartbeat, electrically charged atoms
(ions) move in and out of the heart muscle cells
in a complicated pattern.The most important ions
are Na, K, Ca, and Cl-. The movement of the
ions cause the muscle fibers to contract and the
heart to pump. Normal electrical activity is
important for the pumping function of the
heart. The electrical activity of the heart may
be recorded on the surface of the body. The
recording is called an electrocardiogram (ECG).
K
K
K
K
K
Na
Na
Na
Na
Na
4
ECG
1943 The lead positions are standardized.
1887 The first ECG is recorded in London on
Augustus Wallers dog Jimmy.
2000 Worldwide, approximately 1 million ECGs are
recorded every day.
1911 A commercial ECG machine is constructed by
Willem Einthoven.
5
Heart infarctions
20 of deaths in the western world are due to
hearth infarctions and consequences
thereof. About 50 of patients admitted to
surveillance units with acute chest pain suffer
from a heart infarction. ECG is the most
important tool to diagnose heart infarction and
heart cramp (angina). In some parts of the
heart, the sensitivity of ECG is as low as
60. Normal interpretation of ECG only gives
crude estimates large infarction , small
infarction.
6

7
Computational domain sketch in 2D
8
Cardiac electrical activity

• The heart consists of billions of electrically
  charged cells.
• During a heartbeat these cells leak ions, which
  changes the polarity of the cells
  (depolarization).
• The electrocardiogram records this process on the
  body surface.
• Pathological conditions in the heart can be
  diagnosed from abnormal deviations in the ECG.

9

• Electrical potential in the torso
• J current destiny
• u electrical potential
• M conductivity tensor
• Conservation of current,
• and Ohms law
• gives
• and, in addition

10
Electrical potential in the heart

• Intracellular (space within the cells)
• Ji current density
• ui electrical potential
• Mi conductivity tensor
• qi electrical charge
• Extracellular (space outside the cells)
• Je current density
• ue electrical potential
• Me conductivity tensor
• qe electrical charge
• Transmembrane potential
• vui - ue

11

• Conservation laws for intra and extra cellular
  domains
• where Iion models the ionic current.

12

• From these two equations, the BiDomain model
  is derived
• The BiDomain model was developed by
  Gesolowitz, Miller, Schmilt, and Tung in the
  early 70s. The equations have been studied by a
  series of researchers (Colli Franzone et al,
  Henriques, Trayanova et al, Huang et al, )

13
Model for ionic current (Iion)

• The appropriate model for the ionic current
  will depend upon the application.
• In the simplest type is is just a function of
  v
• In more realistic models it depends on
  several factors
• where s is a vector of variables including
  concentrations of ions and permeability of ion
  channels. These variables are typically governed
  by an ODE system

14
The BiDomain Model
15
Boundary conditions
16
Cell models of ionic currents

• Hodgkin, huxley
• Noble
• Beeler, Reuter
• DiFrancesco, Noble
• Luo, Rudy
• Winslow et al

17
(No Transcript)
18
The BiDomain Model
19
Operator splitting

• Split in two parts
• with appropriate boundary conditions.

20

• Consider the linear problem
• with appropriate boundary conditions.

21
(No Transcript)
22

• The discretization leads to a linear system of
  equations on the form
• where A is a matrix, uh is the FEM-solution,
  and fh is the discretized version of f. Let N be
  the number of grid points.
• We solve this system using the CG-method and
  the MG-method. The stopping criterion is
  in a relative l2-norm on the residual.

23

• Numerical experiments on an IBM-RS6000 show
  that the CPU-time needed by the CG and UG methods
  are
• and
• Suppose we want to solve (2) on a grid
  relevant for the heart.

24

• Numerical experiments indicate that we need
• In the heart. So for 1cm, we need 50 points.
  In 1cm, we need
• Suppose the heart is 300cm, then we need

25

• Solving one linear system of the form (2)
  using the CG method requires
• Which is about 7.5 hours. The MG method
  requires about
• or about 7 minutes.
• We need order optimal methods!

26
Preconditioners

• The number of iterations for the CG method is
  bounded by
• Since, for elliptic problems,
• the convergence is very slow as seen above.
  However, we may solve
• instead of

27

• The goal of preconditioning is to find B such
  that
• is independent of h and where
• can be solved in O(N) operations. Then the
  number of iterations needed to solve (5) is
  independent of h and thus the entire solution
  process uses only O(N) operations.

28

• It can be shown that if there exist constants
  c0, c1 independent of h, such that
• then the condition number of is
  bounded independently of h.

29
Analysis of a BiDomain preconditioner

• Define
• The system (7) becomes

30
The Crank-Nicholsons scheme

• Weak form
• Define

31
(No Transcript)
32

• Linear operators

33
(No Transcript)
34
Spectral equivalence

• Two operators, Tn and Sn are spectrally
  equivalent if there are constants c1 and c2,
  independent of h, such that
• If Tn and Sn (VngVn) are -symmetric and
  positive definite, then

35
Preconditioner

• Define (from (12))
• and, the preconditioner
• Our aim is to prove that
• where c is independent of h and rt.

36
Assumptions

• We assume that there are constants
  independent of h and rt such that
• The last assumption simply state that the
  matrix generated by the Laplace operator is
  spectrally equivalent to the matrices generated
  by the variable coefficient problems. Note that,
  since
• we have from (18)

37
Upper bound

• We start by proving that
• for a suitable constant c.

38
(No Transcript)
39
(No Transcript)
40
(No Transcript)
41

• Lower bound

42
(No Transcript)
43

• From (22) and (23), we have
• so

44
Numerical experiments

• We want to compare the numerical solution of
  the Bidomain equations using CG without
  preconditioning, and with the optimal
  preconditioner.
• Matrix
• Preconditioner

45

• The table shows average number of iterations
  and CPU-time using 10 time-steps (fixed rt).

46

• Models of CPU-time as a function of the number
  of computational nodes, N
• With we expect

47
Parallel solution using Domain Decomposition

• Number of nodes
• Number of unknowns per timestep
• CPU-measurements for one time-step

48

• How fast can we now solve the entire problem?

49
(No Transcript)