Impact Factor 3.394
2017 JCR, Clarivate Analytics 2018

The world's 3rd most-cited Physiology journal

Original Research ARTICLE

Front. Physiol., 03 April 2019 | https://doi.org/10.3389/fphys.2019.00337

Extended Bidomain Modeling of Defibrillation: Quantifying Virtual Electrode Strengths in Fibrotic Myocardium

  • 1Department of Physics and Applied Mathematics, University of Navarra, Pamplona, Spain
  • 2Nora Eccles Harrison Cardiovascular Research and Training Institute, University of Utah, Salt Lake City, UT, United States
  • 3Department of Biomedical Engineering, University of Utah, Salt Lake City, UT, United States

Defibrillation is a well-established therapy for atrial and ventricular arrhythmia. Here, we shed light on defibrillation in the fibrotic heart. Using the extended bidomain model of electrical conduction in cardiac tissue, we assessed the influence of fibrosis on the strength of virtual electrodes caused by extracellular electrical current. We created one-dimensional models of rabbit ventricular tissue with a central patch of fibrosis. The fibrosis was incorporated by altering volume fractions for extracellular, myocyte and fibroblast domains. In our prior work, we calculated these volume fractions from microscopic images at the infarct border zone of rabbit hearts. An average and a large degree of fibrosis were modeled. We simulated defibrillation by application of an extracellular current for a short duration (5 ms). We explored the effects of myocyte-fibroblast coupling, intra-fibroblast conductivity and patch length on the strength of the virtual electrodes present at the borders of the normal and fibrotic tissue. We discriminated between effects on myocyte and fibroblast membranes at both borders of the patch. Similarly, we studied defibrillation in two-dimensional models of fibrotic tissue. Square and disk-like patches of fibrotic tissue were embedded in control tissue. We quantified the influence of the geometry and fibrosis composition on virtual electrode strength. We compared the results obtained with a square and disk shape of the fibrotic patch with results from the one-dimensional simulations. Both, one- and two-dimensional simulations indicate that extracellular current application causes virtual electrodes at boundaries of fibrotic patches. A higher degree of fibrosis and larger patch size were associated with an increased strength of the virtual electrodes. Also, patch geometry affected the strength of the virtual electrodes. Our simulations suggest that increased fibroblast-myocyte coupling and intra-fibroblast conductivity reduce virtual electrode strength. However, experimental data to constrain these modeling parameters are limited and thus pinpointing the magnitude of the reduction will require further understanding of electrical coupling of fibroblasts in native cardiac tissues. We propose that the findings from our computational studies are important for development of patient-specific protocols for internal defibrillators.

Introduction

Several types of cardiac arrhythmias are treated with external or implanted defibrillators, which are devices for application of electrical current to a patient's thorax or heart (Al-Khatib et al., 2017). The electrical currents create an electrical field in the heart, which is determined by the location and geometry of the electrodes, the magnitude and waveform of the applied current, and the distribution, composition and electrical properties of tissues. Targets of defibrillation are constituents of cardiac muscle tissues, i.e., the myocytes. These cells respond to an extracellular electrical field caused by a defibrillator with changes of their transmembrane voltage. Modulated by the electrophysiological state of the myocytes, their transmembrane voltage will exhibit negative or positive shifts due to the electrical field. Thus, defibrillation might cause hyperpolarization and depolarization of myocytes. Also, defibrillation might trigger or modulate action potentials.

Several theories for mechanisms of defibrillation of the heart have been developed (Dosdall et al., 2010). The critical mass theory suggests that defibrillation success requires myocyte depolarization in a sufficient mass of tissue. The upper limit of vulnerability theory proposes the existence of a minimal electrical stimulus strength above which ventricular fibrillation cannot be induced even for stimuli occurring during the vulnerable period of the cardiac cycle. The sawtooth hypothesis and the syncytial heterogeneity hypothesis explains defibrillation effects on the transmembrane voltage of myocytes distally from the electrodes based on microscopic discontinuities of the intracellular space (Fishler, 1998). Similarly, the virtual electrode theory explains defibrillation effects based on structural heterogeneities such as endo- and epicardial surfaces as well as blood vessels. Virtual electrodes are thought to underlie local shifts of transmembrane voltages and thus trigger action potentials and elicit wave fronts that interfere with fibrillation. Recent work using computational modeling focused on understanding virtual electrodes caused by vessels (Connolly et al., 2017a,b) and related curvature of the surface of tissue heterogeneities to myocyte depolarization and initiation of wave fronts (Bittihn et al., 2012).

It is well-established that patients with myocardial fibrosis profit from defibrillator implantation (Iles et al., 2011), but it has also been reported that patients exhibit an increased mortality if they receive defibrillation shocks (Cevik et al., 2009). Several cardiac diseases including fibrosis are associated with discontinuity of the intracellular space and microstructural heterogeneity. For instance, myocardial infarction leads to scar regions with a decreased volume ratio of myocytes, but increased volume ratios of the extracellular space and fibroblasts. Also, our recent studies on a rabbit model of myocardial infarction revealed patches of fibrotic tissue interspersed within working myocardium distal from the scar (Seidel et al., 2017). Furthermore, several types of fibrosis of cardiac tissue are characterized by local increases of extracellular space and fibroblasts. Currently, our knowledge on defibrillation effects in fibrotic tissues is sparse and we do not even know if understanding these effects can help with parameterization of defibrillation protocols. Based on virtual electrode theory, it is conceivable that fibrotic tissue will exhibit substantial effects in response to defibrillation. However, experimental and computational studies on this subject have not been performed.

Here, we investigated the effects of defibrillation on transmembrane voltages in fibrotic tissues. We applied an extended bidomain model of electrical conduction in cardiac tissue. In contrast to the conventional bidomain model, which describes myocyte and extracellular domains only, the extended model comprises a description of a fibroblast domain. We performed simulations in one-dimensional (1D) and two-dimensional (2D) models of cardiac tissues with fibrotic patches of variable sizes, tissue composition and myocyte-fibroblast electrical coupling. We measured shifts of the myocyte membrane voltage caused by application of strong extracellular currents to assess the strength of virtual electrodes.

Methods

Extended Bidomain Modeling of Cardiac Conduction

We utilized the extended bidomain model for computational simulations of defibrillation in 1D and 2D domains (Sachse et al., 2009). The model was developed based on the established bidomain model (Tung, 1978; Henriquez, 1993) and allows description of cardiac tissues comprising more than one cell species. The original extended bidomain model comprised two cell species, i.e., cardiac myocytes and fibroblasts. Intracellular domains for these two species together with the extracellular domain constitute three domains of interest.

The mathematical formulation of the extended bidomain model applies Kirchhoff's electrical current conservation law to myocardium (Plonsey and Barr, 2007). Here we followed the exposition in Sachse et al. (2009), which described the electrical dynamics of the myocardium in three “homogenized” continuous domains. In short, the expression of current conservation at any point in space mathematically translates into a Poisson equation for each of the three domains:

·(σmyoϕmyo)= -fs,myo+ βmyo,fib Imyo,fib+ βmyo Imyo,e    (1)
·(σfibϕfib)= -fs,fib- βmyo,fib Imyo,fib+ βfib Ifib,e    (2)
·(σeϕe)= -fs,e- βmyo Imyo,e- βfib Ifib,e    (3)

