Prior to investigating the SI curves for the vessel configurations shown in Figure 1, we first show how the VE patterns change when the low-conductivity vessel wall is included in the bidomain model. This was done by solving the bidomain equations with a passive membrane (I_ion = V_m/R_m) until steady state is achieved. The steady-state VE patterns for a vessel with and without vessel walls are shown in Figure 3. Figure 3 clearly illustrates that the VE pattern is changed by the presence of the low-conductivity vessel wall. The most significant difference is the swapping of the polarity at the upper and lower surfaces of the vessel, an effect observed previously (Bishop et al., 2010a, 2012). A secondary effect of the presence of a vessel wall is to increase the magnitude of VEs which form on either sides of the vessel due to the anisotropic conductivity ratios. The physical reason for these two effects is that: (a) in the case of no vessel wall, boundary VEs from current transiting the extracellular-to-bath-to-extracellular spaces dominate, as current takes the path of least resistance through the vessel cavity, and (b) in the case of a vessel wall, the insulating vessel wall changes the path of least resistance from through the vessel cavity to around the vessel cavity, thus VEs from dissimilar anisotropy ratios dominate, and boundary VEs are diminished.

The different VE patterns shown above in Figure 3 act to change the active membrane response, as shown in the SI curves in Figure 4. Figure 4 shows that for late diastolic intervals (DIs) of greater than ≈315 ms, the SI curve is approximately linear and near the asymptotic, minimum, value. At such late DIs, the mechanism of tissue capture is via make excitation, with wave propagation initiating from the vessel almost immediately as the S₂ stimulus is applied. The shock strength at late DIs for vessels without vessel walls is approximately half that compared to the case with a vessel wall, and the shock strength required to elicit excitation from larger vessels is lower than that required for smaller vessels, consistent with the literature (Pumir and Krinsky, 1999; Hörning et al., 2010; Luther et al., 2012). For all vessels (in Figure 4), the shock strength monotonically increases with decreasing DI. At early DIs, the shock strength required to elicit propagation for the larger vessel with a vessel wall is lower than that without a vessel wall, a direct consequence of the altered VE pattern around the vessel from the influence of the insulating vessel wall. We include a vessel wall in the rest of the results in this manuscript, as the presence of the vessel wall significantly alters the SI curves.

Figure 5 shows the types of propagation patterns created by the VEs around the large (a = 2.0 mm) blood vessel, with a vessel wall, at different DIs and shock strengths causing make and break excitation patterns. Break excitation (right panel in Figure 5) shows how the depolarized regions rapidly diffuse into the now-excitable hyperpolarized regions, causing waves to propagate through these previously hyperpolarized regions. Without the close proximity of de- and hyperpolarized VE regions, the mechanism of break excitation would not occur.

It should be noted that we investigated the SI response (results not shown) in the absence of the electroporation currents (Ashihara and Trayanova, 2004) and found similar behavior as reported in Bray and Roth (1997); the addition of the electroporation currents acts to reduce the shock strength required to elicit propagation from VEs.

The vessels comprising the coronary vasculature are often colocated, with veins following the same tracts as arteries in vessel pairs. The proximity and orientation of two blood vessels alters the resultant VE pattern produced because the fiber architecture is different for two proximal vessels, compared to that for one vessel, and also because the potential fields in the intracellular and extracellular spaces are linear and obey the principle of superposition (Hörning et al., 2010). In the case where the vessels are aligned perpendicular to the applied field (left panel, Figure 7), the positive and negative VEs between the two vessels are amplified. When the vessels are oriented at π/4 from one another (middle panel, Figure 7), the positive VE at the top of the left vessel and the negative VE at the bottom of the right vessel are amplified. When the vessels are arranged in the direction of the applied field (right panel, Figure 7), the positive and negative VEs between the two vessels are reduced in magnitude.

The complex non-linear dynamics of the membrane response to these VEs of different patterns results in different SI curves for these vessels, as shown in Figure 8. For each of the vessel alignments, the SI curve is similar for late DIs of between 325 and 350 ms; this is because make excitation dominates at late DIs, which is relatively independent of the complex VE patterns produced by the different vessel configurations. At early DIs of less then 300 ms, the different VE patterns result in different SI curves; however, the differences are minimal and each vessel alignment gives a similar response to the single vessel (with a vessel wall), shown in Figure 4. The θ = π/4 alignment requires the lowest shock strength to elicit propagation overall due to its configuration leading to the proximity of the strongest VE polarizations.

The ventricular slice geometry was pre-paced at the tissue level in a similar manner as for the isolated blood vessels; however, the SI curve was not computed. Instead, shocks were applied from anodal and cathodal electrodes in the bath space in the middle of the left ventricular cavity to a ground in the bath space surrounding the epicardium, resulting in an approximately transmural field. Shocks of 10 V/cm (computed via the minimum Euclidean distance from points on the surface of the cathode to the ground) were applied at DIs of between 290 and 350 ms after the previous stimulus. Regardless of the polarity of the electrode, the same type of make excitations were observed at late DIs (between 300 and 350 ms) from the transmural vessels and the same type of break excitations were observed at early DIs (between 290 and 300 ms), as with the idealized vessel geometries considered previously. As with the idealized vessel geometries, break excitations either continued to propagate or decay (as shown in Figure 9) depending on the DI at which the shock was applied.

However, a different mode of break excitation in the form of propagated graded responses (Trayanova and Rantner, 2003; Trayanova et al., 2003) attached to the subepicardial layer was observed around vessels proximal to the epicardium. Figure 9 shows snapshots of the transmembrane potential in the ventricular slice following the application of the shock. These boundary waves were caused by VE depolarizations around subepicardial vessels propagating into VE hyperpolarized regions which intersected the epicardium. The boundary propagation waves had peak transmembrane potentials of approximately −10 mV and action potential durations of approximately 35 ms; values significantly below the normal electrophysiological characteristics of the (ten Tusscher and Panfilov, 2006) cell model and similar to the propagating unstable wavelets described in Boyle et al. (2012). Approximately 120 ms after the shock was applied, these boundary propagation wavefronts re-entered and excited the entire tissue. By saving the cell model state variables (DeBruin and Krassowska, 1998; Ashihara and Trayanova, 2004; ten Tusscher and Panfilov, 2006) from nodes inside the boundary propagation waves and continuing the integration at the single-cell level, we determined that the short action potential durations and low transmembrane amplitudes of these waves were a direct consequence of electrotonic coupling and not an electrophysiological effect of the cell model itself. In contrast to the boundary propagation waves, single-cell action potentials had durations of approximately 200 ms (short in contrast to the normal human ventricular duration) and peak transmembrane potentials of approximately 20 mV.

Note that the transmembrane potential at the epicardium had equalized (from its shock-induced hyperpolarization) with the surrounding transmembrane potential when the boundary propagation was observed. At higher shock strengths, however, break excitations from vessel VEs reached the epicardium earlier and resulted in more rapid boundary propagation; this was due to the latent hyperpolarization from the shock-induced VEs. During anodal shocks, the inverted VE pattern gave rise to break excitations next to the endocardial surface from adjacent hyperpolarized areas; however, these boundary propagation waves rarely gave rise to bulk activation as they were often extinguished after collision with regions of prolonged refractoriness (due to depolarized VE regions) or failed to propagate along the relatively more curved endocardial boundary (trabeculation).

Some insight on the origins of the VEs elicited by blood vessels in isotropic tissue may be gained by using an analytical approach. Unlike anisotropic tissue where VEs may form in the tissue in response to conductivity anisotropy, VEs in isotropic tissue may originate only at tissue surfaces, where current enters or exits the bidomain. At steady state, the governing equations for isotropic tissue simplify to (5)∇2Vm−Vm∕λ2=0,∈Ωt(σi+σe)∇2ϕe=0,∈Ωtσb∇2ϕb=0,∈Ωb, where λ2=Rmσiσe∕βσi+σe is the space constant (here, R_m is the membrane resistance). Equation (5) assume that perturbations to ∇ϕ_e from the induced transmembrane potential (in response to field stimulation) are small, such that they may be neglected, and assumes a parallel combination of the intra- and extracellular isotropic conductivities for Laplace’s equation in the extracellular space (Sobie et al., 1997; Jolley et al., 2008). Equation (5) may be solved by separation of variables as in Pumir and Krinsky (1999), except here instead of assuming a constant extracellular field for ∇ϕ_e (Pumir and Krinsky, 1999); we first solve Laplace’s equation for ϕ_e to satisfy the electrostatic boundary conditions at the interfaces of the vessel cavity and vessel wall, and vessel wall and tissue. Assuming the far-field stimulus, of magnitude E₀, is aligned with the y-axis (E=E0y→) gives, in polar coordinates: (6)Vm(r,θ)=CE0λK1r∕λK′1a∕λsin(θ), where K_ν is the modified Bessel equation of order ν, of the second kind, and C is a scaling factor which comes from the solution of Laplace’s equation for the particular geometry considered: (7)C=2σbσwa2+a′2+σw2a2−a′2  σi+σeσwa2+a′2+σba2−a′2+σbσwa2+a′2+σw2a2−a′2 , where a′ = a − t. Assuming that the electric field is constant around the vessel implies that C = 1, giving the solution in Pumir and Krinsky (1999).

Substituting the transverse values for the intra- and extracellular conductivities (σ_it = 0.019, σ_et = 0.24 S/m) (Clerc, 1976) (approximating an out-of-plane fiber field parallel to the vessel axis), the blood (σ_b = 1.0 S/m) and vessel wall conductivities (σ_w = 0.01 S/m), into C and using the analytical expression for the wall thickness [equation (3) from Podesser et al. (1998)] gives C ≈ 0.15 and C ≈ 0.25 for vessels of radii a = 0.5 and a = 2.0 mm, respectively. In other words, the low-conductivity blood vessel wall reduces the current flux through the vessel cavity, significantly lowering the VE induced by the vessel in response to field stimulation. In the case of no vessel wall, the scaling factor becomes independent of the vessel size and simplifies to C = 2σ_b/(σ_b + σ_i + σ_e) ≈ 1.59.

The SI curve for a vessel of radius a = 2.0 mm in isotropic tissue is shown in Figure 10. The ratio in the field strengths required to elicit excitation, from the vessels with and without a vessel wall, is approximately 6.44 at an interval of 350 ms, a value which is approximately constant within the interval range of validity (290–350 ms). This is close to the ratio of the scaling factors, for the a = 2.0 mm vessel with and without a vessel wall, of 1.59/0.25 ≈ 6.36, suggesting that, in the case of isotropic tissue, the scaling factor is proportional to the minimum field strength required for VE induced excitation from the blood vessels. In isotropic tissue, the VE pattern is identical for vessels with and without vessel walls (albeit with different magnitudes)—these patterns are shown and discussed in detail in Pumir and Krinsky (1999) and Bittihn et al. (2012).

SI curves have been previously computed for unipolar stimulation both experimentally (Dekker, 1970; Sidorov et al., 2005) and in numerical bidomain simulations (Kandel and Roth, 2013, 2014). In general, the same mechanisms of make and break excitations, and the same overall morphology of the curves were observed from field stimulation of blood vessels, as have been observed (numerically and experimentally) in unipolar stimulation of myocardium (Sidorov et al., 2005; Kandel and Roth, 2013, 2014). Direct experimental validation of the VE patterns predicted by our simulations is challenging, due to the fact that vessels are inherently intramural. Thus, the close agreement in overall excitations mechanisms and SI curve morphologies with respect to previous experimental unipolar stimulation studies represents an important means of comparison. It should be noted that in both cases, break excitation occurs due to the co-location of de- and hyperpolarized VEs induced by the stimulus.

VEs may form around any heterogeneities, and the non-dimensional magnitude (|Vm∕(E0λ)|) of such VEs is dependent on the conductivity ratio of the heterogeneity and the surrounding tissue, and the characteristic length scale of the heterogeneity with respect to the local field direction E. For example, it is shown (Pumir et al., 2007) that, in isotropic tissue, elliptical heterogeneities oriented with their semi-major axes parallel to the local field require a stronger stimulus to elicit wave propagation than those oriented with their semi-minor axis parallel to the applied field. However, regardless of the geometry of the heterogeneity, or the direction of the local electric field, in isotropic tissue, the maximum non-dimensional VE is scaled by the conductivity ratio [discussed in equations (6) and (7)]. This does not hold, however, in anisotropic tissue, rendering it difficult to compare the magnitudes of VEs created by cylindrical blood vessels with other heterogeneities, such as fibrosis or sheets. In general, however, for heterogeneities in isotropic tissue with a non-zero conductivity ratio, the magnitude of the non-dimensional VE is approximately constant (and approximately equal to the conductivity ratio), provided the characteristic length scale is much bigger than the space constant (Pumir and Krinsky, 1999). The conductivity of regions of fibrosis is known to be much smaller than for healthy myocardium, so it is likely that the VEs elicited by these fibrosis heterogeneities will be dominated by anisotropy effects.

Previous experimental works highlighting the potential utility of low-energy defibrillation protocols (Fenton et al., 2009; Luther et al., 2012) have suggested the importance of the VEs formed around the coronary vasculature in providing a distributed source of wavefronts that facilitate arrhythmia termination. A consideration not addressed by these earlier works is the specific mechanism by which wavefront propagation is elicited from the vessel by the applied field with respect to the degree of refractoriness of the tissue.

Here, we have shown that break excitations may occur around blood vessels in otherwise refractory tissue. The ability to elicit wavefronts from relatively refractory tissue and not just from recovered diastolic tissue could therefore play a crucial role in optimizing low-energy defibrillation methods (Fenton et al., 2009; Janardhan et al., 2012; Luther et al., 2012; Rantner et al., 2013b) for which removing both excitable and soon to be excitable tissues is the main underlying goal. In the context of low-energy defibrillation, this may be especially relevant considering that the timing of the defibrillation shock, with respect to the time evolution of the excitable volume (Rantner et al., 2013b), is known to influence the defibrillation efficacy. However, an important finding from our study is that the field strength at which break excitations occur appears to be above that of low-voltage regimes, which typically quote field strengths <1 V/cm (Fenton et al., 2009; Luther et al., 2012; Rantner et al., 2013b). Therefore, our analysis suggests that such a mechanism may only be important in protocols that involve relatively higher strength shocks.

In this work, it was assumed that the tissue around the vessel repolarized simultaneously; however, in reality, a repolarization gradient, aligned with the propagation vector of the previous wave, will exist around the vessel. In sinus rhythm, the excitation generally propagates from the endocardium to the epicardium, producing a typical repolarization gradient of approximately 4.5 ms/mm—a value known to be sufficient to cause unidirection propagation from extrasystoles (Laurita and Rosenbaum, 2000) in guinea-pig hearts. It is however difficult to draw any conclusions regarding the arrhythmogenicity of extrasystoles originating from VEs around blood vessels; the induced VEs themselves will perturb the repolarization field and, in context, shocks will typically only be applied during fibrillation when the wave dynamics are chaotic and, by definition, not in sinus rhythm. The results from this work are mainly applicable to regions of tissue with small repolarization gradients, which may exist in regions around the myocardium during fibrillation.

In this work, and in much of the literature on defibrillation, the shock strength is quoted as the potential difference (Δϕ_e) between the shock electrodes divided by the distance between them (l). A single measure of shock strength is a crude approximation when the conductivity in the domain is not constant, as is generally the case, and when one of the shocking electrodes is a catheter. The magnitude of the stimulating electric field, |E|=|∇ϕe|, which measures the shock strength as a function of position, may vary significantly even in simple or idealized geometries such as considered in this work. To highlight this point, see Figure 11 which shows the heterogeneous nature of |∇ϕe| for the (Δϕe∕l) ≈ 10 V∕cm shock, which elicited the boundary wave propagation in the realistic slice geometry. Figure 11 shows how the field strength decays rapidly (approximately as 1/r^d − 1, where d is the number of dimensions) away from the point cathode and concentrates around the endocardial grooves [which, incidentally, is the reason for the preferential wave propagation from these regions as discussed in Rantner et al. (2013b)] and low-conductivity vessel walls. Thus, we highlight that quoted shock strengths for specific protocols may differ significantly from the localized field strengths sensed in specific regions of myocardium.

Subepicardial boundary propagation was observed to occur at early DIs in the image-derived ventricular slice geometry (see Figure 9). The phenomenon, observed after shocks in the relatively refractory period may represent a mode of shock-induced arrhythmogenesis, or an additional mechanism responsible for the earliest post-shock activations after the isoelectric window (Trayanova, 2007; Ashihara et al., 2008; Constantino et al., 2010). We hypothesize that the low transmembrane potential amplitude and short action potential duration of the propagating boundary wave, in addition to the shallow depth in which it propagates, may have prohibited its in vitro observation using voltage sensitive dyes in optical mapping experiments due to depth-averaging effects (Janks and Roth, 2002a,b; Bishop et al., 2007). However, it could be possible to observe this proposed mechanism using optrodes which record fluorescent signals directly from the intramural space.

