A SECOND ORDER SCHEME FOR SOLVING THE COUPLED BIDOMAIN AND FORWARD PROBLEM

J. Sundnes, G. T. Lines, A. Tveito
Simula Research Laboratory,
P.O. Box 134 Lysaker, N-1325, Norway

Abstract: Using the bidomain model to simulate ECG recordings requires the model to be coupled to an equation describing the propagation of the electrical signal in the torso. Realistic simulations also require ODE-based models for the electrical properties of the cardiac cells. The result is a complex non-linear system of PDEs and ODEs, for which it is challenging to construct efficient solution algorithms. In this paper we combine an operator splitting technique with a fully coupled solution of the PDEs, to obtain a second order accurate discretization of the complete system.

INTRODUCTION

The electrical activity in the myocardium may be modeled by the bidomain model[1]. The main variables of the model are the transmembrane potential v and the extracellular potential u_e. To simulate the electrical potentials on the body surface, the model is coupled to an equation modeling the electrical potential in the torso, denoted by u_o. Solving this equation to determine the body surface potential from the surface potential on the heart is referred to as the forward problem. If we denote the heart muscle by H and the surrounding torso by T, the coupled bidomain and forward problem may be formulated as follows:

Here, M_i and M_e are the conductivity tensors for the intracellular and extracellular domain, respectively, and M_o is the conductivity tensor in the torso. The term I_ion(v,s) is a non-linear function describing the transmembrane flow of ions. The variable s is a vector characterizing the conductive properties of the cell membrane. The variations in s are described by (1), which is a non-linear system of ordinary differential equations (ODEs).

To perform accurate simulations of the electrical activity in the heart nd the torso, (1)-(8) should ideally be solved fully coupled. However, because of the complexity of the system, most solution approaches have utilized some form of splitting of the equations, which limits the accuracy of the solution.

In [2] it is described how a Strang splitting technique[3] may be applied to a simplified monodomain model for the electrical activity. The idea is to separate the ODEs from the PDEs and still obtain second order accuracy of the numerical scheme. The purpose of the present paper is to extend these ideas to the complete, coupled bidomain and forward problem. We combine the Strang splitting with a second order accurate fully coupled solution of the PDEs, to obtain a second order accurate time discretization of (1)-(8).

METHOD

By applying a Strang splitting technique to (2), we obtain a sequential solution algorithm for the complex system. Assuming that   vⁿ = v(t_n)   and sⁿ = s(t_n) are known, one time step involves the following operations.

1. Solve the system

   

   for  t_n < t £ t_n+1/2, to obtain s^n+1/2 and the temporary value v^*.
   
   

2. Solve the coupled PDE system

   

   for  t_n < t £ t_n+1, with v(t_n)=v^*. We obtain u_e and u_o, as well as a new temporary value v^** corresponding to time t_n+1/2 for the ODE operator and t_n+1 for the PDE operator.
   
   
3. Compute sⁿ⁺¹ and vⁿ⁺¹ by solving the system

   

   for t_n+1/2 < t £ t_n+1.
   
   

If the sub-problems in steps 1, 2, and 3 are solved sufficiently accurate, this splitting algorithm will give second order accuracy. We have chosen to use a 3rd order semi-implicit Runge-Kutta method for the ODE systems, and a Crank-Nicholson scheme for the PDEs. The Runge-Kutta method chooses the time step adaptively, so that steps 1 and 3 may consist of several steps of the ODE solver. We do not use adaptive time stepping for the PDEs, and Step 2 is hence one Crank-Nicholson step, approximating v at time steps t_n, t_n+1, ¼, and u_e and u_o at steps t_n+1/2, t_n+3/2, ¼.

The PDE system may be simplified by using the continuity conditions (5) and (6) to combine (3) and (7) to one equation. We introduce the new field u, defined over HÈT, which is given by

u = ì í î ue for x Î H uo for x Î T ,

and introduce the inner products

ai(v,f) = ó õ H  MiÑv ·ÑfdH, ai+e(v,f) = ó õ H  (Mi+Me)Ñv ·ÑfdH, ao(v,f) = ó õ T  MTÑv·ÑfdT .

For appropriate function spaces V(H) and V(HÈT), a weakly formulated Crank-Nicholson discretization of (9) - (11) is to find v Î V(H) and u Î V(HÈT) such that

A standard finite element discretization of this system gives a linear system of the form

 

A é ê ë v u ù ú û = é ê ë A B BT C ù ú û é ê ë v u ù ú û = é ê ë a b ù ú û ,

with blocks corresponding to the terms in (12) - (13). This system is solved with the conjugate gradient (CG) method and a block multigrid preconditioner, as described in [4].

RESULTS

The solution technique outlined above has been used to perform experiments in two and three space dimensions. Results from simulations on a 2D domain similar to a cross section of the torso and the heart are shown in Figure 1. The left panel shows the depolarized region of the heart after 40 ms, while the right panel shows a simulated ECG recording. Although the amplitude of the signal is too large, the simulator captures the essential properties of the ECG signal.

Experiments have shown that the algorithm for solving the discretized PDEs is close to optimal, in that the required number of CG iterations increases only slightly as the number of unknowns is increased.

Figure 1: Left: The Transmembrane potential after 40 ms. The gray regions represent depolarized tissue. Right: ECG recording from the 2D simulation.

DISCUSSION

The main purpose of the present work was to introduce a second order accurate time discretization for the coupled bidomain and forward problem. Although we obtain reasonable results, the accuracy of the scheme remains to be verified by numerical experiments. Fot the complex problem at hand such a verification is a non-trivial task, because no analytical solution is available. However, by considering circular or spherical geometries with rotationally invariant solutions, it should be possible to produce numerical reference solutions of sufficient accuracy. Such a reference solution may also be used to compare the performance of the simulator to other solution algorithms. If implemented correctly the increased accuracy of the Strang splitting comes at a very low computational cost, but numerical expreiments are still required to determine if the approach presented here is in fact better than other approaches, in terms of computational work compared to accuracy.

REFERENCES

[1] L. Tung, A Bi-domain model for describing ischemic myocardial D-C potentials. PhD thesis, MIT, Cambridge, MA, 1978.

[2] Z. Qu and A. Garfinkel, "An Advanced algorithm for solving partial differential equation in cardiac conduction", IEEE Transactions on Biomedical Engineering, vol 46(9), pp. 1166-1168, 1999. [3] G. Strang, "On the construction and comparison of difference schemes", SIAM Journal of Numerical Analysis, vol 5, pp. 506-517, 1968. [4] J. Sundnes, G.T. Lines, K.-A. Mardal and A. Tveito, Multigrid block preconditioning for a coupled system of partial differential equations modeling the electrical activity of the heart. Research report, Simula Research Laboratory, 2002.