where σmyo, σfib and σe are the electrical conductivity tensors of the myocyte, fibroblast and extracellular domain, respectively, ϕmyo, ϕfib , and ϕe are the intracellular potential of the myocyte, fibroblast and extracellular domain, respectively, and fs,myo, fs,fib, and fs,e (A/m3), are the stimulus current source densities for the myocyte, fibroblast and extracellular domain, respectively. Membrane voltage of myocytes Vmyo and fibroblasts Vfib was defined as the difference of their intracellular and extracellular potential. The number of myocytes per unit volume βmyo and fibroblasts per unit volume βfib (1/m3) were defined as:

βmyo= Volmyo Volmyo,single     (4)
βfib= Volfib Volfib,single     (5)

with the volume fraction of myocytes Volmyo and fibroblasts Volfib (dimensionless) as well as the individual volume for a myocyte Volmyo, single and fibroblast Volfib,single (m3). Currents between the myocyte and extracellular, fibroblast and extracellular, and myocyte and fibroblast domains are identified by Imyo,e, Ifib,e, and Imyo,fib(A), respectively.

The current between the myocyte and fibroblast domain was defined as:

Imyo,fib=ϕmyo- ϕfib Rmyo,fib    (6)

where Rmyo,fib is the resistance of gap junction channels (Ω). We defined the effective volume density βmyo,fib (1/m3) of electrical connections between fibroblasts and myocytes as:

βmyo,fib= Volfib Volfib,single     (7)

We assumed that the conductivity tensors (S/m) for each domain have a linear relationship with their respective volume fractions:

σmyo=Volmyo σ¯myo    (8)
σfib=Volfib σ¯fib    (9)
σe=Vole σ¯e    (10)

where the tensors σ¯myo, σ¯fib and σ¯e described conductivity for volume ratios of 100%. In the subsequently described simulations, we vary σ¯fib over a wide range of conductivities. We note that the range of σfib is smaller, due to scaling by the fibroblast volume fraction.

We applied a mathematical model of a rabbit ventricular myocyte (Mahajan et al., 2008) and a cardiac fibroblast (Sachse et al., 2008) to calculate the currents Imyo,e and Ifib, e, respectively. Temperature of the fibroblast model was set to the same temperature as the myocyte model (308 K). We applied a temperature coefficient (Q10) of 2 to adjust rate coefficients of the time- and voltage dependent outward current of the fibroblast model.

Implementation of Extended Bidomain Model

The extended bidomain model was implemented using the finite volume method and the programming language Fortran. The Portable, Extensible Toolkit for Scientific Computation (PETSc) (Balay et al., 2018) was used for parallelization and solving the Poisson equation. The myocyte (Cellml, 2008b) and fibroblast (Cellml, 2008a) models were implemented from CellML (Garny et al., 2008).

A simple forward Euler scheme with a fixed time step of 0.5 or 1 μs was applied to solve the time evolution of the model. The Poisson equation for the cellular domain was solved at each time step with a relative tolerance of 10−12. We applied the GMRES method, a member of the iterative Krylov methods that are implemented in PETSc. We used the additive Schwarz preconditioner in order to speed up the convergence for the iterative solver. The internal number of iterations for the Poisson solver for the 2D simulations was in the range of 50–100. This constituted the main burden of the numerical calculation. We checked for current conservation at each link of the discretized domain up to this precision (no numerical artifact of creation or annihilation of currents). In order to validate our numerical implementation, we also checked that the change in grid size discretization does not affect the results. We have compiled measured parameters in Tables S1S3 to show that variation of spatial resolution in the range of 25–100 μm has only minor effects. We further evaluated this new implementation of the extended bidomain model using our prior implementation (Sachse et al., 2009; Seemann et al., 2010). Modifications of the prior implementation were performed to account for discontinuity in the material properties of the cardiac tissue at the interface between the control and fibrotic regions.

Setup of 1D Simulations

We created 1D models of normal rabbit ventricular myocardium with an embedded central fibrotic patch (Figure 1A). The model was discretized with a spatial resolution of 50 or 100 μm. The patch exhibited differences in volume fractions of extracellular space, myocytes and fibroblasts measured in our previous work of fibrotic tissue (Greiner et al., 2018). Two different types of patches representing average and large fibrosis were configured (Table 1). We varied Rmyo,fib and σ¯fib in simulations. The simulations started with applying 9 depolarizing myocyte membrane currents at the location 32.5–33 mm of the domain with an amplitude of 20 μA/cm2 and a duration of 5 ms at a frequency of 4 Hz. At 250 ms after the last stimulus, extracellular currents were applied at 1–1.5 mm and 32.5–33 mm of the domain with an amplitude of 0.7 × 105 and −0.7 × 105 A/m3, respectively. These injected currents created a nearly uniform electric field between the two sites of magnitude E = 0.546 V/cm (Figure S1). This electric field, in turn, was responsible for the generation of the virtual electrodes at the left (referred as point 1) and right (point 2) boundary of the fibrotic patch.

FIGURE 1
www.frontiersin.org

Figure 1. Geometry of (A) 1D and (B,C) 2D models of rabbit ventricular tissue with central fibrotic patch. The sites of extracellular current application are marked in black. Fibrotic patches are shown in grey.

TABLE 1
www.frontiersin.org

Table 1. Parameters of extended bidomain model.

In order to measure virtual electrode strengths, we monitored the difference in the membrane voltage of the myocytes Vmyo and fibroblasts Vfib at patch boundaries during the application of extracellular current with respect to the membrane voltages before current application. Vmyo was not a monotonic function of time during the current application (Figure S2). Therefore, we reported the largest absolute value (maintaining its sign) of the monitored difference during the current application as a measure of the local virtual electrode strength. In addition, we quantified space constants λmyo and λfib by fitting the spatial profile of Vmyo and Vfib, respectively, at and proximal to the boundaries of the patch.

Setup of 2D Simulations

We generated 2D models of normal rabbit ventricular myocardium with an embedded fibrotic patch of square and disk-like shapes (Figures 1B,C). The square side length and the disk diameter were both set to 1.5 mm. We discretized the model with a spatial resolution of 50 μm in x and y direction. The computational domain was 1 cm x 1 cm in size. We applied four intracellular stimuli at 9–9.5 mm along the length and entire width with an amplitude of 20 μA/cm2 for 5 ms at a frequency of 4 Hz. Two hundred and fifty milliseconds after the last stimulation, extracellular currents were applied at 0.5–1 mm and 9–9.5 mm along the length and throughout the width with an amplitude of 0.7 × 105 and −0.7 × 105 A/m3, respectively. These injected currents created a nearly uniform electric field between the two electrodes of magnitude 0.530 V/cm.

Data Analysis

Effects of patch length Wp on ΔVmyo were analyzed using the fitting function:

ΔVmyo=a(1-e-Wpb)    (11)

with the parameters a and b. We used a different fitting function for ΔVfib, because it was not always crossing the origin of axes:

ΔVfib=c-ae-Wpb    (12)

with the parameters a, b and c. We measured the fit quality using the adjusted R2.

Results

1D Simulations

