A convenient scheme for coupling a finite element curvilinear mesh to a finite element voxel mesh: application to the heart

Outline

The present scheme will be described with reference to a particular problem that involves coupling a previously meshed heart[1] to a torso, as shown in the left hand panel in Figure 1. The goal is to solve for the extracellular potential (V _e) throughout the heart and the potential outside of the heart (V _o) in the elliptical equation[7]

∇•(D _i + D _e )∇V _e = -∇•(D _i )∇V _t ; within heart muscle

∇•(σ _o )V _o = 0; in the torso and blood pools     (1)

where D _i and D _e are the intracellular and extracellular conductivity tensors within heart muscle, σ_o represents the conductivity in the torso and blood pools respectively, and V _t are the transmembrane potentials throughout heart muscle. At the interface between the heart and torso, the potential is continuous (V_e = V_o), the normal component of the extracellular current is continuous with the normal component of the current through the torso/blood pools (n^T D _eV_e = n^T σ_o V_o, where n is the unit normal vector), and there is no flow of intracellular current into the torso/blood pools (n^T D _iV_i = 0).

Figure 2 is a flow chart of the coupling scheme. The first step of the method involves embedding the heart in a preliminary voxel mesh that is slightly larger than the heart, as indicated in Figure 3. The next step involves locating preliminary voxel mesh nodes ("boundary nodes") that are outside of the heart but close to the heart, as indicated by the red circles in Figure 3. A surface ("boundary surface") is fitted to these boundary nodes. The boundary surface is essentially a projection of the outer heart surface, as shown in Figure 3. A temporary mesh ("scaffolding mesh") is generated between the outer heart surface and the boundary surface and the stiffness matrix for this mesh is computed. The scaffolding mesh is indicated by the purple area in Figure 3. The scaffolding mesh includes extraneous temporary nodes that are then removed from the stiffness matrix, leaving a coupling matrix that couples the pre-existing outer heart surface nodes to the boundary surface voxel nodes. A voxel mesh is generated for the entire torso by adding voxels to the preliminary mesh. The resulting stiffness matrix for the torso nodes is then computed. Finally, a complete system matrix is assembled from the pre-existing heart stiffness matrix, the heart/torso coupling matrix and the torso stiffness matrix.

Preliminary mesh; boundary extraction

Given the heart surface, a preliminary voxel mesh is generated that is slightly larger than the heart. In the two dimensional example illustrated in Figure 3, the preliminary voxel mesh encloses all of the squares within the blue rectangle. The nodes of this mesh are designated as X_p = (x_p, y_p, z_p). It will be assumed that the voxels are isotropic (side length = s) and chosen to be defined with respect to some convenient origin. The preliminary voxel mesh should extend for a distance of at least 1s in the x, y and z directions beyond any point on the outer heart surface. The voxel size s is preferably the desired voxel size for the entire torso mesh.

The next step involves locating boundary nodes X_b ⊆ X_p that are outside of the heart but closest to it, as illustrated in Figures 3 and 4. The boundary nodes are located by finding boundary voxels V_b, which are defined as those voxels that have at least one node within the heart and at least one node outside of the heart. The boundary nodes X_b ⊆ X_p are defined as all of the nodes within these voxels V_b that are outside of the heart.

There are a number of methods for locating the boundary voxels. One convenient method, which was used in the present example, is applicable to any surface (e.g. the heart surface) that may be characterized by spherical harmonics[1]. In this case, the spherical harmonics provide a continuous representation of a surface, and therefore may be used to determine whether any given point is inside or outside of this surface.

More particularly, the outer heart surface S_h is approximated by R_h = F(θ,Φ), where R_h is the radius of the approximate outer heart surface, F is a spherical harmonic polynomial as will be further described below, and θ and Φ are azimuthal and polar angles, respectively. Figure 4 shows a 2-D representation of R_h. R_h, θ, and Φ are defined with respect to a convenient origin which is at the point (x_h, y_h, z_h). This origin is preferably chosen such that it results in a good spherical harmonic expansion. The grid points X_p are translated so their origin is (x_h, y_h, z_h), which results in a set of translated node points X'_p = (x_p- x_h, y_p - y_h, z_p - z_h). The azimuthal and polar angles θ_p and Φ_p (with respect to the (x_h, y_h, z_h) origin) are computed for X'_p, and then the value F(θ_p, Φ_p) is computed, which is the radius of the outer heart surface at the angles θ_p and Φ_p. If this radius is greater than the radius $\sqrt{{(x_p-x_h)}^2+{(y_p-y_h)}^2+{(z_p-z_h)}^2}$ of a translated node point, then that point is within the heart. The set of node points within the heart will be denoted by the set X_h ⊆ X_p.

