Composite Backward Differentiation Formula for the Bidomain Equations

The bidomain equations have been widely used to model the electrical activity of cardiac tissue. While it is well-known that implicit methods have much better stability than explicit methods, implicit methods usually require the solution of a very large nonlinear system of equations at each timestep which is computationally prohibitive. In this work, we present two fully implicit time integration methods for the bidomain equations: the backward Euler method and a second-order one-step two-stage composite backward differentiation formula (CBDF2) which is an L-stable time integration method. Using the backward Euler method as fundamental building blocks, the CBDF2 scheme is easily implementable. After solving the nonlinear system resulting from application of the above two fully implicit schemes by a nonlinear elimination method, the obtained nonlinear global system has a much smaller size, whose Jacobian is symmetric and possibly positive definite. Thus, the residual equation of the approximate Newton approach for the global system can be efficiently solved by standard optimal solvers. As an alternative, we point out that the above two implicit methods combined with operator splittings can also efficiently solve the bidomain equations. Numerical results show that the CBDF2 scheme is an efficient time integration method while achieving high stability and accuracy.


INTRODUCTION
The monodomain equations or bidomain equations, consisting of a coupled system of partial differential equations (PDEs) and ordinary differential equations (ODEs), are often used to mathematically model the electrical activity of the heart. The PDEs describe the propagation of the electrical signal through the cardiac tissue and the ODEs describe electrochemical reactions in the cells which are usually nonlinear. Modern myocyte models such as the Ohara-Rudy model and anatomically realistic spatial models enable us to achieve a quantitative understanding of the relationship between molecular function and the integrated behavior of the cardiac myocyte in health and disease, e.g., to predict the clinical risk of drug-induced arrhythmias. It is notable that the cardiac electrical activity typically involves multiple and widely varying scales, which makes the monodomain equations and bidomain equations stiff. As a result, it is computationally expensive to get an accurate solution of the electrical activity of the heart (Tung, 1978;Keener and Sneyd, 1998;Sundnes et al., 2007).
In detailed numerical simulations, the finite element method (Sundnes et al., 2005), finite difference method (Pollard et al., 1992), and finite volume method (Johnston, 2010) have all been used for spatial discretization of the domain. While for temporal discretization, the easiest way is to use explicit methods such as the Runge Kutta method and the forward Euler method since they do not require matrix inversion. However, the timesteps allowed in these methods are generally very small in order to satisfy the requirement of stability (dos Santos et al., 2004;Potse et al., 2006). To overcome the drawbacks of explicit methods, the semi-implicit and Crank-Nicolson methods for solving the PDE component of the bidomain equations are proposed (Keener and Bogar, 1998). But the nonlinear term of the bidomain equations still has restriction on the timesteps for stability (Keener and Bogar, 1998). Although fully implicit methods (Hooke, 1992;Pollard et al., 1992;Dal et al., 2012) preferred for stiff systems can guarantee the stability of the numerical schemes and can use large timesteps, they often have the significant drawback of requiring the solution of a large-scale nonlinear system for diffusion and reaction parts of the bidomain equations at each timestep. For example, a fully implicit parallel Newton-Krylov-Schwarz method for solving the bidomain equations has been proposed in Murillo and Cai (2004). They considered a simple case that the conductivity only varies with the co-ordinate system, i.e., the conductivity tensor is a diagonal matrix in the twodimensional case. The methods proposed in Murillo and Cai (2004) are difficult to implement for more complex geometries and membrane models.
Operator splitting (Qu and Garfinkel, 1999;Trangenstein and Kim, 2004;Sundnes et al., 2005) is a popular technique that can avoid solving the large nonlinear system of the bidomain equations by splitting the equations into more manageable parts, e.g., splitting the bidomain equations into the linear diffusion and nonlinear reaction parts. An operator splitting method combined with a semi-implicit approximation or a Crank-Nicolson approximation to update the transmembrane and extracellular potentials of the bidomain equations has been proposed in Whiteley (2006) which allows for much larger timesteps.
In this work, we present a backward Euler (BE) and a twostage composite backward differentiation formula (CBDF2) for solving the cardiac electrical dynamics models. The finite element method is applied for the space discretization. The CBDF2 scheme combines the second-order backward differentiation formula and the BE method with the latter used as fundamental building blocks thus it is fully implicit and easily implementable. We point out that the CBDF2 scheme is unconditionally stable and stiffly accurate for stiff problems, i.e., L-stable. After applying the BE method for the bidomain equations, there are two coupled nonlinear systems: space dependent PDEs and space independent ODEs. The nonlinear elimination method (Ying et al., 2008) is used in this work to solve the space dependent PDEs arising from the fully implicit methods for the bidomain equations. The obtained nonlinear system after elimination has a much smaller size and has a symmetric and positive definite Jacobian when the timestep is small. We solve the nonlinear system by an approximate Newton approach where the residual equations in each Newton iteration can be efficiently solved by standard optimal solvers. As an alternative, we also combine the BE and CBDF2 schemes with Godunov splitting and Strang splitting, respectively, to efficiently solve the bidomain equations. From the numerical results, the CBDF2 scheme allows for timesteps four times larger than that of the BE method while achieving the same high numerical accuracy.