We performed numerical simulations with the 1D extended bidomain model of control tissue with a central patch configured with parameters for average fibrosis (Figure 1A). The model was discretized with a spatial resolution of 100 μm. In Figures 2A–C, we present the results of a simulation of a propagating wave initiated by intracellular stimulation followed by application of extracellular current. The model parameters were Rmyo,fib = 4 × 104 MΩ and σ¯fib=0.1 S/m. After scaling with the fibroblast volume fraction (Equation 9) in control tissue (3%) and patch (10%), σfib amounted to 0.003 and 0.01 S/m, respectively. In this case, the propagation of action potentials was only marginally affected by the fibrotic patch. Extracellular current altered ϕe as well as Vmyo and Vfib at the application sites and the borders of the fibrotic patch.

FIGURE 2
www.frontiersin.org

Figure 2. Representative plots of (A) ϕe, (B) Vmyo, and (C) Vfibfor intracellular pacing and subsequent extracellular current application in the 1D model with the average fibrosis patch, σ¯fib of 0.1 S/m and Rmyo,fib of 40 GΩ. Spatial distribution of (D) ϕe, (E) Vmyo, and (F) Vfib at the end of application of extracellular currents. The extracellular current led to alterations of Vmyo and Vfib not only at the application site, but also at the left and right border of the fibrotic patch.

We present the spatial distribution of ϕe, Vmyo, and Vfib at the end (t = 2.255 s) of current application in Figures 2D–F, respectively. Due to current application, Vmyo at the left border of the patch was increased from −84.95 mV before current application to −75.59 mV at the end of the current application. Vfib was reduced from −54.05 mV before current application to −58.97 mV at the end of the current application. At the right border of the patch Vmyo was decreased from −85.00 to −93.44 mV. Vfib was increased from −54.21 to −51.18 mV. We determined the electric field for average and large fibrosis during the current application (Figures S1A,B, respectively). The spatial variation outside of the region for current application was small.

We applied the 1D model in a series of simulations varying Rmyo,fib and σ¯fib. Example simulated ϕe, Vmyo, and Vfib are presented in Figure S3. Before current application, the biphasic relationship of Vmyo with Rmyo,fib was nearly identical on both boundaries of the patch (Figures 3A,C). Also, variation of σ¯fib did not affect Vmyo. After current application, the left border of the patch exhibited increased Vmyo, measured as ΔVmyo(1), dependent on Rmyo,fib and σ¯fib (Figure 3B). At the right border of the patch, the reduction of Vmyo, measured as ΔVmyo(2), was dependent on Rmyo,fib and σ¯fib (Figure 3D). At both borders of the patch, effects on Vmyo were of similar magnitude for all σ¯fib for large Rmyo,fib.

FIGURE 3
www.frontiersin.org

Figure 3. Measurement of effects of extracellular current application on Vmyo using the 1D model with the average fibrosis patch. Before application of extracellular currents, Vmyo at the (A) left and (C) right border of the fibrotic patch was marginally affected by Rmyo,fib and σ¯fib. (B) At the end of current application, the left patch border exhibited a positive ΔVmyo. (D) In contrast, the right patch border exhibited a negative ΔVmyo.

We present corresponding measures on Vfib in Figure 4. The biphasic relationship of Vfib with Rmyo,fib before application of extracellular current (Figures 4A,C) was similar to the relationship of Vmyo with Rmyo,fib. For low Rmyo,fib (high myocyte-fibroblast coupling), Vfib was similar to the resting Vmyo (Figures 3A,C). For high Rmyo,fib (low myocyte-fibroblast coupling), Vfibwas close to the membrane voltage of isolated fibroblasts, i.e., ~58 mV (Shibukawa et al., 2005; Sachse et al., 2008). Similar as myocytes, fibroblasts at the borders of the patch exhibit alterations of their membrane voltage, ΔVfib(1) and ΔVfib(2), dependent on Rmyo,fib and σ¯fib in response to current application (Figures 4B,D). For low Rmyo,fib, ΔVfib was roughly similar to the ΔVmyo at both borders. In contrast, for high Rmyo,fib, ΔVfib exhibited a reverse sign vs. the corresponding ΔVmyo at both borders. This results in the discontinuity in ΔVfib for Rmyo,fib of 100–1,000 MΩ (Figures 4B,D).

FIGURE 4
www.frontiersin.org

Figure 4. Measurement of effects of extracellular current application on Vfib using the 1D model with the average fibrosis patch. Before application of extracellular currents, Vfib at the (A) left and (C) right border of the fibrotic patch was strongly affected by Rmyo,fib and σ¯fib. At the end of current application, the (B) left and (D) right border exhibited a ΔVfib with magnitude and sign modulated by Rmyo,fib and σ¯fib.

Using the same geometrical model and experimental protocol, we also performed simulations using a patch with large fibrosis (Figures 5, 6). While the results were qualitatively similar as for the average case of fibrosis described above, we note important differences. In particular, the higher degree of fibrosis was associated with a higher magnitude of ΔVmyo and ΔVfib. We observed conduction block for Rmyo,fib of 100–1,000 MΩ (Figures 5A,C). We confirmed that the magnitude ΔVmyo is a monotonously increasing function of the degree of fibrosis (Figures 7A,B). The parameter α described the variation between the low degree of fibrosis (α = 0) to the high degree of fibrosis (α = 1). In the range of settings that we have studied, the variation of ΔVmyo with α was roughly linear. Similarly, for low Rmyo,fib the magnitude ΔVfib was a monotonously increasing function of α (Figures 7C,D). For high Rmyo,fib, the magnitude ΔVfib was approximately constant for α between 0 and 1.

FIGURE 5
www.frontiersin.org

Figure 5. Measurement of effects of extracellular current application on Vmyo using the 1D model with the large fibrosis patch. Before application of extracellular currents, Vmyo at the (A) left and (C) right border of the fibrotic patch was marginally affected by Rmyo,fib and σ¯fib. (B) At the end of current application, the left patch border exhibited a positive ΔVmyo. (D) In contrast, the right patch border exhibited a negative ΔVmyo.

FIGURE 6
www.frontiersin.org

Figure 6. Measurement of effects of extracellular current application on Vfib using the 1D model with the large fibrosis patch. Before application of extracellular currents, Vfib at the (A) left and (C) right border of the fibrotic patch was strongly affected by Rmyo,fib and σ¯fib. At the end of current application, the (B) left and (D) right border exhibited a ΔVfib with magnitude and sign modulated by Rmyo,fib and σ¯fib.

FIGURE 7
www.frontiersin.org

Figure 7. Assessment of effects of extracellular current application on Vmyo and Vfib simulated in 1D models with a patch of varying degrees of fibrosis. The degree of fibrosis ranged from average (α = 0) to large (α = 1). We also varied Rmyo,fib. σ¯fib was set to 0.1 S/m. The effects on the magnitude of ΔVmyo at the (A) left and (B) right boundary increased with the degree of fibrosis. Similarly, effects of ΔVfib at the (C) left and (D) right boundary increased with the degree of fibrosis for small Rmyo,fib. Marginal effects were present for larger Rmyo,fib.