The above procedure locates the nodes X_h ⊆ X_p within the heart. The boundary voxels V_b are then defined as those voxels that have at least one node within the heart (X_h) and at least one node outside of the heart. The boundary points X_b are those nodes that are within the boundary voxels but outside of the heart.

Surface formation

A surface (the boundary surface S_b) is generated by fitting a spherical harmonic expansion to the points X_b about the origin (x_h, y_h, z_h). The radii (R_b) and azimuthal and polar angles θ_b and Φ_b are computed for the boundary points, as indicated in Figure 4, which shows only the azimuthal angle. The spherical harmonic expansion of the radii R_b may be expressed as

$\begin{array}{l}{R}_b={\displaystyle \sum _{i=0}^N\left\{a_i{P}_i(\cos {\varphi }_b)+{\displaystyle \sum _{m=1}^i(a_{im}\cos (m{\theta }_b)+b_{im}\sin (m{\theta }_b))P_i^m(\cos {\varphi }_b)}\right\}}\hfill \\ R_b=P\ast C\hfill \end{array}\left(2\right)$

where N is the number of terms in the spherical harmonic expansion, P_i is the Legendre function of degree i of the first kind, P_i ^m is the associated Legendre function of degree i and order m, a_i, a_im, b_im are coefficients (to be determined), P and C are matrix and vector representations of the Legendre function values and coefficients, respectively, and R _b is a vector representation of the radii R_b. A least squares fit (C = P\R _b in Matlab notation) results in a set of coefficients C_bs that define the boundary surface S_b characterized by R_bs = F(θ,Φ), where F is the spherical harmonic expansion of order N with coefficients C_bs, and R_bs is the radius of a point on the surface S_b at a given pair of azimuthal and polar angles (θ, Φ).

The appropriate value of N will depend on the geometry of the problem and may be determined empirically. For the heart shown in Figure 1, N = 8 produced good results.

Scaffolding mesh generation

The next step involves forming a scaffolding mesh that couples the original heart nodes X_h to the boundary nodes X_b. The scaffolding mesh, which is indicated by the purple region in Figure 3, joins the heart surface S_h to the boundary surface S_b. With reference to Figure 5, which again shows a two dimensional example of the present scheme, the scaffolding mesh is shown in light purple and one of the scaffolding mesh elements is indicated by a dotted black quadrilateral. The scaffolding mesh provides an approximation of the electrical coupling between the original heart mesh nodes (X_h, shown as black filled circles) and the boundary nodes (X_b, shown as red filled circles). The scaffolding mesh is only an estimate of this electrical coupling because the mesh does not include the actual X_b nodes, but approximations (dotted red circles) of the locations of the boundary nodes (X_b).

In three dimensions, the scaffolding mesh comprises prismatic elements. Each prismatic element consists of two similar triangles, one on the heart surface S_h and one on the boundary surface S_b, whose corresponding vertices are joined to one another_. The triangles_, in turn, are formed from combinations of the original heart nodes X_h and boundary nodes X_b, and projections of these nodes_.

More particularly, the heart nodes X_h are projected onto the surface S_b by a line that passes through the heart nodes and the heart origin (x_h, y_h, z_h), as illustrated in Figure 4. Similarly, the boundary nodes X_b are projected onto both the surface S_b and the surface S_h by a line that passes through these nodes and the heart origin (x_h, y_h, z_h). The projection of the boundary nodes on to S_b will be referred to as X_b'. From an implementation standpoint, the projections may be accomplished by calculating the (θ,Φ) pair for any point (about the (x_h, y_h, z_h) origin), and then computing the radius R = F(θ,Φ) for either the heart surface S_h or boundary surface S_b from its corresponding spherical harmonic expansion F(θ,Φ).

The result of the mutual projections is two sets of points, one on S_h and the other on S_b, such that each point at a given (θ,Φ) on a surface has a corresponding point on the other surface at the same (θ,Φ). Thus, a triangular tessellation of one surface may be used to generate an identical triangular tessellation of the other surface.