The Model
The bidomain model describing the electrical activity of the heart or a cardiac tissue can be formulated as a system of nonlinear ordinary and partial differential equations (Keener and Sneyd, 1998). Let ⊂ R d (d > 0) be the bounded computational domain; t and x ∈ be the temporal and spatial independent variables, respectively; Φ i = Φ i (t, x) and Φ e = Φ e (t, x) be the intracellular and extracellular potentials, respectively; x) be the transmembrane potential; q represents a set of state variables which can include both gating and ion concentration variables. We consider the bidomain equations of the following form for t > 0 and x ∈ , subject to the homogeneous Neumann boundary condition n · [D i ∇Φ i ] = 0 and n · [D e ∇Φ e ] = 0 on ∂ and one constraint on the intra-and extra-cellular potentials Here, C m is the membrane capacitance per unit area; β is the surface to volume ratio of the cardiac cells; D i and D e are the space dependent intracellular and extracellular conductivity tensors, respectively; n is the unit outward normal on ∂ ; I ion (V m , q) and M(V m , q) are two known functions, which are typically nonlinear and describe the electrical dynamics of a cardiac myocyte; I stim = I stim (t, x) is an extracellular stimulus current. Provided some appropriate initial conditions on the intracellular potential Φ i , extracellular potential Φ e , and the vector q of state variables, the bidomain equations can be uniquely solved. We remark that the bidomain equations for the electrical activity of the heart may take different but equivalent forms (Sundnes et al., 2005;Whiteley, 2006;Vigmond et al., 2008).