We assessed the influence of patch length on virtual electrode strengths (Figure 8). The simulations were performed with a spatial resolution of 50 μm in order to explore smaller patch lengths. We present fit parameters of ΔVmyo and ΔVfib in Table 2. The patch length was a major modulator of virtual electrode strengths. The larger the patch, the higher the effect on membrane voltages. However, there is a limit to this phenomenon, with the effect tapering off for longer patches. We assessed the length scale, which depends on the coupling between the myocyte and fibroblast domain through the coupling parameter Rmyo,fib. For large Rmyo,fib, the spatial length scales for the myocytes and fibroblasts correspond to their respective space constants λmyo and λfib. In contrast, for low Rmyo,fib the spatial length scales were the same for the two species due to the coupling and settled at an intermediate value between λmyo and λfib. Indeed, larger difference of the ΔVmyo and ΔVfib were at the sub-millimeter scale (Figure 8), which is in agreement with the estimated values for of λmyo and λfib (Tables S1S3).

FIGURE 8
www.frontiersin.org

Figure 8. Effects of patch length Wp on (A,B) Vmyo and (C,D) Vfib at the boundary of the patch in response to application of the extracellular current. We varied Rmyo,fib, but σ¯fib was set on 0.1 S/m. For small Wp, Vmyo, and Vfib were greatly reduced.

TABLE 2
www.frontiersin.org

Table 2. Relationships of membrane voltage changes and Wp according to the fitting functions.

2D Simulations

We investigated virtual electrodes in a 2D spatial domain of control tissue with central patches configured with parameters for average fibrosis (Figures 1B,C). Because we know from the 1D simulations that effects of the patch on ΔVmyo and ΔVfib are larger at the millimeter than sub-millimeter scale, we set the patches (square and disk) to a size of 1.5 mm.

Figure 9 displays the 3 main fields at the end of the extracellular current application in a model with a square-shaped patch, σ¯fib of 0.1 S/m and Rmyo,fib of 1 GΩ. The injected currents at the electrode generated a fairly linear distribution for ϕe (Figures 9A,D), which implies a quasi-uniform electric field. Only close to the cathode (x≈9 mm) the field was distorted by a nascent wave initiated at the cathode location. Closer inspection revealed that the field is really uniform close to the patch (Figure 9D). The corresponding spatial distribution of Vmyo for the same time shows the virtual electrodes at the two edges of the square that are perpendicular to the applied field (Figures 9B,E). The edge corresponding to the small value of x exhibited a depolarization pattern (equivalent to point 1 in the 1D setting). On the contrary the other edge exhibited a hyperpolarization pattern (equivalent to point 2 in the 1D setting). Figures 9C,F display the field for Vfib. In this example, Rmyo, fib was set to 1 GΩ, which is the value for which the coupling between the myocytes and the fibroblast cancels out the effect of the conductivity heterogeneities at the patch boundaries. According to the theory (see section Discussion), we expect to observe hyperpolarization for Vfib at point 1 and depolarization at point 2, which is the opposite of what happens to Vmyo. However, the coupling term between the two domains caused that the Vmyo is pulling the Vfib in the same direction. At the very edge of the patch (Figure 9F) the two effects cancel out.

FIGURE 9
www.frontiersin.org

Figure 9. Spatial distribution of (A) ϕe, (B) Vmyo, and (C) Vfib at the end of extracellular current application in the 2D model with the square-shaped average fibrosis patch, σ¯fib of 0.1 S/m and Rmyo,fib of 1 GΩ. Close-ups of the spatial distribution for (D) ϕe, (E) Vmyo, and (F) Vfib. The extracellular current led to alterations of Vmyo and Vfib not only at the application site, but also at the left and right border of the patch.

We repeated the simulation using a disk-shaped patch with Rmyo,fib = 1 TΩ, which represents negligible coupling between Vfib and Vmyo (Figure 10). The simulation yielded a linear distribution for ϕe and the simulated electric field was almost identical as the field obtained with the square patch (Figures 10A,D). In Figures 10B,E we present the corresponding spatial distribution of Vmyo. As before, we observed the two virtual electrodes at the edges of the patch. Due to the disk-shaped patch, the virtual electrodes were not so sharp and appeared diffuse also in the parallel direction of the electric field. Figures 10C,F confirmed that polarities of ΔVfib are indeed opposite to the ΔVmyo at the patch boundary. Also, the spatial scale for the virtual electrode was much smaller than the virtual electrodes for the Vmyo field, which agrees with the corresponding values for the space constants associated with the two fields.

FIGURE 10
www.frontiersin.org

Figure 10. Spatial distribution of (A) ϕe, (B) Vmyo, and (C) Vfib at the end of extracellular current application in the 2D model with the disk-shaped average fibrosis patch, σ¯fib of 0.1 S/m and Rmyo,fib of 1 TΩ. Close-ups of the spatial distribution for (D) ϕe, (E) Vmyo, and (F) Vfib at the end of application of extracellular currents. The extracellular current led to alterations of Vmyo and Vfib not only at the application site, but also at the left and right border of the patch.

Figure 11 and Table S4 present quantitative measures for the virtual electrode strengths for the two patch geometries. We report the results for ΔVmyo and ΔVfib at both locations of the patch boundaries as a function of Rmyo, fib. We compared them with the corresponding values obtained for the 1D case. Using bidomain models it has been demonstrated that convex shape of boundaries leads to lower activation of the membrane (Entcheva et al., 1998, 1999; Pumir and Krinsky, 1999; Bittihn et al., 2012). This was confirmed in our simulation for ΔVmyo where the virtual electrode strengths were lower for the disk than for the square. However, the difference between the two geometries is rather minute. We measured an approximate 4% difference of ΔVmyo for low values of the Rmyo,fib and an approximate 5% difference for high values of the Rmyo,fib. Also, ΔVmyo for the square-shaped were higher than for the disk-shaped patch, and close to the 1D patch geometry. ΔVfib for the 1D and the two 2D geometries were similar (Figures 11C,D). For ΔVfib, the influence of patch geometry was less pronounced.

FIGURE 11
www.frontiersin.org

Figure 11. Comparison of (A,B) Vmyo and (C,D) Vfib from the 1D simulations and the 2D simulations with disk and square patch geometries. The square and circle symbols refer to the results for the respective square and disk patch geometry.

Discussion

Our computational study provided insights into the relationship between parameters of fibrotic patches and the strength of virtual electrodes caused by defibrillation. We applied the extended bidomain model to reflect that cardiac tissues comprise fibroblasts beyond myocytes. Our study showed that an increased degree of fibrosis and an increased size of the fibrotic patch cause an increased strength of virtual electrodes. Also, we found that the composition of tissue heterogeneities is a modulator of virtual electrode strength. Increased electrical coupling of myocytes with fibroblasts reduced the virtual electrode strength. Intra-fibroblast coupling reduced virtual electrodes in case of high myocyte-fibroblast coupling.

We explain our findings of depolarization and hyperpolarization of membranes during the application of the extracellular currents at the boundaries of the fibrotic patch based on reformulation of the equations for the extended bidomain. The equations can be rewritten in the following simple form for numerical solution see Equations (17–20) in Sachse et al. (2009):

Vmyot=1Cmyo(·(σmyoVmyo)+·(σmyoϕe)+)    (13)
Vfibt=1Cfib(·(σfibVfib)+·(σfibϕe)+)    (14)
·{(σe+σmyo+σfibϕe)}= -fe-·(σmyoVmyo)-·(σfibVfib)    (15)

with the injected current fe due to extracellular current application, the myocyte membrane capacitance Cmyo and the fibroblast membrane capacitance Cfib. Equations (13) and (14) are parabolic equations for the time evolution of the myocytes and fibroblast membrane voltage. Equation (15) is the Poisson equation for calculating the extracellular potential. In these equations, ⋯  indicates terms that are not of interest for our initial discussion and hence omitted.

