- Email: [email protected]
Optimal monodomain approximations of the bidomain equations
Applied Mathematics and Computation 184 (2007) 276–290 www.elsevier.com/locate/amc
Optimal monodomain approximations of the bidomain equations Bjørn Fredrik Nielsen *, Tomas Syrstad Ruud, Glenn Terje Lines, Aslak Tveito Simula Research Laboratory, P.O. Box 134, 1325 Lysaker, Norway
Abstract The bidomain equations are widely accepted to model the spatial distribution of the electrical potential in the heart. Although order optimal methods have been devised for discrete versions of these equations, it is still very CPU demanding to solve the equations numerically on a sufficiently fine mesh in 3D. Furthermore, the equations are hard to analyze; from a mathematical point of view, very little is known about the qualitative behavior of the solutions generated by these equations. It is well known that upon appealing to a certain relation between the extracellular and intracellular conductivities, the bidomain model can be rewritten in terms of a scalar reaction diffusion equation referred to as the monodomain model. This model is of course much easier to solve, and also the qualitative properties of the solutions are well known; such equations have been studied intensively for decades. It is the purpose of the present paper to show how the bidomain equations can be approximated in an optimal manner by the solution of a monodomain model. The key feature here is that this optimal solution can be computed without solving the bidomain model itself. The solution is obtained by putting the problem into a framework of parameter identification problems. Our results are illuminated by a series of numerical experiments. 2006 Elsevier Inc. All rights reserved. Keywords: Electrical behavior of cardiac tissue; Bidomain and monodomain models; Output least squares
1. Introduction The electrical activity in the heart is governed by a complex, fully coupled system of ordinary and partial differential equations: os ¼ F ðs; vÞ in H ; ot vt þ Iðv; sÞ ¼ r ðM i rvÞ þ r ðM i ruÞ in H ;
ð2Þ
r ðM i rvÞ þ r ððM i þ M e ÞruÞ ¼ 0
ð3Þ
in H :
*
ð1Þ
Corresponding author. E-mail addresses: [email protected] (B.F. Nielsen), [email protected] (T.S. Ruud), [email protected] (G.T. Lines), aslak@ simula.no (A. Tveito). 0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.05.158
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
277
In this model, v and u represent the transmembrane and extracellular potentials, respectively, and H the physical domain occupied by the heart. The intracellular and extracellular conductivity tensors1 Mi and Me are typically defined in terms of symmetric and positive definite matrices depending on the spatial position x. The state s of the ionic currents is incorporated into the model by the system (1) of ordinary differential equations (ODEs) and the possibly nonlinear function I, see [25] or [22] for further details concerning this model. In models of this kind, the partial differential equations (PDEs) (2) and (3) are referred to as the bidomain equations, and (1) denotes a (generic) system of ordinary differential equations. More precisely, s is a state vector for the ions, and (1) incorporates the effect of the flow of ions across the cell membrane. Throughout the last 30 years, a wide variety of such cell models have been proposed and analyzed, see Beeler and Reuter [1], DiFrancesco [3], Luo and Rudy [14,15] and Winslow et al. [27] – see [2] for a very nice overview of these models. For further information about the bidomain equations, we would like to refer to [11] and [22]. Eqs. (1)–(3) define a nonlinear, fully coupled system of differential equations for v, u and s. Analytical solutions are not available for such complex models. They must hence be solved by some kind of numerical procedure. Ideally we would of course like to solve the system (1)–(3) in a fully implicit manner, see [16]. However, due to the complexity of the involved ODEs, see e.g. [26], it is common to use some sort of operator splitting scheme [19,23]. Schemes of this kind are typically defined in terms of successive solutions of the two subsystems os ¼ F ðs; vÞ in H ; ot ov ¼ Iðv; sÞ in H ot
ð4Þ ð5Þ
and vt ¼ r ðM i rvÞ þ r ðM i ruÞ in H ;
ð6Þ
r ðM i rvÞ þ r ððM i þ M e ÞruÞ ¼ 0 in H :
ð7Þ
Further information about such methods can be found in [23,22]. A complete simulator for the electrical activity in the heart consists of solving (4), (5) and (6), (7) at each time step. Computer experiments indicate that a very fine computational mesh is needed in order to solve the subsystem (6), (7) with any reasonable accuracy. Typically, 20–40 millions grid-points are needed in 3D, see [13,28]. This leads to an extremely CPU demanding task, and in many cases a need for considering simplified models. In this paper we will consider one of this simplified approaches – namely, the monodomain model. If M e ¼ kM i
ð8Þ
for some scalar k, then (6) and (7) can be written on the form k r ðM i rvÞ in H ; 1þk ð1 þ kÞr ðM i ruÞ ¼ r ðM i rvÞ in H ; vt ¼
ð9Þ ð10Þ
which commonly is referred to as the monodomain model.2 Note that (9) and (10) can be solved in a sequentially manner; first solve (9) for v and thereafter compute u from (10). In fact, if only the transmembrane potential is of interest, then it is sufficient to solve (9). This is of course a simpler task than determining v and u, at each time step, from (6) and (7). For a further discussion of this issue, we refer to [24]. Although the condition (8) is incorrect from a biological point of view, the simplified model (9) is nevertheless frequently used to simulate the electrical activity in the heart. In order to perform such simulations, a suitable value for k must be determined. The purpose of this paper is to shed some light onto this problem, and to design an algorithm for computing optimal, or at least close to optimal, values for k. 1
More precisely, Mi and Me are scaled conductivities, see [22] for details. Usually (9), without (10), is referred to as the monodomain model. For the sake of completeness, we have also included Eq. (10). In the present context, this turns out to be very useful. 2
278
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
One possible approach towards this problem is to minimize the difference between kMi and Me, measured in proper norms. This leads to a simple procedure for determining a ‘‘feasible’’ value ~k for k. However, by applying such a scheme, it is difficult to control the error introduced by solving the monodomain equation instead of the bidomain model. We can simply not guarantee that ~k actually provides a good choice for this parameter. In this paper we will attack this problem by the Output Least Squares (OLS) method. Inspired by the theory for inverse problems, we propose to determine k in such a way that the difference between the solutions of the bidomain and the monodomain equations is approximately minimized. In general, such a procedure would of course require the solution of the bidomain equations. Thus making such a strategy useless - there is no point in solving (9) if the solution of (6), (7) is known! However, by measuring the differences between the involved functions in suitable norms, it turns out that an approximately optimal value for k can be computed without actually solving the bidomain equations. We will see below, in Section 3, that this naturally leads to a timedependent parameter k, i.e. k = k(t). The bidomain model was originally introduced by Tung [25]. It has been thoroughly studied by several scientists, see [6,8,17,12,7,21,22]. Some researchers have also compared the performance and accuracy of monodomain and bidomain simulations, see [9,24]. In [20] a ‘‘reduced bidomain’’, to quote Skouibine and Krassowska, type of scheme was proposed, and some quite impressive simulation results were obtained. However, we are not aware of any results regarding optimal choices of k. This paper is organized as follows: In Section 2 we consider the relation between the monodomain and bidomain models in the case of one space dimension. Section 3 contains our theoretical findings for multidimensional problems, and in Section 4 we derive a modified and optimal monodomain type of scheme.Finally, our numerical experiments and conclusions are presented in Sections 5 and 6. 1.1. Boundary conditions Throughout this text we will, for the sake of simplicity, consider an insulated heart. This leads to the following set of boundary conditions for the bidomain equations (6) and (7): ðM i rv þ M i ruÞ n ¼ 0 ðM e ruÞ n ¼ 0
along oH ;
along oH :
ð11Þ ð12Þ
Furthermore, the transmembrane potential v, as well as the state of the cells s, at time t = 0 is assumed to be known, i.e. vðx; 0Þ ¼ v0 ðxÞ
for x 2 H ;
sðx; 0Þ ¼ s0 ðxÞ
for x 2 H :
Note that, if (8) is satisfied for some k, then (11) and (12) imply that ðM i rvÞ n ¼ 0 along oH ; ðM i ruÞ n ¼ 0 along oH :
ð13Þ ð14Þ
Hence, for the monodomain model it seems to be reasonable to apply homogeneous Neumann boundary conditions both for the transmembrane and extracellular potentials. 2. The bidomain model in 1D In order to shed some light onto the problem outlined above, let us consider the bidomain equations in one space dimension (1D): vt ¼ ðM i vx Þx þ ðM i ux Þx ;
ð15Þ
0 ¼ ðM i vx Þx þ ððM i þ M e Þux Þx ;
ð16Þ
where Mi = Mi(x) and Me = Me(x) are the conductivities. We assume that both u and v have zero derivatives at the boundaries. Obviously, if there exists a constant k such that (8) holds, then (15), (16) can be reduced to
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
279
the monodomain model as explained above. It is, however, worthwhile to note that in 1D a reduction from a bidomain model to a monodomain model is possible also when this condition is not fulfilled. In order to see this, we integrate (16) and get the relation ux ¼
Mi vx ; Mi þ Me
where we have used the assumption that u and v have zero derivatives at the boundaries. Inserting this expression in (15), yields the following scalar equation for v M iM e vt ¼ vx : ð17Þ Mi þ Me x We conclude that the bidomain model can always be reduced to a monodomain type of model when we are addressing a 1D problem. Unfortunately, this seems to be a pure 1D effect. 3. A parameter estimation problem Recall the splitting procedure briefly described in connection with Eqs. (4), (5) and (6), (7). At each time step, both of these subsystems must be solved. We will assume that an efficient solver for the system (4), (5) is available and focus on how to approximate the bidomain equations (6) and (7) with the monodomain model (9) and (10). In view of the splitting approach, it is natural to consider this problem in a time interval [tm1, tm], i.e. in a typical time interval used in the splitting. As mentioned in Section 1, if (8) does not hold, then the solution (v, u) of (6) and (7) will in general not satisfy Eqs. (9) and (10). In order to derive an algorithm for computing appropriate values for k, let us introduce the following notation, cf. (9) and (10), for the monodomain model:3 k r ðM i rwÞ in H ; 1þk 1 r ðM i rwÞ in H ; r ðM i rqÞ ¼ ð1 þ kÞ ðM i rwÞ n ¼ 0 along oH ; wt ¼
ðM i rqÞ n ¼ 0
along oH ;
ð18Þ ð19Þ ð20Þ ð21Þ
where we recall that the boundary conditions (20) and (21) are motivated by the discussion leading to Eqs. (13) and (14). That is, w and q denote the transmembrane and extracellular potentials generated by the monodomain model. Obviously, the solution (w, q) of (18), (21) depends on k, i.e. w ¼ wðx; t; kÞ; q ¼ qðx; t; kÞ: With this notation at hand, our challenge can be formulated as follows: Find k such that the solution (w, q) of (18) and (19) defines an optimal, or close to optimal, approximation of the solution (v, u) of (6) and (7) in the time interval [tm-1, tm]. Expressed in terms of mathematical symbols, we would thus ideally like to solve a problem on the form 2
min kðv; uÞ ðwðkÞ; qðkÞÞk ; k
ð22Þ
where kÆk denotes an appropriate norm.
3
Strictly speaking, an initial condition for the transmembrane potential w is needed, i.e. w(x, tm1) for x 2 H is required. We will return to this issue below.
280
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
In general, applying such an output least squares approach would of course require the solution (v, u) of the bidomain equations in the time interval [tm1, tm]. However, if we are content with considering this problem in terms of a tail-ordered ‘‘energy’’ norm and a semi-discrete version of the PDEs, then it turns out that no explicit knowledge of (v, u) is required in order to solve (22). More precisely, we will discretize the involved equations in time, and show that a suitable output least squares problem for determining k can be easily solved at each time step. Recall that each time interval [tm1, tm] is treated separately, and hence this approach will lead to a time dependent coefficient k = k(t). Consequently, from a global perspective, our algorithm will not necessarily compute an optimal representation for k. Only time-local optimal values for this parameter are obtained by our method. 3.1. Theoretical considerations This section is devoted to a theoretical study of the problem discussed above. Our goal is to investigate whether it is possible to compute optimal, or at least close to optimal, monodomain approximations of (6) and (7) in the time interval [tm1, tm] without actually solving these equations. 3.1.1. Semi-discrete models A semi-discrete approximation of (18) and (19) may be written on the form km r ðM i rwm Þ ¼ wm1 in H ; 1 þ km 1 r ðM i rwm Þ in H ; r ðM i rqm Þ ¼ ð1 þ kÞ
wm Dt
ð23Þ ð24Þ
where km, a constant, denotes k restricted to the time interval [tm1, tm], and Dt ¼ tm tm1 ; wm ðxÞ wðx; tm Þ
for x 2 H ;
m1
w ðxÞ wðx; tm1 Þ for x 2 H ; qm ðxÞ qðx; tm Þ for x 2 H : If km is given, and only the transmembrane potential w is of interest, then it is sufficient to solve (23). However, in the process of determining an appropriate value for km, it turns out that it is convenient to consider both Eq. (23) and (24). We will return to this issue in Section 3.1.4. Let us now focus on the process of advancing the approximate potentials from time step tm1 to tm. Assume that wm1 and qm1 are given,4 and that we want to determine an optimal value for km, to use in (23), (24), for the computation of wm and qm. To this end, we consider the bidomain equations in the time interval t 2 [tm1, tm], i.e. vt ¼ r ðM i rvÞ þ r ðM i ruÞ
r ðM i rvÞ þ r ððM i þ M e ÞruÞ ¼ 0 m1
vðx; tm1 Þ ¼ w m1
ðxÞ
m1
ð25Þ
in H ;
ð26Þ
in H ;
for x 2 H : m1
ð27Þ m1
Since only (w ,q ), and not (v ,u ), are available to us as data, we apply an initial condition expressed in terms of wm1, see (27). The PDEs in the bidomain equations (25) and (26) can be discretized in a similar fashion to the monodomain model, leading to the system vm Dtr ðM i rvm Þ Dtr ðM i rum Þ ¼ wm1 r ðM i rvm Þ þ r ððM i þ M e Þrum Þ ¼ 0
4
That is, wm1 and qm1 have already been computed.
in H ;
in H ;
ð28Þ ð29Þ
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
281
where vm ðxÞ vðx; tm Þ for x 2 H ; um ðxÞ uðx; tm Þ for x 2 H : Note that the initial condition (27) has been incorporated into the right hand side of Eq. (28). In this case, a coupled linear system of two PDEs appears as a results of the discretization procedure. More precisely, Eqs. (28) and (29) defines an elliptic system consisting of two PDEs with two unknowns vm and um. Our goal is to compare the solutions of the systems (28), (29) and (23), (24) and thereby obtain a method for computing an appropriate value for km. More precisely, we will derive a technique for doing so that does not require the solution of the bidomain equations (28) and (29). Prior to presenting our analysis of this problem, we must however introduce some notation. 3.1.2. Function spaces and variational forms Recall that we are considering systems of elliptic equations, cf. (23), (24) and (28), (29). It turns out that the suitable Sobolev spaces in the current situation are o/ o/ V ¼ H 1 ðH Þ ¼ / 2 L2 ðH Þ; ; 2 L2 ðH Þ ; ox oy ð30Þ Z U ¼ w 2 V ; w dx ¼ 0 : H
The variational form of (23) and (24) can now be written on the form: Find (wm, qm) 2 V · U such that Z Z Z km /wm dx þ Dt r/M i rwm dx ¼ /wm1 dx; ð31Þ 1 þ k m H H H Z Z 1 rwM i rqm dx ¼ rwM i rwm dx ð32Þ ð1 þ km Þ H H for all (/, w) 2 V · U. Similarly, the weak form of (28) and (29) reads: Find (vm, um) 2 V · U such that Z Z Z Z m m m /v dx þ Dt r/M i rv dx þ Dt r/ M i ru dx ¼ /wm1 dx; H H H H Z Z rwM i rvm dx þ rwðM i þ M e Þrum dx ¼ 0; H
ð33Þ ð34Þ
H
for all (/, w) 2 V · U. Throughout this paper we will assume that the domain H, the conductivities Mi and Me and the function wm1 are ‘‘well-behaved’’ such that there exist unique solutions (wm, qm) and (vm, um) of (31), (32) and (33), (34), respectively. Further details about this issue can, e.g., be found in [4]. 3.1.3. An inner-product Let us define the following inner-product [Æ, Æ]E on V · U Z Z Z Z ½ð/; wÞ; ðs; nÞE ¼ /s dx þ Dt r/M i rs dx þ Dt r/M i rn dx þ Dt rwM i rs dx H H H H Z þ Dt rwðM i þ M e Þrn dx
ð35Þ
H
for all (/, w), (s, n) 2 V · U, with the associated ‘‘energy’’ norm 2
kð/; wÞkE ¼ ½ð/; wÞ; ð/; wÞE : Note that, for small values of Dt, Z kð/; wÞk2E /2 dx ¼ k/kL2 ðH Þ ; H
ð36Þ
282
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
provided that /, w, Mi and Me are ‘well-behaved’. Consequently, there is a close relation between the kÆkE-norm of (vm, um) and the L2-norm of the transmembrane potential vm. Furthermore, the weak form (33) and (34) of the semi-discrete bidomain equations can now compactly be expressed as: Find (vm, um) 2 V · U such that Z m m ½ð/; wÞ; ðv ; u ÞE ¼ /wm1 dx ð37Þ H
for all (/, w) 2 V · U. 3.1.4. A minimum error analysis Let us now consider the error, measured in the kÆkE-norm, introduced by solving the monodomain model (23) and (24) instead of the more physically accurate bidomain equations (28) and (29) in the time interval [tm1, tm]. To this end, consider the cost-functional 2
Q ¼ Qðkm Þ ¼ kðvm ; um Þ ðwm ; qm ÞkE ¼ ½ðvm ; um Þ ðwm ; qm Þ; ðvm ; um Þ ðwm ; qm ÞE 2
2
¼ kðvm ; um ÞkE 2½ðvm ; um Þ; ðwm ; qm ÞE þ kðwm ; qm ÞkE ; where (wm, qm) and (vm, um) denote the solutions of (23), (24) and (28), (29), respectively. Throughout this section it is also important to keep in mind that both wm and qm depend on km, i.e. wm ¼ wm ðkm Þ; qm ¼ qm ðkm Þ: From Eq. (37), choosing (/, w) = (wm, qm) 2 V · U, we find that Z 2 2 Qðkm Þ ¼ kðvm ; um ÞkE 2 wm wm1 dx þ kðwm ; qm ÞkE : H
m
m
Clearly, (v , u ) is independent of km. Therefore, instead of minimizing Q with respect to km we can consider a cost-functional J on the form Z m m 2 J ðkm Þ ¼ kðw ; q ÞkE 2 wm wm1 dx: ð38Þ H
Note that, in order to minimize J, it is not necessary to solve the semi-discrete bidomain equations (28) and (29)! Thus, we have reached our goal; it is actually possible to compute time-local optimal values km for k without solving the more physically accurate bidomain model. By defining km ¼ max arg min kðvm ; um Þ ðwm ; qm Þk2E ; ð39Þ km 2½;K
our main results may be formulated as follows: Theorem 3.1. Let (wm, qm) = (wm(km), qm(km)) and (vm, um) denote the solutions of Eqs. (31), (32) and (33), (34), respectively. Then km ¼ max arg min J ðkm Þ km 2½;K Z m m 2 m m1 ¼ max arg min kðw ; q ÞkE 2 w w dx ; km 2½;K
H
where kÆkE is the norm defined in (35), (36), and ½; K Rþ represents5 a suitable bound on km. Consequently, it is not necessary to solve the semi-discrete bidomain model (33), (34) in order to compute an optimal value km for km. 5
Appropriate values for e and K will typical depend on the problem under consideration.
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
283
The kÆkE-norm of the solution (wm(km), qm(km)) of (31) and (32) depends continuously on km. Therefore, J(km) will attain its minimum on the closed interval [, K]. The minimum may not be unique, and we have, 2 for the sake of convenience, defined km to be the largest value in [, K] at which kðvm ; um Þ ðwm ; qm ÞkE attains its minimum. By minimizing J with respect to km at each time interval [tm1, tm], we obviously obtain a time dependent parameter k in (18) and (19), i.e. kðtÞ ¼ km
for t 2 ½tm1 ; tm and m ¼ 1; . . . ; M;
where M denotes the number of time intervals used in the discrete model. This approach will clearly yield a piecewise constant function k(t) which depends on the temporal discretization parameter Dt, i.e. k ¼ kðt; DtÞ: From a theoretical point of view, it would certainly be interesting to try to analyze the limit of k(t; Dt) as Dt ! 0. That topic is, however, beyond the scope of the present text. Instead, we will use our analysis to derive a suitable algorithm for solving a modified monodomain type of model. 4. Practical considerations We now turn our attention towards the practical aspects of the problem discussed above. More specifically, our goal is to use Theorem 3.1 to propose a modified monodomain type of scheme suitable for simulation the electrical activity in the heart. Eqs. (30) and (32) imply that Z 1 1 m m q ¼ w w dx ; ð1 þ km Þ jH j H m
ð40Þ
R where jH j ¼ H 1 dx (i.e. the volume (area in 2D) of the domain occupied by the heart H). Thus, the cost-functional J, defined in (38), takes the form Z 2 J ðkm Þ ¼ kðwm ; qm ÞkE 2 wm wm1 dx H 2 Z Z m 1 1 m m ¼ w ; ð1 þ k Þ w w dx 2 wm wm1 dx: m jH j H H E Here 2 Z Z Z m 2Dt m m 2 w ; ð1 þ km Þ1 wm 1 w dx ¼ ðw Þ dx þ Dt rwm M i rwm dx jH j H 1 þ km H H E Z Z Dt rwm M i rwm dx þ rwm ðM i þ M e Þrwm dx 2 ð1 þ km Þ H H 2 Z Z 1 2 ¼ ðwm Þ dx þ Dt 1 rwm M i rwm dx 1 þ km H H Z Dt þ rwm M e rwm dx; 2 ð1 þ km Þ H cf. (35), and consequently Z J ðkm Þ ¼ ðwm Þ2 dx þ Dt H
þ Dt
2
ðkm Þ
Z
1 2
ð1 þ km Þ
ð1 þ km Þ2
H
Z
rwm M i rwm dx
H
rwm M e rwm dx 2
Z
wm wm1 dx; H
ð41Þ
284
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
where it is important to remember that wm depends on km, i.e. wm ¼ wm ðkm Þ: Next, define a¼
km 2 ½aa ; ab 1 þ km
for km 2 ½; K;
where aa ¼
; 1þ
ab ¼
K 1þK
and note that the cost-functional J may conveniently be written on the form Z Z Z Z 2 2 J ðaÞ ¼ ðwm Þ dx þ Dta2 rwm M i rwm dx þ Dtð1 aÞ rwm M e rwm dx 2 wm wm1 dx H
H
H
ð42Þ
H
for a 2 (aa, ab). 4.1. A formula for wm(km) In order to minimize J, more information about how wm depends on km is needed. More precisely, we will now derive a formula for wm(km) = wm(a) and use it to construct an algorithm suitable for estimating the optimal value km for km. In fact, it turns out that km can be approximately determined by minimizing a fourth order polynomial. The details are as follows: Let w denote the continuous solution of the monodomain equation in the time interval [tm1, tm]: wt ¼ ar ðM i rwÞ for x 2 H ; t 2 ðtm1 ; tm ; aðM i rwÞ n ¼ 0 for x 2 oH ; t 2 ðtm1 ; tm ;
ð43Þ ð44Þ
wðx; tm1 Þ ¼ wm1 ðxÞ for x 2 H ;
ð45Þ
with the discrete counterpart (23). Our goal is to characterize how the solution w of this problem depends on km the parameter a ¼ 1þk . It turns out that this can be accomplished by considering a simple auxiliary problem: m rt ¼ r ðM i rrÞ for x 2 H ; t 2 ðtm1 ; tm ; ðM i rrÞ n ¼ 0 for x 2 oH ; t 2 ðtm1 ; tm ;
ð46Þ ð47Þ
rðx; tm1 Þ ¼ wm1 ðxÞ for x 2 H :
ð48Þ
We will now show that the solution w of (43)–(45) is given by the formula wðx; tÞ ¼ rðx; a½t tm1 þ tm1 Þ:
ð49Þ
Clearly, defining l(t) = a[t tm1] + tm1, ½rðx; lðtÞÞt ¼ art ðx; lðtÞÞ ¼ ar ½M i rrðx; lðtÞÞ and we conclude that the function w defined in (49) satisfies the differential equation (43). Moreover, it is easy to verify that this function also fulfills the boundary and initial conditions expressed in Eqs. (44) and (45), respectively. The key point of this analysis is that the solution r of (46)–(48) does not depend on the parameter a. Consequently, the solution w = w(a) of (43)–(45) can be determined for all values of a by simply solving this auxiliary problem. This feature can of course be used to compute an optimal value a* for a (and consequently an optimal value km for km). 4.1.1. Discretization Recall the semi-discrete monodomain equation (23). In a similar fashion we get the semi-discrete form of (46)–(48):
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
rm Dtr ðM i rrm Þ ¼ wm1 m
ðM i rr Þ n ¼ 0
for x 2 H ;
for x 2 oH :
285
ð50Þ ð51Þ
From formula (49) we now obtain the approximation wm ðxÞ wðx; tm Þ ¼ rðx; aDt þ tm1 Þ ð1 aÞrm1 þ arm ¼ ð1 aÞwm1 þ arm ; where we have applied linear interpolation to the function r.
ð52Þ
4.2. A fourth order polynomial Substituting the approximation (52) for wm into formula (42) yields a fourth order polynomial J ðaÞ ¼ c0 þ c1 a þ c2 a2 þ c3 a3 þ c4 a4
for a 2 ½aa ; ab :
ð53Þ
The coefficients c0, c1, c2, c3, c4 can obviously be expressed in terms of integrals involving wm1, rm and their gradients, see Appendix A. km Since a ¼ 1þk , it follows that it is sufficient to minimize a fourth order polynomial in order to determine an m approximation6 of the optimal value km of km to use in the monodomain model in the time interval [tm1, tm]. Recall that our analysis of the monodomain model was presented in terms of an operator splitting procedure for simulating the electrical activity in the heart, cf. the discussion of Eqs. (4)–(7) and Section 3. If a Godunov splitting is applied, our findings lead to the following scheme: Algorithm 4.1. For given positive numbers < K and initial conditions v0 and s0 for v and s, we define the following algorithm: 1. 2. 3. 4. 5.
K aa ¼ 1þ , ab ¼ 1þK Dt = T/M tm = m Æ Dt for m = 0, 1, . . . , M w0 = v0 and s0 = s0 for m = 1 . . . M do (a) Solve
os ¼ F ðs; vÞ; sðtm1 Þ ¼ sm1 ; ot ov ¼ Iðv; sÞ; vðtm1 Þ ¼ wm1 ot (b) (c) (d) (e)
ð54Þ ð55Þ
for x 2 H and t 2 (tm1, tm] by, e.g., a Runge–Kutta method. wm1 = v(tm) Solve (50), (51) for rm Compute the coefficients c0, c1, c2, c3, c4, see Appendix A Compute a ¼ max arg min J ðaÞ ¼ c0 þ c1 a þ c2 a2 þ c3 a3 þ c4 a4 : a2½aa ;ab
m
m1
(f) w = (1 a*)w + a*rm 1 1 R m (g) q ¼ ð1 þ km Þ ðwm jH j wm dxÞ. (This step can be omitted if only the transmembrane potential is H of interest.) a (h) km ¼ 1a ð2 ½; KÞ Step 5(f) can be omitted. It has only been included in the algorithm to make it possible to monitor the time dependency of the optimal choice of k in the monodomain model (18) and (19). This scheme can be modified in a straightforward manner in order to apply our results to a Strang splitting of the system (1)–(3). 6
Keep in mind that linear interpolation was used in (52).
286
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
Note that the computational needs of this algorithm are only slightly larger than that of the standard monodomain approach. More specifically, the additional CPU time required is that of steps 5(d) and (e), and is thus only linked to the construction and minimization of a fourth order polynomial. Such tasks are ‘‘cheap’’ compared to the computational demands of solvers of ODEs and PDEs (steps 5(a) and (c)). 5. Numerical experiments This section is devoted to an experimental study of the scheme proposed above. We will compare the performance of Algorithm 4.1 with that of the traditional monodomain method. The numerical experiments were carried out on a two dimensional domain representing an axial slice through the heart, midway between the apex and the basal plane. The fiber directions were defined to run circumferentially to the left ventricle, see Fig. 1 and [13] for details. This approach leads to the following formulas for the intracellular and extracellular conductivities Mi and Me: M i ¼ rti I þ ðrli rti Þal aTl ; Me ¼
rte I
þ
ðrle
rte Þal aTl ;
ð56Þ ð57Þ
where rti ; rli ; rte and rle are scalars and al = al(x) is a vector function incorporating the fibrous structure of the heart into the PDE models.7 In the traditional monodomain method a suitable value for k must somehow be determined, cf. the discussion of Eqs. (8)–(10). Based on expressions (56) and (57), a fair choice seems to be 2 ~ k ¼ arg min kðrte ; rle Þ kðrti ; rli Þkl2 ¼ 0:80; k
or alternatively ~ k ¼ arg min k
Z
kM e ðxÞ kM i ðxÞk2F dx; ¼ 0:71;
ð58Þ
ð59Þ
H
where kÆkF denotes the Frobenius norm. In all the examples presented below, ~k is equal to 0.80; the results obtained with ~ k ¼ 0:71 were very similar. In order to compare the performance of Algorithm 4.1 with that of the standard monodomain scheme, we also solved the bidomain equations and computed the relative L2-errors of the two methods: kvm wm kL2 for m ¼ 0; 1; . . . ; M; kvm kL2 kvm zm kL2 ¼ for m ¼ 0; 1; . . . ; M; kvm kL2
Eopt m ¼
ð60Þ
Emono m
ð61Þ
where vm, wm and zm denote the transmembrane potentials generated by the bidomain equations, Algorithm 4.1 and the traditional monodomain scheme (with ~k ¼ 0:80), respectively. Throughout this section it is important to keep in mind that the value km for k, see (39) and Theorem 3.1, is optimal with respect to the ‘‘energy’’ norm kÆkE, and not with respect to the classical L2-norm used in (60) and (61). In addition, due to the use of linear interpolation in (52), step 5(h)) in Algorithm 4.1 only provides an approximation of km . All the tests were carried out with a time step Dt = 0.125 ms and 4059 mesh-points in the 2D domain occupied by the ‘‘heart’’ H. 7
In 3D these formulas take the form M i ¼ rti I þ ðrli rti Þal aTl þ ðrni rti Þan aTn ; M e ¼ rte I þ ðrle rte Þal aTl þ ðrne rte Þan aTn ;
where al and an are vector functions used to model the fibrous structure of the heart.
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
287
Fig. 1. An illustration of the fiber directions in the heart.
5.1. Case I Recall that the function I = I(v, s) implements the coupling between the ODEs (1), modeling the cell dynamics, and the bidomain equations (2), (3). For the sake of academic interest, let us now consider the problem obtained by defining I = 0. In this case, the bidomain equations are decoupled from the system (1), and we can thus not expect (1)–(3) to define an accurate model for the electrical activity in the heart. Nevertheless, the resulting equations (2) and (3) for u and v are linear and well-posed, and consequently interesting to consider from both a mathematical and an algorithmic point of view. If I = 0, then the operator splitting approach discussed above is not needed, and steps 5 (a) and (b) in Algomono rithm 4.1 can be omitted. Fig. 2 shows the behavior of the relative L2-errors Eopt for this problem, cf. m and Em the discussion of definitions (60) and (61). Clearly, Algorithm 4.1 delivers far more accurate transmembrane potentials than the traditional monodomain scheme. Furthermore, the optimal value km for k is significantly larger than the value ~ k ¼ 0:80 computed by the ‘‘brute force’’ technique (58). 5.2. Case II In simple heart models, I is typically defined as a cubic polynomial in the transmembrane potential v: IðvÞ ¼ A2 ðv vrest Þðv vthreshold Þðv vpeak Þ;
4.5
x 10
1.7
MonoDomain λ~=0.8 Optimal
4
1.65
3.5
1.6
3
1.55
* λm
Relative error (L2 −norm)
* for I=0 λm
I=0
3
2.5
1.5
2
1.45
1.5
1.4
1
1.35
0.5
0
2
4
6
8
10
time (ms)
12
14
16
18
20
1.3
0
2
4
6
8
10
12
14
16
18
20
time (ms)
mono Fig. 2. The results obtained in Case I. In the left figure we have plotted the relative errors Eopt (solid line) m (dashed line) and Em associated with Algorithm 4.1 and the traditional monodomain scheme, see the discussion of Eqs. (60) and (61). The right figure shows how the value of km , generated in step 5(h)) in Algorithm 4.1, changes with time.
288
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290 * for Cubic λm
Cubic 0.4
~ MonoDomain λ =0.8 Optimal
0.35
1.9
0.3 1.8
0.25
m
1.7 λ*
Relative error (L2−norm)
2
0.2
1.6
0.15 1.5
0.1
1.4
0.05
0
0
2
4
6
8
10
time (ms)
12
14
16
18
20
1.3
0
2
4
6
8
10
12
14
16
18
20
time (ms)
mono Fig. 3. The results obtained in Case II. In the left figure we have plotted the relative errors Eopt (solid line) m (dashed line) and Em associated with Algorithm 4.1 and the traditional monodomain scheme. The right figure shows how the value of km , generated in step 5(h)) in Algorithm 4.1, changes with time.
where we use the parameters A ¼ 0:04; vrest ¼ 85 mV; vthreshold ¼ 65 mV; vpeak ¼ 40 mV; see e.g. [10] for further details. This means that step 5(a) in Algorithm 4.1 takes the form: Solve ov ¼ IðvÞ; vðtm1 Þ ¼ wm1 ; ot for x 2 H and t 2 (tm1, tm]. There is thus no model for the cell dynamics involved, and the operator splitting is only applied in order to handle the nonlinearities introduced into equations (2), (3) by the function I. Fig. 3 shows the relative errors Emono and Eopt m m , along with the value of km , at each time step for this problem. The qualitative behavior of these results are similar to the observations made in Case I; the error associated with Algorithm 4.1 is significantly smaller than that of the standard monodomain scheme, and the value ~ k ¼ 0:80 defined in (58) is far too small. 5.3. Case III In the final test we use the FitzHugh–Nagumo cell model, see [5] and [18]. The results obtained in this case are shown in Fig. 4. During the first 4ms of the simulation, the traditional monodomain method produces results that are slightly more accurate than those obtained with our method.8 For t > 7 ms, the performance of Algorithm 4.1 is considerably better than that of the standard monodomain scheme. 5.4. A remark In the three cases presented above we have seen that the optimal value km for k varies between 1.33 and 2.05. Based on this observation we also performed a series of experiments with k = 1.75 in the standard monodomain scheme. These simulations generated results that were almost as good as those obtained with Algorithm 4.1, indicating that the classical approach might work very well, provided that a suitable value for k is Please keep in mind that the value km for k is optimal with respect to the energy kÆkE-norm, see (39) and Theorem 3.1. Fig. 4 shows the relative L2 errors. 8
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290 λ*m for FitzHugh
FitzHugh 0.35
2.1
~
MonoDomain λ =0.8 Optimal
2
0.3
1.9
0.25 1.8
2
Relative error (L −norm)
289
λ*
m
0.2 1.7
0.15 1.6
0.1 1.5
0.05
0
1.4
0
2
4
6
8
10
12
14
16
18
20
1.3
0
2
4
6
time (ms)
8
10
12
14
16
18
20
time (ms)
mono Fig. 4. The results obtained in Case III. In the left figure we have plotted the relative errors Eopt (solid line) m (dashed line) and Em associated with Algorithm 4.1 and the traditional monodomain scheme. The right figure shows how the value of km , generated in step 5(h) in Algorithm 4.1, changes with time.
used. However, in order to determine such a number for k, we must run Algorithm 4.1; the ‘‘brute force’’ techniques (58) and (59) seem to significantly underestimate the appropriate values. 6. Summary In this paper we have proposed a modified monodomain scheme intended for fast simulations of the electrical activity in the heart. More precisely, a method for computing optimal monodomain approximations of the more physically correct bidomain equations has been introduced and analyzed. As expressed in Theorem 3.1, the key feature of our approach is that the optimal results can be obtained without actually solving the bidomain model, provided that optimality is defined in terms of a suitable ‘‘energy’’ norm. Our theoretical findings were illuminated by a series of numerical experiments. In each of these examples, our method produced significantly more accurate results than the traditional monodomain approach. Furthermore, the new algorithm is only slightly more CPU demanding than the standard scheme. Appendix A The formulas for the coefficients of the polynomial J(a) in Eq. (53) are Z Z m1 2 c0 ¼ w dx þ Dt rwm1 M e rwm1 dx; H Z H Z m1 m1 c1 ¼ 4Dt rw M e rw dx þ 2Dt rwm1 M e rrm dx; H Z Z Z H Z m1 2 m m1 m 2 c2 ¼ w dx 2 r w dx þ ðr Þ dx þ Dt rwm1 M i rwm1 dx H H HZ H Z Z m1 m1 m1 m þ 6Dt rw M e rw dx 6Dt rw M e rr dx þ Dt rrm M e rrm dx; ZH Z H ZH m1 m1 m1 m c3 ¼ 2Dt rw M i rw dx þ 2Dt rw M i rr dx þ 6Dt rwm1 M e rrm dx H ZH ZH m1 m1 m m 4Dt rw M e rw dx 2Dt rr M e rr dx; H H Z Z Z m1 m1 c4 ¼ Dt rw M i rw dx 2Dt rwm1 M i rrm dx þ Dt rrm M i rrm dx HZ HZ H Z m1 m1 m1 m þ Dt rw M e rw dx 2Dt rw M e rr dx þ Dt rrm M e rrm dx: H
H
H
290
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
References [1] G.W. Beeler, H. Reuter, Reconstruction of the action potential of ventricular myocardial fibres, J. Physiol. 268 (1977) 177–210. [2] Cellml. Available from:. [3] D. DiFrancesco, D. Noble, A model of cardiac electrical activity incoportating ionic pumps and concentration changes, Philos. Trans. Roy. Soc. Lond. Biol. 307 (1985) 353–398. [4] L.C. Evans, Partial Differential Equations, American Mathematics Society, 1998. [5] R.A. FitzHugh, Impulses and physiological states in theoretical models of nerve membrane, Biophys. J. 1 (1961) 445–466. [6] P. Colli Franzone, L. Guerri, S. Tentoni, Mathematical modeling of the excitation process in myocardial tissue: influence of fiber rotations on wavefront propagation and potential field, Math. Biosci. 101 (1990) 155–235. [7] I.J. Le Grice, P.J. Hunter, B.H. Smaill, Laminar structure of the heart: a mathematical model, Am. J. Physiol. Heart. Circ. Physiol. 272 (1997) H2466–H2476. [8] C.S. Henriquez, Simulating the electrical behavior of cardiac tissue using the bidomain model, Crit. Rev. Biomed. Eng. 21 (1993) 1–77. [9] Q. Huang, J.C. Eason, F.J. Claydon, Membrane polarization induced in the myocardium by defibrillation fields: an idealized 3-d finite element bidomain/monodomain torso model, IEEE Trans. Biomed. Eng. 46 (1) (1999) 26–34. [10] P.J. Hunter, P.A. McNaughton, D. Noble, Analytical models of propagation in the exitable cells, Prog. Biophys. Mol. Biol. 30 (1975) 99–144. [11] J. Keener, J. Sneyd, Mathematical Physiology, Springer-Verlag, 1998. [12] W. Krassowska, J.C. Neu, Effective boundary conditions for syncytial tissue, IEEE Trans. Biomed. Eng. 41 (1994) 143–150. [13] G.T. Lines, P. Grøttum, A. Tveito, Modeling the electrical activity of the heart, a bidomain model of the ventricles embedded in a torso, Comput. Visual. Sci. 5 (4) (2002) 195–213. [14] C.H. Luo, Y. Rudy, A model of the ventricular cardiac action potential: depolarisation, repolarisation, and their interaction, Circul. Res. 68 (1991) 1501–1526. [15] C.H. Luo, Y. Rudy, A dynamic model of the cardiac ventricular action potential, Circul. Res. 74 (1994) 1071–1096. [16] M. Murillo, X.C. Cai, A fully implicit parallel algorithm for simulating the non-linear electrical activity of the heart, Numer. Linear Algebra Appl. 11 (2004) 261–277. [17] A.L. Muzikant, C.S. Henriquez, Bipolar stimulation of a three-dimensional bidomain incorporating rotational anisotropy, IEEE Trans. Biomed. Eng. 45 (1998) 449–462. [18] J. Nagumo, S. Arimoto, S. Yoshizawa, An active pulse transmission line simulating 1214-nerve axons, Proc. IRL 50 (1960) 2061– 2070. [19] Z. Qu, A. Garfinkel, An advanced algorithm for solving partial differential equation in cardiac conduction, IEEE Trans. Biomed. Eng. 46 (9) (1999) 1166–1168. [20] K. Skouibine, W. Krassowska, Increasing the computational efficiency of a bidomain model of defibrillation using a time-dependent activating function, Biomed. Eng. 28 (2000) 772–780. [21] K. Skouibine, N. Trayanova, P. Moore, A numerically efficient model for simulation of defibrillation in an active bidomain sheet of myocardium, Math. Biosci. 166 (2000) 85–100. [22] J. Sundnes, G.T. Lines, X. Cai, B.F. Nielsen, K.A. Mardal, A. Tveito, Computing the Electrial Activity in the Heart, Springer-Verlag, 2006. [23] J. Sundnes, G.T. Lines, A. Tveito, An operator splitting method for solving the bidomain equations coupled to a volume conductor model for the torso, Math. Biosci. 194 (2) (2005) 233–248. [24] J. Sundnes, B.F. Nielsen, K.A. Mardal, X. Cai, A. Tveito, On the computational complexity of the bidomain and the monodomain models of electrophysiology, Ann. Biomed. Eng., in press. [25] L. Tung, A Bi-domain model for describing ischemic myocardial D-C potentials, Ph.D. thesis, MIT, Cambridge, MA, 1978. [26] R.L. Winslow, A.L. Kimball, A. Varghese, D. Noble, Simulating cardiac sinus and atrial network dynamics on the connection machine, Physica D 64 (1993) 281–298. [27] R.L. Winslow, J. Rice, S. Jafri, E. Marban, B. O’Rourke, Mechanisms of altered excitation-contraction coupling in canine tachycardia-induced heart failure, II, model studies, Circul. Res. 84 (1999) 571–586. [28] R.L. Winslow, D.F. Scollan, J.L. Greenstein, C.K. Yung, W. Baumgartner, G. Bhanot, D.L. Gresh, B.E. Rogowitz, Mapping, modelling, and visual exploration of structure-function relationships in the heart, IBM Syst. J. 40 (2001) 342–359.
Optimal monodomain approximations of the bidomain equations Bjørn Fredrik Nielsen *, Tomas Syrstad Ruud, Glenn Terje Lines, Aslak Tveito Simula Research Laboratory, P.O. Box 134, 1325 Lysaker, Norway
Abstract The bidomain equations are widely accepted to model the spatial distribution of the electrical potential in the heart. Although order optimal methods have been devised for discrete versions of these equations, it is still very CPU demanding to solve the equations numerically on a sufficiently fine mesh in 3D. Furthermore, the equations are hard to analyze; from a mathematical point of view, very little is known about the qualitative behavior of the solutions generated by these equations. It is well known that upon appealing to a certain relation between the extracellular and intracellular conductivities, the bidomain model can be rewritten in terms of a scalar reaction diffusion equation referred to as the monodomain model. This model is of course much easier to solve, and also the qualitative properties of the solutions are well known; such equations have been studied intensively for decades. It is the purpose of the present paper to show how the bidomain equations can be approximated in an optimal manner by the solution of a monodomain model. The key feature here is that this optimal solution can be computed without solving the bidomain model itself. The solution is obtained by putting the problem into a framework of parameter identification problems. Our results are illuminated by a series of numerical experiments. 2006 Elsevier Inc. All rights reserved. Keywords: Electrical behavior of cardiac tissue; Bidomain and monodomain models; Output least squares
1. Introduction The electrical activity in the heart is governed by a complex, fully coupled system of ordinary and partial differential equations: os ¼ F ðs; vÞ in H ; ot vt þ Iðv; sÞ ¼ r ðM i rvÞ þ r ðM i ruÞ in H ;
ð2Þ
r ðM i rvÞ þ r ððM i þ M e ÞruÞ ¼ 0
ð3Þ
in H :
*
ð1Þ
Corresponding author. E-mail addresses: [email protected] (B.F. Nielsen), [email protected] (T.S. Ruud), [email protected] (G.T. Lines), aslak@ simula.no (A. Tveito). 0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.05.158
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
277
In this model, v and u represent the transmembrane and extracellular potentials, respectively, and H the physical domain occupied by the heart. The intracellular and extracellular conductivity tensors1 Mi and Me are typically defined in terms of symmetric and positive definite matrices depending on the spatial position x. The state s of the ionic currents is incorporated into the model by the system (1) of ordinary differential equations (ODEs) and the possibly nonlinear function I, see [25] or [22] for further details concerning this model. In models of this kind, the partial differential equations (PDEs) (2) and (3) are referred to as the bidomain equations, and (1) denotes a (generic) system of ordinary differential equations. More precisely, s is a state vector for the ions, and (1) incorporates the effect of the flow of ions across the cell membrane. Throughout the last 30 years, a wide variety of such cell models have been proposed and analyzed, see Beeler and Reuter [1], DiFrancesco [3], Luo and Rudy [14,15] and Winslow et al. [27] – see [2] for a very nice overview of these models. For further information about the bidomain equations, we would like to refer to [11] and [22]. Eqs. (1)–(3) define a nonlinear, fully coupled system of differential equations for v, u and s. Analytical solutions are not available for such complex models. They must hence be solved by some kind of numerical procedure. Ideally we would of course like to solve the system (1)–(3) in a fully implicit manner, see [16]. However, due to the complexity of the involved ODEs, see e.g. [26], it is common to use some sort of operator splitting scheme [19,23]. Schemes of this kind are typically defined in terms of successive solutions of the two subsystems os ¼ F ðs; vÞ in H ; ot ov ¼ Iðv; sÞ in H ot
ð4Þ ð5Þ
and vt ¼ r ðM i rvÞ þ r ðM i ruÞ in H ;
ð6Þ
r ðM i rvÞ þ r ððM i þ M e ÞruÞ ¼ 0 in H :
ð7Þ
Further information about such methods can be found in [23,22]. A complete simulator for the electrical activity in the heart consists of solving (4), (5) and (6), (7) at each time step. Computer experiments indicate that a very fine computational mesh is needed in order to solve the subsystem (6), (7) with any reasonable accuracy. Typically, 20–40 millions grid-points are needed in 3D, see [13,28]. This leads to an extremely CPU demanding task, and in many cases a need for considering simplified models. In this paper we will consider one of this simplified approaches – namely, the monodomain model. If M e ¼ kM i
ð8Þ
for some scalar k, then (6) and (7) can be written on the form k r ðM i rvÞ in H ; 1þk ð1 þ kÞr ðM i ruÞ ¼ r ðM i rvÞ in H ; vt ¼
ð9Þ ð10Þ
which commonly is referred to as the monodomain model.2 Note that (9) and (10) can be solved in a sequentially manner; first solve (9) for v and thereafter compute u from (10). In fact, if only the transmembrane potential is of interest, then it is sufficient to solve (9). This is of course a simpler task than determining v and u, at each time step, from (6) and (7). For a further discussion of this issue, we refer to [24]. Although the condition (8) is incorrect from a biological point of view, the simplified model (9) is nevertheless frequently used to simulate the electrical activity in the heart. In order to perform such simulations, a suitable value for k must be determined. The purpose of this paper is to shed some light onto this problem, and to design an algorithm for computing optimal, or at least close to optimal, values for k. 1
More precisely, Mi and Me are scaled conductivities, see [22] for details. Usually (9), without (10), is referred to as the monodomain model. For the sake of completeness, we have also included Eq. (10). In the present context, this turns out to be very useful. 2
278
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
One possible approach towards this problem is to minimize the difference between kMi and Me, measured in proper norms. This leads to a simple procedure for determining a ‘‘feasible’’ value ~k for k. However, by applying such a scheme, it is difficult to control the error introduced by solving the monodomain equation instead of the bidomain model. We can simply not guarantee that ~k actually provides a good choice for this parameter. In this paper we will attack this problem by the Output Least Squares (OLS) method. Inspired by the theory for inverse problems, we propose to determine k in such a way that the difference between the solutions of the bidomain and the monodomain equations is approximately minimized. In general, such a procedure would of course require the solution of the bidomain equations. Thus making such a strategy useless - there is no point in solving (9) if the solution of (6), (7) is known! However, by measuring the differences between the involved functions in suitable norms, it turns out that an approximately optimal value for k can be computed without actually solving the bidomain equations. We will see below, in Section 3, that this naturally leads to a timedependent parameter k, i.e. k = k(t). The bidomain model was originally introduced by Tung [25]. It has been thoroughly studied by several scientists, see [6,8,17,12,7,21,22]. Some researchers have also compared the performance and accuracy of monodomain and bidomain simulations, see [9,24]. In [20] a ‘‘reduced bidomain’’, to quote Skouibine and Krassowska, type of scheme was proposed, and some quite impressive simulation results were obtained. However, we are not aware of any results regarding optimal choices of k. This paper is organized as follows: In Section 2 we consider the relation between the monodomain and bidomain models in the case of one space dimension. Section 3 contains our theoretical findings for multidimensional problems, and in Section 4 we derive a modified and optimal monodomain type of scheme.Finally, our numerical experiments and conclusions are presented in Sections 5 and 6. 1.1. Boundary conditions Throughout this text we will, for the sake of simplicity, consider an insulated heart. This leads to the following set of boundary conditions for the bidomain equations (6) and (7): ðM i rv þ M i ruÞ n ¼ 0 ðM e ruÞ n ¼ 0
along oH ;
along oH :
ð11Þ ð12Þ
Furthermore, the transmembrane potential v, as well as the state of the cells s, at time t = 0 is assumed to be known, i.e. vðx; 0Þ ¼ v0 ðxÞ
for x 2 H ;
sðx; 0Þ ¼ s0 ðxÞ
for x 2 H :
Note that, if (8) is satisfied for some k, then (11) and (12) imply that ðM i rvÞ n ¼ 0 along oH ; ðM i ruÞ n ¼ 0 along oH :
ð13Þ ð14Þ
Hence, for the monodomain model it seems to be reasonable to apply homogeneous Neumann boundary conditions both for the transmembrane and extracellular potentials. 2. The bidomain model in 1D In order to shed some light onto the problem outlined above, let us consider the bidomain equations in one space dimension (1D): vt ¼ ðM i vx Þx þ ðM i ux Þx ;
ð15Þ
0 ¼ ðM i vx Þx þ ððM i þ M e Þux Þx ;
ð16Þ
where Mi = Mi(x) and Me = Me(x) are the conductivities. We assume that both u and v have zero derivatives at the boundaries. Obviously, if there exists a constant k such that (8) holds, then (15), (16) can be reduced to
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
279
the monodomain model as explained above. It is, however, worthwhile to note that in 1D a reduction from a bidomain model to a monodomain model is possible also when this condition is not fulfilled. In order to see this, we integrate (16) and get the relation ux ¼
Mi vx ; Mi þ Me
where we have used the assumption that u and v have zero derivatives at the boundaries. Inserting this expression in (15), yields the following scalar equation for v M iM e vt ¼ vx : ð17Þ Mi þ Me x We conclude that the bidomain model can always be reduced to a monodomain type of model when we are addressing a 1D problem. Unfortunately, this seems to be a pure 1D effect. 3. A parameter estimation problem Recall the splitting procedure briefly described in connection with Eqs. (4), (5) and (6), (7). At each time step, both of these subsystems must be solved. We will assume that an efficient solver for the system (4), (5) is available and focus on how to approximate the bidomain equations (6) and (7) with the monodomain model (9) and (10). In view of the splitting approach, it is natural to consider this problem in a time interval [tm1, tm], i.e. in a typical time interval used in the splitting. As mentioned in Section 1, if (8) does not hold, then the solution (v, u) of (6) and (7) will in general not satisfy Eqs. (9) and (10). In order to derive an algorithm for computing appropriate values for k, let us introduce the following notation, cf. (9) and (10), for the monodomain model:3 k r ðM i rwÞ in H ; 1þk 1 r ðM i rwÞ in H ; r ðM i rqÞ ¼ ð1 þ kÞ ðM i rwÞ n ¼ 0 along oH ; wt ¼
ðM i rqÞ n ¼ 0
along oH ;
ð18Þ ð19Þ ð20Þ ð21Þ
where we recall that the boundary conditions (20) and (21) are motivated by the discussion leading to Eqs. (13) and (14). That is, w and q denote the transmembrane and extracellular potentials generated by the monodomain model. Obviously, the solution (w, q) of (18), (21) depends on k, i.e. w ¼ wðx; t; kÞ; q ¼ qðx; t; kÞ: With this notation at hand, our challenge can be formulated as follows: Find k such that the solution (w, q) of (18) and (19) defines an optimal, or close to optimal, approximation of the solution (v, u) of (6) and (7) in the time interval [tm-1, tm]. Expressed in terms of mathematical symbols, we would thus ideally like to solve a problem on the form 2
min kðv; uÞ ðwðkÞ; qðkÞÞk ; k
ð22Þ
where kÆk denotes an appropriate norm.
3
Strictly speaking, an initial condition for the transmembrane potential w is needed, i.e. w(x, tm1) for x 2 H is required. We will return to this issue below.
280
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
In general, applying such an output least squares approach would of course require the solution (v, u) of the bidomain equations in the time interval [tm1, tm]. However, if we are content with considering this problem in terms of a tail-ordered ‘‘energy’’ norm and a semi-discrete version of the PDEs, then it turns out that no explicit knowledge of (v, u) is required in order to solve (22). More precisely, we will discretize the involved equations in time, and show that a suitable output least squares problem for determining k can be easily solved at each time step. Recall that each time interval [tm1, tm] is treated separately, and hence this approach will lead to a time dependent coefficient k = k(t). Consequently, from a global perspective, our algorithm will not necessarily compute an optimal representation for k. Only time-local optimal values for this parameter are obtained by our method. 3.1. Theoretical considerations This section is devoted to a theoretical study of the problem discussed above. Our goal is to investigate whether it is possible to compute optimal, or at least close to optimal, monodomain approximations of (6) and (7) in the time interval [tm1, tm] without actually solving these equations. 3.1.1. Semi-discrete models A semi-discrete approximation of (18) and (19) may be written on the form km r ðM i rwm Þ ¼ wm1 in H ; 1 þ km 1 r ðM i rwm Þ in H ; r ðM i rqm Þ ¼ ð1 þ kÞ
wm Dt
ð23Þ ð24Þ
where km, a constant, denotes k restricted to the time interval [tm1, tm], and Dt ¼ tm tm1 ; wm ðxÞ wðx; tm Þ
for x 2 H ;
m1
w ðxÞ wðx; tm1 Þ for x 2 H ; qm ðxÞ qðx; tm Þ for x 2 H : If km is given, and only the transmembrane potential w is of interest, then it is sufficient to solve (23). However, in the process of determining an appropriate value for km, it turns out that it is convenient to consider both Eq. (23) and (24). We will return to this issue in Section 3.1.4. Let us now focus on the process of advancing the approximate potentials from time step tm1 to tm. Assume that wm1 and qm1 are given,4 and that we want to determine an optimal value for km, to use in (23), (24), for the computation of wm and qm. To this end, we consider the bidomain equations in the time interval t 2 [tm1, tm], i.e. vt ¼ r ðM i rvÞ þ r ðM i ruÞ
r ðM i rvÞ þ r ððM i þ M e ÞruÞ ¼ 0 m1
vðx; tm1 Þ ¼ w m1
ðxÞ
m1
ð25Þ
in H ;
ð26Þ
in H ;
for x 2 H : m1
ð27Þ m1
Since only (w ,q ), and not (v ,u ), are available to us as data, we apply an initial condition expressed in terms of wm1, see (27). The PDEs in the bidomain equations (25) and (26) can be discretized in a similar fashion to the monodomain model, leading to the system vm Dtr ðM i rvm Þ Dtr ðM i rum Þ ¼ wm1 r ðM i rvm Þ þ r ððM i þ M e Þrum Þ ¼ 0
4
That is, wm1 and qm1 have already been computed.
in H ;
in H ;
ð28Þ ð29Þ
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
281
where vm ðxÞ vðx; tm Þ for x 2 H ; um ðxÞ uðx; tm Þ for x 2 H : Note that the initial condition (27) has been incorporated into the right hand side of Eq. (28). In this case, a coupled linear system of two PDEs appears as a results of the discretization procedure. More precisely, Eqs. (28) and (29) defines an elliptic system consisting of two PDEs with two unknowns vm and um. Our goal is to compare the solutions of the systems (28), (29) and (23), (24) and thereby obtain a method for computing an appropriate value for km. More precisely, we will derive a technique for doing so that does not require the solution of the bidomain equations (28) and (29). Prior to presenting our analysis of this problem, we must however introduce some notation. 3.1.2. Function spaces and variational forms Recall that we are considering systems of elliptic equations, cf. (23), (24) and (28), (29). It turns out that the suitable Sobolev spaces in the current situation are o/ o/ V ¼ H 1 ðH Þ ¼ / 2 L2 ðH Þ; ; 2 L2 ðH Þ ; ox oy ð30Þ Z U ¼ w 2 V ; w dx ¼ 0 : H
The variational form of (23) and (24) can now be written on the form: Find (wm, qm) 2 V · U such that Z Z Z km /wm dx þ Dt r/M i rwm dx ¼ /wm1 dx; ð31Þ 1 þ k m H H H Z Z 1 rwM i rqm dx ¼ rwM i rwm dx ð32Þ ð1 þ km Þ H H for all (/, w) 2 V · U. Similarly, the weak form of (28) and (29) reads: Find (vm, um) 2 V · U such that Z Z Z Z m m m /v dx þ Dt r/M i rv dx þ Dt r/ M i ru dx ¼ /wm1 dx; H H H H Z Z rwM i rvm dx þ rwðM i þ M e Þrum dx ¼ 0; H
ð33Þ ð34Þ
H
for all (/, w) 2 V · U. Throughout this paper we will assume that the domain H, the conductivities Mi and Me and the function wm1 are ‘‘well-behaved’’ such that there exist unique solutions (wm, qm) and (vm, um) of (31), (32) and (33), (34), respectively. Further details about this issue can, e.g., be found in [4]. 3.1.3. An inner-product Let us define the following inner-product [Æ, Æ]E on V · U Z Z Z Z ½ð/; wÞ; ðs; nÞE ¼ /s dx þ Dt r/M i rs dx þ Dt r/M i rn dx þ Dt rwM i rs dx H H H H Z þ Dt rwðM i þ M e Þrn dx
ð35Þ
H
for all (/, w), (s, n) 2 V · U, with the associated ‘‘energy’’ norm 2
kð/; wÞkE ¼ ½ð/; wÞ; ð/; wÞE : Note that, for small values of Dt, Z kð/; wÞk2E /2 dx ¼ k/kL2 ðH Þ ; H
ð36Þ
282
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
provided that /, w, Mi and Me are ‘well-behaved’. Consequently, there is a close relation between the kÆkE-norm of (vm, um) and the L2-norm of the transmembrane potential vm. Furthermore, the weak form (33) and (34) of the semi-discrete bidomain equations can now compactly be expressed as: Find (vm, um) 2 V · U such that Z m m ½ð/; wÞ; ðv ; u ÞE ¼ /wm1 dx ð37Þ H
for all (/, w) 2 V · U. 3.1.4. A minimum error analysis Let us now consider the error, measured in the kÆkE-norm, introduced by solving the monodomain model (23) and (24) instead of the more physically accurate bidomain equations (28) and (29) in the time interval [tm1, tm]. To this end, consider the cost-functional 2
Q ¼ Qðkm Þ ¼ kðvm ; um Þ ðwm ; qm ÞkE ¼ ½ðvm ; um Þ ðwm ; qm Þ; ðvm ; um Þ ðwm ; qm ÞE 2
2
¼ kðvm ; um ÞkE 2½ðvm ; um Þ; ðwm ; qm ÞE þ kðwm ; qm ÞkE ; where (wm, qm) and (vm, um) denote the solutions of (23), (24) and (28), (29), respectively. Throughout this section it is also important to keep in mind that both wm and qm depend on km, i.e. wm ¼ wm ðkm Þ; qm ¼ qm ðkm Þ: From Eq. (37), choosing (/, w) = (wm, qm) 2 V · U, we find that Z 2 2 Qðkm Þ ¼ kðvm ; um ÞkE 2 wm wm1 dx þ kðwm ; qm ÞkE : H
m
m
Clearly, (v , u ) is independent of km. Therefore, instead of minimizing Q with respect to km we can consider a cost-functional J on the form Z m m 2 J ðkm Þ ¼ kðw ; q ÞkE 2 wm wm1 dx: ð38Þ H
Note that, in order to minimize J, it is not necessary to solve the semi-discrete bidomain equations (28) and (29)! Thus, we have reached our goal; it is actually possible to compute time-local optimal values km for k without solving the more physically accurate bidomain model. By defining km ¼ max arg min kðvm ; um Þ ðwm ; qm Þk2E ; ð39Þ km 2½;K
our main results may be formulated as follows: Theorem 3.1. Let (wm, qm) = (wm(km), qm(km)) and (vm, um) denote the solutions of Eqs. (31), (32) and (33), (34), respectively. Then km ¼ max arg min J ðkm Þ km 2½;K Z m m 2 m m1 ¼ max arg min kðw ; q ÞkE 2 w w dx ; km 2½;K
H
where kÆkE is the norm defined in (35), (36), and ½; K Rþ represents5 a suitable bound on km. Consequently, it is not necessary to solve the semi-discrete bidomain model (33), (34) in order to compute an optimal value km for km. 5
Appropriate values for e and K will typical depend on the problem under consideration.
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
283
The kÆkE-norm of the solution (wm(km), qm(km)) of (31) and (32) depends continuously on km. Therefore, J(km) will attain its minimum on the closed interval [, K]. The minimum may not be unique, and we have, 2 for the sake of convenience, defined km to be the largest value in [, K] at which kðvm ; um Þ ðwm ; qm ÞkE attains its minimum. By minimizing J with respect to km at each time interval [tm1, tm], we obviously obtain a time dependent parameter k in (18) and (19), i.e. kðtÞ ¼ km
for t 2 ½tm1 ; tm and m ¼ 1; . . . ; M;
where M denotes the number of time intervals used in the discrete model. This approach will clearly yield a piecewise constant function k(t) which depends on the temporal discretization parameter Dt, i.e. k ¼ kðt; DtÞ: From a theoretical point of view, it would certainly be interesting to try to analyze the limit of k(t; Dt) as Dt ! 0. That topic is, however, beyond the scope of the present text. Instead, we will use our analysis to derive a suitable algorithm for solving a modified monodomain type of model. 4. Practical considerations We now turn our attention towards the practical aspects of the problem discussed above. More specifically, our goal is to use Theorem 3.1 to propose a modified monodomain type of scheme suitable for simulation the electrical activity in the heart. Eqs. (30) and (32) imply that Z 1 1 m m q ¼ w w dx ; ð1 þ km Þ jH j H m
ð40Þ
R where jH j ¼ H 1 dx (i.e. the volume (area in 2D) of the domain occupied by the heart H). Thus, the cost-functional J, defined in (38), takes the form Z 2 J ðkm Þ ¼ kðwm ; qm ÞkE 2 wm wm1 dx H 2 Z Z m 1 1 m m ¼ w ; ð1 þ k Þ w w dx 2 wm wm1 dx: m jH j H H E Here 2 Z Z Z m 2Dt m m 2 w ; ð1 þ km Þ1 wm 1 w dx ¼ ðw Þ dx þ Dt rwm M i rwm dx jH j H 1 þ km H H E Z Z Dt rwm M i rwm dx þ rwm ðM i þ M e Þrwm dx 2 ð1 þ km Þ H H 2 Z Z 1 2 ¼ ðwm Þ dx þ Dt 1 rwm M i rwm dx 1 þ km H H Z Dt þ rwm M e rwm dx; 2 ð1 þ km Þ H cf. (35), and consequently Z J ðkm Þ ¼ ðwm Þ2 dx þ Dt H
þ Dt
2
ðkm Þ
Z
1 2
ð1 þ km Þ
ð1 þ km Þ2
H
Z
rwm M i rwm dx
H
rwm M e rwm dx 2
Z
wm wm1 dx; H
ð41Þ
284
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
where it is important to remember that wm depends on km, i.e. wm ¼ wm ðkm Þ: Next, define a¼
km 2 ½aa ; ab 1 þ km
for km 2 ½; K;
where aa ¼
; 1þ
ab ¼
K 1þK
and note that the cost-functional J may conveniently be written on the form Z Z Z Z 2 2 J ðaÞ ¼ ðwm Þ dx þ Dta2 rwm M i rwm dx þ Dtð1 aÞ rwm M e rwm dx 2 wm wm1 dx H
H
H
ð42Þ
H
for a 2 (aa, ab). 4.1. A formula for wm(km) In order to minimize J, more information about how wm depends on km is needed. More precisely, we will now derive a formula for wm(km) = wm(a) and use it to construct an algorithm suitable for estimating the optimal value km for km. In fact, it turns out that km can be approximately determined by minimizing a fourth order polynomial. The details are as follows: Let w denote the continuous solution of the monodomain equation in the time interval [tm1, tm]: wt ¼ ar ðM i rwÞ for x 2 H ; t 2 ðtm1 ; tm ; aðM i rwÞ n ¼ 0 for x 2 oH ; t 2 ðtm1 ; tm ;
ð43Þ ð44Þ
wðx; tm1 Þ ¼ wm1 ðxÞ for x 2 H ;
ð45Þ
with the discrete counterpart (23). Our goal is to characterize how the solution w of this problem depends on km the parameter a ¼ 1þk . It turns out that this can be accomplished by considering a simple auxiliary problem: m rt ¼ r ðM i rrÞ for x 2 H ; t 2 ðtm1 ; tm ; ðM i rrÞ n ¼ 0 for x 2 oH ; t 2 ðtm1 ; tm ;
ð46Þ ð47Þ
rðx; tm1 Þ ¼ wm1 ðxÞ for x 2 H :
ð48Þ
We will now show that the solution w of (43)–(45) is given by the formula wðx; tÞ ¼ rðx; a½t tm1 þ tm1 Þ:
ð49Þ
Clearly, defining l(t) = a[t tm1] + tm1, ½rðx; lðtÞÞt ¼ art ðx; lðtÞÞ ¼ ar ½M i rrðx; lðtÞÞ and we conclude that the function w defined in (49) satisfies the differential equation (43). Moreover, it is easy to verify that this function also fulfills the boundary and initial conditions expressed in Eqs. (44) and (45), respectively. The key point of this analysis is that the solution r of (46)–(48) does not depend on the parameter a. Consequently, the solution w = w(a) of (43)–(45) can be determined for all values of a by simply solving this auxiliary problem. This feature can of course be used to compute an optimal value a* for a (and consequently an optimal value km for km). 4.1.1. Discretization Recall the semi-discrete monodomain equation (23). In a similar fashion we get the semi-discrete form of (46)–(48):
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
rm Dtr ðM i rrm Þ ¼ wm1 m
ðM i rr Þ n ¼ 0
for x 2 H ;
for x 2 oH :
285
ð50Þ ð51Þ
From formula (49) we now obtain the approximation wm ðxÞ wðx; tm Þ ¼ rðx; aDt þ tm1 Þ ð1 aÞrm1 þ arm ¼ ð1 aÞwm1 þ arm ; where we have applied linear interpolation to the function r.
ð52Þ
4.2. A fourth order polynomial Substituting the approximation (52) for wm into formula (42) yields a fourth order polynomial J ðaÞ ¼ c0 þ c1 a þ c2 a2 þ c3 a3 þ c4 a4
for a 2 ½aa ; ab :
ð53Þ
The coefficients c0, c1, c2, c3, c4 can obviously be expressed in terms of integrals involving wm1, rm and their gradients, see Appendix A. km Since a ¼ 1þk , it follows that it is sufficient to minimize a fourth order polynomial in order to determine an m approximation6 of the optimal value km of km to use in the monodomain model in the time interval [tm1, tm]. Recall that our analysis of the monodomain model was presented in terms of an operator splitting procedure for simulating the electrical activity in the heart, cf. the discussion of Eqs. (4)–(7) and Section 3. If a Godunov splitting is applied, our findings lead to the following scheme: Algorithm 4.1. For given positive numbers < K and initial conditions v0 and s0 for v and s, we define the following algorithm: 1. 2. 3. 4. 5.
K aa ¼ 1þ , ab ¼ 1þK Dt = T/M tm = m Æ Dt for m = 0, 1, . . . , M w0 = v0 and s0 = s0 for m = 1 . . . M do (a) Solve
os ¼ F ðs; vÞ; sðtm1 Þ ¼ sm1 ; ot ov ¼ Iðv; sÞ; vðtm1 Þ ¼ wm1 ot (b) (c) (d) (e)
ð54Þ ð55Þ
for x 2 H and t 2 (tm1, tm] by, e.g., a Runge–Kutta method. wm1 = v(tm) Solve (50), (51) for rm Compute the coefficients c0, c1, c2, c3, c4, see Appendix A Compute a ¼ max arg min J ðaÞ ¼ c0 þ c1 a þ c2 a2 þ c3 a3 þ c4 a4 : a2½aa ;ab
m
m1
(f) w = (1 a*)w + a*rm 1 1 R m (g) q ¼ ð1 þ km Þ ðwm jH j wm dxÞ. (This step can be omitted if only the transmembrane potential is H of interest.) a (h) km ¼ 1a ð2 ½; KÞ Step 5(f) can be omitted. It has only been included in the algorithm to make it possible to monitor the time dependency of the optimal choice of k in the monodomain model (18) and (19). This scheme can be modified in a straightforward manner in order to apply our results to a Strang splitting of the system (1)–(3). 6
Keep in mind that linear interpolation was used in (52).
286
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
Note that the computational needs of this algorithm are only slightly larger than that of the standard monodomain approach. More specifically, the additional CPU time required is that of steps 5(d) and (e), and is thus only linked to the construction and minimization of a fourth order polynomial. Such tasks are ‘‘cheap’’ compared to the computational demands of solvers of ODEs and PDEs (steps 5(a) and (c)). 5. Numerical experiments This section is devoted to an experimental study of the scheme proposed above. We will compare the performance of Algorithm 4.1 with that of the traditional monodomain method. The numerical experiments were carried out on a two dimensional domain representing an axial slice through the heart, midway between the apex and the basal plane. The fiber directions were defined to run circumferentially to the left ventricle, see Fig. 1 and [13] for details. This approach leads to the following formulas for the intracellular and extracellular conductivities Mi and Me: M i ¼ rti I þ ðrli rti Þal aTl ; Me ¼
rte I
þ
ðrle
rte Þal aTl ;
ð56Þ ð57Þ
where rti ; rli ; rte and rle are scalars and al = al(x) is a vector function incorporating the fibrous structure of the heart into the PDE models.7 In the traditional monodomain method a suitable value for k must somehow be determined, cf. the discussion of Eqs. (8)–(10). Based on expressions (56) and (57), a fair choice seems to be 2 ~ k ¼ arg min kðrte ; rle Þ kðrti ; rli Þkl2 ¼ 0:80; k
or alternatively ~ k ¼ arg min k
Z
kM e ðxÞ kM i ðxÞk2F dx; ¼ 0:71;
ð58Þ
ð59Þ
H
where kÆkF denotes the Frobenius norm. In all the examples presented below, ~k is equal to 0.80; the results obtained with ~ k ¼ 0:71 were very similar. In order to compare the performance of Algorithm 4.1 with that of the standard monodomain scheme, we also solved the bidomain equations and computed the relative L2-errors of the two methods: kvm wm kL2 for m ¼ 0; 1; . . . ; M; kvm kL2 kvm zm kL2 ¼ for m ¼ 0; 1; . . . ; M; kvm kL2
Eopt m ¼
ð60Þ
Emono m
ð61Þ
where vm, wm and zm denote the transmembrane potentials generated by the bidomain equations, Algorithm 4.1 and the traditional monodomain scheme (with ~k ¼ 0:80), respectively. Throughout this section it is important to keep in mind that the value km for k, see (39) and Theorem 3.1, is optimal with respect to the ‘‘energy’’ norm kÆkE, and not with respect to the classical L2-norm used in (60) and (61). In addition, due to the use of linear interpolation in (52), step 5(h)) in Algorithm 4.1 only provides an approximation of km . All the tests were carried out with a time step Dt = 0.125 ms and 4059 mesh-points in the 2D domain occupied by the ‘‘heart’’ H. 7
In 3D these formulas take the form M i ¼ rti I þ ðrli rti Þal aTl þ ðrni rti Þan aTn ; M e ¼ rte I þ ðrle rte Þal aTl þ ðrne rte Þan aTn ;
where al and an are vector functions used to model the fibrous structure of the heart.
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
287
Fig. 1. An illustration of the fiber directions in the heart.
5.1. Case I Recall that the function I = I(v, s) implements the coupling between the ODEs (1), modeling the cell dynamics, and the bidomain equations (2), (3). For the sake of academic interest, let us now consider the problem obtained by defining I = 0. In this case, the bidomain equations are decoupled from the system (1), and we can thus not expect (1)–(3) to define an accurate model for the electrical activity in the heart. Nevertheless, the resulting equations (2) and (3) for u and v are linear and well-posed, and consequently interesting to consider from both a mathematical and an algorithmic point of view. If I = 0, then the operator splitting approach discussed above is not needed, and steps 5 (a) and (b) in Algomono rithm 4.1 can be omitted. Fig. 2 shows the behavior of the relative L2-errors Eopt for this problem, cf. m and Em the discussion of definitions (60) and (61). Clearly, Algorithm 4.1 delivers far more accurate transmembrane potentials than the traditional monodomain scheme. Furthermore, the optimal value km for k is significantly larger than the value ~ k ¼ 0:80 computed by the ‘‘brute force’’ technique (58). 5.2. Case II In simple heart models, I is typically defined as a cubic polynomial in the transmembrane potential v: IðvÞ ¼ A2 ðv vrest Þðv vthreshold Þðv vpeak Þ;
4.5
x 10
1.7
MonoDomain λ~=0.8 Optimal
4
1.65
3.5
1.6
3
1.55
* λm
Relative error (L2 −norm)
* for I=0 λm
I=0
3
2.5
1.5
2
1.45
1.5
1.4
1
1.35
0.5
0
2
4
6
8
10
time (ms)
12
14
16
18
20
1.3
0
2
4
6
8
10
12
14
16
18
20
time (ms)
mono Fig. 2. The results obtained in Case I. In the left figure we have plotted the relative errors Eopt (solid line) m (dashed line) and Em associated with Algorithm 4.1 and the traditional monodomain scheme, see the discussion of Eqs. (60) and (61). The right figure shows how the value of km , generated in step 5(h)) in Algorithm 4.1, changes with time.
288
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290 * for Cubic λm
Cubic 0.4
~ MonoDomain λ =0.8 Optimal
0.35
1.9
0.3 1.8
0.25
m
1.7 λ*
Relative error (L2−norm)
2
0.2
1.6
0.15 1.5
0.1
1.4
0.05
0
0
2
4
6
8
10
time (ms)
12
14
16
18
20
1.3
0
2
4
6
8
10
12
14
16
18
20
time (ms)
mono Fig. 3. The results obtained in Case II. In the left figure we have plotted the relative errors Eopt (solid line) m (dashed line) and Em associated with Algorithm 4.1 and the traditional monodomain scheme. The right figure shows how the value of km , generated in step 5(h)) in Algorithm 4.1, changes with time.
where we use the parameters A ¼ 0:04; vrest ¼ 85 mV; vthreshold ¼ 65 mV; vpeak ¼ 40 mV; see e.g. [10] for further details. This means that step 5(a) in Algorithm 4.1 takes the form: Solve ov ¼ IðvÞ; vðtm1 Þ ¼ wm1 ; ot for x 2 H and t 2 (tm1, tm]. There is thus no model for the cell dynamics involved, and the operator splitting is only applied in order to handle the nonlinearities introduced into equations (2), (3) by the function I. Fig. 3 shows the relative errors Emono and Eopt m m , along with the value of km , at each time step for this problem. The qualitative behavior of these results are similar to the observations made in Case I; the error associated with Algorithm 4.1 is significantly smaller than that of the standard monodomain scheme, and the value ~ k ¼ 0:80 defined in (58) is far too small. 5.3. Case III In the final test we use the FitzHugh–Nagumo cell model, see [5] and [18]. The results obtained in this case are shown in Fig. 4. During the first 4ms of the simulation, the traditional monodomain method produces results that are slightly more accurate than those obtained with our method.8 For t > 7 ms, the performance of Algorithm 4.1 is considerably better than that of the standard monodomain scheme. 5.4. A remark In the three cases presented above we have seen that the optimal value km for k varies between 1.33 and 2.05. Based on this observation we also performed a series of experiments with k = 1.75 in the standard monodomain scheme. These simulations generated results that were almost as good as those obtained with Algorithm 4.1, indicating that the classical approach might work very well, provided that a suitable value for k is Please keep in mind that the value km for k is optimal with respect to the energy kÆkE-norm, see (39) and Theorem 3.1. Fig. 4 shows the relative L2 errors. 8
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290 λ*m for FitzHugh
FitzHugh 0.35
2.1
~
MonoDomain λ =0.8 Optimal
2
0.3
1.9
0.25 1.8
2
Relative error (L −norm)
289
λ*
m
0.2 1.7
0.15 1.6
0.1 1.5
0.05
0
1.4
0
2
4
6
8
10
12
14
16
18
20
1.3
0
2
4
6
time (ms)
8
10
12
14
16
18
20
time (ms)
mono Fig. 4. The results obtained in Case III. In the left figure we have plotted the relative errors Eopt (solid line) m (dashed line) and Em associated with Algorithm 4.1 and the traditional monodomain scheme. The right figure shows how the value of km , generated in step 5(h) in Algorithm 4.1, changes with time.
used. However, in order to determine such a number for k, we must run Algorithm 4.1; the ‘‘brute force’’ techniques (58) and (59) seem to significantly underestimate the appropriate values. 6. Summary In this paper we have proposed a modified monodomain scheme intended for fast simulations of the electrical activity in the heart. More precisely, a method for computing optimal monodomain approximations of the more physically correct bidomain equations has been introduced and analyzed. As expressed in Theorem 3.1, the key feature of our approach is that the optimal results can be obtained without actually solving the bidomain model, provided that optimality is defined in terms of a suitable ‘‘energy’’ norm. Our theoretical findings were illuminated by a series of numerical experiments. In each of these examples, our method produced significantly more accurate results than the traditional monodomain approach. Furthermore, the new algorithm is only slightly more CPU demanding than the standard scheme. Appendix A The formulas for the coefficients of the polynomial J(a) in Eq. (53) are Z Z m1 2 c0 ¼ w dx þ Dt rwm1 M e rwm1 dx; H Z H Z m1 m1 c1 ¼ 4Dt rw M e rw dx þ 2Dt rwm1 M e rrm dx; H Z Z Z H Z m1 2 m m1 m 2 c2 ¼ w dx 2 r w dx þ ðr Þ dx þ Dt rwm1 M i rwm1 dx H H HZ H Z Z m1 m1 m1 m þ 6Dt rw M e rw dx 6Dt rw M e rr dx þ Dt rrm M e rrm dx; ZH Z H ZH m1 m1 m1 m c3 ¼ 2Dt rw M i rw dx þ 2Dt rw M i rr dx þ 6Dt rwm1 M e rrm dx H ZH ZH m1 m1 m m 4Dt rw M e rw dx 2Dt rr M e rr dx; H H Z Z Z m1 m1 c4 ¼ Dt rw M i rw dx 2Dt rwm1 M i rrm dx þ Dt rrm M i rrm dx HZ HZ H Z m1 m1 m1 m þ Dt rw M e rw dx 2Dt rw M e rr dx þ Dt rrm M e rrm dx: H
H
H
290
B.F. Nielsen et al. / Applied Mathematics and Computation 184 (2007) 276–290
References [1] G.W. Beeler, H. Reuter, Reconstruction of the action potential of ventricular myocardial fibres, J. Physiol. 268 (1977) 177–210. [2] Cellml. Available from: