PDE-induced connection of moving frames for the Atlas of the cardiac electric propagation on 2D atrium

03/11/2020 ∙ by Sehun Chun, et al. ∙ Yonsei University Imperial College London 0

As another critical implementation of moving frames for partial differential equations, this paper proposes a novel numerical scheme by aligning one of three orthogonal unit vectors at each grid point along the direction of a wave propagation to construct an organized set of frames, called a connection. This connection characterizes the geometry of wave propagation depending on (1) the initial point, (2) type of wave, and (3) shape of the domain with conduction properties. The constructed connection is differentiated again to derive the Riemann curvature tensor of orthonormal bases corresponding to important physical and biological meanings in wave propagation. As a practical application, the proposed scheme is applied to diffusion-reaction equations to obtain the Atlas, or a geometric map with connections, of an atrium with cardiac fibers, for the quantitative and qualitative analysis of cardiac action potential propagation, which could contribute to the clinical and surgical planning of atrial fibrillation.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 17

page 23

page 24

page 31

page 32

page 35

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In the previous works on adapting moving frames to obrain the numerical solution of PDEs on curved surfaces such as conservational laws MMF1, diffusion equation MMF2, shallow water equations MMF3, and Maxwell’s equations MMF4, the direction of the moving frames is given at random on the tangent plane in the absence of anisotropy. Because the first moving frames, denoted as , should be differentiable within each element, the direction of moving frames cannot be independent of the direction of neighboring moving frames to ensure differentiability. However, their average direction is still random. Let us call this type of distribution as moving frames of order zero according to the Cartan’s original definitions Cartan4. In this case, moving frames are nothing more than an orthonormal local reference reflecting the geometry of the tangent plane.

The moving frames of first order are naturally obtained by aligning in the direction of the tangent vector along the trajectory of a wave. In Maxwell’s equations, for example, the first moving frame is in the same direction as the electric or magnetic field, and in the diffusion-reaction equation, it is in the same direction as the gradient of the variable that represents an electric potential or material concentration. The moving frames of second order and moving frames of third order are obtained by additionally aligning and in the binormal and normal directions, respectively. In many cases, the first-order moving frames are naturally the second and third-order moving frames. Thus, we may refer to such moving frames of non-zero order collectively as the Darboux frame HCartan.

(a) Moving frames of order zero
(b) Moving frames of first order
Figure 1: Illustrations of moving frames of different kinds. The dashed line is the trajectory. and are the first and second moving frames, respectively.

1.1 Riemann curvature tensor of wave propagation

The main objective of aligning the first moving frame in the direction of propagation is to analyze the trajectory of wave propagation qualitatively. This can be achieved by computing the Riemann curvature tensor of trajectory. This curvature is related to the shape of the domain, but the shape is not the sole factor of the Riemann curvature. For example, the plane is flat with zero curvature. However, the trajectory of wave propagation can have non-zero curvature; a planar propagation with flat wavefront has zero curvature because the trajectory is equally spaced and parallel. However, a point-initiated wave propagation has a trajectory that is aligned along the radial direction in the polar coordinate axis. Thus, the curvature of the trajectory is the same as the polar coordinate axis with non-zero curvature. The trajectory can have various curvature depending on many factors, particularly (1) initiation type and location, (2) type of wave, (3) shape of the domain, and (4) property of the medium.

The Riemann curvature tensor has twenty-one components in three dimensions to indicate how moving frames change in each direction. One of the most critical components in this tensor is the component of in the direction of the wavefront when is aligned along the propagational direction, as shown in Fig. 2. If is parameterized with , then indicates the acceleration of along the wavefront. Note that is the normalization of the separation vector Misner; thus, corresponds to the degree that the trajectory diverges or converges depending on its sign and magnitude.

A significant magnitude of the Riemann curvature component always has critical meaning in physics and biology. For example, in spacetime physics, the non-trivial magnitude of this curvature component for electromagnetic wave propagation can suggest the presence of mass Misner; Dray. In the heart, the non-trivial magnitude and the positive of this curvature component may indicate the stopping condition for the biological electric flow, as suggested in ref. MyBIOP. The objective of creating an Atlas of the geometry map of the electric current flow on the atrium is to qualitatively and quantitatively analyze the stopping condition of the cardiac electric propagation on a patient-specific atrium structure to reveal a unidirectional block in order to cause atrial reentry and consequently fibrillation. Moreover, the Atlas may explain the role of the cardiac fiber in the heart, which has been unknown to the cardiology community and may promote more efficient surgical planning.

Figure 2: Illustrations of the separation vector of trajectory and its parameter.

1.2 Motivations for cardiac electric propagation analysis

This paper describes the general numerical scheme to derive the Atlas of wave propagation with the subsequent connection form and Riemann curvature tensor of orthonormal bases while also providing an exemplary use for the analysis of the cardiac electric signal propagation. The motivation of this practical application will be described briefly using relevant literature in this subsection.

Cardiac electric signal propagation is often identified by a series of wavefronts in silico or in vivo observations. Then, the analysis of the stopping conditions or conditions for spiral wave or fibrillation is conducted by studying the shape of the wavefront, otherwise known as kinematics analysis Mikhailov; Zykov; alternatively, spatiotemporal analysis by phase mapping after introduced in the 1980s by Winfree, can be used Winfree; Gray; Clayton. However, both methods seem to fail to reflect the multidimensional anisotropic features of cardiac electric signal propagation in the atrium and ventricles completely. This is because the kinematic analysis cannot incorporate both the three-dimensional shape and anisotropy into the analysis MyBIOP, whereas the phase mapping analysis has other limitations Umapathy.