During the current application, we injected positive charge at the electrode located at 1–1.5 mm and negative charge of the same magnitude at the electrode located at 32.5–33 mm of the extracellular domain. This current injection produced a rather uniform gradient of ϕe (Figure 2). This corresponds to a quasi-uniform E in the system (Figure S1):

E= - ϕe    (16)

With our computational set-up, E is oriented in the same direction as the x-axis. E together with the discontinuity of the conductivity at the patch boundary is responsible for the initial changes of the membrane voltage during the current application. This can be shown by decomposing the second right hand term of Equation (13) into:

Vmyot -(σmyo)·E+σmyo2ϕe    (17)

This represents the generalized activating function previously described to predict the distribution of virtual electrodes (Sobie et al., 1997). The crucial term in determining the direction of the membrane voltage change during the current application is the first right hand side term. If ∇σmyo is oriented in the same direction as E one observes a decrease of Vmyo (hyperpolarization). If ∇σmyo points in the opposite direction of E, the current application leads to an increase of Vmyo (depolarization). We extend this description for Vfib using a decomposition of the second right hand term of Equation (14):

Vfibt -(σfib)·E+σfib2ϕe    (18)

In our setup, due to the variations of the volume fraction of the myocytes and fibroblasts, ∇σmyo is pointing in the opposite direction of E at point 1 of the domain, which results in a depolarization of Vmyo. On the contrary, ∇σfib at point 1 is aligned with E, which results in a hyperpolarization of Vfib at point 1.

However, this simple mechanistic explanation ignores the coupling term between the fibroblast and the myocyte domain. The current Imyo, fib, which is determined by the parameter Rmyo, fib and the voltage between the domains, explains why in some cases, Vfib follows Vmyo. In particular, at high coupling between myocyte and fibroblast domain (i.e., low values of Rmyo,fib), Vmyo is very similar to Vfib. Interestingly, at intermediate values of Rmyo,fib, the simulations revealed a competition of the two antagonistic effects for Vfib (Figures 4, 6). These two effects can even cancel each other out in a sort of tug of war situation, leaving Vfib at the patch boundary apparently unaffected by the application of the current (Figure 9F).

We note that the arguments presented here are valid in any spatial dimension. Based on the same arguments we also explain the mechanism behind the decrease of the virtual electrode strengths for large fibroblast-myocyte coupling. Our studies showed that due to heterogeneity of conductivity at the patch boundary, Vfib and Vmyo tend to change in opposite directions. One membrane is depolarizing and the other one is hyperpolarizing, but the effect of the heterogeneities can be canceled out and even be inverted for Vfib because of current flow between the fibroblast and myocyte domain. This flow is lowering the effect of the conductivity heterogeneities and the virtual electrode strength. Theoretically, the current flow can go in either direction. If the fibroblast is the dominant constituent of the cardiac tissue at some heterogeneities, Vfib could pull Vmyo to follow its electric depolarization or hyperpolarization. With the parameters applied in this study, myocytes were volume-wise dominant in the tissue. Thus, we did not observe such behavior.

For the 1D simulations, we reported extrema of the time evolution of ΔVfib during the extracellular current application. In Figures 4B,D, the reported ΔVfib is the largest absolute value (keeping its sign) of the monitored difference. These indicators exhibited a discontinuous behavior for intermediate values for the Rmyo,fib as exemplified in Figures 4B,D. The discontinuity is caused by the transition between dominant minima and maxima of Vfib as illustrated in Figure S2B. The transition occurred at decreased Rmyo,fib, when increasing σ¯fib. This can be explained by the fact that the pull exerted by the myocytes on the fibroblasts diffuses faster when the diffusion is increased and is therefore less effective.

The 1D simulations for the large fibrosis revealed a complex behavior of Vmyo. In Figures 5A,C, discontinuities occurred at intermediate values of Rmyo, fib. This was caused by block of the propagating waves prior to the extracellular current application at the fibrotic patch, which modified the values of Vmyo. Though this effect is important for the dynamics of the system, in the context of our study it did affect only marginally ΔVfib and ΔVmyo.

In our 2D simulations, the influence of patch geometry on ΔVfib was small (Figures 11C,D). This can be explained by the fact that the space constant for the fibroblast model is very small (computed previously around 0.15 mm) vs. the patch size of 1.5 mm. We refer the reader to Bittihn et al. (2012) for a more detailed description of the underlying mechanism.

A straightforward conclusion from our study is that fibrotic patches, similar as endo- and epicardial surfaces as well as blood vessels, constitute a potential site for wave initiation in response to defibrillation. This finding could help in parameterization of defibrillators. We suggest that further studies will lead to defibrillator protocols specifically for patients with known myocardial fibrosis to ensure defibrillation success. Also, based on our findings we propose that knowledge in the spatial distribution of fibrosis is valuable for optimization of positions for defibrillation electrodes. Potential clinical benefits of accounting for cardiac fibrosis are reduction in myocardial injury due to optimized defibrillation energy.

Prior studies on modeling of effects of fibrosis aimed at characterization of conduction and conduction defects as well as arrhythmogenesis and maintenance of arrhythmia in the atria and ventricles. Fibrosis was introduced in the models with various approaches including embedding of conduction barriers of various sizes and geometries, reduction of inter-myocyte coupling, regional replacement of myocytes with fibroblasts or myofibroblasts, and incorporation of myocyte-fibroblast coupling (Fishler, 1998; Fishler and Vepa, 1998; Mcdowell et al., 2012; Costa et al., 2014; Zeigler et al., 2016). Here, we applied the extended bidomain model to integrate various aspects of fibrotic remodeling.

A difficulty of applying the extended bidomain model is that we have only vague knowledge on parameters related to the fibroblast domain, in particular, the parameters Rmyo, fib and σfib. Quantitative measurements of these parameters have not been performed in normal and diseased myocardium (Kohl and Gourdie, 2014). However, a small number of in vitro measurements provided quantitative information to constrain these parameters. Beyond those, qualitative and indirect evidence suggests that etiologies and the progression of fibrotic remodeling affect the degree of myocyte-fibroblast and fibroblast-fibroblast electrical coupling in diseased myocardium. For these reasons, computational studies with the extended bidomain model and other computational modeling approaches commonly applied measurements from in vitro studies to justify parameter settings or varied parameters to reconstruct experimental findings. In this study, we chose to vary parameters within a wide range to go beyond the in vitro measurements and comprehensively explore potential effects of coupling.