The Backward Euler Method
We first introduce the BE method which is used as the fundamental building block of the CBDF2 scheme. Set t > 0 as the timestep. Time integration of the bidomain equations by the BE method from time t n to time t n+1 leads to the semi-discrete bidomain equations where Φ n i , Φ n e , V n m , and q n are the finite difference approximations of the variables Φ i , Φ e , V m , and q at time t n , respectively.
To eliminate the state variables q n+1 in Equations (5) and (6), we express the ionic current I ion (V n+1 m , q n+1 ) as a function of V n+1 m denoted by i.e., the vector of state variable q n+1 in Equation (7) can be regarded as a vector value function of the transmembrane potential V n+1 m . In our implementation in each timestep, we first give V n+1 m an initial guess, e.g., extrapolated from V n m and V n−1 m , and numerically solve the nonlinear system (7) to obtain the corresponding state variables q n+1 which is substituted into the ionic current function I ion (V n+1 m , q n+1 ) to estimate I ion (V n+1 m ). The value of I ion (V n+1 m ) is used in the Newton-type iteration formula for solving Equations (5) and (6) to update V n+1 m . The above steps are performed iteratively until the residual is less than the tolerance.
We rewrite the Equations (5) and (6) as or where µ = C m / t. Assuming that the domain can be partitioned into a quasiuniform grid if d = 2 or a quasi-uniform mesh if d = 3, the semidiscrete bidomain Equations (10) can be further discretized by using the continuous piecewise linear finite element method on the grid. Let M be the finite element mass matrix; K i and K e be the finite element stiffness (non-negative) matrices corresponding to the intracellular and extracellular diffusion operators (−β −1 ∇ · D i ∇) and (−β −1 ∇·D e ∇), respectively. Let Φ n i , Φ n e , V n m , I ion (V n m ), and I stim be the column vectors whose entries are the discrete values of the potential and current variables at the corresponding grid nodes. The dimensions of the matrices and vectors are the same and both are equal to the number of nodes in the triangular grid. Using the matrix-vector notation, the discrete finite element bidomain equations now read as which is nonlinear and can be solved with an approximate Newton method. We denote J ion (V m ) as the derivative of the current function I ion (V m ) with respect to V m , i.e., be the diagonal matrix with the diagonal entries being the values of J ion (V n+1 m ). Then we have the Jacobian matrix of the nonlinear system (11) Note that the evaluation of the current function I ion (V n+1 m ) involves the solution of a local nonlinear system for the vector of state variables q n+1 . Since the derivative function J ion (V n+1 m ) of I ion (V n+1 m ) may not be directly available, for the convenience of implementation, it is computed approximately by the centered difference where δ is a small positive parameter, typically on the order of 10 −6 .
are the piecewise linear basic functions for the finite element method and N is the total number of nodes in the quasi-uniform grid or mesh of . Since the conductivity tensors D i and D e are symmetric and positive definite, for a non-zero vector c, the quadratic forms of the stiffness matrices are non-negative and vanish if and only if ∇y ≡ 0 for any x ∈ , which implies that c is parallel to (1, 1, . . . , 1) T . Therefore, the stiffness matrices K i and K e are symmetric positive semi-definite and have a rank of N − 1. We then discuss the property of the Jacobian matrix A in Equation (12). Since µ = C m / t, the first part of matrix A in Equation (12) will dominate as t approaches zero which is positive semi-definite and has a rank of 2N − 1. The second part of the matrix A is symmetric and negligible compared with the first part when t approaches zero. Therefore, with a rank of 2N − 1, the Jacobian matrix A is irreversible, which is consistent with the fact that the solution of the bidomain equations is not unique if the constraint (4) is absent.
To obtain a unique solution, we take the following discrete formula of the constraint (4) as where W = (1, 1, . . . , 1) T / √ N has a dimension of N. Then in each approximate Newton iteration, we use the following matrix to solve the bidomain equations. Note that the column vector W 1 = (1, 1, . . . , 1) T is an eigenvector corresponding to the eigenvalue 0 of matrix A. Denote the other eigenvectors corresponding to the non-zero eigenvalues of matrix A by W 2 , W 3 , . . . , W 2N . Then we have Therefore, the modified Jacobian matrix B is symmetric and positive definite and the system has a unique solution after adding condition (14). In principle, we can solve the stabilized system which has a coefficient matrix B by a direct method such as the Gauss elimination method or an iterative method such as the Gauss-Seidel iteration and the (preconditioned) conjugate gradient iteration. Because the bidomain equations often need to be discretized on a fine grid with a large number of nodes, the computational work involved with the solution of the linear equations by the above methods may be large and not optimal. To achieve optimal performance, we can apply the V-cycle geometric multigrid method or multigrid preconditioned conjugate gradient method (Wesseling, 1995;Saad, 2003;Ying, 2005) to solve the stabilized linear system in each Newton iteration.
To numerically solve the bidomain equations by the BE method, we give two nonlinear vector-valued functions from Equations (7) and (11) as and . (17) The system and from time t n to t n+1 can be solved by the following detailed Algorithm 1.
e ) ν and update the state variables q n+1 of each space point by solving the local nonlinear system (18) using V n+1 m .

5
Compute I ion (V n+1 m ) by substituting q n+1 into the ionic current I ion (V n+1 m , q n+1 ) and compute J ion (V n+1 m ) by Equation (13). 6 Solve the linear system by a standard V-cycle multigrid method with Gauss-Seidel smoother such that Update the iteration number ν ← ν + 1.
It is worth noting that the updated membrane potential e ) has a corresponding updated q n+1 , which is computed by solving the nonlinear system of Equation (18), e.g., by the Newton iteration or the damped Newton iteration method.