The proposed novel scheme uses the trajectory, instead of the wavefront, to analyze the spatial-temporal propagational pattern of cardiac electric signal propagation. One advantage of this approach is that the geometrical factors, such as the three-dimensional shape, direction, and strength of the cardiac fiber, can be conveniently incorporated into the Atlas in the propagation. Thus, the direct causal relation between the geometry and the propagational pattern can be more easily identified. Another advantage is that, similar to space-time physics, all different types of cardiac physiological phenomena that affect the electric signal propagation, such as cardiac restitution, cardiac memory, or the presence of infarction or fibrosis, can be incorporated into this geometric model to be interpreted considering geometric deformation.

1.3 Goals

This paper may provide a more significant result that can be immediately used for the cardiac electric signal propagation in the complex multidimensional and anisotropic heart. However, as the first step, this paper focuses on the description of the numerical scheme, validations, and exemplary use of this numerical scheme for cardiac electrophysiological analysis. The goals of this paper are as follows: (1) to develop an algorithm to efficiently align the first moving frames to the gradient of a variable in order to achieve the Darboux frames; (2) to compute the connection form and Riemann curvature tensor using Cartan’s structure equation for diverse propagation on two-dimensional curved surfaces; (3) to apply the scheme to an actual atrium with cardiac fiber to understand the role of cardiac fiber in the propagation of the cardiac action potential.

The remainder of this paper is organized as follows. Chapter 2 describes the formulation to compute the connection form, and Chapter 3 explains the importance of a component of the connection form,

, with respect to the Riemann curvature tensor. Chapter 4 elaborates this importance in cardiac electric propagation. Chapter 5 describes the numerical scheme to derive the propagational direction from a diffusion-reaction equation. In Chapter 6, the numerical scheme to solve a diffusion-reaction equation is briefly described. Chapters 7 and 8 provide the exponential convergence of moving frames, connection form, and Riemann curvature for the test in the plane and on the sphere, respectively. Chapter 9 illustrates the meaning of anisotropy in trajectory tracking for cardiac electric propagation. Chapter 10 describes the pre-processing of raw fiber data. In Chapter 11, the two-dimensional cardiac fiber is analyzed to provide a geometric meaning to the propagation. In Chapter 12, the connection from the numerical simulation of the cardiac electric propagation is presented and explained for future clinical applications. The discussion is provided in Chapter 13.

2 Connection Form

Consider a set of orthonormal vectors, known as

moving frames at every grid point in a curved surface . Let two frames and lie on the tangent plane of the surface and be aligned along the surface normal direction. Because moving frames are orthonormal, they satisfy the following conditions

where is the Kronecker delta. The only freedom that the moving frames have is with respect to the direction of the first moving frame on the tangent plane. The directions of other frames are subsequently determined because of the orthogonality of the frames. However, refs. MMF1 and MMF2 proved that the direction of can affect the integration and differentiation error in differential operators such as divergence and Laplacian operators.

In this paper, is aligned such that the connectivity between the moving frames is considered; this is termed as a specifically connected and ordered set of frames. We may refer to this configuration of frames in the domain as a connection. There are infinitely many ways to align moving frames. However, we choose to be aligned along the direction of the wave propagation, which is the same as the direction of a vector or the gradient of a variable. Then, the pattern of the trajectory, i.e., the convergence or divergence of the bundle of trajectories, can be expressed by the connection form.

At every point in the domain, the set of moving frames is expressed in the following matrix form: for unit vectors along the Cartesian coordinate axes , respectively, we obtain

(1)

Or,

(2)

where we introduced a new tensor , . The matrix represents the orientation of moving frames; thus, it is known as the orientation matrix or attitude matrix ONeil. By using the fact that is fixed everywhere and by computing the infinitesimal displacement of on a surface, the 1-form is expressed as follows.

where we used Eq. (2) and the fact that is an orthonormal matrix. Briefly speaking, the 1-form is the general type of derivative such that it has a scalar value when a vector is chosen. For example, let be a scalar; then, the exterior derivative of , or , is the 1-form. Consequently, is a scalar value corresponding to the directional derivative of in the direction of , or Misner.

Then, is a 1-form with nine components and is called the connection form Cartan1; Cartan2. In words, represents the amount of rotation of with respect to when it moves in the direction of Piuze2015. However, the orthonormality of moving frames yields the following equalities.

(3)
(4)

Eq. (3) yields for and Eq. (4) yields . Thus, only has three independent components, namely

Or, we have

Additionally, 1-form is linear such that , it would be sufficient to compute for . For simplicity, the following notation is used for connection form such that the magnitude of is the same as that of , but possibly with opposite signs, and . Then, the connection form can be obtained as

Thus,

(5)

where the matrix is the Jacobian matrix for . There are nine components of the connection form, three for each of , , and . Out of the three, and are related to the Gaussian curvature () and mean curvature () because they can be also expressed in terms of connections as follows Treibergs

In a two-dimensional plane, and are zero along all the directions. However, various connections of moving frames can create a non-zero along a certain direction. In this paper, we are particularly interested in the 1-form along the direction of the second frame because the corresponding quantity is related to an important component of the Riemann curvature tensor as will be explained in the following section.

3 Importance of

One of the benefits of deriving the connection form of moving frames is obtaining a Riemann curvature tensor representing important features of a physical or biological phenomena. However, computing the Riemann curvature tensor requires the metric tensor of the curved domain, or the length of the tangent vector of the curved axis of the surface. We refer to this tensor as the Riemann curvature tensor of coordinate basis (the Riemann curvature tensor of holonomic basis). However, computing this Riemann curvature tensor of coordinate basis is computationally challenging for the general curved surface, and most of all, is not always required.

