Conduction in the Heart Wall: Helicoidal Fibers Minimize Diffusion Bias

The mammalian heart must function as an efficient pump while simultaneously conducting electrical signals to drive the contraction process. In the ventricles, electrical activation begins at the insertion points of the Purkinje network in the endocardium. How does the diffusion component of the subsequent excitation wave propagate from the endocardium in a healthy heart wall without creating directional biases? We show that this is a consequence of the particular geometric organization of myocytes in the heart wall. Using a generalized helicoid to model fiber orientation, we treat the myocardium as a curved space via Riemannian geometry, and then use stochastic calculus to model local signal diffusion. Our analysis shows that the helicoidal arrangement of myocytes minimizes the directional biases that could lead to aberrant propagation, thereby explaining how electrophysiological principles are consistent with local measurements of cardiac fiber geometry. We discuss our results in the context of the need to balance electrical and mechanical requirements for heart function.

. Left: a rat left ventricle, with a slice of fiber data from diffusion MRI (color denotes transmural position). Middle: displays the fiber field (from the white cube of the left inset), as well as a local GHM fit (shown as the larger, thicker curves). Right: shows the diffusion behavior in the planes of corresponding border color from the central inset (where columns represent times t = 0.02, 0.3, 1.0) using the GHM parameters from the local fit (see Fig. 5 for further details).  54,55 . Left: v f = 3, v t = 1. Right: v f = 1, v t = 1. The line segments on the left represent the local fiber directions along the z-axis. Notice that the anisotropic case (left) leads to spatially skewed diffusion, in comparison with the behavior of the isotropic case (right). This is particularly noticeable along the x axis, which coincides with the fiber orientation at the origin, where each trajectory begins.
SCIENtIfIC RePoRTS | (2018) 8:7165 | DOI: 10.1038/s41598-018-25334-7 the z-axis to lie along the transmural direction, and define an orientation (or fiber angle) function θ(x, y, z), which describes the angle of rotation of the local myofiber at the given position.
The full expression for the orientation function of the GHM may be written as 4 : where k B , k T , k N ∈  are constants (see Fig. 3 for an illustration). The GHM is a local model, which assumes the fibers have no component out of the local tangent plane to the wall (i.e. no imbrication angle). It describes a series of "planes" of fibers along the transmural axis, corresponding to the local heart wall. Thus, the orientation function above allows defining the GHM as a vector field 33 v : 3 x y z x y z ( , , ) (cos( ( , , )), sin( ( , , )), 0) (2) which is illustrated for several parameter values in Fig. 4 using 3D streamlines. Recent works using diffusion imaging suggest that |k N | and |k T | are quite small across species, whereas k B is always distinctly non-zero 4,33 . In the present article, therefore, we also analyze a GHM where the in-plane curvatures vanish (i.e. k T = k N = 0). Stochastic Calculus on the Cardiac Manifold. Previous work by Young and Panfilov 17 has modeled diffusion propagation of the electrophysiological signal responsible for cardiomyocyte contraction using Riemannian geometry. The essence of this approach is to consider distance measures in the heart wall to be warped by the local anisotropy of the fibers.
Let  g ( , ) 3 be a Riemannian manifold. Consider the following moving frame field to be a basis for the tangent space:ê e e cos( ) sin( ) 0 sin( ) cos( ) 0 0 Then, in this dynamic frame, the metric tensor is given by: which is the inverse of the diffusivity tensor D ∼ in local coordinates. In keeping our notation consistent with prior work 17 , we denote the diffusivity entries along the local fiber direction and in the plane orthogonal to the local fiber as v f 2 and v t 2 , since they are proportional to the squares of the conduction velocities (e.g. 34 ). We further note that previous work utilized an orthotropic metric tensor, with a third diffusion speed constant accounting for the propagation rate within the laminar sheets of the heart 17 . Here, we have simplified our analysis by considering only one speed in directions transverse to the fiber. Empirical measurements 35 showing that the difference between v f and transverse propagation speeds were much larger than the difference in speeds within the transverse plane suggest that this simplification is reasonable.
Intuitively, the metric tensor describes infinitesimal lengths and distances on a manifold. In this case, g describes the curved space in which the signal moves, where the warping is caused by the presence of the fibers. In other words, instead of explicitly modeling the fibers, we model space itself as being intrinsically curved, such that the same movement in different directions leads to different distances being covered (with movement faster along fibers than orthogonal to them).
On the cardiac manifold, it is possible to describe diffusion as a continuous space-time stochastic process using the Ito calculus, which generalizes classical calculus to handle random variables. This is in contrast to the standard method of analyzing diffusion in the heart via partial differential equations (PDEs), as it allows probabilistic interpretations of propagation behavior that deterministic approaches do not. Informally, the PDEs approach describes the macroscopic behavior of diffusion, while the stochastic representation describes the nature of a single infinitesimal member of the ensemble that combines to generate the macroscopic average behavior. Another advantage is that the stochastic description of diffusion can permit intuitive and analytic understanding, whereas the PDEs model is generally analytically intractable.
First, the Laplace-Beltrami operator on a Riemannian manifold (the analogue of the Laplacian in Euclidean space) can be defined in local coordinates by where g ij are elements of the inverse metric tensor g −1 . Note that the Einstein summation convention is used to sum over repeated indices. Then, diffusion on a Riemannian manifold is a stochastic process with an infinitesimal generator given by = ∆ /2 g  36,37 . Further, define the drift coefficient = Γ b g i j k jk i and diffusion coefficient where Γ jk i are the Christoffel symbols. Then the following system of stochastic differential equations (SDEs) describes a diffusion process on the manifold 37 . The Ito process defined by this system of SDEs describes the behavior of diffusion on the cardiac manifold.
Diffusion on the Cardiac Manifold. Classically, the propagation of the contraction signal in the heart can be modeled as a reaction-diffusion equation, with a spatially varying diffusion tensor D (with components D ij ) depending on the local fiber orientation 27,28 , via where Ψ → , Φ are non-linear reaction functions describing the cardiomyocyte electrophysiology, while u and v → track cellular activation state variables, and the first term is generated by the standard anisotropic diffusion gen- Such studies tend to convert between the diffusion PDE and Riemannian geometric representations by defining the metric to be the inverse diffusivity tensor, with components D ij . However, from a stochastic perspective, it can be shown that a more natural approach is to consider the metric tensor to be the adjugate matrix of the diffusion tensor 39,40 , as it expresses isotropic diffusion in the curved space as anisotropic diffusion in Euclidean space. There is a deep duality between the heat (isotropic diffusion) equation and BM on the cardiac manifold: the heat kernel, i.e. the fundamental solution to the heat equation ( , is exactly the transition density function of BM on the manifold 38,41 . More intuitively, we can formalize the notion that the heat equation describes the expected behavior of an ensemble of BM processes by noting that solves the manifold heat equation with initial conditions u(0, x) = f(x), where X t,x is a BM on the manifold starting from x 37,42 . Thus, analysis of the Ito diffusion process above corresponds to understanding the heat equation on the manifold, which comprises the diffusion term in the reaction-diffusion equation above. Figure 5 illustrates how the curved GHM space affects diffusion from a point in a helicoidal medium, such as from an excitation site of the Purkinje network in the endocardium 2 .  Herein, we focus on the behavior of diffusion on the manifold from a stochastic perspective, to understand how the GHM affects contraction signal propagation.

Results
The GHM Manifold Possesses Negative Ricci Curvature. As noted above, the cardiac tissue can be considered a Riemannian manifold with a metric tensor defining spatial distances on the manifold g, warped by the presence of the fibers. For any such manifold, the curvature of the space can be measured via the Ricci curvature scalar R, which in the case of the GHM is given by (see Supplemental Information (Derivations): Ricci Curvature of the Cardiac Manifold) in agreement with previous calculations by Young and Panfilov 17 . Since the scalar curvature is spatially constant and always negative, it accelerates diffusion propagation on the manifold. As also noted in previous work 17,43 , this acceleration can be seen in the fact that, for small t, the volume of activated tissue grows as where B t is a metric ball at time t. We note another connection to the speed of diffusion lies in considering the short-time asymptotic expansion of the heat kernel p(t, q 1 , q 2 ), the fundamental solution of the diffusion equation  . For the GHM manifold, it can be shown that the autodiffusion function h t (q) = p(t, q, q) can be expanded as 42,44 when t is small, showing that the local Ricci scalar curvature dominates the local behavior of diffusion on such time scales. Recalling that p is the transition density of the Ito diffusion defining BM on the manifold, we again see the duality between the stochastic and PDE formulations of the signal diffusing. One can thus interpret h t (q) as either the flow of a quantity from a point to itself or the probability of BM staying at q over time. Informally, if the probability of staying near a start point q (i.e. h t (q)) is lower, then the diffusion is flowing away from its start point q more quickly. As such, in the case that R < 0, notice that the curvature term decreases the value of h t (q), meaning that the term is encouraging the diffusion flow to move faster from its start point. Hence, negative R relates directly to faster diffusion from a point.
A natural question is whether the GHM optimizes the Ricci curvature in the variational sense, which would imply a maximization or minimization of the curvature-derived acceleration of the diffusion rate. Using the calculus of variations, we treat the Ricci curvature as a functional on the space of orientation functions. Akin to classical calculus, the first variation shows whether an input is an extremal point in function space, while the second variation classifies its type. We consider two different function spaces: the set of "planar" fiber fields (with no imbrication angle) and the space of general 3D unit vector fields.
Using variational techniques, it can be shown that the GHM lies at a stationary point of the Ricci curvature scalar, in the planar case (see Supplemental Information (Derivations): Variational Analysis of the GHM Ricci Curvature). In contrast, in the general space, only the GHM with k N = k T = 0 has vanishing first variation; the full GHM is not a stationary point for other parameter values. However, the stationary point is neither a minimum nor a maximum, as it fails to satisfy the Legendre condition, which is a necessary requirement for the function to lie at an extremum 45 .
We note that the variational status of the GHM as a stationary point with respect to the Ricci curvature in the full space requires the vanishing of the in-plane curvatures (k N and k T ). Recent imaging studies with GHM fits confirm that these curvatures are small in magnitude 4,33 . Another example of the importance of the vanishing of these curvatures will be shown in the next section, with respect to the drift vector of diffusion on the manifold.
In summary, corroborating previous analyses 17 , the GHM has a constant negative Ricci curvature, which has an accelerating effect on the diffusion process. Our variational analysis of the GHM shows that it always lies at a stationary point in the planar space of fiber angle functions, but it is only stationary in the full space when the in-plane curvatures (k N , k T ) vanish. It is interesting that the GHM lies only at a stationary point, rather than a minimum or maximum, perhaps indicating a balance between competing factors, such as wave speed and stability. We note that previous studies 22,23,26 suggest that increased rate of transmural rotation (determined by k B ) may be linked to cardio-electrophysiological instability. Hence, although increasing |k B | would accelerate the signal propagation, this may not be beneficial to the organism overall.
The Diffusion Drift is Minimized by the Helicoidal Architecture. The Ito diffusion process of BM on a manifold is described by a system of SDEs that is parametrized by a diffusion coefficient σ = For the GHM, the former is given by  Thus, we have shown two results linking the in-plane curvatures of the GHM to the drift vector. The first shows that the drift vector is identically zero across the cardiac manifold generated by the GHM if and only if the in-plane curvatures vanish. The second is a complementary result, showing that the full GHM (i.e. with arbitrary in-plane curvatures) is not generally a variational minimizer of the drift magnitude, unless k N and k T vanish. As such, the empirical findings 4,33 that the curvature parameters k T , k N in the tangential plane to the mammalian heart wall are small relates also to the minimization of the Ito drift. Although there is no deterministic drift (or convection) component to the GHM encountered empirically (i.e. that with vanishing in-plane curvatures, or k N = k T = 0), there is still a stochastic bias in the components of the variance vector so that diffusion in some directions is favored over others. We show that for the GHM, as long as k B ≠ 0, in the heart wall tangent (xy) plane, this bias is removed exponentially quickly. We also show that the case of k B = 0, corresponding to no transmural rotation, is biologically problematic, because then the variance vector has a strong constant bias in the fiber direction x.

In-Plane Variance Becomes Isotropic Exponentially
For the GHM with k N = k T = 0, we obtain the following system of SDEs:  Fig. 2 (left) for sample trajectories obtained by numerically integrating the SDE above. One can see the influence of fibers on the diffusion processes, analogous to the acceleration along fibers seen in Fig. 1 (right).
It can be shown (see Supplemental Information (Derivations): Moments of the GHM Stochastic Diffusion Process) that the variance of the stochastic process described by these SDEs is given by: where = c k v 2 B t 2 2 and that the expected value is /constant at zero. See Fig. 6 for plots of the numerical convergence of the moments of the process. Notice that the first term of the variance vector has identical components in the two in-plane directions x, y in the heart wall tangent plane, with a smaller component in the transmural direction z. The second term has components with identical magnitudes in the fiber direction x and the in-plane direction y, but with the latter having a minus sign. This causes an increase in variance in the x direction with a corresponding decrease in the y direction. Since the first term in the variance equation above has no dependence on k B , it is in fact the second term that relates to the directional bias to the variance vector caused by the transmural rotation of fibers. Now, recall that the variance of a diffusion process can be interpreted as being related to the expected growth (movement) of the diffusion in that direction. For instance, the variance of BM in 1D isotropic Euclidean space ; i.e. linear in time. Hence, we can gain intuition about the behavior of diffusion on the GHM by examining the asymptotics of the variance as we vary some of the parameters of the GHM diffusion model. As a simple case, suppose that v f = v t . Then, as one would expect, , meaning the directional anisotropy in the variance vector disappears. In particular, if v f = v t = 1, we get the classical result for Brownian diffusion in  3 . Next, consider the case where k B → 0, so that there is no transmural rotation of fibers in the heart wall.

It can be shown (see Supplemental Information (Derivations): Moments of the GHM Stochastic Diffusion
, meaning that growth in the variance vector is accelerated in the fiber direction x, whereas it is not in y and z.
We next consider the asymptotic behavior of the terms of the variance vector in time. As t increases, the first variance term grows linearly, while the second converges exponentially fast to a constant. This can be interpreted as an isotropic linear growth of the variance vector in the heart wall tangent plane (xy) in time, with a reshaping to increase the variance of the diffusion process along fibers in the direction x while decreasing it across them in the direction y. However, the effect of this reshaping becomes quickly negligible as time increases. We also observe that this effect is identical to that of increasing k B , the helicoidal fiber rotation rate. Intuitively, a larger k B means greater homogeneity in x and y, converging to a growth averaging the diffusion rates in the fiber and non-fiber directions. Thus: ; i.e. growth expands at the same rate in the x and y directions. Intuitively, the fibers create a local directional bias initially, but this vanishes exponentially fast into an isotropically accelerated propagation pattern in the heart wall tangent plane. We also showed that the bias of the variance vector grows linearly in time when k B → 0, which is the case where all fibers lie along the x-axis (i.e. no transmural turning). There is strong empirical evidence that k B is non-zero (i.e. that there is always transmural turning of fibers) 4,31-33 . Thus, similar to how the empirical minimization of k T and k N can be interpreted as minimization of the drift vector b, the presence of non-zero k B can be viewed as reducing the bias of the variance vector within the heart wall plane.
Our expression for the variance function allows us characterize the timescale of the exponential convergence. We can consider the difference between the x and y components of the variance as the in-plane stochastic bias: × − m 2 /s. As a caveat, we note that measurement of the conductivities is itself a difficult problem, and there is variability in the values reported in the literature 47,48 as well as spatially in the heart wall 46 . With these values, we can estimate that c ∞ = 5.19 × 10 −5 m 2 and τ = 0.23 s. This result suggests that, while the shrinking of the in-plane stochastic bias is exponential, it seems plausible that at least some effect on the wave front due to the anisotropy can persist for a non-trivial fraction of the heartbeat. Uncertainty in the measurement of diffusivities as well as their spatial variation preclude surety in our conclusion. After approximately τ = 0.23 seconds have passed, the ratio of the magnitude of the x (or y) component of the first term in the variance equation to the magnitude of the x (or y) component of its second term, is 2.63, meaning that the first (isotropic) term dominates.
In summary, based on our results concerning the drift vector, we restricted our attention to the GHM with vanishing in-plane curvatures, and considered the variance of the diffusion process, which can serve as a measure of its growth rate along each dimension. We showed the stochastic bias of this process in the heart wall tangent plane vanished exponentially quickly, as long as k B ≠ 0. In combination with our results showing the vanishing of the deterministic bias (i.e. Ito drift), the shrinking of the anisotropy of the variance vector in the heart wall tangent plane suggests that the helicoidal fiber geometry found in mammalian hearts 4,33 , where k N , k T ≈ 0 and k B ≠ 0, minimizes the diffusion bias.

Discussion
The contraction of myocytes in the heart wall is controlled by current flowing through an elaborate Purkinje network, which raises the question of how current diffuses through the wall around the insertion points. Classically, the propagation of the contraction signal is represented using a reaction-diffusion equation with the inverse diffusivity tensor as its metric 17,27,28 . We combined this Riemannian approach with a minimal surface model of the cardiac fiber geometry, the generalized helicoid model (GHM) 4 , and exploited the Ito calculus to describe and analyze diffusion on it.
Among our results, we showed that the helicoidal structure of cardiac fibers minimizes local biases in the diffusion, both deterministic (embodied by the Ito drift vector) and stochastic (derived from differing x and y components of the second moment of the Ito diffusion SDEs). Since the GHM is a local approximation anywhere in the left ventricular wall 4 , one would expect it to not favor biased flow in a particular direction. In line with these SCIENtIfIC RePoRTS | (2018) 8:7165 | DOI:10.1038/s41598-018-25334-7 theoretical expectations, we find that the GHM with k N = k T = 0 minimizes the Ito drift, providing a natural interpretation for the recent imaging studies that suggest |k N | and |k T | are relatively small across species 4,33 . Moreover, using the variance of the Ito diffusion on the cardiac manifold with vanishing in-plane curvatures, we also show that stochastic bias in the heart wall tangent plane shrinks exponentially fast in time, becoming isotropic with a rate given by averaging the squares of the speeds along and across fibers.
This solution is also good in a functional analytic sense, for both the drift and Ricci curvature. Using the calculus of variations, we showed that the GHM lies at a stationary point of the Ricci curvature in the space of orientation functions, provided the in-plane curvatures are zero. However, it is neither a minimum nor a maximum. One might speculate that this is because increased negativity of the Ricci scalar, equivalent to increased transmural rotation rate (or |k B | for the GHM) could lead to cardiomyopathic disturbances in the signal wavefront as it propagates 26 , such as scroll wave turbulence 22,23 . On the other hand, there are a number of potential mechanical reasons for the presence of the helicoidal structure, including shear wave filtering 49,50 , energy dissipation 6 , and mechanical efficiency 11,29 . Thus, the competing effects of changing k B must be balanced to allow safe and fast wave propagation, as well as mechanical efficacy. Physical constraints, such as the myocyte size and necessary accompanying extracellular matrix structure, which idealized models often do not capture, also affect the structure.
While we have considered a number of properties conferred upon the heart by its helicoidal geometry, there are a few limitations inherent to our approach. Currently, we consider only local, static fiber geometry within the left ventricle. Our analysis also only considers the diffusion term in the reaction-diffusion contraction wave propagation equation, since this is the term that is affected by myocyte orientation and it largely governs the local spreading of the signal. In addition, we do not consider the effect of cardiomyocyte laminar sheets, which may contribute to both the electrophysiology 17,51 and mechanics 52 of the heart. Our approach provides useful insights, despite these approximations, because the effect of the fibers on diffusion of the contraction signal has been empirically shown to be much greater than that of the sheets 35 . Lastly, the GHM assumes that fibers lie in planes tangent to the heart wall. Although this is an acceptable local approximation 4,33 , some studies suggest that the presence of fibers with an out-of-plane component may play a role in cardiac function 53 . Finally, exploring the link between local diffusion bias and long-term wave propagation in the heart is a promising avenue for future work.
In the end, our analysis reveals that the GHM arrangement minimizes any in-plane directional biases in the diffusion process, provided that there is transmural rotation and small in-plane curvature of the fibers in the heart wall. Normally, one might have thought about the orientation of fibers from a purely mechanical perspective. However, given the ubiquity of helicoidal arrangements in biology and artificial composites [5][6][7][8][9] , it should perhaps not be a surprise that nature has used this geometry to satisfy both the mechanical and electrical requirements of the heart.