We varied Rmyo, fib from 1 MΩ to 1 TΩ to describe strong and negligible electrical coupling, respectively, between fibroblasts and myocytes. The lower limit is close to the junctional resistance of a myocyte pair isolated from rabbit ventricular tissue (0.81 MΩ) (Kieval et al., 1992). The upper limit led to negligible currents between the fibroblast and myocyte domain, and effectively represents the case of electrical uncoupling of the domains. Experimental measurements on myocyte-fibroblast pairs in 24 h-old culture of neonatal rat hearts cells suggest that Rmyo, fib is variable and within a range of 125 MΩ and 3.23 GΩ (Rook et al., 1992). The variability was explained by the variable number of gap junction channels of the myocyte-fibroblast pairs. Prior modeling work applied measurements of gap junction channel conductances in the same preparation and assumed 10 to 30 gap junction channels per cell pair to vary Rmyo, fib from 1.11 to 3.33 GΩ for simulation of fibroblast effects in the sinoatrial node (Rook et al., 1989; Kohl et al., 1994). Assuming 0 to 75 gap junctions per cell pair, Rmyo, fib was varied from 0.44 to ∞ GΩ in a modeling study of fibroblast effects on conduction (Jacquemet and Henriquez, 2007). Based on the prior experimental and computational studies, the lower Rmyo, fib limit (1 MΩ) in our studies appears small. We note that effects of Rmyo, fib between 1 and 10 MΩ on virtual electrodes were in general similar. This range represents the extreme case of strong myocyte-fibroblast coupling. Increasing the lower limit of the parameter range to e.g., 1 GΩ does not change our conclusions on the relationship between Rmyo, fib and virtual electrodes.

We varied σ¯fib from 0 to 0.5 S/m to describe absence of electrical coupling and strong coupling, respectively, within the fibroblast domain. The lower limit represents the case of electrical uncoupling of fibroblasts. The upper limit is identical to our setting of σ¯myo. Dependent on Volfib (Equation 9), this led to a range of σfib of 0 and up to 0.1 S/m in the region with high fibrosis with Volfib of 20%. An experimental framework for estimating σfib comprises measurements of junctional conductance in pairs of cultured rat cardiac cells. Conductances of fibroblast-fibroblasts pairs were in the range of 175 pS and 6 nS (Rook et al., 1992). Measurements of conductance in myocyte-myocyte pairs were up to 42 nS (Rook et al., 1990). Estimates of conductance based on gap-junctional area in myocyte pairs range between 0.9 and 216.6 nS. These measurements and estimates of conductance between cell pairs are not sufficient to determine σfib in normal and diseased myocardium. Further information on fibroblast geometry, gap junction distribution and cytosolic conductivity would be required to estimate σfib using computational approaches introduced by us previously (Bauer et al., 2013; Greiner et al., 2018). Considering that conductances of freshly isolated ventricular myocyte-myocyte pairs from control rabbit hearts were much higher (1.24 ± 0.25 μS) (Kieval et al., 1992), it is also unclear how the in vitro measurements relate to conductance in native myocardium. However, based on this information and the explored volume fractions, we argue that σfib cannot be larger than σmyo even in the region of high fibrosis with identical volume fractions for the fibroblast and myocyte domain. In this region, our calculations (Equations 8, 9) using the maximal σ¯fib (0.5 S/m) led to identical σfib and σmyo (0.1 S/m). Prior modeling studies of electrical conduction in human atrial tissues applied σfib in the range of 0.02–0.1 S/m (Greisas and Zlochiver, 2016). Also, σfib was set to 0.06 S/m in a study on rabbit heart (Corrias et al., 2012). Based on these arguments and the prior work, we limited the exploration of σfib to values up to 0.1 S/m in our simulations. We note that this value is the maximum in our simulations and it is likely that smaller values in the explored range are more realistic. For small Volfib and thus sparse inter-fibroblast coupling, we point the reader at our simulations with σfib equal to zero.

We note that, despite the large flexibility of the extended bidomain model to describe tissue remodeling, the computational demands of the model are similar as for conventional bidomain models. A minor increase of computational demands is caused by the integration of the fibroblast domain.

We suggest that the extended bidomain model provides a more comprehensive framework for modeling of conduction in control and diseased tissue than the conventional bidomain model. While the additional parameters enable already more detailed parameterization, various refinements can be envisioned. For instance, an alternative method to define the effective volume density (1/m3) of electrical connections between fibroblasts and myocytes is:

βmyo,fib=β(r)  Volfib Volfib,single     (19)

Because the electrical coupling between the fibroblasts and myocytes is mainly unknown and potentially varies from one situation to another we can define a parameter, β(r) (dimensionless), for characterizing this coupling. A value of β(r) = 1 corresponds to a situation where each of the fibroblasts has an electrical connection with the myocyte domain. A larger (smaller) value of β(r) allows changing this situation. For β(r) < 1 not all the fibroblasts are coupled to a myocyte and conversely for β(r) > 1 each fibroblast is coupled to more than one myocyte. Both, Rmyo, fib and β(r) are factors that contribute to the coupling between myocytes and fibroblasts and hence their effective influence can be combined in a single variable, β(fract) defined as:

β(fract) = β(r) Rmyo,fib    (20)

Limitations

We acknowledge several limitations related to the presented work. The studies were performed in 1D and 2D domains with isotropic conductivities. We did not consider unequal anisotropy ratios of conductivities. We applied simple descriptions of fibrotic remodeling derived from a rabbit model of cardiac infarction. While fibrotic remodeling is known to be a multi-faceted process comprising e.g., activation of fibroblasts, their differentiation into myofibroblasts, and subcellular remodeling of structural und electrophysiological properties, we primarily considered increased volume ratios of cells and the extracellular space in our study. Also, we applied mathematical models of normal myocytes and fibroblasts, and did not adjust these models to account for fibrosis and associated remodeling.

A limitation of our work is related to incomplete understanding of myocyte-fibroblast and fibroblast-fibroblast coupling in the heart. As discussed above, experimental measurements exist only for in vitro models and it is difficult to apply these measurements for parameterization of models of fibrotic remodeling in native myocardium. Also, we applied a simple linear model for scaling of σfib by Volfib (Equation 9). Non-linear models might be more realistic. In particular, for small Volfib, it is likely that fibroblasts are not electrically coupled with other fibroblasts, which can be reflected by a model using piece-wise defined functions. We note the controversy on myocyte-fibroblast coupling and the contribution of fibroblasts in physiological and pathophysiological conduction (Kohl and Gourdie, 2014). Based on our simulations, myocyte-fibroblast coupling is a major determinant of virtual electrode strength. However, for the entire simulated range of myocyte-fibroblast coupling, boundaries of fibrotic patches were associated with virtual electrodes.

In our comparison of software for extended bidomain modeling, we confirmed that the new and our prior implementation yielded very similar results when simulating a homogenous tissue. Peak Vmyo and Vfib differed only by 5.08% (13.98 vs. 13.3 mV) and 5.23% (13.69 vs. 13.03 mV), respectively. However, we noticed that the approach for discretization of the patch boundary affects the magnitude of the altered membrane voltages. Conventional finite difference discretization (Sachse, 2004) led to significant differences vs. the finite volume discretization used in our studies. We reduced differences to the finite volume method by implementation of a finite difference approximation based on Kirchhoff's law (Witwer et al., 1972). However, the appropriate modeling approach should be dictated by the spatial distribution of the structural remodeling. While our finite volume approach described structural remodeling as an abrupt change from control tissue to the fibrotic patch, a gradual approach might be more appropriate (Seidel et al., 2017).

A related limitation of our study is the geometry and size of patches explored in the 2D simulations. We focused on disk and square shaped patches with a single size. Larger effects of geometry are expected for smaller sizes, i.e., below the millimeter scale, but we leave this for a future exploration. Also, we suggest that microscopic imaging can provide insights into the detailed distribution of remodeling.

Data Availability

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

JB and FS designed the study. All authors developed software for this project, analyzed, and interpreted the data, drafted and critically revised the manuscript, and approved the version to be published.

Funding