Instead, the length of the tangent vector of the trajectory is normalized and the Riemann curvature tensor is computed for the normalized trajectory. This yield the Riemann curvature tensor of orthonormal basis. Geometrically, this Riemann curvature describes the same qualitative behavior for the divergence and convergence of the trajectory in the orthogonal direction of propagation, or in the direction of the separation vector, but its magnitude may be changed proportionally. Note that the Riemann curvature tensor is exactly same as this Riemann curvature tensor of orthonormal basis if the metric tensor is constant.

To derive the Riemann curvature tensor of orthonormal basis, the second Cartan structure equation is recalled as follows Cartan1; Cartan2.

where is a 2-form, is the wedge product, and is known as the curvature 2-form that is defined by Dray. The Riemann curvature tensor is defined in relation with the curvature 2-form as follows.

where the 1-form representing an orthonormal dual frame satisfies where is the Kronecker delta. The overline notation is introduced to indicate the corresponding Riemann curvature computed with orthonormal bases. For the connection component , can be computed in the following manner.

(6)

where we used the property and . The numerical evaluation of can be obtained as follows. Consider expressed on a two-dimensional surface as follows.

Then, the derivative of this 1-form is obtained in the form of directional derivative by using the variable from Eq. (5) as follows.

(7)

The exact calculation of should be obtained based on , the axis aligned along the same direction as but with the non-unit magnitude of the tangent vector. Here, is differentiated based on the unit axis rather than the underlying curved axis. This yields the Riemann curvature tensor of orthonormal bases, which differs from that of coordinate bases.

The second component in the right hand side of Eq. (6) is known as the Gaussian torsion Cartan1 and can be computed as follows.

(8)

However, there exists a property that reduces this Gaussian torsion to null. Substituting Eqs. (7) and (8) into Eq. (6), the following formula are obtained for the nontrivial components of the Riemann curvature tensor of orthonormal basis.

(9)

Note that is negligible because the rotation of with respect to when moving toward remains fixed orthogonally. Then, is approximately the same as the change of in the direction of on a surface.

(a) ,
(b) .
(c) .
Figure 3: Intuitive relationship between wavefront curvature and Riemann curvature. Solid line = trajectory, Dashed line = wavefront.

4 Stopping criteria in Riemann curvature

As a large Riemann curvature implies the presence of an object with a large mass in the universe such as a black hole, a large Riemann curvature in the cardiac excitation propagation corresponds to a similar unique feature in the cardiac tissue. In the cardiac excitation propagation, one of the most important analyses of the propagation is related to whether it continues to propagate or stops. The stopping conditions are the most critical knowledge for understanding many cardiac pathologies. However, reliable analytical or computational methodologies remain unknown considering the strong anisotropy and complex shape of the heart. Initial attempts to identify the general stopping conditions were related to the shape of the wavefront, known as the kinematic analysis; the linearized version of the propagation velocity () can be obtained as Davydov1991

(10)

where is the conducting velocity determined by the cardiac tissue, and is the diffusivity constant, and is the curvature of the wavefront. Eq. (10) implies that the propagation velocity of the cardiac propagation depends on the curvature of the wavefront. If the conducting properties of cardiac tissue are homogeneous, then the only important factor for slowing down the propagation is the shape of the wavefront. If the curvature is too large or above a certain threshold that is determined by the conducting properties of the cardiac tissue, then the propagation velocity dramatically reduces, even to zero.

On the cellular scale, the large curvature of the wavefront is related to a phenomenon called impedence mismatch or sink-source mismatch Rohr; Xie, which is a higher ratio between the number of excited cells and the number of excitable cells than the threshold ratio at which propagation can occur. In this paper, it is assumed that this sink-source mismatch can also be measured by the changes in the connection form or the Riemann curvature tensor of orthonormal basis as follows.

Hypothesis 0.

If the connection component changes abruptly in the direction of , or equivalently the positive is sufficiently large, then the cardiac action potential propagation slows down dramatically, even to zero.

The rigorous validation of this hypothesis is beyond the scope of this paper, but the complementary relationship between the wavefront and trajectory can be intuitively given in Fig. 3. Moreover, Fig. 4 demonstrates the corresponding large quantity of and for the classical gap propagation test case to demonstrate a relation between strong wavefront curvature and large positive . From the geometric perspective, the Riemann curvature tensor measures the convergence or divergence of the geodesics of propagational rays Misner; Dray. To measure whether trajectories are pushed closer or pulled apart, the first moving frame () is aligned along the propagational direction and the second moving frame () is aligned orthogonal to and the surface normal vector. This implies that is aligned along the wavefront of the wave, as shown in Fig. 2. Therefore, the second covariant differentiation of along as it moves along , which is mathematically the same as , is equivalent to the Riemann curvature of the orthonormal basis representing the convergence or divergence of trajectories such as Misner

where we used the notation that . This is coincident with the analysis of relative acceleration that was used for the reentry around the pulmonary vein MyBIOP.

(a) Moving frames
(b)
(c)
Figure 4: Propagation of the cardiac electric propagation through a small gap.

5 Aligning along the diffusion propagation

The diffusion-reaction model is one popular mathematical model of the cardiac action potential propagation for the simulation and analysis of cardiac action potential propagation on an anisotropic atrium . Since the original FitzHugh-Nagumo model in the 1960’s FitzHugh as the macroscopic version of the famous Hodgkin-Huxley model Hodgkin, several accurate models have been proposed. The two variable Aliev-Panfilov model Aliev and twenty-one variable Courtemanche-Ramirez-Nattel (CRN) model Coutemanche are two examples. They exhibit similar wave patterns, called the action potentials, as shown in Fig. 5.

In the propagation, the diffusion-reaction process is often represented by the wavefront of the action potential Zykov. The trajectory in this process is often considered the time series of particle traces in the diffusion process Tejedor. However, its counterpart may not correspond to the wavefront on the same scale. In this paper, the trajectory in the diffusion-reaction process is defined as follows.

Definition 1.