The Composite Backward Differentiation Formula
We then introduce the CBDF2 scheme (Ying et al., 2008(Ying et al., , 2009) to solve the bidomain equations. The CBDF2 scheme has secondorder accuracy and is L-stable. It uses the BE method as a fundamental building block. As we have presented how to discretize and solve the bidomain equations with the BE method, the CBDF2 scheme can be readily and easily implemented.
For ease of illustration, we first give the CBDF2 scheme for the abstract form which can then be straightforwardly extended to the bidomain equations. The CBDF2 scheme for Equation (20) is given by (22) where the parameter γ takes the value of 1 − √ 2/2 to obtain a second-order numerical accuracy (Ying et al., 2009). Here, u n is an approximation of u(t n ) at time t n = n t and u n+γ is an approximation of u(t n+γ ) at the intermediate time t n+γ = t n + γ t. The right hand side of Equation (22) is an approximation of u(t n+1−γ ) at time t n+1−γ = t n+1 − γ t which is extrapolated from the solutions at t n and t n+γ . Equation (21) is called the first stage and can be regarded as the BE discretization from time t n to time t n+γ , while Equation (22) is called the second stage and can be regarded as the BE discretization from time t n+1−γ to time t n+1 . The form of the two stages of the CBDF2 scheme are the same as that in the BE temporal discretization except that the timestep is replaced by γ t. Thus, u n+γ and u n+1 in Equations (21) and (22) can be computed iteratively as previously discussed in Algorithm 1, i.e., the BE method is the building block of the CBDF2 scheme. To extend the CBDF2 scheme to the bidomain equations, we only need to rewrite the system (1)-(3) to a vector form of Equation (20). The detailed CBDF2 scheme for the bidomain equations is given in the following Algorithm 2.

Stability Analysis of The CBDF2 Scheme
We first analyze the stability function of the CBDF2 scheme of a simple ODE which is The stability function S(z) has been proven to be bounded by one for all numbers on the left half complex plane and converge to zero as the real part of the complex number z tends to negative infinity with a particular value of the characteristic constant γ = 1 − √ 2/2 (Ying et al., 2008), i.e., the CBDF2 scheme is L-stable.
Algorithm 2: the CBDF2 method 1 Denote the discrete initial condition by V 0 m and q 0 . 2 Denote the simulating time as T and let K = T/ t. 3 for i = 0 to K do 4 Get the concrete form of the nonlinear system by substituting V i m and q i into the Equations (18) and (19).

5
Solve the system f(q i+γ ) = 0 and i+γ e ) T = 0 from time t i to t i+γ by the Newton iteration method as shown in Algorithm 1 but using timestep of γ t.

6
Compute the intermediate variable Get the concrete form of the nonlinear system by substituting the variable V i+1−γ m and q i+1−γ into the equations (18) and (19).

8
Solve the system f(q i+1 ) = 0 and F (Φ i+1 i , Φ i+1 e ) T = 0 from time t i+1−γ to t i+1 by the Newton iteration method. 9 end 10 Accept the vector (Φ K i , Φ K e ) T as an approximation of (Φ i , Φ e ) T at the simulating time T.
We then derive the global error of the CBDF2 scheme for Equation (20). Assume the function f (u) in Equation (20) is Lipschitz continuous in u over some feasible domain D, i.e., there exists some constant L ≥ 0 such that |f (u)−f (u * )| ≤ L|u−u * | for all u and u * in D. The local truncation error (LTE) of the CBDF2 scheme (Ying, 2005) has been proven to be which is on the order of ( t) 3 . From the form of the second stage of the CBDF2 scheme (22), we have where u(t n+γ ) satisfies Let Equation (25) subtract Equation (21) and we can obtain (1 − Lγ t)| e n+γ | ≤ |e n | where e n = u(t n ) − u n and e n+γ = u(t n+γ ) − u n+γ . By subtracting Equation (24) from Equation (22), we have the following inequality i.e., ≤ 1 + 2L t and we can further obtain the global error of the CBDF2 scheme at simulation time T as Similarly, the LTE of the BE method for the abstract form (20) is t 2 f ′ /2, while its global error satisfies and has the form at simulation time T as when L t ≤ 0.5. Since the restriction for the CBDF2 scheme is L t ≤ 1.12, the CBDF2 scheme allows a larger timestep than the BE method for the abstract form (20).

