Negative Curvature Boundaries as Wave Emitting Sites for the Control of Biological Excitable Media

Philip Bittihn,* Marcel Hörning, and Stefan Luther Max Planck Institute for Dynamics and Self-Organization, Göttingen, Germany Institute for Nonlinear Dynamics, Georg-August-Universität, Göttingen, Germany German Center for Cardiovascular Research (DZHK), partner site Göttingen, and Heart Research Center Göttingen, D-37077 Göttingen, Germany Department of Physics, Graduate School of Science, Kyoto University, Kyoto, Japan Laboratory for Physical Biology, RIKEN Center for Developmental Biology, Kobe, Japan Department of Biomedical Sciences, Cornell University, Ithaca, New York, USA (Received 27 March 2012; published 14 September 2012)

Self-organized dynamics is ubiquitous in extended biological, chemical, and technical systems, and their control remains a major scientific challenge [1].In biological excitable media, such as cardiac and neural tissue, the complexity of spatial-temporal patterns is often associated with pathological states [2,3].In the brain, excitation patterns have been observed, e.g., during spreading depression [4] and in experiments with disinhibited cortical tissue [5].In the heart, complex spatial-temporal dynamics underlies life-threatening arrhythmias [6,7].
Electric fields are commonly used to control biological excitable media.For example, the termination of cardiac fibrillation relies on high-energy electric shocks that depolarize every cell but may have severe side effects such as tissue damage and pain.In contrast, low-energy pacing strategies use weak electric-field pulses to recruit heterogeneities in electrical conductance [8,9] as endogenous control sites and synchronize the heart.Low-energy antifibrillation pacing (LEAP) has been shown to require at least 80% less energy compared to conventional defibrillation in vitro and in vivo in canine dogs [10,11].A similar energy reduction was found in Refs.[12,13].
In order to develop and optimize such low-energy pacing strategies, it is essential to understand the mechanisms that determine the location of electric-field-induced wave emission.The role of tissue heterogeneities [14] has already been subject to a number of studies, which have been focused on various kinds of heterogeneities, including changes in fiber direction with respect to the electric field [15][16][17], nonconducting inclusions in the tissue [8,[18][19][20][21][22], and the geometry and tissue-bath interface conditions [23,24].In this Letter, we provide a generalized theory for the depolarization of biological excitable media near curved boundaries of arbitrary shapes.
When cardiac tissue is stimulated by an electric field, the activations for the lowest field strengths originate at the boundaries of the tissue that are perpendicular to the electric field.This effect can be observed in experiments on cardiac myocyte cultures with artificial boundaries but an otherwise homogeneous cell monolayer.However, the experiments illustrated in Fig. 1 indicate that the field strength necessary to induce excitation is significantly different for tissue boundaries of convex and flat shape.Specific geometrical features thus enable the recruitment of additional wave emitting sites in the low field-strength range-an effect that has been neglected before.
The interaction of an external electric field with a substrate of interconnected cells can be modeled by the bidomain equations [25].These equations describe the dynamics of the extracellular and intracellular potentials e and i in a tissue domain D. For simplicity, we consider isotropic tissue (a zeroth order approximation for the case of unequal anisotropy ratios; cf.[19]), in which case the bidomain equations reduce to a system of partial differential equations for the membrane potential V m ¼ i À e , the so-called monodomain model: @h @t ¼ HðV m ; hÞ: Here, D is the (scalar) diffusion coefficient, hðx; tÞ describes the local ion channel states, g is proportional to the transmembrane currents, and H describes the ion channel dynamics.The appropriate boundary conditions of the original bidomain equations are where n is the local unit vector perpendicular to the boundary @D and the additional potential o is defined outside the tissue domain D. Equation (3) ensures the continuity of currents across the tissue border, Eq. ( 4) ensures the continuity of the corresponding potentials, and Eq. ( 5) is a no-flux boundary condition for the intracellular space.In the isotropic case, we obtain a boundary condition for the membrane potential V m by subtracting Eq. (3) from Eq. ( 5): We assume a homogeneous electric field o ðxÞ ¼ ÀE Á x for the potential outside the tissue and linearize Eq. ( 1) around the resting state.If the function H describing the ion channel dynamics is independent of V m below some finite threshold V thresh m > V resting m , then this linearization corresponds to the passive dynamics of the membrane [15,16,23].Neglecting the conductivity ratio o = e (which only scales the electric field) yields for the deviation e from the resting potential, where is the electrotonic constant and Eq. ( 8) is the boundary condition of Eq. ( 6).In the following, we will solve Eqs. ( 7) and ( 8) on domains D of different prototypical shapes, determining the dependence of the maximum depolarization in D on the curvature of the boundary @D, where @D can be up to two-dimensional.A suitable measure for comparing the curvature of the boundary of two-and three-dimensional tissue domains D is the mean curvature , where 1 and 2 are the two principal curvatures of the surface (with 2 ¼ 0 for two-dimensional D).It should be noted that, with the assumed form of o , Eq. ( 8) is the only place in the derivation which depends on the applied electric field.Thus, monodomain simulations can account for the effect of an applied electric field by replacing the usual no-flux boundary condition with Eq. ( 8).
All numerical simulations carried out in this Letter use the Fenton-Karma model [26].The parameters are d ¼ C m =g fi (with g fi ¼ 8:7 mS and w ¼ 1000 ms, À w ¼ 65 ms, u c ¼0:13, u v ¼0:025, u si c ¼0:85, k¼10, and D ¼ 1 cm 2 =s.can be calculated analytically as Þg (see Supplemental Material [27] for a detailed derivation).The above parameters result in ¼ 1:12 mm.Circular boundaries are analytically tractable, yet important, geometrical shapes, as they represent cross sections of blood vessels in the heart [11].The convex case, not considered by Pumir et al. [8], can be thought of as a cross section through, e.g., papillary or trabecular muscles.These are cylindrical isolated muscle strands connected to the bulk of tissue only at their ends.For these two cases, depicted in Figs.2(a) and 2(b), two analytical solutions e 1 and e 2 of Eqs. ( 7) and ( 8) can be calculated, respectively (see Supplemental Material [27] for details): K 1 and I 1 are modified Bessel functions of the second and first kind, respectively.In Fig. 3, the maximum depolarization for these two cases is plotted against the normalized mean curvature ¼ 1 2 ðR=Þ À1 of the boundary, where positive curvature (gray line) refers to a concave medium and negative curvature (green line) refers to a convex medium.There are two interesting results to note here: First, the nonmonotonicity for negative curvature means that there is an ''optimal'' curvature for which the effect of the electric field is largest.Second, the maximum depolarization can exceed the theoretical limit of jEj, as calculated in Ref. [18], meaning that some boundaries with negative curvature will induce excitation for a smaller electric field strength than all boundaries with positive curvature.The decline of the maximum depolarization towards large negative curvatures can be attributed to the diffusive interaction with the hyperpolarized area at the opposite boundary.The shape shown in Fig. 2(c) eliminates this effect by cutting away the hyperpolarized half of the boundary.In Fig. 3, the maximum depolarization for this shape is plotted against the curvature of the protuberance as a black dashed line.The curve features a maximum and declines towards jEj for !À1, because features much smaller than the space constant are not ''seen'' by the dynamics due to diffusion.Quantitatively, the values of the depolarization occurring for this geometry are always larger than for circular boundaries-in agreement with the absence of the compensatory effect between adjacent regions of de-and hyperpolarization described above.
The full effect of boundary curvature is obtained in the parabolic medium D ¼ fðx; yÞ 2 R 2 jx < ay 2 g, which has mean curvature ¼ a at the origin [compare Fig. 2(d)].The smooth nature of the boundary makes it possible to calculate an exact solution for Eq. ( 7) which satisfies Eq. ( 8) with an accuracy of $10 À10 by composing it as a linear combination of a few terms of the form e i ðxÞ ¼ expð x;i xÞ coshð y;i yÞ with 2 x;i þ 2 y;i ¼ 2 (see Supplemental Material [27] for a detailed derivation).The resulting maximum depolarization is plotted as a red line in Fig. 3.The curve has the same asymptotic behavior close to 0 as the dashed line for the semicircular geometry, but it increases monotonically instead of approaching a smaller limiting value.From the finite part of the curve shown here, we cannot conclude whether there will be a limiting value for large !À1.However, the curve seems to level off towards larger negative curvatures.Thus, the most apparent difference in the reaction of this kind of boundary to electric fields should be noted up to curvatures of about 1= or 2=.Results for two three-dimensional tissue domains are plotted as blue and cyan dashed lines in Fig. 3.For the blue line, D ¼ fðx; y; zÞ 2 R 3 jx < ay 2 g; i.e., D is a threedimensional extension of Fig. 2(d) with zero curvature in the z direction.An exemplary depolarization pattern is given in Fig. 2(e).The two principal curvatures 1 ¼ 2a and 2 ¼ 0 lead to a mean curvature of ¼ a as in the two-dimensional case.The blue curve falls exactly on top of the red curve for the two-dimensional parabolic case, since any solution eðx; yÞ of Eqs. ( 7) and ( 8) on a twodimensional domain D 2d is also a solution on the threedimensional domain D 3d ¼ fðx; y; zÞ 2 R 3 jðx; yÞ 2 D 2d g.This is because @ 2 e=@z 2 and the z components of all normal vectors vanish, reducing Eqs. ( 7) and ( 8) to their two-dimensional originals.In this way, all results presented for two dimensions above are also valid for three-dimensional shapes, with cross sections in every z ¼ const plane corresponding to Figs. 2(a)-2(d).
If the added dimension introduces additional curvature, this can enhance (or diminish) the observed depolarization.The tissue domain depicted in Fig. 2(f) is defined by D ¼ fðx; y; zÞ 2 R 3 jx < aðy 2 þ z 2 Þg, which corresponds to the parabolic shape rotated around the x axis.Its principal curvatures at the origin are k 1 ¼ 2a and k 2 ¼ 2a, resulting in a mean curvature of ¼ 2a.The resulting maximum depolarization at the tip of the paraboloid is shown as a cyan curve in Fig. 3.For larger negative curvatures, the  The mean curvature at the location of maximum depolarization is identical for all six cases and % AE0:34 in units of 1= (cf.Fig. 3).All length units on the axes are given in units of ; the depolarization e Ã is in units of jEj.
PRL 109, 118106 ( 2012) week ending 14 SEPTEMBER 2012 118106-3 additional curved dimension increases the induced depolarization.This is an inherent effect of the threedimensional topology, since all curves are plotted against the mean curvature.In contrast, towards zero curvature, the maximum depolarization asymptotically behaves like the one for two-dimensional tissue domains, meaning that for this regime, the mean curvature alone determines the strength of depolarization.
Figure 3 was obtained by solving Eqs. ( 7) and ( 8), which is a low-field-strength, monodomain approximation.This approximation is validated by using numerical simulations on the complex geometry of the left ventricular wall obtained from microcomputed tomography (micro-CT; Fig. 4).For this purpose, the full, nonlinear, time-dependent Fenton-Karma model was simulated, with Eq. ( 8) instead of the usual no-flux boundary condition.We now use field strengths that are capable of evoking action potentials; i.e., we leave the scope of the linear (passive membrane) approximation introduced at the beginning.The results in Fig. 4 qualitatively reproduce the mechanisms described above: A large protuberance on the endocardium is activated at a lower field strength than the flat part of the endocardium.This finding agrees with the cell culture experiments shown in Fig. 1.Inner boundaries such as blood vessels are being activated at higher field strengths, since they correspond to the regime of > 0 in Fig. 3, for which the depolarization is generally lower.Good agreement with the theoretical results shown in Fig. 3 suggests that other effects which might depend on the geometry such as the liminal area [28] play only a secondary role.
Our investigations have revealed an interplay of different mechanisms governing the response of cardiac tissue to electric-field stimulation: curvature effects, interaction of depolarized and hyperpolarized regions from close opposite boundaries, and the scale of anatomical features with respect to the electrotonic space constant .The findings summarized in Fig. 3 indicate that the depolarization of tissue boundaries can well exceed the value jEj for a straight boundary.Specifically for quasicircular shapes, the results reported here indicate that, if these structures have an optimal radius (in the order of ), they will respond to particularly low electric-field strengths.The nonmonotonic shape of the green curve in Fig. 3 therefore singles out a particular size from the size distribution of papillary, trabecular, or pectinate muscles.The importance of structures that are linked to the bulk of cardiac tissue in distant places for the activation dynamics has already been shown for the Purkinje system, which was reported to have an increased excitability compared to the bulk of tissue due to it being a one-dimensional geometric structure [29].The results shown here now suggest that there are even more structures in the heart that, only due to their geometry, can be excited very easily.More realistic experimental studies on the complex endocardial surface could verify this spatially heterogeneous response to electric fields.A next step is the application of this theory for the recruitment of control sites during ongoing activity in the medium.
In the presence of an electric field, the arborizing tubular network ( > 0) of the coronary vasculature provides an anatomical substrate for an adjustable intramural wave source density, which can be employed for multisite pacing and low-energy termination of fibrillation [11].This Letter confirms the strong correlation of wave emission locations with anatomical features of the substrate and generalizes the underlying theory towards boundaries with arbitrary shapes.With respect to low-energy control of fibrillation, our results indicate that further significant energy reduction may be achieved by control strategies that aim at inducing wave emission from structures with negative curvature, such as papillary, trabecular, and pectinate muscles or protrusions.However, located at the endocardium, they provide control sites at or near the inner surface of the heart only.
Although we have used a cardiac model for the local dynamics in this study, it should be emphasized that the described curvature effects are solely based on the different extracellular and intracellular boundary conditions of Eqs.(3) and ( 5) and the electrically excitable nature of the system.Both facts apply to any substrate composed of cells which are capable of action potentials.Promising applications include nerve fibers and neuronal tissue, which may have a complex shape and are stimulated

FIG. 1 (
FIG.1(color online).Boundary effects in two-dimensional cardiac cell culture experiments (homogeneous cell distribution).Cultured neonatal rat cardiac cells were stimulated with vertical electric fields (white arrows) of (a) 0.91 and (b) 1:37 V=cm, and activity was imaged by using the fluorescent Ca 2þ indicator Fluo-8 (for details, see Supplemental Material[27]).Colors indicate time from pacing to activation; gray color marks regions without tissue.At higher field strengths, the flat parts of the boundary are excited as well as the tip of the V-shaped boundary.For the lower electric field strength, excitation can be induced only at the latter location.A series of N ¼ 9 experiments yielded a field-strength ratio of 1:49 AE 0:26 between the two types of boundaries.

FIG. 3 (
FIG. 3 (color online).Maximum depolarizations occurring for the geometries shown in Fig.2(see the legend).Solid lines refer to analytical or semianalytical solutions of Eqs.(7) and (8), dashed lines denote purely numerical solutions, and the dots are numerical values to confirm the analytical calculations.

FIG. 2 (
FIG. 2 (color online).Stationary membrane potential distributions for infinitesimal electric fields.(a)-(d) Two-dimensional geometries.White color indicates regions that are not part of the active tissue domain D. The arrows indicate the location of maximum depolarization.(e),(f) Three-dimensional geometries.The surfaces indicate the boundary of the excitable media.The mean curvature at the location of maximum depolarization is identical for all six cases and % AE0:34 in units of 1= (cf.Fig.3).All length units on the axes are given in units of ; the depolarization e Ã is in units of jEj.

FIG. 4 (
FIG. 4 (color online).Numerical simulation of the Fenton-Karma model on the complex anatomical geometry of a twodimensional cross section of the left ventricle (rabbit, micro-CT).An electric field (field strength E, duration 5 ms) results in curvature-dependent activation.(a) E ¼ 0:2 V=cm: Wave emission is observed at a protuberance with large negative curvature, while depolarization at the boundary remains below threshold.(b) E ¼ 0:4 V=cm: Excitation of almost the entire tissue boundary results in a quasiplanar wave.(c) E ¼ 1:0 V=cm: Internal boundaries with positive curvature (blood vessels indicated by arrows) emit waves.