The tangent unit vector of the trajectory for a diffusion-reaction process is the same as the normalized gradient of the action potential.

In an isotropic medium, the gradient of the action potential is the same as the orthogonal to the wavefront. Thus, the above definition implies that the trajectory is orthogonal to the wavefront. However, in an anisotropic medium, the gradient of the action potential is not necessarily orthogonal to the wavefront but depends on the charge distribution of the action potential, and its gradient determines the advancement of propagation. This definition represents the exact mechanism of a diffusion-reaction process but may have additional requirements in computational modeling due to discretization errors. This is called the validity of the direction, as will be explained in the next subsection.

(a) AP model
(b) CRN model
Figure 5: AP model: is the solid line, is the dashed line. CRN model: is the solid line, is the dashed line. The wave travels from left to right.

5.1 Validity

For the validity of the direction of the gradient, the independence of the coordinate systems is introduced for the solution of partial differential equations. Consider the following formulation of the diffusion-reaction equation:

(11)

where is the diffusivity tensor and is the reaction function. Then, the constructed moving frames are the local reference frames, for which the solution obtained using one set of moving frames should not be swapped with the solution using another set of moving frames. This property is defined as the validity of moving frames, described as follows.

Definition 2.

Let be the initially constructed moving frames that yield the stable and efficient numerical solution of Eq. (11). Then, the newly adapted moving frames are called valid if they satisfy

(12)

for a sufficiently small tolerance or, approximately

(13)

where the tolerance is a function of and , i.e., and is independent of .

The above condition implies that if the difference of with the new moving frames is sufficiently lower than the discretization error, then adapting the new moving frames does not affect the numerical accuracy and stability during time advances. Therefore, the corresponding moving frames can be adapted as a valid connection because the direction of should not affect the computed value of as long as they are orthogonal and lie on the tangent plane. This should always be geometrically true and inscribed in the equation. On the other hand, any invalid moving frames quickly yield unstable and inaccurate computations.

In actual computations, the condition Eq. (13) discards any non-differentiable directions of the propagation within an element due to numerical noise or discretization error. For example, suppose that the direction of the gradient at a point is abruptly changed from a smooth distribution of the directions at the other points in an element. Then, the non-differentiable direction at the point is not likely constructed from the diffusion-reaction process because of its dissipative and smooth property. However, the condition Eq. (13) is also not satisfied when the curvature of the flow is large, and the spatial resolution is not sufficient in this region. This occurs because the current numerical scheme to compute the covariant derivative requires more grid points in more highly curved regions. Therefore, the parameter must be adjusted to include all the curvatures of the propagation while excluding the gradients of non-realistic propagation.

5.2 Position of the gradient

In constructing the trajectory, the unique shape of the wave should be considered. Let wavefront be the area where the membrane potential () increases and waveback be the area where decreases. The AP and CRN models share the generic features of action potential, such as stiffer wavefront, as shown in Fig. (5), but the increasing or decreasing rates of the wavefront and waveback are different. Thus, integrating the gradient over a time interval should be sign sensitive, depending on whether the gradient is on the wavefront or waveback. The location of time integration is chosen through the following methods.

  1. Wavefront where is positive. is considered as the direction of the tangent vector. (WF)

  2. Waveback where is negative. is considered as the direction of the tangent vector. (WB)

5.3 Extracting the gradient

Moreover, Definition 2 indicates that the trajectory of the diffusion-reaction model is constructed using the unit tangent vector, which is constructed through one of the following methods.

  1. To choose the unit tangent vector when the magnitude of the gradient is the largest (Str).

  2. To integrate all the unit tangent vectors over a time interval when the magnitude of the gradient is non-trivial (Int).

  3. To perform weighted integration over all the tangent vectors during a time interval proportional to the magnitude of the gradient (Wint).

The three methods are compared, and one superior method is chosen. Combining the procedures mentioned above for accurate and stable moving frames for Eq. (11), we obtain an implementable algorithm for aligning moving frames along the propagational direction, as shown in Algorithm 1.

Data: Gradient vector (GV) and moving frames (MF) Result: Construct a connection - Initialization; while Magnitude of GV is not negligible do        - Align MF along GV;        if Newly aligned MF is derived from the following procedure        (1) Validity Check (Eq.(13)),        (2) Choose specific position (WF, WB),        (3) Extracting the gradient (Str, Int, Wint)        then              - Accept the newly aligned MF;                    else              - Retrieve the original MF;                     end if        end while Algorithm 1 Align moving frames along the gradient vector

6 Numerical scheme for diffusion-reaction equation

To apply the alignment scheme along the propagation direction, as explained in the previous section, a stable and accurate numerical scheme to solve diffusion-reaction equations on curved surfaces is required. The numerical scheme to solve Eqs. (17) and (18), called the MMF-diffusion

scheme, on an anisotropic and curved surface is achieved in the context of discontinuous Galerkin methods and is implemented at an open-sourced spectral/hp library, called Nektar++

Nektar++. The MMF-diffusion scheme is obtained by expanding the gradient on moving frames as follows.

(14)

Note that the new moving frame is orthogonal, but it is not necessarily of unit length depending on the diffusion coefficient, particularly in the presence of an anisotropic medium such as a cardiac fiber. Then, we have the following MMF-diffusion scheme MMF2

(15)
(16)

where the tilde sign indicates that the corresponding quantity is the numerical flux at the interface of elements. Owing to the strong restriction of time marching interval due to spatial discretization, an implicit time marching scheme is adapted by using the Helmholtz solver CGHDG.

