Verification of computational models of cardiac electro-physiology.

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN BIOMEDICAL ENGINEERING Int. J. Numer. Meth. Biomed. Engng. 2014; 30:525–544 Published online 20 November 2013 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/cnm.2615

Verification of computational models of cardiac electro-physiology Pras Pathmanathan 1, * ,† and Richard A. Gray 2 1 Center

for Devices and Radiological Health, U.S. Food and Drug Administration, Silver Spring, Maryland, USA; Computational Biology Group, Oxford University, UK 2 Center for Devices and Radiological Health, U.S. Food and Drug Administration, Silver Spring, Maryland, USA

SUMMARY For computational models of cardiac activity to be used in safety-critical clinical decision-making, thorough and rigorous testing of the accuracy of predictions is required. The field of ‘verification, validation and uncertainty quantification’ has been developed to evaluate the credibility of computational predictions. The first stage, verification, is the evaluation of how well computational software correctly solves the underlying mathematical equations. The aim of this paper is to introduce novel methods for verifying multi-cellular electro-physiological solvers, a crucial first stage for solvers to be used with confidence in clinical applications. We define 1D–3D model problems with exact solutions for each of the monodomain, bidomain, and bidomain-with-perfusing-bath formulations of cardiac electro-physiology, which allow for the first time the testing of cardiac solvers against exact errors on fully coupled problems in all dimensions. These problems are carefully constructed so that they can be easily run using a general solver and can be used to greatly increase confidence that an implementation is correct, which we illustrate by testing one major solver, ‘Chaste’, on the problems. We then perform case studies on calculation verification (also known as solution verification) for two specific applications. We conclude by making several recommendations regarding verification in cardiac modelling. Copyright © 2013 John Wiley & Sons, Ltd. Received 26 April 2013; Revised 23 July 2013; Accepted 20 October 2013 KEY WORDS:

cardiac modelling; VVUQ; software

1. INTRODUCTION The first computational model of the electrical activity of a cardiac myocyte was developed more than half a century ago [1], and modelling of cardiac electro-physiology (EP) has progressed to the point where whole-organ simulations using highly anatomically detailed geometrical data and extremely complex representations of sub-cellular processes have been performed [2]. Computational cardiac models are a powerful tool for investigating mechanisms of heart function, and there is now an increasing interest in using these models for clinical applications, such as the prediction of ablation sites to terminate arrhythmias [3], optimising lead placement in cardiac resynchronization therapy [4] and optimising defibrillation therapy [5]. For simulation results to be used in clinical settings, the accuracy and robustness of model predictions need to be thoroughly investigated. In order to determine the credibility of computational model predictions, the field of verification, validation and uncertainty quantification (VVUQ)—also known as verification and validation (V&V)—has been developed by the engineering and physical sciences communities, and provides formalism, methodologies and best practices for evaluating the reliability of computational models [6]. To present precise definitions of each of these terms and describe the difference between verification and validation, we must first distinguish between a mathematical model and a computational model. *Correspondence to: Pras Pathmanathan, Center for Devices and Radiological Health, U.S. Food and Drug Administration, Silver Spring, Maryland, USA; Computational Biology Group, Oxford University, UK. † E-mail: [email protected] Copyright © 2013 John Wiley & Sons, Ltd.

526

P. PATHMANATHAN AND R. A. GRAY

We define a mathematical model as the underlying equations used to model a (bio-)physical process, generally derived based on various simplifying assumptions. Mathematical models are most often too complex to solve analytically, and software (code) is created (written) to solve the mathematical model numerically, which is the computational model. Verification, validation and uncertainty quantification are then defined as follows [7] (see also Figure 1): verification is the process of ensuring that the computational model accurately solves the mathematical model, validation is the process of using data to evaluate the extent that the computational model accurately represents the real-world process which it attempts to simulate and uncertainty quantification (UQ) is the process of determining the extent that uncertainty in inputs to the computational model (such as parameters and initial conditions) affects the results of the model. Verification, validation and uncertainty quantification has been used extensively in certain fields to evaluate computational models. For example, the American Society of Mechanical Engineers has produced a guide for V&V in computational solid mechanics [8, 9] and a standard for V&V in computational fluid dynamics and heat transfer [10]. VVUQ is also used by the National Aeronautics and Space Administration [11], the nuclear modelling community [12, 13] and others. For detailed introductions to V&V and VVUQ, see, for example, [6, 7, 14, 15] and references therein. In contrast, rigorous VVUQ-based evaluation of computational cardiac models has rarely been performed. This can be attributed to the complexity of biological models, a lack of knowledge about the field of VVUQ and the fact that currently these models are mostly used as tools for generating hypotheses related to physiological mechanisms. However, stringent testing is clearly required for model use in the clinical setting, and VVUQ provides a formal way to do this [16]. This article is focused on the verification stage for cardiac electro-physiological modelling. Such models can either be of a single isolated cardio-myocyte only, in which case they are referred to as cell models or they can be multicellular and used to represent tissue or the whole-organ. Single-cell models are almost always systems of ODEs; tissue models are composed of a cell model coupled to PDEs. An illustration of the components of a multi-cellular cardiac EP solver is given in Figure 1. There are also two separate parts to verification, named code and calculation verification, which we begin by defining in the succeeding texts. We then briefly state why both types of verification are relatively straightforward in the single-cell case, before reviewing verification of tissue-level models, which will motivate the aims of this paper. 1.1. Code verification and calculation verification The first part of verification is code verification, which is the process of confirming that the software correctly implements the required algorithms [6]. For software which solves ODEs or PDEs, a major aspect of this is confirming that the numerical solution converges to the true solution of

Figure 1. Left: stages in developing a computational model imitating reality; assumptions and analysis are used to develop a mathematical model, which is then solved numerically after choosing a numerical method and developing computational software. Verification compares computational model results with the underlying mathematical model, whereas validation compares with the real-world process. Right: components of a multi-cellular cardiac electro-physiology solver, with major parts that require verification highlighted in grey. Copyright © 2013 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Biomed. Engng. 2014; 30:525–544 DOI: 10.1002/cnm

VERIFICATION OF COMPUTATIONAL MODELS OF CARDIAC ELECTRO-PHYSIOLOGY

527

the ODEs/PDEs as discretisation parameters (spacestep and/or timestep) are decreased and ideally checking that convergence is at the rate predicted by theory. The main approaches to code verification are to solve problems with known analytic solutions or to compare several independently written codes on a common problem. The former approach allows exact errors and convergence rates to be determined but is likely only possible for simplified (non-physical/non-physiological) problems. The latter approach has the advantage of being feasible on more complex problems but suffers from a lack of a gold-standard solution. Comparisons of codes have been performed on astronomy codes [17], meteorology codes [18], seismic data processing codes [19] (where results are described as ‘deeply disturbing’) and on cardiac EP codes in [20] (discussed later). The two approaches are complementary, and it is important that fields develop a set of well-defined community-wide benchmark problems, with and without exact solutions, to aid code verification. The second stage is calculation (or solution) verification, which is the process of estimating or bounding the numerical error in the particular application simulation, that is, when the solver is used to make predictions for the application of interest, and with a given final choice of parameters, computational mesh resolution, timestep, and so on [6]. One major recommendation in [6], a white paper on VVUQ, is that calculation verification (as well as validation and uncertainty quantification) is performed explicitly for the quantity of interest (QOI) to be analysed in regard to the application. Whilst the application dictates the QOI choice, the nature of the QOI, in particular the sensitivity of the QOI to potential coding errors and to numerical parameters, will significantly impact the ability to perform verification. (Similarly, the sensitivity of the QOI to modelling parameters significantly impacts validation and uncertainty quantification (UQ), and overall the nature of QOI is a major factor in the ability to perform VVUQ and develop a credible model). Examples of methods for estimating or bounding the numerical error in a QOI are the following: computing the QOI on several meshes of increasing resolution, including using Richardson extrapolation to approximate the true solution [14]; and analytical error estimation methods [21], including adaptive mesh refinement methods [22]. For more information on verification in general, see [7] or [6]. 1.2. Verification of single-cell models As mentioned previously, single-cell cardiac EP models are almost always systems of ODEs. Dozens of models exist, and they can be extremely complex, involving dozens of equations and over a hundred parameters. See [23] for a discussion. If a general ODE solver is implemented to solve a cell model, code verification of the solver is relatively straightforward—one may simply test the solver on (preferably stiff) ODEs with known solutions. If a purpose-built solver is used, which takes into account the structure of a particular cell model for increased efficiency (for example, an implementation of [24]), direct verification is more difficult. However, one straightforward and computationally cheap option in this case is to check that results agree when compared with a general ODE solver, such as a verified commercial ODE solver. Whilst it is relatively easy to verify the solver, it is important to also be confident that the cell model has been implemented correctly. Because the models are so complex, such implementations are not trivial, and any attempt to translate equations into code manually is likely to lead to errors. For this reason, the CellML language [25] has been developed to allow cell models to be stored in an XML format, and automatically and reliably converted to other forms (e.g. C/Fortran/Matlab code). There are several hundred cell models stored in the CellML repository‡ . Tools are available for automatic and reliable conversion to human-readable output or various coding languages [26], and significant work falling broadly under the field of verification has been carried out by the CellML community, for example, on model curation [27] or units and dimensional consistency [25, 28]. Calculation verification is also relatively straightforward, compared to tissue-level simulations. The use of adaptive ODE solvers§ , for example, allows explicit control of numerical error and