There are many ways to tessellate a surface. The approach that was adopted took advantage of the regularity of the (θ,Φ) coordinates of the heart nodes X_h. These nodes formed a regular grid in the (θ,Φ) coordinates by choice; it is easy to form a such a regular grid when the surface at issue (S_b) is characterized by spherical harmonics. Because of the cyclicity of the azimuthal angle θ, the original (θ,Φ) grid (-π<=θ<=π) points were augmented such that those points (θ_a, Φ_a) with θ_a<=-4π/5 were identified, and new points with (θ_a+2*π, Φ_a) added to the grid. This "wrapped around" grid was triangulated, and the triangles characterized by all three points having θ>2*π were removed. This type of tessellation is generally acceptable only in cases when there are points at both poles (Φ = +/-pi) and there is a reasonably smooth gradation in the Φ values of the points. This condition was satisfied in the present case because the original heart points X_h were generated via a spherical harmonic expansion with an even division in the Φ direction.

The left panel in Figure 6 shows the surface of the example problem heart with both the original heart surface points (blue) and projected grid boundary points X_b (red) triangulated in the manner described above. The boundary surface (right hand panel in Figure) has shape very similar to that of the heart surface and the triangulation is identical to that shown for the heart surface.

Because each triangle on the heart surface has a corresponding triangle on the boundary surface, connecting the corresponding triangle nodes by radial line segments defines a prismatic element. All of the prismatic elements define the scaffolding mesh.

Heart to torso transfer matrix

The finite element stiffness matrix (A^k) for the prismatic element (scaffolding) mesh describes the electrical coupling between the heart surface and the boundary surface. However, this matrix contains extraneous elements, i.e. the projected heart nodes on the boundary surface (X_hp), and the projected boundary nodes on the heart surface (X_bp). These extraneous nodes must be removed, resulting in a transfer matrix between the heart nodes and the boundary nodes.