(a) T=
(b) T=
Figure 6: Distribution of the membrane potential (contour) and aligned moving frames (arrow) from the point initialization at in the plane. Mesh of and . WaveBack. Gradient = Str. The domain is .
(a) Moving frames ()
(b) Connection ()
(c) Curvature ()
Figure 7: Error -convergence of moving frames, connection, and Riemann curvature for the wave from the point-initiation at in the plane. = . Str = Square, Int = Triangle, Wint = Circle. WB = solid line, WF = Dashed line. Measured at and for WB and WF, respectively. Area within the radius of from the initiation point is ignored. .

7 Test in the plane

For the first test in the plane, consider the following Aliev-Panfilov model for the cardiac electric propagation.

(17)
(18)

where , , . The variable represents the membrane potential and represents the collective ion channel openness, including the refractory region. All the variables are dimensionless but can be converted into real measure such as , .

Let the wave be point-initialized and propagate in the plane to form a connection, which is the same as the polar coordinate system centered at the initiation point. The moving frames are aligned by Algorithm 1, and the corresponding connection and Riemann curvature tensor of the moving frames are compared to the exact value. Fig. 6 displays the propagation of the point-initialized wave with an initial radius of and strength of and the creation of moving frames aligned along the wave trajectory. The gradient is obtained at the waveback when the magnitude of the gradient is the largest.

h () 4.08 (1946) 2.97 (4260) 2.21 (7664) 1.58 (14074)
WB 0.00484769 0.000724808 0.000153335 1.21913e-05
order - 5.98 5.26 7.55
0.0271494 0.00587723 0.00167592 0.00019338
order - 4.82 4.25 6.44
0.0394072 0.00406901 0.0016778 0.000226476
order - 7.15 3.00 5.97
WF 0.0284854 0.00895313 0.00228055 0.000230573
order - 3.6449 4.6270 6.83
0.177675 0.0853012 0.0257 0.00416253
order - 2.3108 4.0589 5.42
1.12586 0.24467 0.0431835 0.00866062
order - 4.8070 5.8682 4.79
Table 1: Plane Test: Error -convergence for the wave from a point-initiation in the plane. = maximum length of edge. = number of elements. Gradient Computation = Wint. dt=0.01. Measured at and for WB and WF, respectively. The area within the radius of from the initiation point is ignored. .
h () 2.48 (1152) 1.99 (1152) 1.70 (2304) 1.40 (3072)
WB 0.00148048 0.000783078 0.000272058 7.6965e-05
order - 2.9160 6.7042 6.4565
0.00589095 0.00319238 0.00109853 0.000619431
order - 2.8050 6.7649 2.9296
0.0222421 0.00799031 0.00242781 0.000995709
order - 4.6873 7.5540 4.5575
WF 0.0115381 0.00717967 0.00169167 0.000454362
order - 2.1721 9.1666 6.7217
0.0588998 0.0234089 0.00466264 0.00199622
order - 4.2247 10.2319 4.3377
0.492037 0.329433 0.0301431 0.0114059
order - 1.8368 15.1648 4.9691
Table 2: Sphere test: Error -convergence for a point-initialized wave from the north pole on the sphere. = . = maximum length of the edge. = number of elements. Gradient Computation = Wint. dt=0.01. Measured at and for WB and WF, respectively. Area within the radius of from the initiation point is ignored. .

For this wave propagation, the propagational direction of the wave is the same as the radial direction and is expressed as follows.

where . The connection should have the following values

and the Riemann curvature tensor of orthonormal bases has the following non-trivial component.

The Riemann curvature tensor is the same as the above tensor because the component of the metric tensor for the axis is the same along the axis. The gradient is measured at the wavefront when and at the waveback when , respectively.

Fig. 7 displays the exponential -convergence of moving frames (left), connection (center), and Riemann curvature (right). Three different methods are used for computing the propagation direction from the gradient: Str (to choose the gradient when the gradient is the largest), Int (to average the integration of the gradients which is larger than a threshold value), Wint (to perform a weighted integration proportion to the magnitude of the gradient). The result shows that Wint shows slightly better accuracy, however, overall, the difference is very small for this test case. Moreover, the accuracy at the waveback (solid line) is considered better than the wavefront (dashed line), possibly because the waveback is much smoother than the wavefront. In the wavefront, the membrane potential increases dramatically; thus, the differentiation of the insufficiently resolved membrane potential may not yield an accurate derivative of the membrane potential.

Table 1 shows the exponential -convergence of error for moving frames (), connection (), and Riemann curvature () by Wint. The convergence at = on meshes with various confirms the convergence order of for the waveback and , and for the wavefront. The exact convergence order should be , , . The convergence order of connection is one order less than the order of moving frames because the computation of connection requires additional differentiation in space. Similar order reduction from connection occurs for the Riemann curvature.

(a) T=
(b) T=
Figure 8: Distribution of the membrane potential () and moving frames for a point-initialized wave at the north pole. The mesh of and is used. WaveBack. Gradient = Wint. The domain is the spherical shell with a radius of .
(a) Moving frames ()
(b) Connection ()
(c) Curvature ()
Figure 9: Error -convergence for the wave from a point-initiation on the sphere. = . Str = Square, Int = Triangle, Wint = Circle. WB = Solid line, WF = Dashed line. Measured at and for WB and WF, respectively, when the moving frames of of the area are aligned along the propagational direction. The area within the radius of from the initiation point is ignored. .

8 Test on the sphere

The second test is to initiate a point-excitation at the north pole on the sphere of radius . Fig. 8 displays the propagation of the wave from the north pole and the alignment of the moving frames by computing the gradient in the waveback. Then, the sphere is parametrized as follows.

The moving frame should be aligned along the latitudinal direction of the sphere, or , as follows.

Note that is a unit vector. The tangent vector of the axis has the magnitude of . The corresponding connection should have the following values.

The Riemann curvature of orthonormal basis has the following non-trivial component.