‡ §

http://www.cellml.org Such as the sundials adaptive ordinary differential equation solver (https://www.llnl.gov/casc/sundials), which has been used extensively in cardiac modelling.

Copyright © 2013 John Wiley & Sons, Ltd.


528


therefore facilitates calculation verification. We reiterate though that calculation verification ought be performed for the specific QOI for the application being studied. 1.3. Verification of tissue-level models In contrast to single-cell models, verification of tissue-level or organ-level cardiac EP is highly challenging and the focus of this paper. We will concentrate solely on continuous representations of cardiac EP. Such models are created by coupling a chosen cell model with PDEs governing the spatio-temporal behaviour of the electrical activity. There are three commonly used formulations of tissue-level cardiac EP: the monodomain formulation, which provides the trans-membrane voltage V throughout the tissue as a function of space and time; the more sophisticated bidomain formulation, which provides voltage V and extracellular potential e ; and bidomain case where an extra-cardiac electrically conductive bath (e.g. torso) is also modelled, which we will refer to as the bidomain-with-bath formulation. The governing equations for each of these will be given in Section 2.1. Many computational models (software) have been developed to solve some or all of these mathematical models, including CARP (Cardiac Arrhythmia Research Package) [29], Chaste (Cancer, Heart and Soft Tissue Environment) [30, 31], CMISS/openCMISS (Continuum Mechanics, Image analysis, Signal processing and System Identification) [32], Continuity¶ , and Fenics [33]. These computational models use the finite difference method or FEM to solve the equations, together with a variety of numerical methods. In all cases, code verification is required to confirm that the solver has been implemented correctly. As mentioned previously, one approach to verify code is comparison with other solvers. There has only been one attempt, to the best of our knowledge, at comparing the results of several large-scale cardiac EP solvers on a common problem to confirm that they agree. In Niederer et al. 2011 [20], 11 solvers were compared on an unambiguously defined benchmark problem, using the monodomain formulation, a regular geometry and a complex cell model. Surprisingly, large differences in the results were observed between the solvers, later explained in [34]. The alternative, and complementary, approach is to use problems with analytic solutions. As far as we are aware, there are no known cardiac EP problems in two or three dimensions with realistic boundary conditions and physiological cell models that can be solved analytically. Some simplified physiological problems may be solved analytically for travelling wave solutions in one-dimensional (1D), for example, [35, 36]; however, because these problems are defined in 1D and on the entire real line, they provide limited utility for testing general codes. Approximate analytical solutions are possible in higher dimensions if asymptotic analysis is performed (e.g. [37]), but such non-exact solutions are also of limited use for verifying solvers. One major contribution of the present paper is to define a set of model problems with exact solutions for each of the monodomain, bidomain and bidomain-with-bath settings, in 1D, 2D and 3D. These problems use a simplified non-physiological cell model (it is highly unlikely that a problem with a physiological cell model may be solved exactly). The model problems are carefully constructed so that they can be used with general codes, and are, to the best of our knowledge, the only fully coupled 2D and 3D cardiac EP problems with natural boundary conditions and non-trivial cell models for which the exact solutions are known. As such, the model problems will be of great value for verifying solvers. However, they should not be considered as the only method for verifying a cardiac EP solver. Rather, they should be considered as part of a suite of tests for verifying a tissue-level solver. Such a suite of tests might comprise the following: 1. Decouple the cell model solver and test the ODE solver in isolation (Section 1.2). 2. Confirm that the complex cell model has been implemented correctly (Section 1.2). 3. Decouple the PDEs from the cell model (set Iion in Section 2.1, see later, to zero) and test the PDE solver against exact solutions in all dimensions. This is straightforward for the monodomain equation because the decoupled system is just the heat equation, for which problems with exact solutions are easy to construct. However, this is less easy for the bidomain ¶

http://continuity.ucsd.edu




529

and bidomain-with-bath settings; also, there are as far as we are aware no established benchmark problems for any bidomain cases. 4. Solve the model problems defined in this paper, fully coupled problems with a simplified cell model, in all dimensions, for each of the monodomain, bidomain and bidomain-with-bath settings, and confirm that the error converges to zero at expected rates. 5. Solve the Niederer benchmark problem [20], a fully coupled problem with realistic cell model (but 3D monodomain only), and compare results with other solvers. 6. Carefully test all aspects of the code that are not covered in the aforementioned tests, such as the use of spatially varying conductivities, post-processing methods and so on. The first three of these involve testing individual components of the entire solver, illustrated in Figure 1. Stages 4 and 5 however involve the entire solver, for a particular choice of cell model. The importance of the model problems defined in this paper is easily seen if one compares the information provided by these stages with and without stage 4. The previous texts also illustrates that there is currently a lack of well-defined benchmark problems with which to test cardiac EP solvers, which prohibits both verification and the ability to provide evidence of strong code verification— both of which are of clear importance if solvers are to be used in future clinical decision-making. The Niederer benchmark test [20] (physiological, monodomain-only and no exact solution) began to address this deficiency, and the model problems in this paper (non-physiological and exact solutions) complement the Niederer benchmark and constitute a substantial step forward. 1.4. Structure of the paper The remainder of this paper is devoted to verification of tissue-level and organ-level cardiac EP models. Section 2 deals with code verification. Here, we state the governing equations of the three formulations considered and then provide novel model problems with exact solutions for all three formulations and on 1D, 2D and 3D domains. These problems are carefully constructed so that they can be easily run by most cardiac solvers and can be used to vastly increase confidence that a general cardiac EP solver is implemented correctly. We will test the cardiac EP solver ‘Chaste’ [30, 31] on the model problems. Section 3 is concerned with calculation verification. Here, we investigate the feasibility of performing rigorous QOI-based calculation verification in cardiac applications, using two case studies. Finally, we discuss our results in Section 4, where we also make a number of general recommendations. 2. CODE VERIFICATION USING MODEL PROBLEMS WITH EXACT SOLUTIONS 2.1. Governing equations Let represent the domain, with boundary @. We start with the bidomain equations, because the monodomain equations are a special case of the bidomain equations. The bidomain equations are derived by performing a mathematical homogenisation of equations representing conservation of current within the intra-cellular and extra-cellular spaces; in the resultant formulation, the intracellular and extra-cellular electrical potentials, i and e , are defined at all points in . Defining the transmembrane voltage V D i e , the governing equations are [35] @V Cm (1) C Iion .u, V / r .i r .V C e // D Ii.stim/ , @t .stim/ r ..i C e / r e C i r V / D Itotal ,

(2)

@u D f .u, V / , @t

(3)

where Cm is the specific capacitance of the cell membrane, is the membrane surface-area-tovolume ratio, i and e are intra-cellular and extra-cellular conductivity tensors, Ii.stim/ is an Copyright © 2013 John Wiley & Sons, Ltd.


530


.stim/ intra-cellular volume stimulus current per unit volume and Itotal D Ii.stim/ CIe.stim/ , where Ie.stim/ .stim/ is an extra-cellular volume stimulus (usually (implicitly) chosen so that Itotal D 0). u u.t , x/ is a vector of state variables representing the current state of the cell at that location, and Iion and f are prescribed functions, which together make up the cell model. Typical boundary conditions for (1) and (2) are the specification of zero current across the boundary:

n .i r .V C e // D 0 n .e r e / D 0

on @ on @

(4) (5)

where n is the outward-pointing unit normal vector. The system of equations (1)–(5) is then completed by specifying suitable initial conditions for V and u. e does not require initial conditions and is only specified up to a constant function of time. The monodomain equations can be derived from the bidomain equations when the intra-cellular and extra-cellular conductivities are proportional: @V Cm (6) C Iion .u, V / r . r V / D I .stim/ , @t @u D f .u, V / , @t