The concept of extraneous nodes will be further described with reference to Figure 5. The desired final matrix should reflect only electrical coupling between the original heart mesh nodes (X_h) shown as filled black circles and the boundary nodes (X_b) shown as filled red circles. The positions of the boundary nodes (X_b) are approximated by the dotted red circles (X_b'), so in actuality the final desired stiffness matrix will reflect the coupling between the filled black circles and the dotted red circles. However, the scaffolding matrix also contains electrical coupling between the desired nodes (filled black and dotted red circles) and extraneous nodes (circles with black and red diagonal lines), which are the projected heart nodes on the boundary surface (X_hp), and the projected boundary nodes on the heart surface (X_bp). These extraneous nodes must not exist in the desired final system matrix.

First, to assemble the finite element stiffness matrix (A^k) for the prismatic element mesh, the local stiffness matrix for each prism element is constructed. This may be done in a number of ways. For a Galerkin type FEM, a closed form local matrix based on the prismatic geometry may be used. Alternatively, each prism may be decomposed into tetrahedrons, and the stiffness matrix, based on linear basis functions, may be computed for each tetrahedron.

To remove the extraneous nodes in A^k, it was partitioned as follows:

$A^k=\left[\begin{array}{ccc}{A}_{hh}^k& A_{he}^k& A_{h{b}^{\prime }}^k\\ A_{eh}^k& A_{ee}^k& A_{e{b}^{\prime }}^k\\ A_{b^{\prime }h}^k& A_{b^{\prime }e}^k& A_{b^{\prime }{b}^{\prime }}^k\end{array}\right]\left(3\right)$

where the subscript h denotes the heart surface nodes X_h, the subscript e denotes the extraneous nodes (X_hp ∪ X_bp) to be removed, and the subscript b' denotes the projected boundary nodes X_b'. The desired coupling matrix (A^c) between X_h and X_b' is found by solving for the extraneous nodes in the set of equations A^k = 0:

$A^c=\left[\begin{array}{cc}{A}_{hh}^c& A_{h{b}^{\prime }}^c\\ A_{b^{\prime }h}^c& A_{b^{\prime }{b}^{\prime }}^c\end{array}\right]\left(4\right)$

where

$\begin{array}{l}{A}_{hh}^c=A_{hh}^k-A_{he}^k{(A_{ee}^k)}^{-1}{A}_{eh}^k\hfill \\ A_{hb}^c=A_{h{b}^{\prime }}^k-A_{he}^k{(A_{ee}^k)}^{-1}{A}_{e{b}^{\prime }}^k\hfill \\ A_{b^{\prime }h}^c=A_{b^{\prime }h}^k-A_{b^{\prime }e}^k{(A_{ee}^k)}^{-1}{A}_{eh}^k\hfill \\ A_{b^{\prime }{b}^{\prime }}^c=A_{b^{\prime }{b}^{\prime }}^k-A_{b^{\prime }e}^k{(A_{ee}^k)}^{-1}{A}_{e{b}^{\prime }}^k\hfill \end{array}$

The set of equations A^k = 0 is somewhat of an oversimplification of the process of eliminating extraneous nodes. More precisely, the right hand side corresponding to each extraneous node is set to 0. To take a simple example, assume it is desired to eliminate a variable x_2. If the row in the matrix A^k corresponding to x₂ is equal to -3x₁ + 4x₂ - 2x₃, then setting this row equal to 0 and solving for x₂ yields x₂ = 1/4(3x₁ + 2x₃). 1/4(3x₁ + 2x₃) is then substituted for x₂ wherever it appears in any other equation, effectively eliminating x₂ from the system of equations.

Total system matrix

The total system matrix comprises a combination of the preexisting FEM heart matrix (A^o), the FEM heart to torso coupling matrix A^c, and the FEM torso matrix A^t.

The FEM torso matrix A^t is based on a voxel based mesh of the torso. If the voxel torso mesh did not exist a priori, it may be generated by extending the preliminary mesh until it approaches the outer torso, as indicated by the green filled circles shown in Figure 3. The outer torso nodes may be found in a number of ways. In the approach adopted, the original torso surface nodes were multiply interpolated with two dimensional Legendre polynomials, resulting in a set of surface indicator functions analogous to the spherical harmonic surface indicator function described with reference to generation of the boundary surface.

The FEM stiffness matrix A^t for the voxel matrix was computed once again with a Galerkin type scheme with linear basis functions. For voxel meshes, this type of FEM matrix may be generated very quickly, on the order of a few seconds assuming reasonable computer processing speed and memory.

The FEM matrix (A^s) for the entire system is:

$A^s=\left[\begin{array}{cc}{A}_{ii}^o& A_{ih}^o\\ A_{hi}^o& A_{hh}^c+A_{hh}^o& A_{h{b}^{\prime }}^c\\ A_{b^{\prime }h}^c& A_{b^{\prime }{b}^{\prime }}^c+A_{b^{\prime }{b}^{\prime }}^t& A_{b^{\prime }g}^t\\ A_{g{b}^{\prime }}^t& A_{gg}^t\end{array}\right]\left(5\right)$

where the subscript g denotes the set of torso mesh nodes excluding the boundary nodes, and the subscript i denotes the set of original heart nodes excluding the outer surface heart nodes X_h. In this equation, the preexisting FEM heart matrix A^o has been partitioned into components: (i) that connect interior heart nodes to other interior heart nodes ($A_{ii}^o$); (ii) that connect interior heart nodes to nodes on the heart surface ($A_{hi}^o$ and $A_{ih}^o$); and (iii) that connect outer heart nodes to other outer heart nodes ($A_{hh}^o$). Similarly, A^c and A^t have been partitioned into portions: (i) that connect "internal nodes" (i.e. nodes that are not a part of any other matrix) to other "internal nodes"; (ii) that connect "internal nodes" to "overlapping nodes" (i.e. nodes that appear in two matrices); and (iii) that that connect "overlapping nodes" to "overlapping nodes."

For example, $A_{hh}^o$ represents the coupling between the outer heart surface nodes, which are overlapping nodes, to themselves in the matrix A ^o . Similarly, $A_{hh}^c$ represents the coupling between the outer heart surface nodes to themselves in the matrix A ^c . The matrix entry in A^s for the coupling of the outer heart nodes to themselves is thus the combined entry $A_{hh}^o$ + $A_{hh}^c$.

The discretized version of equation (1) is A^s V _e = A^{o, rhs} V _t, where V _e and V _t are the vectors corresponding to the continuous variables in equation (1), and A^{o, rhs} is the original, pre-computed stiffness matrix corresponding to the right hand side of equation (1). A^{o, rhs} and V _t must be appropriately augmented with 0's to match the size of A^s V _e. The boundary conditions at the heart/torso boundary, that V_e is continuous and the normal component of the current is continuous[3], are necessarily enforced by the finite element formulation used to create A^s.

Test

The above described scheme was tested as follows. A spherical shell, analogous to the outer heart surface, was coupled to a surrounding spherical "torso". The potentials on the inner spherical shell were set equal to those which would have occurred had there been an ideal current dipole, normalized to unit strength, located at the sphere's origin and oriented along the z axis. Given this potential distribution, the potentials throughout the spherical "torso" were computed and compared with the analytic solution[8]:

$V(r,\varphi )=\cos (\varphi )\ast \left(\frac{1}{r^2}+\frac{2R}{r^3}\right)\left(6\right)$

where Φ is the polar angle (as before), R is the radius of the outer spherical torso surface, and r is the radius to an observation point. For the test problem, R was set at 50. The radius of the spherical shell/heart surface was set at 40, roughly the same size as the heart shown in Figure 1. To create the potential distribution equivalent to a dipole oriented along the z axis, the potentials on this sphere/heart surface were thus set at:

$V(40,\varphi )=\cos (\varphi )\ast \left(\frac{1}{{40}^2}+\frac{2\ast 50}{{40}^3}\right).\left(7\right)$

A voxel mesh, akin to the torso mesh shown in Figure 3, was generated to fill in the space between the inner and outer spherical surfaces, and the corresponding FEM matrix A^t was computed. Next, a coupling matrix A^c to link the torso to the nodes on the spherical shell/heart surface was created. The system of equations A^s *V _c = 0 was constructed, where A^s is the total system FEM matrix and V _c is the solution. The spherical shell/heart surface nodes were set equal to their corresponding potential V(40, Φ), resulting in an augmented system A'*V' = b, where V(i)' = 1 where i represents the spherical shell/heart surface nodes. This system was solved using the biconjugate gradient stable method, after preconditioning by scaling each row in the matrix so that its diagonal element was equal to 1.

Relative error was defined as:

$RE=\sqrt{\frac{{\displaystyle \sum {(V_c-V_a)}^2}}{{\displaystyle \sum V_a^2}}}\left(8\right)$

where V_c is the computed solution (excluding the spherical shell/heart nodes) and V_a is the analytical solution. Both the number of nodes on the spherical shell/heart surface and the number of nodes in the "torso" were varied.

Heart beat

The heart in Figure 1, which is based on the Auckland canine model[2], was meshed[1] and oriented within a human shaped torso, as shown. A cellular automata like scheme was used to generate an activation and repolarization sequence throughout the heart, resulting in a corresponding sequence of transmembrane potentials V _t . Maximum V _t was 100 mV, with the resting potential set at 0 mV. Initial left ventricular and right ventricular activation sites were chosen roughly in accord with Selvester et al.[9] but were adjusted to generate better matches of epicardial potentials with recordings made by Spach et al.[10] in intact chimpanzees.

A propagation sequence was generated by a simple eikonal like scheme, according to which propagation velocity (v) was assumed to be proportional to the square root of bulk cardiac tissue conductivity[11]. The eikonal like scheme was implemented by calculating a propagation time between connected nodes in the heart mesh, selecting initial activation sites, and then stepping through time to activate the remaining heart nodes. The activation time of a heart node was set equal to the activation time of its nearest neighbor plus the approximate propagation time from the heart node to this nearest neighbor. The propagation time between neighboring nodes was set equal to v/d, where v is velocity (proportional to the square root of the conductance between nodes) and d is the distance between nodes.

Repolarization was initiated from the right ventricular breakthrough area[10] and allowed to spread like a wave across the epicardium in the manner in which activation was simulated. Repolarization was then allowed to proceed inward transmurally toward the epicardium based simply on relative transmural depth. In other words, repolarization proceeded from the epicardium inward. Although this oversimplified repolarization sequence is not necessarily strictly physiologically correct[12], it was chosen because of convenience and because it produced reasonable intra-epicardial gradients and epicardial/endocardial transmembrane gradients that resulted in realistic T waves.

Within the heart muscle: (i) intracellular conductivity along and transverse to fiber direction was set equal to 1 and 1/6, respectively; (ii). extracellular conductivity along and transverse to fiber direction was set equal to 1 and 1/3, respectively. Ventricular blood and torso conductivity were set at 5 and 2.5, respectively.

The heart mesh had approximately 270000 nodes and the torso mesh had approximately 120000 nodes. The system of equations was solved with Matlab's preconditioned conjugate gradient routine. The preconditioner was an incomplete Cholesky (0 fill in) factorization of the system matrix.