Note that is not the same as , that is, equal to . Fig. 9 shows a similar exponential convergence of moving frames (), connection (), and Riemann curvature () as the exponential convergence for the plane. For error convergence, WB shows much better accuracy than WF. Out of the three gradient computation methods, Str shows the worst accuracy and Wint shows the best accuracy. As shown in Fig. (b)b, error stagnation for the connection component for is possibly due to the geometric approximation error of the spherical mesh GAerror as it appears again at a finer mesh for -convergence. However, this error stagnation in does not affect exponential convergence in as shown in Fig. (c)c.

Table 2 demonstrates the -convergence for this test problem on the sphere. Similar to the test case in the plane, , , and shows , , and order or higher, respectively. As observed in the -convergence, a low order of between and by WB is possibly due to the geometric approximation error of the spherical mesh. The order of between and by WF is , which is close to the exact value of .

(a)
(b)
Figure 10: Domain = [-20,20] [-20,20]. Anisotropy is aligned along the circular arc centered at (a) and (a). The plane wave is initiated from the left wall at .

9 Test in anisotropy medium

In an isotropic medium, the direction of propagation depends on many factors, especially on the initiation type and point. Even with the same geometry and wave type, the change of initiation point would lead to a different direction of propagation and consequently, a different connection map. However, in the anisotropic medium where one direction is dominantly faster than the other direction in terms of propagation speed, the location of initiation is considered much less significant to the connection map of the propagation. Let us verify this.

In the weak formulation of the MMF-diffusion scheme in the anisotropy medium, the moving frames are aligned along the direction of the anisotropy, and the length of the moving frame is adjusted according to the strength of the anisotropy. Then, the coefficient of the gradient in the mixed formulation in Eq. (16) is modified with the new expression of Eq. (14). A non-uniform length of moving frames yields the same effect as multiplying the diffusivity tensor to the gradient.

(a) Initialized at
(b) Initialized at
Figure 11: Domain = [-20,20] [-20,20]. Point-initialized. Anisotropy is aligned along within the circle.

Fig. 10 presents the map of aligned moving frames along the propagational direction on two different anisotropies within the center circle in the domain of [-20,20] [-20,20]: For the left plot, the anisotropy is aligned along the circular arc centered at , i.e., the anisotropy is on the line of . For the right plot, the anisotropy is aligned along the circular arc centered at , i.e., the anisotropy is on the line of . The plane wave from the left wall follows the direction of the anisotropy within the circle. Observe the coincidence (blue color) of the propagational direction and anisotropy within the circle in Fig. 10.

The main difference between cardiac electric propagation and electromagnetic wave is that cardiac electric propagation does not always follow the direction of anisotropy. Fig. 11 displays the propagation after the point initialization near the circle where the anisotropy is aligned along the line . In the left plot, the wave initiates in the low left corner, i.e., at . In the right plot, the wave initiates in the top left corner, i.e., at . To follow the anisotropy within the circle, the propagation should take a very sharp turn or the propagational velocity along the orthogonal to the fiber should be almost zero. However, in some areas within the circle, the wave fails to follow the anisotropy. In Fig. 11, non-blue color indicates a non-trivial angle between the anisotropy and the direction of propagation. The deviation is more significant for the initiation, which started from the top-left corner because anisotropy is curved toward the bottom of the domain. The directional difference between the propagational direction and anisotropy could be the key feature in analyzing and predicting the abnormal propagation of the cardiac electric flow in the heart.

(a) Original fiber map
(b) Projected fiber map
Figure 12: In (a), the contour represents the magnitude of the cardiac fiber. In (b), the contour represents for the angle between the cardiac fiber and the surface normal vector ().

10 Applications to Atlas of atrium

This section describes how the proposed algorithm can be applied to the generation of an Atlas of a 2D atrium, or the connection map of the cardiac electric signal propagation. The most significant benefit of the Atlas is to show the spatial-temporal map of the divergence and convergence of the propagation, possibly to show the location where the propagation stops according to Hypothesis 1. In contrast to the contour plots of the wavefront, this Atlas is unique and valid in the presence of anisotropy; thus, it can be used to analyze the pattern of propagation in a complex cardiac tissue.

In the following sections, we demonstrate the several possible uses of the Atlas for the analysis of the cardiac electric signal propagation in the atrium. For this purpose, the three-dimensional shape and cardiac fiber of the atrium are constructed from magnetic resonance angiography images and histological examination of ex-vivo atria, respectively Cantwell, as shown in Fig. (a)a. These mesh and fiber data are not the actual data of the human heart but serves the purpose as an exemplary use of the proposed scheme to the realistic atrium.

10.1 Preprocessing

The length of the cardiac fiber representing the ratio of the strength of conductivity is not constant, but slightly variable between and . The cardiac tissue of the atrium is anatomically known to be relatively thin. It can be modeled as a curved surface, but still has the features of a three-dimensional structure. Thus, the direction of the cardiac fiber may not be orthogonal to the surface normal vector. In other words, if the cardiac fiber is represented as a vector , then in a two-dimensional model of atrium, the vector should be orthogonal to the surface normal vector such as . For this reason, some of the original fiber that has an oblique angle with is projected on the tangent plane, but the conductivity velocity was also adjusted by the factor of because the velocity parallel to the tangent plane is only considered in the 2D model such as

Fig. (b)b displays the projected fiber (arrow) and (contour) between the 3D cardiac fiber and the surface normal vector. To simplify notation, let mean in the remainder of this paper.

11 Analysis of 2D cardiac fiber

The following analysis of a 2D cardiac fiber yields the two fundamental features of the cardiac fiber on controlling the propagation: The first is to block the propagation, and the second is to remove the sources of self-initiation.

11.1 To block the propagation