NUMERICAL RESULTS
In this section, we compare the numerical results from the BE and CBDF2 schemes with different spatial step sizes and timestep sizes. The computational domain is partitioned into a quasiuniform triangular grid in the two-dimensional space and a hexahedral mesh in the three-dimensional space, respectively. A linear approximation is used to obtain the solution on each element. The detailed refinement of a coarser grid to get a finer grid and the numbering principle of the triangular grid are shown in Supplementary Figure 1 The transmembrane polarization is driven by an extracellular stimulus which is modeled as a virtual battery with an anode and a cathode in the center of the domain (Vigmond et al., 2008) and is given by where V m stim is a constant, r is the distance between x and the anode and s is the distance between x and the cathode.
The operator splitting technique (Qu and Garfinkel, 1999;Trangenstein and Kim, 2004;Sundnes et al., 2005) is also used to efficiently solve the bidomain equations where the linear diffusion part and the nonlinear reaction part are solved separately. To be specific, the Godunov splitting combined with the BE method and the Strang splitting combined with the CBDF2 scheme are used to solve the bidomain equations (see Supplementary Material for details). The linear system obtained from discretizing the diffusion part is solved by the standard multigrid V-cycle with Gauss-Seidel smoother. Notice that the coefficient matrix of the discrete linear system obtained from the operator splitting does not contain the evaluation of the derivative function J ion (V n+1 m ), i.e., the second part of matrix A. In the numerical simulation of the above methods, the membrane capacitance per unit area and the surface-to-volume ratio are set as C m = 1 (units: µF/cm 2 ) and β = 1, 000 (units: cm −1 ), respectively. Other constant parameters are set as tol Newton = 10 −7 , tol Multigrid = 10 −7 , and the number of smoothing iterations ν Multigrid = 6. The absolute tolerance in the Newton iteration method for solving the nonlinear system (18) is chosen to be 10 −10 . The parameter δ in Equation (13) is selected to be 10 −6 for all simulations. The units for the time steps, spatial steps, and voltage are milliseconds (ms), centimeters (cm), and millivolts (mV), respectively. The unit of conductivity is millisiemens per centimeter (mS/cm −1 ). The fully implicit integration methods are implemented with custom codes written in C++ and the numerical simulations are all performed on a 3.6 GHz computer with an Intel Core i3-4160 CPU.
We should mention that the computational domain is not limited to the regular areas. For example, circular regions can be divided into a quasi-uniform triangular grid by the same principle as shown in Supplementary Figure 8 and a solid sphere can be partitioned into a quasi-uniform tetrahedral mesh (Liu and Joe, 1996;Everett, 2010). The additional numerical results solved on a circle and sphere are shown in Supplementary Figures 9-13.