It is thought that the boundary propagation waves are a form of “propagating unstable wavelets,” as described in Boyle et al. (2012). It is thought that due to the sealed boundary, current flows along the edge and accumulates longitudinally. This is just enough current to excite longitudinally but not transversely, with a partial activation of the fast sodium channel (Boyle et al., 2012). The electrotonic loading is also lower near the boundary due to the boundary condition, so a lower transmembrane current density is required to excite tissue down the concentration gradient (Kelly et al., 2013; Bishop et al., 2014). In addition, the high extracellular conductivity adjacent to the boundary (in the bath space) significantly increases the conduction velocity of the wave close to the surface (Henriquez et al., 1996, 2007; Bishop and Plank, 2011; Bishop et al., 2011). As the largest blood vessels in the coronary vasculature reside in the subepicardial layer and the magnitude of the induced VE depends on the size of the vessel, it is more likely that shocks involving an intracardiac cathode will lead to boundary propagation along the subepicardium, as this shock polarity leaves the epicardium hyperpolarized and excitable post-shock.

As shown in Figure 3 and in other works (Bishop et al., 2010a, 2012), specifically representing the insulating blood vessel wall within the computational model significantly alters the VE pattern created around the vessel, upon application of a field stimulus. The effect of resolving the blood vessel wall on the electrophysiological response of the tissue to field stimulus, as a function of the surrounding tissue refractoriness (DI) and shock strength, is shown to be significant (see Figure 4). Thus, it is important to resolve the coronary vasculature in computational bidomain models of field stimulation and to take account of the vessel wall, which was experimentally measured to have a conductivity of approximately two orders of magnitude lower than blood (Bishop et al., 2010a). The physical resolution of the blood vessel wall may however be computationally prohibitive, due to the necessarily small dimensions of the vessel wall. We propose that in more anatomically detailed 3D models, it may be possible to represent the macroscopic effect of the vessel wall by artificially lowering the conductivity of the medium inside the vessel cavity, such that the induced VE pattern is similar.

We demonstrate this in Figure 12 for the a = 2.0 mm vessel, where we applied a homogeneous equivalent conductivity σ_eq ∈ [σ_w, σ_b] to the vessel cavity. The equivalent conductivity may be computed by solving Laplace’s equation for the potential field in an (piecewise) isotropic medium, for the vessel geometry, ensuring the current density at the vessel surface, with the equivalent homogeneous conductivity, equals that when the vessel wall and blood are resolved (i.e., at the surface of the vessel, σwn→⋅∇ϕ=σeqn→⋅∇ϕ where n→ is the unit normal vector). The equivalent conductivity is stated as (8)σeq=t2σwσbw+2atσwσbw+2a2σbσw2atσbw−t2σbw+2a2σw, where σ_bw = σ_b − σ_w. Equation (8) provides a good approximation in the case of isotropic tissue (results not shown), but under-predicts the VE strength in anisotropic tissue, as shown in Figure 12. It is worth noting that, when the conductivities and dimensions are substituted into equation (8), the equivalent conductivity σ_eq is small (approximately two orders of magnitude lower than the minimum extracellular conductivity, for both the small and large vessels), implying that a zero-flux boundary condition for the extracellular potential at the vessel surface may provide an adequate approximation. These approaches require more investigation, however. An alternative approximation may be to apply a Robin-type boundary condition for the extracellular field across element edges/faces which constitute the blood vessel, in order to modify the extracellular current flux in a similar manner to that done for the intracellular current in Costa et al. (2014). Each of these proposed approximations requires the blood vessel surface to be spatially resolved in the computational models, however.

An additional implication of resolving the coronary vasculature is that it naturally distributes VE sources in the myocardium, causing break excitations in relatively refractory tissue which may result in boundary propagation along the hyperpolarized surface as shown in Figure 9. Without the proximity of VE depolarizations to the hyperpolarized boundary, this type of excitation may not occur. However, this mode of excitation does not require the vessel wall to be resolved per se.

The models used in this study were two dimensional, and thus, no out-of-plane effects on the VE patterns were investigated. However, using symmetry, we have investigated the effects of a homogeneous out-of-plane fiber field in the section on isotropic tissue. The accuracy of the experimentally measured conductivity for the blood vessel wall has some uncertainty, as this has only been measured once (Bishop et al., 2010a); the vessel wall conductivity is shown to significantly influence the shape of the SI curve (see Figure 6). The S₁ stimulus procedure resulted in no repolarization gradients around the vessels, which, although unphysiological, allowed us to generate SI curves which were not dependent upon initial propagation patterns. The propagating unstable wavelets have been, to the authors’ best knowledge, hitherto experimentally unobserved—thus, there is some uncertainty regarding the accuracy of these results. The effects of biphasic shocks were also not investigated in this work as here we sought to relate more closely with lower energy protocols which use monophasic shocks (Fenton et al., 2009; Janardhan et al., 2012; Luther et al., 2012; Rantner et al., 2013b); such analysis could provide an avenue for future work.