JB acknowledges the support of a Fulbright Fellowship for his stay at the Nora Eccles Harrison Cardiovascular Research and Training Institute at the University of Utah and also partial support through a grant project by the Spanish Ministerio de Economía y Competitividad (MINECO) under number SAF2014-58286-C2-2-R. FS acknowledges financial support from the Nora Eccles Treadwell Foundation and the National Institutes of Health (R01 HL 135077 and R01 HL 132067).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Mr. Marcus Blackburn, Mr. Joachim Greiner, and Dr. Gunnar Seemann for their help and discussions. We also thank Mr. Haonan Yang and Mrs. Nancy Allen for technical support.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2019.00337/full#supplementary-material

References

Al-Khatib, S. M., Stevenson, W. G., Ackerman, M. J., Bryant, W. J., Callans, D. J., Curtis, A. B., et al. (2017). 2017 AHA/ACC/HRS Guideline for Management of Patients With Ventricular Arrhythmias and the Prevention of Sudden Cardiac Death: A Report of the American College of Cardiology/American Heart Association Task Force on Clinical Practice Guidelines and the Heart Rhythm Society. Circulation 138:e272–e391. doi: 10.1161/CIR.0000000000000549

CrossRef Full Text | Google Scholar

Balay, S., Abhyankar, S., Adams, M. F., Brown, J., Brune, P., Buschelman, K., et al. (2018). PETSc Users Manual. Argonne National Laboratory.

Bauer, S., Edelmann, J. C., Seemann, G., Sachse, F. B., and Dossel, O. (2013). Estimating intracellular conductivity tensors from confocal microscopy of rabbit ventricular tissue. Biomed. Tech. 58(Suppl. 1). doi: 10.1515/bmt-2013-4333. Available online at: https://www.degruyter.com/view/j/bmte.2013.58.issue-s1-N/bmt-2013-4333/bmt-2013-4333.xml

PubMed Abstract | CrossRef Full Text | Google Scholar

Bittihn, P., Horning, M., and Luther, S. (2012). Negative curvature boundaries as wave emitting sites for the control of biological excitable media. Phys. Rev. Lett. 109:118106. doi: 10.1103/PhysRevLett.109.118106

PubMed Abstract | CrossRef Full Text | Google Scholar

Cellml (2008a). Electrophysiological Modeling of Fibroblasts and Their Interaction with Myocytes [Online]. Available online at: https://models.cellml.org/exposure/4e478093480e0d7e08a2f24177ec4906. (Accessed October 1, 2018).

Cellml (2008b). A Rabbit Ventricular Action Potential Model Replicating Cardiac Dynamics at Rapid Heart Rates [Online]. Available online at: https://models.cellml.org/exposure/a5586b72d07ce03fc40fc98ee846d7a5 (Accessed October 1, 2018).

Cevik, C., Perez-Verdia, A., and Nugent, K. (2009). Implantable cardioverter defibrillators and their role in heart failure progression. Europace 11, 710–715. doi: 10.1093/europace/eup091

PubMed Abstract | CrossRef Full Text | Google Scholar

Connolly, A., Vigmond, E., and Bishop, M. (2017a). Virtual electrodes around anatomical structures and their roles in defibrillation. PLoS ONE 12:e0173324. doi: 10.1371/journal.pone.0173324

PubMed Abstract | CrossRef Full Text | Google Scholar

Connolly, A. J., Vigmond, E., and Bishop, M. J. (2017b). Bidomain predictions of virtual electrode-induced make and break excitations around blood vessels. Front Bioeng Biotechnol. 5:18. doi: 10.3389/fbioe.2017.00018

PubMed Abstract | CrossRef Full Text | Google Scholar

Corrias, A., Pathmanathan, P., Gavaghan, D. J., and Buist, M. L. (2012). Modelling tissue electrophysiology with multiple cell types: applications of the extended bidomain framework. Integr. Biol. 4, 192–201. doi: 10.1039/c2ib00100d

PubMed Abstract | CrossRef Full Text | Google Scholar

Costa, C. M., Campos, F. O., Prassl, A. J., Dos Santos, R. W., Sanchez-Quintana, D., Ahammer, H., et al. (2014). An efficient finite element approach for modeling fibrotic clefts in the heart. IEEE Trans. Biomed. Eng. 61, 900–910. doi: 10.1109/TBME.2013.2292320

PubMed Abstract | CrossRef Full Text | Google Scholar

Dosdall, D. J., Fast, V. G., and Ideker, R. E. (2010). Mechanisms of defibrillation. Annu. Rev. Biomed. Eng. 12, 233–258. doi: 10.1146/annurev-bioeng-070909-105305

PubMed Abstract | CrossRef Full Text | Google Scholar

Entcheva, E., Eason, J., Efimov, I. R., Cheng, Y., Malkin, R., and Claydon, F. (1998). Virtual electrode effects in transvenous defibrillation-modulation by structure and interface: evidence from bidomain simulations and optical mapping. J. Cardiovasc. Electrophysiol. 9, 949–961. doi: 10.1111/j.1540-8167.1998.tb00135.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Entcheva, E., Trayanova, N. A., and Claydon, F. J. (1999). Patterns of and mechanisms for shock-induced polarization in the heart: a bidomain analysis. IEEE Trans. Biomed. Eng. 46, 260–270. doi: 10.1109/10.748979

PubMed Abstract | CrossRef Full Text | Google Scholar

Fishler, M. G. (1998). Syncytial heterogeneity as a mechanism underlying cardiac far-field stimulation during defibrillation-level shocks. J. Cardiovasc. Electrophysiol. 9, 384–394. doi: 10.1111/j.1540-8167.1998.tb00926.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Fishler, M. G., and Vepa, K. (1998). Spatiotemporal effects of syncytial heterogeneities on cardiac far-field excitations during monophasic and biphasic shocks. J. Cardiovasc. Electrophysiol. 9, 1310–1324. doi: 10.1111/j.1540-8167.1998.tb00107.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Garny, A., Nickerson, D. P., Cooper, J., Weber Dos Santos, R., Miller, A. K., Mckeever, S., et al. (2008). CellML and associated tools and techniques. Philos. Trans. A Math. Phys. Eng. Sci. 366, 3017–3043. doi: 10.1098/rsta.2008.0094

PubMed Abstract | CrossRef Full Text | Google Scholar

Greiner, J., Sankarankutty, A. C., Seemann, G., Seidel, T., and Sachse, F. B. (2018). Confocal microscopy-based estimation of parameters for computational modeling of electrical conduction in the normal and infarcted heart. Front. Physiol. 9:239. doi: 10.3389/fphys.2018.00239

PubMed Abstract | CrossRef Full Text | Google Scholar

Greisas, A., and Zlochiver, S. (2016). The multi-domain fibroblast/myocyte coupling in the cardiac tissue: a theoretical study. Cardiovasc. Eng. Technol. 7, 290–304. doi: 10.1007/s13239-016-0266-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Henriquez, C. S. (1993). Simulating the electrical behavior of cardiac tissue using the bidomain model. Crit. Rev. Biomed. Eng. 21, 1–77.

PubMed Abstract | Google Scholar