Numerical Example in Two-Dimensional Space
For the numerical simulations of the bidomain equations in the two-dimensional space, we model the square consisting of fibers which rotate counterclockwise around the point c = (c 0 , c 1 ) = (−0.5, −0.5). The conductivity tensor is space dependent with unequal anisotropy ratios: in the intracellular space, it is 12.0 along the fiber and 2.0 perpendicular to the fiber; in the extracellular space, it is 8.0 along the fiber and 4.0 perpendicular to the fiber. The conductivity tensors D i (x) and D e (x) with x = (x 0 , x 1 ) can be expressed by D i (x) = σ l i e l e T l + σ n i e n e T n and D e (x) = σ l e e l e T l + σ n e e n e T n , where σ l i = 12, σ n i = 2, σ l e = 8, σ n e = 4, e l = (sinθ , cosθ ) T , and e n = (−cosθ , sinθ ) T with cosθ = The positions of the anode and cathode are set as (0.875, 1.0) and (1.125, 1.0), respectively. In this work, we first use a variant of the FitzHugh-Nagumo (Vfhn) model as the membrane dynamics which is governed by with the parameters set as V peak = 1, G = 20, V th = 0.125, α = 3, and γ = 1. The gating and ion concentration variables q, the intracellular potential i , and the extracellular potential e are set to be at rest at time t = 0.
We first discuss the numerical performance of the CBDF2 scheme. Using the same timestep size t = 1/32 and spatial step size h = 1/64, we compare the CBDF2 scheme with the other different time integration methods. A high precision solution obtained from the BE method using a very small timestep size t = 1/256 and spatial step size h = 1/256 is set as a benchmark. As shown in Figure 1, the CBDF2 scheme using a large timestep can obtain a highly accurate trajectory of the membrane voltage compared with the benchmark while that obtained from the BE method using a large timestep is inaccurate. Besides, the Strang operator splitting (SOS) combined with the CBDF2 scheme is stable but less accurate than the CBDF2 scheme in the membrane voltage trajectory. We point out that the Godunov operator splitting combined with the BE method is unstable using the timestep size t = 1/32 and spatial step size h = 1/64. Figure 2 shows the isopotential lines V m = 0.1 at time T = 8.75 for the bidomain simulations from the above methods, while Figure 3 shows the isopotential lines V m = 0.1 and the Newton iteration number obtained from the CBDF2 scheme using different timesteps.
We next study the convergence of the four different time schemes, i.e., the BE method, the CBDF2 scheme, the Godunov splitting combined with the BE method, and the Strang splitting combined with the CBDF2 scheme with the Vfhn membrane model. Comparing with the high precision numerical solution obtained from the BE method with a quite small timestep size t = 1/512 and spatial step size h = 1/512, we compute the errors of the membrane voltage which is summed over the whole space at the discrete times. To be specific, we use the scaled to compute the error of voltage summed over the whole domain, i.e., e (x) = (e(x 1 ), e(x 2 ), · · · , e(x N )) T . By halving both the timestep and spatial step size, we can obtain a first-order convergence of the BE method and the Godunov splitting combined with the BE method and a second-order convergence of the CBDF2 scheme and the Strang splitting combined with the CBDF2 scheme, respectively, as shown in Table 1.
We further consider two more realistic membrane models, i.e., the DiFrancesco and Noble (DFN) model (DiFrancesco and Noble, 1985;Cabo and Barr, 2006) and the Courtemanche et al. (CRN) model (Courtemanche et al., 1998) which have 15 and 20 state variables other than the transmembrane potential, respectively. In each simulation, the state variables q are all assumed to be at rest at time t = 0 and the initial condition of the membrane potential V m is given by The iso-contours and trajectories of the membrane voltage obtained from the CBDF2 scheme using the more realistic DFN and CRN models in the two-dimensional space is shown in Supplementary Figure 6.
We study the stability of the four different time integration methods with the DFN membrane model. The simulations for the bidomain equation with the DFN model use timestep  sizes of 1/8, 1/16, 1/32, 1/64, 1/128, and spatial step sizes of 1/128. As shown in Table 2, the CBDF2 scheme is stable for the timestep sizes ≤1/8 while the BE method is stable for timestep sizes ≤1/32. The Strang operator splitting together with the CBDF2 scheme is stable for the timestep sizes less than or equal to 1/8 while the Godunov operator splitting together with the BE method is stable for the timestep sizes less than or equal to 1/16. The numerical results show that the CBDF2 scheme has better stability compared with the BE method which is consistent with the stability analysis in section 2.4.
As shown in Table 2, the Godunov splitting combined with the BE method is less stable than the Strang splitting combined with the CBDF2 scheme. A possible reason is that the BE method is less stable than the CBDF2 scheme. To fairly compare the stability of the Godunov splitting and Strang splitting, we further perform numerical experiments using Strang splitting combined with the BE method which allows using a much larger time step (maximum t ≈0.20) than the Godunov splitting combined with the BE method. Our experiments indicate that the Strang splitting has better stability than the Godunov splitting for the bidomain equations.
To compare the numerical accuracy of each scheme, a high precision solution is obtained from the BE method with a very small timestep t = 1/1, 024 and spatial step size h = 1/128. As shown in Table 2, the CBDF2 scheme using a large timestep t = 1/8 is more accurate and efficient compared with the BE method using a small timestep 1/64. Besides, the CBDF2 scheme, combined with the nonlinear elimination method, has better accuracy than the CBDF2 scheme combined with Strang operator splitting although taking a bit more time as shown in Table 2. Note that the errors in Table 2 do not convergence to zero as the timestep size is reduced. This is because the errors  depend on both the timestep size and spatial step size which have a lower bound when using a fixed spatial step size. We emphasize that the computational domain is not limited to the regular areas such as rectangle regions. Numerical results solved on a circular region which is partitioned into a quasiuniform triangular grid (Liu and Joe, 1996;Everett, 2010) are given in Supplementary Figures 9, 10. The performance of the CBDF2 scheme is as well as that solved on the rectangle region.