Applying Eq. (5) to the obtained cardiac fiber reveals the connection map and the Riemann curvature of orthonormal bases, as shown in Fig. (13). Fig. (a)a shows the distribution of , which is overall relatively low, i.e., below , in most of the atrium except the entrance of each vein. Observe that the region with a large magnitude of is approximately overlapped with the region with a large positive as shown in Fig. (b)b. A positive means a diverging fiber, and a negative means a converging fiber. Thus, only the red area is related to the possible stopping area for the cardiac electric propagation. If the cardiac electric signal propagates along the cardiac fiber, then the propagation may stop or dramatically slowdown in the red area. However, if the cardiac signal flows in a different direction than the fiber, then this connection map does not tell anything for certainty, and it is possible that the electric flow passes through the red area in another direction.

Fig. 14 presents the magnified view of the cardiac fiber around the vein. The red and blue highlights roughly correspond to the region with high and . Fig. (b)b shows that the red area for is related to the abrupt change of the cardiac fiber around the vein. The abrupt change of microscopic fiber is anatomically observed at the canine heart Hamabe, sheep posterior left atrium Klos, and the human heart Tan. The connection map confirms the causal relationship between the anatomical structure and propagational pattern for a unidirectional block around the pulmonary vein to initiate the atrial reentry MyBIOP; Nattel2008; Quan1990.

(a)
(b)
Figure 13: Projected cardiac fiber: Distribution of the connection component (a) and (b) .
(a)
(b)
Figure 14: Magnified view for cardiac fiber with high curvature: (a) region of interest, (b) magnified around the vein.

11.2 To remove the sources

Other interesting geometric properties of the cardiac fiber can be better understood by computing the convariant divergence and curl of the cardiac fiber, which can be computed as follows MMFCovariant.

where is the Christoffel symbol of the second kind and is computed using the connection of the connection form because of the equality for orthonormal basis .

Fig. (a)a shows a similar distribution of the divergence of the cardiac fiber except around the veins with a large magnitude of . With regard to electromagnetics, the cardiac electric flow is related to the flow of charged ions to induce the electric field () along the fiber. The divergence of the fiber, or the divergence of the electric field if the electric propagation completely follows the fiber, means the amount of charge density () as shown by Gauss’s law as

Thus, the negligible divergence of the cardiac fiber eliminates the possibility of self-excitation due to a high charge density by the fiber structure. Because many cardiac electric disorders are initiated by self-excitation such as a spiral wave, it is obvious that the cardiac fiber of a healthy heart should avoid a structure with a high covariant divergence.

Another interesting geometric feature arises concerning suppressing the magnetic field under the structure of the cardiac fiber. According to Maxwell’s equations, the magnetic field () is generated by the non-negligible curl of the electric field by Faraday’s law such as

The above equation implies that if is curl-free, the magnetic field remains zero if no magnetic field is initially present. Fig. (b)b displays the negligible amount of that is approximately the same as if the electric propagation completely follows the fiber. The suppression of the magnetic field in the propagation of the cardiac electric flow is crucial because the magnetic field could change the propagational direction of charged ions to deviate from the cardiac fiber Jackson. Observe that the area with a strong magnetic field is again approximately overlaps with the area with a high . This means that an abnormal cardiac electric flow, such as a spiral wave, could generate a relatively significant amount of magnetic field.

(a)
(b)
Figure 15: (a) Divergence of the fiber. (b) Projection of the curl of the fiber onto the surface normal direction.

12 Connection map by trajectory tracing

To derive the connection map of the cardiac electric signal propagation in the obtained atrium, the twenty-one variable Courtemanche-Ramirez-Nattel (CRN) model is solved on the mesh of the atrium with a initialization from a point close to the sinoatrial node (SAN). The CRN model is known to be more accurate in the atrium than the Aliev-Panfilov model. However, similar results would be obtained by the two-variable Aliev-Panfilov model. The CRN model is described as follows Coutemanche: Consider the following diffusion-reaction equation.

where is the total ionic current which is given as

where is the fast inward current, is the inward rectifier current, is the transient outward current, is the ultrarapid delayed rectifier current, is the rapid delayed rectifier current, is the slow delayed rectifier current, is the L-type inward current, is the sarcoplasmic pump current, is the - pump current, is the / exchanger current, is the background current, and is the background current. For exact parameter values, refer to Coutemanche.

For numerical simulation, we used the following specifications: the order of solution polynomial is with regular elements. The diffusion operator is marched in time implicitly by the IMEX third-order time marching scheme. The weak formulation is adapted by the MMF-diffusion scheme in the context of discontinuous Galerkin methods MMF2. The gradient of propagation is computed at WB by Wint. The time step is set as , and the final time is . is given as .

12.1 What is the role of the fiber?

The first simulation is to run the CRP model on the atrium, but without the fiber to reveal the role of the cardiac fiber in the cardiac electric signal propagation. Fig. (a)a displays the aligned moving frames along the propagation direction in the isotropic cardiac tissue. The heart without the cardiac fiber would produce a significantly decreased mechanical efficiency in pumping the oxygenated blood without twisting and squeezing. Electrophysiologically speaking, Fig. (b)b and (c)c reveal the presence of strong curvature (1) around the initiation point and (2) in the middle of the veins. Strong and consequently strong around the initiation point requires a larger magnitude or larger radius of initial excitation to overcome the initial strong curvature restraint. Moreover, the strong in the middle of the veins could be closely linked to atrial reentry because a small change in the direction to the veins can stop or allow the propagation into the vein. The magnified view in Fig. 17 confirms a relatively strong distribution of around the veins (red color) in isotropic atrium.

(a) Aligned moving frames
(b)
(c)
Figure 16: Without the fiber: (a) Aligned moving frames along the propagational direction. Distribution of (b) the connection form and (c) .
(a) Front view
(b) Back view
Figure 17: Magnification of the region of the veins where is relatively high. A corresponding high curvature occurs due to the indented shape of the vein.