Iles, L., Pfluger, H., Lefkovits, L., Butler, M. J., Kistler, P. M., Kaye, D. M., et al. (2011). Myocardial fibrosis predicts appropriate device therapy in patients with implantable cardioverter-defibrillators for primary prevention of sudden cardiac death. J. Am. Coll. Cardiol. 57, 821–828. doi: 10.1016/j.jacc.2010.06.062

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacquemet, V., and Henriquez, C. S. (2007). Modelling cardiac fibroblasts: interactions with myocytes and their impact on impulse propagation. Europace 9(Suppl. 6), vi29–37. doi: 10.1093/europace/eum207

PubMed Abstract | CrossRef Full Text | Google Scholar

Kieval, R. S., Spear, J. F., and Moore, E. N. (1992). Gap junctional conductance in ventricular myocyte pairs isolated from postischemic rabbit myocardium. Circ. Res. 71, 127–136. doi: 10.1161/01.RES.71.1.127

PubMed Abstract | CrossRef Full Text | Google Scholar

Kohl, P., and Gourdie, R. G. (2014). Fibroblast-myocyte electrotonic coupling: does it occur in native cardiac tissue? J. Mol. Cell. Cardiol. 70, 37–46. doi: 10.1016/j.yjmcc.2013.12.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Kohl, P., Kamkin, A. G., Kiseleva, I. S., and Noble, D. (1994). Mechanosensitive fibroblasts in the sino-atrial node region of rat heart: interaction with cardiomyocytes and possible role. Exp. Physiol. 79, 943–956. doi: 10.1113/expphysiol.1994.sp003819

PubMed Abstract | CrossRef Full Text | Google Scholar

Mahajan, A., Shiferaw, Y., Sato, D., Baher, A., Olcese, R., Xie, L. H., et al. (2008). A rabbit ventricular action potential model replicating cardiac dynamics at rapid heart rates. Biophys. J. 94, 392–410. doi: 10.1529/biophysj.106.98160

PubMed Abstract | CrossRef Full Text | Google Scholar

Mcdowell, K. S., Vadakkumpadan, F., Blake, R., Blauer, J., Plank, G., Macleod, R. S., et al. (2012). Methodology for patient-specific modeling of atrial fibrosis as a substrate for atrial fibrillation. J. Electrocardiol. 45, 640–645. doi: 10.1016/j.jelectrocard.2012.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Plonsey, R., and Barr, R. C. (2007). Bioelectricity: A Quantitative Approach. New York, NY: Springer.

Google Scholar

Pumir, A., and Krinsky, V. (1999). Unpinning of a rotating wave in cardiac muscle by an electric field. J. Theor. Biol. 199, 311–319. doi: 10.1006/jtbi.1999.0957

PubMed Abstract | CrossRef Full Text | Google Scholar

Rook, M. B., De Jonge, B., Jongsma, H. J., and Masson-Pevet, M. A. (1990). Gap junction formation and functional interaction between neonatal rat cardiocytes in culture: a correlative physiological and ultrastructural study. J. Membr. Biol. 118, 179–192. doi: 10.1007/BF01868475

PubMed Abstract | CrossRef Full Text | Google Scholar

Rook, M. B., Jongsma, H. J., and De Jonge, B. (1989). Single channel currents of homo- and heterologous gap junctions between cardiac fibroblasts and myocytes. Pflugers Arch. 414, 95–98. doi: 10.1007/BF00585633

PubMed Abstract | CrossRef Full Text | Google Scholar

Rook, M. B., Van Ginneken, A. C., De Jonge, B., El Aoumari, A., Gros, D., and Jongsma, H. J. (1992). Differences in gap junction channels between cardiac myocytes, fibroblasts, and heterologous pairs. Am. J. Physiol. 263, C959–977. doi: 10.1152/ajpcell.1992.263.5.C959

PubMed Abstract | CrossRef Full Text | Google Scholar

Sachse, F. B. (2004). Computational Cardiology: Modeling of Anatomy, Electrophysiology, and Mechanics. Berlin; New York, NY: Springer. doi: 10.1007/b96841

CrossRef Full Text | Google Scholar

Sachse, F. B., Moreno, A. P., and Abildskov, J. A. (2008). Electrophysiological modeling of fibroblasts and their interaction with myocytes. Ann. Biomed. Eng. 36, 41–56. doi: 10.1007/s10439-007-9405-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Sachse, F. B., Moreno, A. P., Seemann, G., and Abildskov, J. A. (2009). A model of electrical conduction in cardiac tissue including fibroblasts. Ann. Biomed. Eng. 37, 874–889. doi: 10.1007/s10439-009-9667-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Seemann, G., Sachse, F. B., Karl, M., Weiss, D. L., Heuveline, V., and Dössel, O. (2010). “Framework for modular, flexible and efficient solving the cardiac bidomain equation using PETSc,” in Progress in Industrial Mathematics at ECMI 2008, eds A. D. Fitt, J. Norbury, H. Ockendon, and E. Wilson (Berlin; Heidelberg: Springer), 363–369. doi: 10.1007/978-3-642-12110-4_55

CrossRef Full Text | Google Scholar

Seidel, T., Sankarankutty, A. C., and Sachse, F. B. (2017). Remodeling of the transverse tubular system after myocardial infarction in rabbit correlates with local fibrosis: a potential role of biomechanics. Prog. Biophys. Mol. Biol. 130, 302–314. doi: 10.1016/j.pbiomolbio.2017.07.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Shibukawa, Y., Chilton, E. L., Maccannell, K. A., Clark, R. B., and Giles, W. R. (2005). K+ currents activated by depolarization in cardiac fibroblasts. Biophys. J. 88, 3924–3935. doi: 10.1529/biophysj.104.054429

PubMed Abstract | CrossRef Full Text | Google Scholar

Sobie, E. A., Susil, R. C., and Tung, L. (1997). A generalized activating function for predicting virtual electrodes in cardiac tissue. Biophys. J. 73, 1410–1423. doi: 10.1016/S0006-3495(97)78173-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Tung, L. (1978). A Bi-Domain Model for Describing Ischemic Myocardial d-c Potentials. Ph.D. Thesis, MIT, Cambridge, MA.

Witwer, J. G., Trezek, G. J., and Jewett, D. L. (1972). The effect of media inhomogeneities upon intracranial electrical fields. IEEE Trans. Biomed. Eng. 19, 352–362. doi: 10.1109/TBME.1972.324138

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeigler, A. C., Richardson, W. J., Holmes, J. W., and Saucerman, J. J. (2016). Computational modeling of cardiac fibroblasts and fibrosis. J. Mol. Cell. Cardiol. 93, 73–83. doi: 10.1016/j.yjmcc.2015.11.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: defibrillation, cardiac tissue, fibrosis, computational modeling, multidomain modeling

Citation: Bragard J, Sankarankutty AC and Sachse FB (2019) Extended Bidomain Modeling of Defibrillation: Quantifying Virtual Electrode Strengths in Fibrotic Myocardium. Front. Physiol. 10:337. doi: 10.3389/fphys.2019.00337

Received: 17 December 2018; Accepted: 13 March 2019;
Published: 03 April 2019.

Edited by:

Javier Saiz, Universitat Politècnica de València, Spain

Reviewed by:

Bradley John Roth, Oakland University, United States
Edward Joseph Vigmond, Université de Bordeaux, France

Copyright © 2019 Bragard, Sankarankutty and Sachse. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Jean Bragard, jbragard@unav.es
Frank B. Sachse, frank.sachse@utah.edu