Numerical Example in Three-Dimensional Space
The numerical methods introduced in this work can be straightforwardly extended to the case of three-dimensional space. For ease of illustration, we first simulate the bidomain equations with the Vfhn membrane on a unit cube x ∈ [0, 1] 3 (x = (x 1 , x 2 , x 3 ) T ). The extracellular stimulus in the three-dimensional space is given by applying a virtual battery with an anode and a cathode located at (0.375, 0.5, 0.5) and (0.625, 0.5, 0.5) in the center of the domain and is given by where r is the distance between x and the anode and s is the distance between x and the cathode. The computational domain in three-dimensional space is partitioned into a quasi-uniform hexahedron mesh. We use a linear approximation to the solution on each element. The parameters, tolerances, and modules are set the same as those noted previously in this work.
In the case of three-dimensional space, the axially isotropic conductivity tensors are formulated as D i = σ n i I + (σ l i − σ n i )e l e T l and D e = σ n e I + (σ l e − σ n e )e l e T l , where I is the 3 × 3 identity matrix, σ l i = 12, σ n i = 2, σ l e = 8, and σ n e = 4. The readers can refer to Neu and Krassowska (1993) and Ying (2005) for more details of the conductivity tensors.
The iso-contour and iso-surface plots of the membrane potential simulated on a unit cube are shown in Supplementary Figures 2-5. Additional numerical results obtained from simulating the bidomain equations in a solid sphere are shown in Supplementary Figures 11-13.
We then consider the DFN and CRN models in threedimensional space. In each simulation, the state variables q are all assumed to be at rest at time t = 0 and the initial condition of the membrane potential V m is given by The iso-contours and trajectories of the membrane voltage obtained from the CBDF2 scheme using the more realistic DFN and CRN models in the three-dimensional space are shown in Supplementary Figure 7.

DISCUSSION
In this work, we have presented two fully implicit methods, the BE and CBDF2 schemes, to numerically simulate the bidomain equations arising from modeling the electrical activity of the heart. The second-order CBDF2 scheme is L-stable for the stiff problems and uses the BE method as building blocks. The CBDF2 scheme has better accuracy and efficiency using a large timestep t = 1/8 compared with the BE method using a small timestep t = 1/64. When the error of membrane voltage for the CBDF2 and BE schemes are approximately the same, the CBDF2 scheme allows for timestep an order of magnitude larger than that used in the BE method. A further advantage of the CBDF2 scheme is that it allows temporal adaptivity to speedup, i.e., it can use a small timestep during the action potential upstroke period (stiff region) and a much larger timestep when the electrical activity is slowly varying (non-stiff region) (Ying, 2005;Whiteley, 2007). The CBDF2 scheme satisfies the requirement of stability, accuracy, and efficiency well which enable us to achieve a quantitative understanding of the relationship between molecular function and the integrated behavior of the cardiac myocyte in health and disease more effectively, e.g., to predict the clinical risk of drug-induced arrhythmias.
When using fully implicit methods to solve the bidomain equations, the obtained nonlinear systems are generally very large. We use a variant of the nonlinear elimination method (Lanzkron and Rose, 1997) to reduce the size of the system which includes the evaluation of the Jacobian matrix for each timestep. The spatially and temporally discretized nonlinear system (11) is solved with an approximate Newton method where the coefficient matrix of the residual system is symmetric and possibly positive definite consisting of the stiffness matrix, mass matrix, and the derivative function J ion . The linear residual equation can be efficiently solved by standard optimal solvers such as the V-cycle multigrid method (Saad, 2003). The derivative function J ion is time dependent and can be numerically evaluated with sufficient accuracy.
The computational domains in this work are square and cube (see Supplementary Figures 9-13 for results simulated on circle and sphere), but we emphasize that the CBDF2 scheme can be applied to more complicated computational domains since the generation rules of the finite element mesh in both two-and three-dimensional space do not depend on the particularity of the square or cube. The generation rules are applicable to more complex computational domains such as the realistic human left ventricle model (Cai et al., 2015) which will be presented in our future work.
In addition, the CBDF2 scheme using an adaptive mesh refinement algorithm can further improve the computational efficiency (Trangenstein and Kim, 2004;Ying, 2005;Whiteley, 2007). For the parallel implementation, the fully implicit methods proposed here leads to a full space independent system which is highly localized. As the inter-processor communication for updating the state variables by solving Equation (7) will be much less, the speedup for the simulations will be impressive when the membrane model including tens of or more state variables. In future works, we will consider more realistic computational domains such as the human left ventricle domain (Cai et al., 2015) combined with more realistic membrane models such as the Tusscher-Noble-Noble-Panfilov model (ten Tusscher et al., 2004) which can potentially help study new drugs and methods of treatment.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.