(7)

K i where K is the proportionality where is the effective conductivity (defined by D 1CK .stim/ is a stimulus current. Typical boundary conditions are constant, e D Ki ) and I

n . r V / D 0 .

(8)

Finally, the bidomain-with-bath equations are used to model the case of cardiac tissue contained in a conductive bath (e.g. the human torso). Let b represents the bath domain, assumed to surround . V is defined only in (i.e. only in the tissue), and (1)–(3) still hold within . However, e is now defined everywhere on [ b (i.e. throughout the tissue and the bath), and outside the tissue satisfies r .b r e / D 0 in b ,

(9)

where b is the conductivity of the bath (usually a scalar). The boundary/interface conditions are the following: (4) (zero flux of i across @), continuity of e across @ and continuity of the extracellular current across @—the extracellular current flowing out of the tissue, e re n, should be equal to that entering the bath, b re n, everywhere on @. The remaining boundary condition on the edge of the bath domain is n .b r e / D IE.surf/

on @b n@ .

(10)

where IE.surf/ is a stimulus current (per unit area) applied to the edge of the bath domain and may be used to represent defibrillating electrodes. See Figure 2, later, for an example of this set of boundary R conditions. (Note that the prescribed function IE.surf/ should satisfy IE.surf/ dS D 0 for a solution to exist (conservation of current), or alternatively a Dirichlet boundary condition on e should be applied somewhere on @b n@, corresponding to a ground electrode). 2.2. Model problems with exact solutions We now provide a method for computing exact errors in fully coupled solvers, using model 1D–3D problems with exact solutions for each of the monodomain, bidomain and bidomain-with-bath settings. We use the ‘Method of Manufactured Solutions’ [7] to generate these problems, which means we choose (non-physiological) V .t , x/, u.t , x/ and, for bidomain, e .t , x/, and insert them into the governing equations, to determine what the forcing terms (Iion , f) should be for that choice of V and Copyright © 2013 John Wiley & Sons, Ltd.



531

Figure 2. For the monodomain and bidomain model problems, the domain is just the unit line, square or cube, depending on dimension. For the bidomain-with-bath problems, the geometries are illustrated here in 1D, 2D and 3D (top). The bath domain (b ) is white/thin, and the tissue domain () is shaded/thick. V is defined in , e in both. Surfaces with different boundary/interface conditions are labelled. The boundary conditions are the following: b r e n D ˛ on region A, ˛ on region B and 0 on region D; i r .V C e / n D 0 on regions C and E. Interface conditions are the following: e is continuous across region C and b r e n D e r e n on region C. Bottom p left: voltage solution of M-2D and B-2D at t D 1, which is also the voltage initial condition scaled by 2. Bottom right: solution of BB-1D at t D 1, which is also the solution of BB-2D* or BB-3D* at any y or ´ (because the solution of BB-2D* or BB-3D* is independent of y and ´).

u to be a solution. The goal, however, is to construct problems which can easily be solved with a typical cardiac EP solver, which means we need to choose V .t , x/ and u.t , x/ carefully. In particular, we would like to define problems that can be solved by cardiac EP solvers without the need for additional implementation code, excluding the implementation of a simple new cell model. We .stim/ thus wish to choose V .t , x/ and u.t , x/ such that (4), (5) and (8) still hold, such that Itotal is zero, and such that the resultant Iion and f can be written as function of the chosen solution V and u. It would be trivial to construct a problem for which Iion Iion .t , x/ and f f.t , x/, but this would likely require solvers to be altered to run the problem. We also wish for Iion and f to be dependent on both V and u, and that I stim and Iistim are zero. Finally, we would like to define problems involving functions that are not overly complex, to mitigate the possibility of errors upon implementation. These constraints make the construction of problems with exact solutions non-trivial. The following non-physiological, three-variable cell model has been constructed for use in all the model problems: u D .u1 , u2 , u3 /, specified by 2 3 .u1 C u3 V /2 u22 C 12 .u1 C u3 V /u22 .V u3 / 5 f.u, V / D 4 (11) .u1 C u3 V /u32 0 Iion .u, V / D

Cm ˇ.V u3 / .u1 C u3 V /u22 .V u3 / C 2

(12)

where ˇ is a free parameter. As required, f and Iion are nonlinear functions of V and u, and not especially complicated. Clearly, this is essentially a two-variable model; u3 is constant in time and just a spatially varying parameter. In fact, it will be zero for the monodomain and bidomain problems but will be convenient for the bidomain-with-bath problem. Copyright © 2013 John Wiley & Sons, Ltd.


532


Two spatially varying functions, F F .x/ and G G.x/, are used in the model problems. Whilst specific choices will be made, in general, we just need the following. Let G be any function which is never zero on . Let O be a symmetric positive-definite matrix and choose F satisfying the eigenfunction problem: r O r F D cF

and

O r F n D 0 on @

(13)

for some c 2 R. An example is D Œ0, 1n (the unit line/square/cube), O D diag.s1 , : : : , sn / and F .x/ D …niD1 cos.mi xi /, where si > 0 and mi are integers, for which (13) is satisfied with c D 2 †niD1 si m2i . In the succeeding texts, we first define a general monodomain model problem with an analytic solution, which holds in all dimensions and any , O and F satisfying (13), including spatially varying O . It may therefore be used for the (future) development of more advanced model problems. We then define three concrete monodomain model problems, referred to as M-1D, M-2D and M-3D, which involve making a specific choice of , F , G and parameters. Similarly, we define a general bidomain problem and then three concrete problems B-1D, B-2D and B-3D, as well as three bidomain-with-bath problems, named BB-1D, BB-2D* BB-3D*. As we shall discuss, it is difficult to construct genuinely 2D and 3D problems for the bidomain-with-bath case. The bidomain-withbath model problems, although posed in all three dimensions, are essentially 1D problems. The stars in BB-2D* and BB-3D* denote this. Note that we do not prescribe a discretisation or solution grid when defining the model problems. The model problems are independent of the choice of numerical method for solving the mathematical models—this allows them to be used to test any numerical scheme proposed for solving the governing equations. 2.2.1. Monodomain problems. The following holds in all spatial dimensions and for any F F .x/ and G G.x/ as defined previously: Let > 0, Cm > 0 and I (stim) D 0. Setting D O , using the cell model (11)–(12) with ˇ D c, and taking the initial conditions V .0, x/ D F .x/ and u.0, x/ D .G.x/ C F .x/, G.x/1=2 , 0/, the exact solution of the monodomain equations (6)–(8) is V .t , x/ D .1 C t /1=2 F .x/

(14)

u1 .t , x/ D .1 C t /G.x/ C .1 C t /1=2 F .x/

(15)

u2 .t , x/ D .1 C t /1 .G.x//1=2

(16)

u3 .t , x/ D 0

(17)

This can be easily confirmed by substituting (14)–(17) in (6)–(8), (11)–(12). As concrete examples, we define the following three monodomain model problems, which we will use to evaluate Chaste in Section 2.3. Defining these problems involves making somewhat arbitrary choices for parameters values and G. See also Figure 2. Problem M-1D: Take the domain to be D Œ0, 1. Let F .x/ D cos.x/ and G.x/ D 1 C x. Using the parameter values given in Table I, solve the monodomain equations (6)–(7) with I .stim/ D 0, using cell model (11)–(12), zero-flux boundary conditions (8), and initial conditions V .0, x/ D F .x/ and u.0, x/ D .G.x/ C F .x/, G.x/1=2 , 0/. The exact solution is (14)–(17). Problem M-2D: As M-1D, except take the domain to be the unit square D Œ0, 1 Œ0, 1, and use F .x/ D cos.x/ cos.2y/ and G.x/ D 1 C xy 2 (and the appropriate column of Table I). Problem M-3D: As M-1D, except take the domain to be the unit cube D Œ0, 1 Œ0, 1 Œ0, 1, and use F .x/ D cos.x/ cos.2y/ cos.3´/ and G.x/ D 1Cxy 2 ´3 (and the appropriate column of Table I). The exact solution in all three cases is (14)–(17). 2.2.2. Bidomain problems. Similarly, we can now define bidomain problems with exact solutions: (stim) Let > 0, Cm > 0 and Ii(stim) D Itotal D 0. Choose k 2 .0, 1/ and let i D O and e D 1k O . k Using the cell model (11)–(12), this time with ˇ D c.1 k/, and again choosing initial condiCopyright © 2013 John Wiley & Sons, Ltd.



k e ˛ sb ˇ

i

3 2

Cm

1.1

1D (see caption) -

M-1D

Parameter

5.9

-

2D

3 2

M-2D

8.6

-

3D

3 2

M-3D

3D p 1= 2 .1 k/i =k 8.6.1 k/

2D p 1= 2 .1 k/i =k 5.9.1 k/

1D (see caption) p 1= 2 .1 k/i =k 1.1.1 k/

3 2 -

B-3D

3 2 -

B-2D

3 2 -

B-1D

p 1= 2 .1 k/i =k 0.01 se =2 1.1.1 k/

1D

3 2 -

BB-1D

p 1= 2 .1 k/i =k 0.01 se =2 1.1.1 k/

2D

3 2 -

BB-2D*

p 1= 2 .1 k/i =k 0.01 se =2 1.1.1 k/

3D

3 2 -

BB-3D*

Table I. Parameter values to be used in model problems, where 1D D 1.1 2 , 2D D 2 diag.1.1, 1.2/ and 3D D 2 diag.1.1, 1.2, 0.3/. As model problems are non-physiological, parameters are dimensionless. Most can be arbitrarily chosen; any choice of > 0, Cm > 0, ˛ 2 R, k 2 .0, 1/ and sb > 0 would work. se is defined in the text. ˇ is dependent on other parameters.


533


534


tions V .0, x/ D F .x/ and u.0, x/ D .G.x/ C F .x/, G.x/1=2 , 0/, the exact solution of the bidomain equations (1)–(5) is (14)–(17) together with e .t , x/ D k.1 C t /1=2 F .x/ C C.t / ,

(18)

where C.t / is an arbitrary function of time. Problem B-1D: Take the domain to be D Œ0, 1. Let F .x/ D cos.x/ and G.x/ D 1 C x. Using .stim/ parameter values given in Table I, solve the bidomain equations (1)–(3) with Ii.stim/ D Itotal D 0, using cell model (11)–(12), with zero-flux boundary conditions (4)–(5), and initial conditions V .0, x/ D F .x/ and u.0, x/ D .G.x/ C F .x/, G.x/1=2 , 0/. Problem B-2D: As B-1D, except take the domain to be the unit square D Œ0, 1 Œ0, 1, and use F .x/ D cos.x/ cos.2y/ and G.x/ D 1 C xy 2 (and the appropriate column of Table I). Problem B-3D: As B-1D, except take the domain to be the unit cube D Œ0, 1 Œ0, 1 Œ0, 1, and use F .x/ D cos.x/ cos.2y/ cos.3´/ and G.x/ D 1 C xy 2 ´3 (and the appropriate column of Table I). The exact solution in all three cases is (14)–(17) and (18). See also Figure 2. 2.2.3. Bidomain-with-bath problems. It is not straightforward to extend the ideas behind the aforementioned problems to the bidomain-with-bath case. One difficulty is that if the tissue region, , is completely contained within the bath region, b , the bath domain is non-convex, making it difficult to construct problems with an exact solution to Laplace’s equation (9). We relax this assumption, choosing a domain which does not have completely contained in b . We will use a region of tissue surrounded on two sides by bath, as shown in Figure 2. Also, the additional complexity of the bidomain-with-bath equations makes it difficult to construct a genuinely 2D or 3D problem. Hence, the problems in the succeeding texts are all essentially 1D problems, in that the solution is dependent on x only, and not y or ´, regardless of dimension. The construction of genuinely 2D/3D model problems for this case is left as an open problem|| . Problem BB-1D: Let al l D Œ1, 2 and define b D fx 2 al l W 1 6 x 6 0 or 1 6 x 6 2g D fx 2 al l W 0 6 x 6 1g that is, tissue surrounded on two sides by bath, as shown in Figure 2. Let F .x/ D cos.x/. Let G.x/ D 1 C x. For ease of notation, write se for .e /11 . Using parameter values in Table I, solve the bidomain-with-bath equations as described in Section 2.1, using cell model (11)–(12), .stim/ a scalar bath conductivity b D sb , setting Ii.stim/ D Itotal D 0, using an electrode stimulus of 8 < ˛ if x D 1 if x D 2 IE.surf/ D ˛ : 0 otherwise and initial conditions V .0, x/ D F .x/ exact solution of this problem is:

˛x se

and u.0, x/ D .G.x/ C F .x/, G.x/1=2 , ˛x /. The se

V .t , x/ D .1 C t /1=2 F .x/

˛ x se

8 ˛ 1=2 ˆ if 1 6 x 6 0 < k.1 C t / C sb x C C.t / ˛ 1=2 if 0 6 x 6 1 e .t , x/ D k.1 C t / cos.x/ C se x C C.t / ˆ : k.1 C t /1=2 cos./ C ˛ C ˛ .x 1/ C C.t / if 1 6 x 6 2 se sb ||

(19)

(20)

We remark that extending the one-dimensional bidomain-with-bath problem (BB-1D) to a cylindrical geometry will likely lead to a solution involving Bessel functions.




535

u1 .t , x/ D .1 C t /G.x/ C .1 C t /1=2 F .x/

(21)

u2 .t , x/ D .1 C t /1 .G.x//1=2

(22)

u3 .t , x/ D

˛ x se

(23)

where C.t / is an arbitrary function of time. Figure 2 includes plots of (19) and (20). Problem BB-2D*: As BB-1D except use al l D Œ1, 2 Œ0, 1, G.x/ D 1 C xy 2 (but do not change F .x/), and the appropriate column of Table I. The solution is (19)–(23). Problem BB-3D*: As BB-1D except use al l D Œ1, 2 Œ0, 1 Œ0, 1, G.x/ D 1 C xy 2 ´3 , and the appropriate column of Table I. The solution is (19)–(23). We remark that if we want to include the case of a ground electrode, this problem can be modified to apply the Dirichlet boundary condition e D 0 on x D 1 and set IE.surf/ D ˛ on x D 2 and zero on the other boundaries. Then the exact solution is just (19)–(23) with C.t / D k.1 C t /1=2 C s˛b . 2.3. Evaluation of the Chaste electro-physiology solver on model problems We now test the Chaste cardiac EP solver [30, 31]** on the model problems defined in Section 2.2. The end time is taken to be T D 1. The code for running these models problems using Chaste can be viewed at, or downloaded from, https://chaste.cs.ox.ac.uk/trac/wiki/PaperTutorials. We can measure the error of the numerical solutions using various norms. Let k RkL2 and k kH 1 represent the standard L2 and Sobolev H 1 spatial norms, that is, kvk2L2 D v 2 dV R and kvk2H 1 D v 2 C jr vj2 dV . For a general function w w.t , x/ defined on Œ0, T , the spatio-temporal L1 .L2 / and L2 .H 1 / norms are defined by†† kwkL1 .L2 / D max kw.t /kL2 ./

(24)

t 2Œ0,T

Z kwkL2 .H 1 / D

! 12

T 0

kw.t /k2H 1 ./

dt

,

(25)

Let h and t represent the spacestep and timestep, respectively. Chaste uses linear basis functions in the finite element implementation of the EP equations and a first-order approximation of the time-derivatives. The error in the L2 .H 1 /-norm, kVcomputed Vexact kL2 .H 1 / , is then expected to be O.t C h/, whereas the error in the L1 .L2 /-norm is expected to be O.t C h2 / [38]. We let h D 0.1, 0.05, 0.025, : : :, and choose the timesteps to be proportional to h2 , in which case O.t C h/ D O.h/ and O.t C h2 / D O.h2 /. Hence, we expect linear convergence in the L2 .H 1 /-norm and quadratic convergence in the L1 .L2 /-norm. These errors can be computed using numerical quadrature for spatial integrals and the trapezoid rule for time integrals. We note, because it is not possible to discern from the following convergence plots, that the numerical solutions are generally visually indistinguishable from the analytical solutions (such as that shown in Figure 2.2) even on the coarser meshes used. The errors when Chaste is used to solve the monodomain model problems M-1D, M-2D and M-3D are plotted in Figure 3. We see that all solutions converge to the exact solution of the problems, and at the expected rates of convergence, in both the L1 .L2 / and L2 .H 1 / norms. Figure 3 also plots the error in both V and e using the bidomain model problems B-1D, B-2D and R B-3D. (Note: for bidomain, the numerical method used in Chaste is such that the solution with e d3 x D 0 is found [39], allowing us to determine that C.t / D 0). Again, the errors converge to zero at the expected rates. ** ††

Version 3.1, revision 18000. More precisely, essential supremum instead of maximum in the definition of the L1 norm, although these are equivalent for our purposes.



536

P. PATHMANATHAN AND R. A. GRAY M-1D, M-2D, M-3D 10

Error

10 10 10 10 10

1 0

-1

3D

-2 -3

1D

2D

-4

10

-2

10-1

B-1D, B-2D, B-3D 10 10 10 10 10

B-1D, B-2D, B-3D

1

10

0

Error in φe

Error in V

10

-1 -2 -3

1D

3D 2D

10 10 10 10

-4

10

-2

10

10

-1

1 0

-1 -2

2D

-4

10

-2

BB-2D* 10 10 10 10

-1

10

BB-2D*

0

0

10

-1

Error in φe

Error in V

10

1D

3D -3

-2 -3 -4

-1

10

-2

10

-3

10

-4

10

-2

10

Spatial stepsize, h

-1

10

10

-2

-1

10

Spatial stepsize, h

Figure 3. Evaluation of Chaste using the model problems. Top: error in voltage, in the L2 .H 1 / (dashed lines) and L1 .L2 / (solid lines) norms, against spatial stepsize, for monodomain model problems, M-1D, M-2D and M-3D. The triangles have gradient exactly 1 and 2, indicating that convergence is O.h/ in the L2 .H 1 /-norm and O.h2 / in the L1 .L2 /-norm, the required rates of convergence. (Note: less points are computed in 3D for reasons of computational cost (recall we use t / h2 )) Middle: similarly, error in voltage (left) and extra-cellular potential (right) for bidomain model problems, B-1D, B-2D, B-3D. Bottom: error in voltage (left) and extra-cellular potential (right) for 2D bidomain-with-bath model problem, BB-2D*.

For the bidomain-with-bath case, we only run BB-2D* to emphasise that, as stated previously, all the bidomain-with-bath problems are essentially 1D problems. The results are also given in Figure 3. R2 (The numerical method in this solver is such that the solution satisfying 1 e dx D 0 is found [39], again allowing us to determine C.t /). The errors in V and e in both norms are displayed; again, we see that the equations are solved correctly and at the expected rates of convergence. It is instructive to consider which aspects of Chaste have been tested by these model problems and how much confidence the results give us. Chaste has numerous unit tests of all aspects of the functionality [30]; we neglect the existence of such tests for the purposes of this discussion. The results give good confidence that the monodomain implementation in Chaste is correct in all dimensions for this cell model, and because Chaste auto-generates code for cell models directly from CellML files, good confidence that the implementation is correct for general cell models. The stimulus is always zero in the model problems, so correct implementation of stimulus currents needs testing separately. Because the model problems all use a simple geometry, functionality related to reading in of complex meshes also needs separate testing, as does functionality related to setting up spatially varying conductivity tensors from spatially varying fibres. The same conclusions, with the same caveats, hold for the bidomain implementation. As for the bidomain-with-bath implementation, whilst in principle the model problem could be passed even if there was a major error in how y and ´ derivatives are computed (because BB-2D* is essentially a 1D problem), the fact that the 2D and 3D bidomain tests were passed implies this is not the case, and we can be reasonably confident the bidomain-with-bath implementation is correct, with the same caveats, and with extra testing required to check non-scalar bath conductivities is implemented correctly. Copyright © 2013 John Wiley & Sons, Ltd.



537

2.4. Evaluation and limitations of the model problems Excluding the fact that the bidomain-with-bath model problems are essentially 1D, there are a few aspects of the model problems which might appear unsatisfactory, which we now discuss. The question to be asked when critiquing the model problems is what kind of errors could be present in a cardiac EP solver that would still lead to the correct solution of these problems being computed, that is, what kind of errors would not be picked up upon by solving the model problems. So, for example, whilst none of the solutions involves any propagation‡‡ or high spatial or temporal gradients, it is not clear that this is a failing in this context. Similarly, although the bidomain model problems use proportional conductivities tensors (which means that analytically the bidomain problems could be reduced to the monodomain problems), this does not mean it does not provide a good test of a bidomain solver. As long as the solver does not use different code for the case of proportional conductivities, correctly solving the bidomain problems provides good confidence that the solver will correctly solve the equations in the case of non-proportional conductivities. On the other hand, one failing in the model problems is that if a finite difference solver§§ incorrectly implemented r V n D 0 rather than (8); the exact solutions will still be found. The models problems only work with the cell model (11)–(12). Therefore, they can only be used to verify solvers which can run general cell models, or solvers for which new cell models can be implemented. However, because there are over 100 cardiac cell models in the literature, and no established ‘best’ model, it is reasonable to expect a large-scale solver to be able to handle a new cell model. Some solvers may use numerical methods for solving the cell model that take into account the structure of the cell model. In this case, solving the model problems will only verify other aspects of the solver (PDE solver and implementation of the coupling). Verification of the part of the solver used to solve the specific cell model would be more difficult; this falls under the discussion of purpose-built solvers in Section 1.2. Finally, we re-iterate that the models problems should be considered part of a suite of tests for solvers, as was discussed in Section 1.3. The most information, and most confidence, is gained if each of the stages listed in Section 1.3 is carried out (maybe excluding stage 3, the need for which is somewhat reduced if stage 4 is carried out). 3. CALCULATION VERIFICATION We now consider calculation verification in the context of cardiac electro-physiological simulations. In calculation verification, the goal is to estimate the errors in simulations used for the application of interest, with a final choice of model and numerical parameters, in particular at a specific mesh resolution. (We focus on mesh resolution rather than timestep because it is well-known that, for typical values used in studies, simulation results in cardiac EP are far more sensitive to mesh resolution than timesteps, see for example [20]). It is highly recommended [6] that calculation verification be carried out explicitly for the QOI for the chosen application. The aim of this section is to investigate how difficult it will be to perform strong calculation verification in cardiac EP modelling. We do not need to choose an application for this; instead, we perform two case studies, one in 1D and one in 3D on a realistic geometry, analysing errors in chosen QOIs. We again use the cardiac EP solver in Chaste for the simulations. Therefore, any conclusions on magnitudes of errors are specific to Chaste, specific to its particular numerical methods and implementation and specific to the method used for computing the QOIs. This is a fundamental aspect of calculation verification in VVUQ; calculation verification should be repeated for any new code and any new QOI. Error estimates given in this section may be used to gain some insight into expected errors when low-order numerical schemes are used, but they should not be used as likely magnitudes of errors when strict calculation verification is required with a different solver.

Note: choosing a standard propagating wave as the solution of the monodomain equations (e.g. V .x, t/ D V.x t/ for some function V) is not possible as the zero Neumann boundary conditions need to be satisfied exactly at all times. §§ In a finite element formulation, (8) is automatically imposed and requires no implementation. ‡‡



538


There are a number of ways in which calculation verification can be performed. The simplest is a brute force approach; by computing the QOI on meshes of differing resolutions, the error at a particular resolution can be estimated from the convergence details for the QOI. Alternatively, one may employ more sophisticated methods for error estimation, including mesh adaptivity. Efforts at error estimation in cardiac EP include the following: a threshold-based method [40], a posteriori error analysis [41], a priori estimates [42], a curvature-based method [43], Richardson extrapolation based adaptivity [44] and a Hessian-based method [45]. We re-iterate that such techniques aid but do not necessarily accomplish calculation verification, because calculation verification should be done for the particular QOI in the application at hand. Advanced error-estimation methods are beyond the scope of this paper. We use the brute force approach, although for the 1D problem, we investigate the efficacy of Richardson extrapolation to improve the estimate and provide error bounds. 3.1. One-dimensional fibre For the first case study, we take a 1D fibre of length 1 cm, stimulate one side and compute the voltage profile over time using the monodomain formulation. We choose the conduction velocity (CV) as the QOI for which we want to estimate the error. We compute the CV by measuring the time of activation (defined as the time that V reaches zero, calculated by finding the first timepoint ti the voltage Vi becomes positive and using the interpolated value t D .Vi ti 1 Vi 1 ti /=.Vi Vi 1 /), at the nodes at x D 0.25 and x D 0.75. We use the Luo-Rudy 1991 [46] cell model and the following parameters: D 1.75 mS/cm [47], Cm D 1 F/cm2 and D 1400 cm1 . This problem is solved using Chaste, using a first-order forward Euler method for the ODEs. The mesh resolutions are h D 0.05, 0.025, 0.0125, : : : (cm), and the ODE and PDE timesteps are set to 0.01 .h=0.05/ (ms). Richardson extrapolation is then applied on the computed conduction velocities. Briefly, this assumes that spatial stepsize and timesteps are small enough that the results are in the asymptotic regime of convergence for the particular QOI q 2 R and that convergence is monotonic [48], in which case the QOI will satisfy qh qexact D C hp C h.o.t. where qh is the computed QOI using a particular mesh resolution h, qexact is the unknown exact value of the QOI, C and p are unknown, and ‘h.o.t.’ stands for higher-order terms, which are neglected. By applying this equation with three different choices of h, we can eliminate the three unknowns C , p and qexact . The computed qexact is the Richardson extrapolated QOI using these three meshes, which we denote qR . For more details on this method, see, for example, [9]. An error band, a region in which the true QOI is highly likely to lie in (although not guaranteed) can also be defined, via the so-called ‘grid convergence index’ (GCI). See [9, 14] for details. The GCI is such that the true solution is expected to lie in the band qh .1 ˙ GCI/. Table II provides Table II. Richardson extrapolation of the conduction velocity in a 1D monodomain fibre. h is the spatial stepsize in the mesh, C Vh is the computed conduction velocity, C VR is the Richardson extrapolated conduction velocity using this mesh and the two coarser meshes (not defined on first two meshes, and not shown on third mesh as non-real result (due to oscillatory C Vh on first three meshes)), pR is the Richardson extrapolated order of convergence and GCI is the ‘grid convergence index’ (see text). h (cm) 0.05 0.025 0.0125 0.00625 0.003125 0.0015625 0.0007813 0.0003906 0.0001953

C Vh (cm/ms)

C VR (cm/ms)

pR

GCI

0.073638 0.079808 0.077012 0.075230 0.074676 0.074505 0.074455 0.074437 0.074429

— — 62 R 0.072103 0.074425 0.074430 0.074434 0.074426 0.074423

— — 62 R 0.651 1.683 1.703 1.761 1.467 1.223

— — — 0.05196 0.00420 0.00127 0.00035 0.00017 0.00010




539

Table III. Richardson extrapolation using monodomain model problem M-1D, with quantity of interest (QOI) q D kV kL1 .L2 / (exact solution: q D 1). h is the spatial stepsize in the mesh, qh is the computed QOI, qR is the Richardson extrapolated QOI using this mesh and the two coarser meshes, pR is the Richardson extrapolated order of convergence and GCI is the ‘grid convergence index’ (see text). h 0.1 0.05 0.025 0.0125 0.00625 0.003125

qh

qR

pR

GCI

0.95326 0.98799 0.99698 0.99924 0.99981 0.99995

— — 1.00010848 1.00000765 1.00000049 1.00000003

— — 1.951 1.986 1.996 1.991

— — 0.00393 0.00096 0.00024 0.00006

the computed CV, C Vh , at various mesh resolutions, together with the corresponding Richardson extrapolated CV, C VR , using the given mesh and the previous two coarser meshes. Also provided are the corresponding computed values of p, denoted pR and the GCI. We have also taken the 1D monodomain model problem M-1D and computed using Richardson extrapolation the QOI q D kV kL1 .L2 / . The exact value of this QOI is 1. The results are given in Table III. Consider the model problem first (Table III), which clearly showcases the effectiveness of Richardson extrapolation. We see that pR quickly converges to 2, the expected true value of p, and that Richardson extrapolated QOIs, qR , improve the estimate by several orders of magnitude on each mesh. Overall, Richardson extrapolation works well on this problem with kV kL1 .L2 / as the QOI. However, for the physiological problem with CV as the QOI (Table II), the values of pR do not obviously settle down to any value for the meshes considered, although they may be slowly converging to 1. This implies that the asymptotic regime of convergence has not been reached—despite the extremely high resolution of the last few meshes—which means that the Richardson extrapolation results should be used with caution, especially the error bands. We may still use the convergence results to estimate the exact solution; to three significant figures, it appears to be 0.0744 cm/ms. We observe that the error in CV is about 7% on the 0.025 cm mesh and about 3.5% on the 0.0125 cm mesh. It appears that the error is only less than 1% when h 0.0031 cm and smaller. Richardson extrapolation only appears to improve the estimate of the CV on the h 0.0031 cm and smaller meshes. It is not quite true that (our estimate of) the exact solution always lies in qh .1˙GCI/. Hence, whilst the GCI is considered a reasonably robust error estimator in computational fluid dynamics and other fields (Oberkampf et al. state that ‘Recent studies have shown that the GCI method is fairly reliable, even for solutions that are not in the asymptotic-convergence region’ [7]), it appears more research is required to identify such estimators for cardiac EP. 3.2. Three-dimensional rabbit heart We now move to a 3D anatomically accurate computational mesh of a rabbit heart. Here, we first look at changes in the computed voltage as mesh resolution is altered and then consider the impact of the findings on potential QOIs. We use the ‘Oxford Rabbit Heart’ [2], a highly detailed MR-based rabbit mesh of 24.1 million elements and an average element edge-length of 0.0125 cm, which makes it one of the highest resolution cardiac meshes in the world. We also make use of two coarsened versions of this mesh, a medium-resolution mesh (average edge-length of 0.0251 cm), and a coarse-resolution mesh (average edge-length of 0.0482 cm). Our goal is to investigate the feasibility of rigorous calculation verification in whole-heart simulations. We suppose the application mesh for a simulation study is the fine, anatomically detailed, mesh. In calculation verification, the aim is to estimate or bound numerical errors on QOIs given this choice of mesh. Computing simulations on an ever finer mesh is computationally intractable; however, we may hope to use results from simulations on the coarser meshes to estimate numerical error. This convergence experiment is somewhat unusual, because as well as being lower resolution, the coarser meshes involve a loss of geometrical Copyright © 2013 John Wiley & Sons, Ltd.



Color Online, B&W in Print

540

Figure 4. Voltage (dark, repolarised; bright, depolarised) at various times after an S1–S2 protocol is used to induce reentry on a highly anatomically detailed realistic rabbit mesh (labelled ‘Fine’, average edgelength: 0.0125cm), and two coarsened versions of this mesh (‘Medium’, average edge-length: 0.0251 cm; and ‘Coarse’, average edge-length: 0.0482cm). Copyright © 2013 John Wiley & Sons, Ltd.



541

detail (trabeculae and papillary muscles, see [2]) compared to the finest mesh. A more conventional approach would be to begin with the coarsest mesh and refine twice to generate three meshes of differing resolution covering exactly the same geometry. However, such as experiment would not involve the application mesh, so would not be appropriate for rigorous calculation verification. An S1–S2 protocol, composed of an initial apical stimulus followed by a second stimulus applied to a large left-ventricular region 170 ms later, is used to induce reentry. The monodomain formulation is used, together with the Mahajan et al. rabbit cell model [49] for membrane kinetics, parameters Cm D 1 F/cm2 , D 1400cm1 and an isotropic conductivity of 0.93P mS/cm. The conductivity is chosen by taking the monodomain average conductivity of an intracellular conductivity of 1.7 mS/cm and extracellular conductivity of 6.2 mS/cm [47] and scaling this conductivity by 70%, chosen arbitrarily so that reentry is sustained for more than 500 ms. The following conclusions are therefore specific to this choice of parameters and to the solver used, Chaste. The code for running this simulation can be viewed at, or downloaded from, https://chaste.cs.ox.ac.uk/trac/wiki/ PaperTutorials. Figure 4 shows several snapshots of the epicardial voltage at various times following the S2 stimulus, on the three meshes. The results are interesting, surprising and somewhat difficult to interpret. First, consider the differences between the medium and coarse meshes. Despite the coarse mesh being very coarse for cardiac EP simulations, and it being known that coarse meshes can lead to inaccurate conduction velocities (e.g. see Table II), the voltage maps for these two meshes are surprisingly similar, agreeing reasonably closely until about 400 ms. After that, numerical differences lead to divergence of these solutions, and activity on the coarse mesh in fact ends 130 ms before the medium mesh. The differences between the fine and medium solutions are actually greater than between the medium and coarse solutions. This can be attributed to loss of geometrical detail from fine to medium mesh. We observe that this has a significant effect on simulation results (in agreement with [2], which compares simulation results on the fine mesh and a mesh of about the same resolution but less geometric detail). These results highlight the inherent difficulties in trying to perform rigorous calculation verification on realistic cardiac meshes. If the application mesh is the fine anatomically-detailed mesh, and the QOI is for example the time it takes for electrical activity to self-terminate under this protocol, it is clear that there is essentially no error band that can be estimated from these results. Performing another convergence study with the coarse mesh and two refined meshes (i.e. the more conventional approach mentioned previously) may provide additional information, but might still be unlikely to lead to confident estimates of numerical error on the application mesh. Often, the QOI will be the difference in a result between this base simulation and a second simulation using a potentially pro-arrhythmic drug. Obtaining an error bound for this case appears to be a challenging problem indeed. 4. DISCUSSION AND RECOMMENDATIONS There has been a great deal of research into developing methods and best practices in the field of VVUQ, which may be used to evaluate the credibility of predictions by computational models. VVUQ has been used extensively in the engineering and physical sciences, but there has been very little overlap between this field and cardiac electro-physiological modelling, due to the complexity of biological models and because the consequences of incorrect predictions by cardiac EP models have, up until now, not been high. Stronger VVUQ of cardiac models will be required before they can be used in safety-critical clinical decision-making. The rigour required in the VVUQ activities will depend on the application and the impact of an incorrect prediction. A framework for deciding how rigorous the VVUQ should be has yet to be created, although one is currently in development by the U.S. Food and Drug Administration for modelling used for medical devices [50]. This paper has been concerned with verification, which can be split into two stages, code verification and calculation verification. For safety-critical simulations, it is important that both types of verification are performed and that it is well-documented. Code verification amounts to building a body of evidence that suggests that the software implementation is correct. Tools for code verification complement rather than supplant each other, and it is important that a suite of benchmark Copyright © 2013 John Wiley & Sons, Ltd.


542


problems be available to test solvers. Prior to the present article, only one benchmark problem, the Niederer monodomain benchmark [20], had been defined for cardiac EP. As far as we are aware, there were no benchmark problems with analytic solutions that are fully coupled and 2D or 3D. In this paper, we have defined 1D–3D monodomain and bidomain problems, and a bidomain-withbath problem (set in 1D–3D, although each essentially 1D), all with known exact solutions. These problems have been carefully constructed so that they can be easily run by a general cardiac solver. Therefore, they could be used to vastly increase confidence that a cardiac EP solver is implemented correctly, as well as confirm that convergence occurs at the expected rate. Chaste [30, 31] has been shown to correctly solve the model problems. Further problems with analytic solutions should be developed in the future to continue to improve confidence in simulation results. These could include monodomain or bidomain problems on other domains or non-constant/non-diagonal conductivities (this could be accomplished by utilising other solutions of (13)), or 2D or 3D bidomain-withbath problems. Following from [20], further physiological benchmark problems for comparisons of solvers, for example, using the bidomain formulation or complex geometries, would also help improve confidence in cardiac solvers. We have also considered calculation verification in this paper. It is important that calculation verification be performed explicitly in the QOI to be used in the application, and with the solver and numerical methods to be used in the application problem. Using one particular solver, we have considered whether it is possible to determine error bounds on QOIs in a 1D problem and a 3D problem on a realistic rabbit geometry. The results suggest that performing reliable calculation verification in cardiac EP applications is a very challenging problem. We finish by making some recommendations, which simply summarise the aforementioned discussion. We believe that 1. Computational cardiac modellers should familiarise themselves with VVUQ. 2. Monodomain, bidomain or bidomain-with-bath software should be run on the model problems presented here (or extensions of these), to identify errors. 3. For verification, further benchmark problems, both with and without exact solutions, should be developed. 4. When running simulations, a researcher should consider how much VVUQ will likely be required, given the application and potential impact of incorrect predictions. 5. If strong evidence of verification is required a. Hand-coding of complex cell models should be avoided. b. Software should be run on the Niederer benchmark [20] and compared to other solvers. c. All code verification steps should be documented in detail. d. For the specific application problem of interest, error bounds in the computed quantities of interest should be estimated, as much as possible.

DISCLAIMER The mention of commercial products, their sources or their use in connection with material reported herein is not to be construed as either an actual or implied endorsement of such products by the Department of Health and Human Services.

ACKNOWLEDGEMENTS

Pras Pathmanathan would like to gratefully acknowledge the input from Chris Arthurs and David Gavaghan (Computational Biology Group, Oxford University) and support from the UK Engineering and Physical Sciences Research Council (EPSRC)-funded ‘2020 Science’ programme (grant number: EP/I017909/1), http://www.2020science.net. The authors would also like to thank Leonardo Angelone and Tina M. Morrison (Center for Devices and Radiological Health (CDRH), Food and Drug Administration (FDA)) and two anonymous reviewers for their feedback. The rabbit heart simulations were carried out using FDA super-computing resources (supported by the Office of Science and Engineering Laboratories in CDRH). Copyright © 2013 John Wiley & Sons, Ltd.



543

REFERENCES 1. Noble D. A modification of the Hodgkin-Huxley equations applicable to Purkinje fibre action and pacemaker potentials. Journal of Physiology 1962; 160:317–352. 2. Bishop M, Plank G, Burton R, Schneider J, Gavaghan D, Grau V, Kohl P. Development of an anatomically detailed MRI-derived rabbit ventricular model and assessment of its impact on simulations of electrophysiological function. American Journal of Physiology - Heart and Circulatory Physiology 2010; 298(2):H699–H718. 3. Trayanova N, O’Hara T, Bayer J, Boyle P, McDowell K, Constantino J, Arevalo H, Hu Y, Vadakkumpadan F. Computational cardiology: how computer simulations could be used to develop new therapies and advance existing ones. Europace 2012; 14:v82–v89. 4. Niederer S, Shetty A, Plank G, Bostock J, Razavi R, Smith N, Rinaldi C. Iophysical modeling to simulate the response to multisite left ventricular stimulation using a quadripolar pacing lead. Pacing and Clinical Electrophysiology 2012; 35(2):204–214. 5. Trayanova N, Plank G, Rodríguez B. What have we learned from mathematical models of defibrillation and postshock arrhythmogenesis? Application of bidomain simulations. Heart rhythm 2006; 3(10):1232–1235. 6. National Research Council. Assessing the Reliability of Complex Models: Mathematical and Statistical Foundations of Verification, Validation, and Uncertainty Quantification. The National Academies Press: Washington, D.C., 2012. 7. Oberkampf W, Trucano T, Hirsch C. Verification, validation, and predictive capability in computational engineering and physics. Applied Mechanics Reviews 2004; 57(5):345–384. 8. American Society of Mechanical Engineers. ASME V&V 10-2006: Guide for Verification and Validation in Computational Solid Mechanics. American Society of Mechanical Engineers: New York City, 2006. 9. American Society of Mechanical Engineers. ASME V&V 10.1-2012: An Illustration of the Concepts of Verification and Validation in Computational Solid Mechanics. American Society of Mechanical Engineers: New York City, 2012. 10. American Society of Mechanical Engineers. ASME V&V 20-2009: Standard for Verification and Validation in Computational Fluid Dynamics and Heat Transfer. American Society of Mechanical Engineers: New York City, 2009. 11. National Aeronautics and Space Administration. NASA-STD-7009: Standard for models and simulation. NASA: Washington DC, 2009. 12. International Atomic Energy Agency. Verification and Validation of Software Related to Nuclear Power Plant Instrumentation and Control, Technical reports series. International Atomic Energy Agency: Vienna, Austria, 1999. 13. Harvego E, Schultz R, Crane R. Development of a standard for verification and validation of software used to calculate nuclear system thermal fluids behavior. American Society of Mechanical Engineers: New York City, 2010. 14. Roache P. Verification and Validation in Computational Science and Engineering. Hermosa: Albuquerque, NM, 1998. 15. Oberkampf W, Roy C. Verification and Validation in Scientific Computing. Cambridge University Press: Cambridge, UK, 2010. 16. Post D, Votta L. Computational science demands a new paradigm. Physics Today 2005; 58(1):35–41. 17. Frenk C, White S, Bode P, Bond J, Bryan G, Cen R, Couchman H, Evrard A, Gnedin N, Jenkins A, Khokhlov A, Klypin A, Navarro J, Norman M, Ostriker J, Owen J, Pearce F, Pen U, Steinmetz M, Thomas P, Villumsen J, Wadsley J, Warren M, Xu G, Yepes G. The Santa Barbara cluster comparison project: a comparison of cosmological hydrodynamics solutions. The Astrophysical Journal 2009; 525(2):554–582. 18. Andren A, Brown A, Mason P, Graf J, Schumann U, Moeng C-H, Nieuwstadt F. Large-eddy simulation of a neutrally stratified boundary layer: a comparison of four computer codes. Quarterly Journal of the Royal Meteorological Society 2006; 120(520):1457–1484. 19. Hatton L, Roberts A. How accurate is scientific software? IEEE Transactions on Software Engineering 1994; 20(10):785–797. 20. Niederer S, Kerfoot A, Benson E, Bernabeu M, Bernus O, Bradley C, Cherry E, Clayton R, Fenton F, Garny A, Heidenreich E, Land S, Maleckar M, Pathmanathan P, Plank G, Rodríguez R, Sachse F, Seemann G, Skavhaug O, Smith N. Verification of cardiac tissue electrophysiology simulators using an N-version benchmark. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 2011; 369(1954):4331–435. 21. Ainsworth M, Oden J. A posteriori error estimation in finite element analysis. Computer Methods in Applied Mechanics and Engineering 1997; 142(1):1–88. 22. Bangerth W, Rannacher R. Adaptive Finite Element Methods for Differential Equations. Birkhäuser: Basel, 2003. 23. Noble D, Garny A, Noble P. How the Hodgkin–Huxley equations inspired the cardiac Physiome Project. J The Journal of Physiology 2012; 590(11):2613–2628. 24. Rush S, Larsen H. A practical algorithm for solving dynamic membrane equations. Biomedical Engineering, IEEE Transactions on 1978; 4:389–392. 25. Miller A, Marsh J, Reeve A, Garny A, Britten R, Halstead M, Cooper J, Nickerson D, Nielsen P. An overview of the CellML API and its implementation. BMC Bioinformatics 2010; 11(1):178–189. 26. Garny A, Nickerson D, Cooper J, Weber dos Santos R, McKeever S, Nielsen P, Hunter P. CellML and associated tools and techniques. Philosophical Transactions of the Royal Society A 2008; 366(1878):3017–3043. 27. Lloyd C, Lawson J, Hunter P, Nielsen P. The CellML model repository. Bioinformatics 2008; 24(18): 2122–2123. 28. Cooper J. Automatic validation and optimisation of biological models. PhD Thesis, University of Oxford, 2009. Copyright © 2013 John Wiley & Sons, Ltd.


544


29. Vigmond E, Hughes M, Plank G, Leon L. Computational tools for modeling electrical activity in cardiac tissue. Journal of electrocardiology 2003; 36:69–74. 30. Pitt-Francis J, Pathmanathan P, Bernabeu M, Bordas R, Cooper J, Fletcher A, Mirams G, Murray P, Osbourne J, Walter A, Chapman S, Garny A, van Leeuwen I, Maini P, Rodriguez B, Waters S, Whiteley J, Byrne H, Gavaghan D. Chaste: a test-driven approach to software development for biological modelling. Computer Physics Communications 2009; 180:2452–2471. 31. Mirams G, Arthurs C, Bernabeu M, Bordas R, Cooper J, Corrias A, Davit Y, Dunn S-J, Fletcher A, Harvey D, Marsh M, Osborne J, Pathmanathan P, Pitt-Francis J, Southern J, Zemzemi N, Gavaghan D. Chaste: an open source C++ library for computational physiology and biology. PLoS Computational Biology 2013; 9(3):e1002790. 32. Bradley C, Bowery A, Britten R, Budelmann V, Camara O, Christie R, Cookson A, Frangi A, Gamage T, Heidlauf T. et al. OpenCMISS: A multi-physics & multi-scale computational infrastructure for the VPH/Physiome project. Progress in Biophysics and Molecular Biology 2011; 107(1):32–47. 33. Mardal K-A, Skavhaug O, Lines G, Staff C, Odegard A. Using Python to solve partial differential equations. Computing in Science & Engineering 2007; 9(3):48–51. 34. Pathmanathan P, Bernabeu M, Niederer S, Gavaghan D, Kay D. Computational modelling of cardiac electrophysiology: explanation of the variability of results from different numerical solvers. International Journal for Numerical Methods in Biomedical Engineering 2011; 28(8):890–903. 35. Keener J, Sneyd J. Mathematical Physiology volume 8 of Interdisciplinary Applied Mathematics. Springer: New York City, 1998. 36. Chernyak Y. Steady state plane wave propagation speed in excitable media. Physical Review E 1997; 56(2):2061–2073. 37. Hakim V, Karma A. Spiral wave meander in excitable media: the large core limit. Physical Review Letters 1997; 79(4):665–668. 38. Thomée V. Galerkin Finite Element Methods for Parabolic Problems, volume 25 of Springer Series in Computational Mathematics. Springer–Verlag: Berlin–Heidelberg–New York, 1997. 39. Pathmanathan P, Bernabeu M, Bordas R, Cooper J, Garny A, Pitt-Francis J, Whiteley J, Gavaghan D. A numerical guide to the solution of the bidomain equations of cardiac electrophysiology. Progress in Biophysics and Molecular Biology 2010; 102:136–155. 40. Whiteley J. Physiology driven adaptivity for the numerical solution of the bidomain equations. Annals of Biomedical Engineering 2007; 35(9):1510–1520. 41. Franzone P, Deuflhard P, Erdmann B, Lang J, Pavarino L. Adaptivity in space and time for reaction-diffusion systems in electrocardiology. SIAM Journal on Scientific Computing 2006; 28(3):942–962. 42. Arthurs C, Bishop M, Kay D. Efficient simulation of cardiac electrical propagation using high order finite elements. Journal of Computational Physics 2012; 231(10):3946–3962. 43. Southern J, Gorman G, Piggott M, Farrell P, Bernabeu M, Pitt-Francis J. Simulating cardiac electrophysiology using anisotropic mesh adaptivity. Journal of Computational Science 2010; 1(2):82–88. 44. Cherry E, Greenside H, Henriquez C. A space-time adaptive method for simulating complex cardiac dynamics. Physical Review Letters 2000; 84(6):1343–1346. 45. Belhamadia Y, Fortin A, Bourgault Y. Towards accurate numerical method for monodomain models using a realistic heart geometry. Mathematical Biosciences 2009; 220(2):89–101. 46. Luo C, Rudy Y. A model of the ventricular cardiac action potential: depolarization, repolarization, and their interaction. Circulation Research 1991; 68:1501–1526. 47. Clerc L. Directional differences of impulse spread in trabecular muscle from mammalian heart. The Journal of Physiology 1976; 255(2):335–346. 48. Li J, Hu G, Shaffer C. Limitations of Richardson extrapolation and some possible remedies. Journal of Fluids Engineering 2005; 127:795–805. 49. Mahajan A, Shiferaw Y, Sato D, Baher A, Olcese R, Xie L, Yang M, Chen P, Restrepo J, Karma A, Garfinkel A, Qu Z, Weiss J. A rabbit ventricular action potential model replicating cardiac dynamics at rapid heart rates. Biophysical Journal 2008; 94(2):392–410. 50. Sainani K. Getting it right: Better validation key to progress in biomedical computing. Biomedical Computation Review Fall 2012; Fall 2012:9–17.

SUPPORTING INFORMATION Additional supporting information may be found in the online version of this article at the publisher’s web site.



Cell-specific cardiac electrophysiology models.

Computational cardiac electrophysiology is moving towards translation medicine.

Computational cardiac electrophysiology is ready for prime time.

A brief history of tissue models for cardiac electrophysiology.

Computational modelling of cardiac electrophysiology: explanation of the variability of results from different numerical solvers.

A quantitative comparison of the behavior of human ventricular cardiac electrophysiology models in tissue.

Uncertainty quantification of fast sodium current steady-state inactivation for multi-scale models of cardiac electrophysiology.

Estimability Analysis and Optimal Design in Dynamic Multi-scale Models of Cardiac Electrophysiology.

Regional segmentation of ventricular models to achieve repolarization dispersion in cardiac electrophysiology modeling.

The Cardiac Electrophysiology Web Lab.

The Impact of Opioids on Cardiac Electrophysiology.

Computational models of planning.

Rabbit-specific computational modelling of ventricular cell electrophysiology: Using populations of models to explore variability in the response to ischemia.

Robotics in invasive cardiac electrophysiology.

A second-generation computational modeling of cardiac electrophysiology: response of action potential to ionic concentration changes and metabolic inhibition.

Human atrial cell models to analyse haemodialysis-related effects on cardiac electrophysiology: work in progress.

Computational models of heart disease.

Automatic Parameterization Strategy for Cardiac Electrophysiology Simulations.

Ensuring reliability of safety-critical clinical applications of computational cardiac models.

National registry on cardiac electrophysiology (2012).

A comparison of solver performance for complex gastric electrophysiology models.

Blood flow in intracranial aneurysms treated with Pipeline embolization devices: computational simulation and verification with Doppler ultrasonography on phantom models.

Variability in cardiac electrophysiology: Using experimentally-calibrated populations of models to move beyond the single virtual physiological human paradigm.

Electrophysiology of human cardiac atrial and ventricular telocytes.