Adding the anisotropic cardiac fiber into the mesh of the atrium dramatically changes the trajectory of the cardiac electric propagation. Fig. (a)a displays the propagational direction that is approximately the same as the fiber map of Fig. (b)b. A quantitative comparison between the propagational direction and the fiber will be shown later. First, the strong curvature around the initiation point disappears due to anisotropy. The initial excitation by SAN could be less intensive than the isotropic medium. Second, there is a strong in a larger area in the veins, including the roots. This broader and stronger at the entrance and the veins prevent the cardiac electric signal from propagating into the veins to prevent atrial reentry and subsequent fibrillation. This new high-curvature distribution can significantly reduce the risk of generating a unidirectional block for atrial reentry by avoiding angle-dependent pathways of the veins. Fig. 19 shows the stronger Riemann curvature throughout the veins and the roots of the veins. The presence of the cardiac fiber dramatically changes the curvature, especially around the veins. However, the critical roles in the prevention of atrial reentry should be studied more closely in the future.

(a) Aligned moving frames
(b)
(c)
Figure 18: With the fiber: (a) Aligned moving frames along the propagational direction. Distribution of (b) the connection form and (c) .
(a)
(b)
Figure 19: Magnification of the region of the veins where is relatively high.

12.2 Cardiac fiber validation

This subsection briefly demonstrates how cardiac fiber data can be validated by the proposed scheme. The cardiac fiber is optimally designed to deliver the electric signal along the fiber direction efficiently. If the cardiac fiber is not along the propagational direction, then many factors are to be considered. A few of those have been mentioned below.

  1. The data of the obtained cardiac fiber contain noise.

  2. The electric propagation models or their parameters for the atrium are not accurate.

  3. The atrium of the patient is in an abnormal condition.

Fig. 20 presents the match between the cardiac fiber and the propagational direction in the atrium to display the distribution of where is the angle between the cardiac fiber and the propagational direction. This figure presents the sufficient coincidence of the two directions except for some regions, especially around the veins. As the anisotropic strength increases as (Fig. (a)a) , (Fig. (b)b), and (Fig. (c)c), the cardiac electric propagation follows the cardiac fiber more closely. This computational scheme can be used to analyze the personalized atrium fiber and shape for clinical studies and prevention of atrial fibrillation.

(a)
(b)
(c)
Figure 20: Distribution of where is the angle between the cardiac fiber and the propagational direction. The contour value of (blue) means parallel (agreement), (red) means orthogonal (disagreement).

12.3 Effect of fibrosis

Another important application of this computational scheme is the analysis of cardiac electrical signal propagation in the atrium with fibrosis, or even scar tissue. From the point of views of electrophysiology, a fibrosis corresponds to low or even zero conductivity in the cardiac tissue. Fig. 21 presents an artificially generated map of the various conductivities from to . A scale implies normal conductivity and means no conductivity. We use this conductivity map to simulate cardiac electrical propagations both with and without the cardiac fibers.

Figure 21: Distribution of the conductivity of cardiac tissue. The scale of corresponds to the normal conductivity and corresponds to no conductivity.

Fig. 22 displays the aligned moving frame along the propagational direction and its corresponding . In comparison with Fig. 16, where the conductivity is all , we observe a notable change in the propagational direction as well as the corresponding and . This implies that the propagation in an isotropic atrium is very sensitive to fibrosis, thereby yielding a considerable change in the propagation and, consequently, its excitation sequence for mechanical pumping. On the other hand, the same fibrosis with the cardiac fiber does not experience a much different propagational direction or the corresponding and . This is clearly observed in Fig. 23, which is only slightly different from Fig. 18.

(a) Moving frames
(b)
(c)
Figure 22: Atrial Isotropy: Propagation in the presence of scar tissue: (a) aligned moving frames (b) connection component , and (c) .
(a) Moving frames
(b)
(c)
Figure 23: Atrial anisotropy: Propagation in the presence of scar tissue: (a) aligned moving frames (b) connection component , and (c) .

13 Discussion

The connection map provides an efficient quantification of the stopping conditions of the cardiac electrical signal propagation. This map includes both spatial and temporal components when the propagation is steady. Thus, all the features of the time-dependent propagation can be viewed in a single plot. The aligned moving frames can be obtained from computational simulations, equivalent mapping, or clinical mapping data. However, the analysis of this map could be more multidimensional and contain geometric information regarding propagation. Therefore, we expect that this analysis would complement the current analysis techniques of cardiac electrical signal propagation in the map format, such as the kinematic analysis (wavefront analysis) or phase map analysis.

The first drawback of this method is that the construction of the aligned moving frames is achieved by solving the two-dimensional diffusion-reaction equations for the atrium. Thus, it is computationally expensive and not available in clinic timeline. If the numerical solution of the diffusion-reaction equations could be obtained from a set of solutions such as the inertial manifold Mallet-Paret, then the construction time could be significantly reduced. Second, certain parameters need to be fitted for every simulation, such as and the minimum magnitude of the gradient, which may affect the final connection map. In the future works, these parameters could be entirely eliminated or be guided for each case individually. Third, the computation of the Riemann curvature tensor is the third derivative of the aligned moving frames, which requires the sufficient smoothness of the geometry and cardiac fiber data.

The future work would therefore involve the (1) use of the anatomically-accurate atrial and ventricular data, (2) study of the representations of spiral waves in terms of the connection form and Riemann curvature tensor, and (3) finding the differences in the connection map for cardiac restitution and cardiac memory.

Acknowledgements

The authors thank Professor Eun-Jae Park (Dept. of Computational Science and Engineering, Yonsei University) for inspirational and encouraging discussion. This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) and funded by the Ministry of Education, Science and Technology (No. 2016R1D1A1A02937255).

References