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

^{1}Department of Physics and Applied Mathematics, University of Navarra, Pamplona, Spain^{2}Nora Eccles Harrison Cardiovascular Research and Training Institute, University of Utah, Salt Lake City, UT, United States^{3}Department 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:

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 *f*_{s,myo}, *f*_{s,fib}, and *f*_{s,e} (*A/m*^{3}), are the stimulus current source densities for the myocyte, fibroblast and extracellular domain, respectively. Membrane voltage of myocytes *V*_{myo} and fibroblasts *V*_{fib} 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/m^{3}) were defined as:

with the volume fraction of myocytes *Vol*_{myo} and fibroblasts *Vol*_{fib} (dimensionless) as well as the individual volume for a myocyte *Vol*_{myo, single} and fibroblast *Vol*_{fib,single} (m^{3}). Currents between the myocyte and extracellular, fibroblast and extracellular, and myocyte and fibroblast domains are identified by *I*_{myo,e}, *I*_{fib,e}, and *I*_{myo,fib}*(A)*, respectively.

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

where *R*_{myo,fib} is the resistance of gap junction channels (Ω). We defined the effective volume density β_{myo,fib} (1/m^{3}) of electrical connections between fibroblasts and myocytes as:

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

where the tensors ${\overline{\sigma}}_{myo}$, ${\overline{\sigma}}_{fib}$ and ${\overline{\sigma}}_{e}$ described conductivity for volume ratios of 100%. In the subsequently described simulations, we vary ${\overline{\sigma}}_{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 *I*_{myo,e} and *I*_{fib, 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 S1–S3 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 *R*_{myo,fib} and ${\overline{\sigma}}_{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/cm^{2} 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 × 10^{5} and −0.7 × 10^{5} A/m^{3}, 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**. 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.

In order to measure virtual electrode strengths, we monitored the difference in the membrane voltage of the myocytes *V*_{myo} and fibroblasts *V*_{fib} at patch boundaries during the application of extracellular current with respect to the membrane voltages before current application. *V*_{myo} 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 *V*_{myo} and *V*_{fib}, 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/cm^{2} 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 × 10^{5} and −0.7 × 10^{5} A/m^{3}, 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 *W*_{p} on Δ*V*_{myo} were analyzed using the fitting function:

with the parameters a and b. We used a different fitting function for Δ*V*_{fib}, because it was not always crossing the origin of axes:

with the parameters a, b and c. We measured the fit quality using the adjusted *R*^{2}.

## 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 *R*_{myo,fib} = 4 × 10^{4} MΩ and ${\overline{\sigma}}_{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 *V*_{myo} and *V*_{fib} at the application sites and the borders of the fibrotic patch.

**Figure 2**. Representative plots of **(A)** ϕ_{e}, **(B)** *V*_{myo}, and **(C)** *V*_{fib}for intracellular pacing and subsequent extracellular current application in the 1D model with the average fibrosis patch, ${\overline{\sigma}}_{fib}$ of 0.1 S/m and *R*_{myo,fib} of 40 GΩ. Spatial distribution of **(D)** ϕ_{e}, **(E)** *V*_{myo}, and **(F)** *V*_{fib} at the end of application of extracellular currents. The extracellular current led to alterations of *V*_{myo} and *V*_{fib} 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}, *V*_{myo}, and *V*_{fib} at the end (*t* = 2.255 s) of current application in Figures 2D–F, respectively. Due to current application, *V*_{myo} 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. *V*_{fib} 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 *V*_{myo} was decreased from −85.00 to −93.44 mV. *V*_{fib} 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 *R*_{myo,fib} and ${\overline{\sigma}}_{fib}$. Example simulated ϕ_{e}, *V*_{myo}, and *V*_{fib} are presented in Figure S3. Before current application, the biphasic relationship of *V*_{myo} with *R*_{myo,fib} was nearly identical on both boundaries of the patch (Figures 3A,C). Also, variation of ${\overline{\sigma}}_{fib}$ did not affect *V*_{myo}. After current application, the left border of the patch exhibited increased *V*_{myo}, measured as $\Delta {V}_{myo}^{(1)}$, dependent on *R*_{myo,fib} and ${\overline{\sigma}}_{fib}$ (Figure 3B). At the right border of the patch, the reduction of *V*_{myo}, measured as $\Delta {V}_{myo}^{(2)}$, was dependent on *R*_{myo,fib} and ${\overline{\sigma}}_{fib}$ (Figure 3D). At both borders of the patch, effects on *V*_{myo} were of similar magnitude for all ${\overline{\sigma}}_{fib}$ for large *R*_{myo,fib}.

**Figure 3**. Measurement of effects of extracellular current application on *V*_{myo} using the 1D model with the average fibrosis patch. Before application of extracellular currents, *V*_{myo} at the **(A)** left and **(C)** right border of the fibrotic patch was marginally affected by *R*_{myo,fib} and ${\overline{\sigma}}_{fib}$. **(B)** At the end of current application, the left patch border exhibited a positive Δ*V*_{myo}. **(D)** In contrast, the right patch border exhibited a negative Δ*V*_{myo}.

We present corresponding measures on *V*_{fib} in Figure 4. The biphasic relationship of *V*_{fib} with *R*_{myo,fib} before application of extracellular current (Figures 4A,C) was similar to the relationship of *V*_{myo} with *R*_{myo,fib}. For low *R*_{myo,fib} (high myocyte-fibroblast coupling), *V*_{fib} was similar to the resting *V*_{myo} (Figures 3A,C). For high *R*_{myo,fib} (low myocyte-fibroblast coupling), *V*_{fib}was 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, $\Delta {V}_{fib}^{(1)}$ and $\Delta {V}_{fib}^{(2)}$, dependent on *R*_{myo,fib} and ${\overline{\sigma}}_{fib}$ in response to current application (Figures 4B,D). For low *R*_{myo,fib}, Δ*V*_{fib} was roughly similar to the Δ*V*_{myo} at both borders. In contrast, for high *R*_{myo,fib}, Δ*V*_{fib} exhibited a reverse sign vs. the corresponding Δ*V*_{myo} at both borders. This results in the discontinuity in Δ*V*_{fib} for *R*_{myo,fib} of 100–1,000 MΩ (Figures 4B,D).

**Figure 4**. Measurement of effects of extracellular current application on *V*_{fib} using the 1D model with the average fibrosis patch. Before application of extracellular currents, *V*_{fib} at the **(A)** left and **(C)** right border of the fibrotic patch was strongly affected by *R*_{myo,fib} and ${\overline{\sigma}}_{fib}$. At the end of current application, the **(B)** left and **(D)** right border exhibited a Δ*V*_{fib} with magnitude and sign modulated by *R*_{myo,fib} and ${\overline{\sigma}}_{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 Δ*V*_{myo} and Δ*V*_{fib}. We observed conduction block for *R*_{myo,fib} of 100–1,000 MΩ (Figures 5A,C). We confirmed that the magnitude Δ*V*_{myo} 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 Δ*V*_{myo} with α was roughly linear. Similarly, for low *R*_{myo,fib} the magnitude Δ*V*_{fib} was a monotonously increasing function of α (Figures 7C,D). For high *R*_{myo,fib}, the magnitude Δ*V*_{fib} was approximately constant for α between 0 and 1.

**Figure 5**. Measurement of effects of extracellular current application on *V*_{myo} using the 1D model with the large fibrosis patch. Before application of extracellular currents, *V*_{myo} at the **(A)** left and **(C)** right border of the fibrotic patch was marginally affected by *R*_{myo,fib} and ${\overline{\sigma}}_{fib}$. **(B)** At the end of current application, the left patch border exhibited a positive Δ*V*_{myo}. **(D)** In contrast, the right patch border exhibited a negative Δ*V*_{myo}.

**Figure 6**. Measurement of effects of extracellular current application on *V*_{fib} using the 1D model with the large fibrosis patch. Before application of extracellular currents, *V*_{fib} at the **(A)** left and **(C)** right border of the fibrotic patch was strongly affected by *R*_{myo,fib} and ${\overline{\sigma}}_{fib}$. At the end of current application, the **(B)** left and **(D)** right border exhibited a Δ*V*_{fib} with magnitude and sign modulated by *R*_{myo,fib} and ${\overline{\sigma}}_{fib}$.

**Figure 7**. Assessment of effects of extracellular current application on *V*_{myo} and *V*_{fib} 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 *R*_{myo,fib}. ${\overline{\sigma}}_{fib}$ was set to 0.1 S/m. The effects on the magnitude of Δ*V*_{myo} at the **(A)** left and **(B)** right boundary increased with the degree of fibrosis. Similarly, effects of Δ*V*_{fib} at the **(C)** left and **(D)** right boundary increased with the degree of fibrosis for small *R*_{myo,fib}. Marginal effects were present for larger *R*_{myo,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 Δ*V*_{myo} and Δ*V*_{fib} 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 *R*_{myo,fib}. For large *R*_{myo,fib}, the spatial length scales for the myocytes and fibroblasts correspond to their respective space constants λ_{myo} and λ_{fib}. In contrast, for low *R*_{myo,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 Δ*V*_{myo} and Δ*V*_{fib} were at the sub-millimeter scale (Figure 8), which is in agreement with the estimated values for of λ_{myo} and λ_{fib} (Tables S1–S3).

**Figure 8**. Effects of patch length W_{p} on **(A,B)** *V*_{myo} and **(C,D)** *V*_{fib} at the boundary of the patch in response to application of the extracellular current. We varied *R*_{myo,fib}, but ${\overline{\sigma}}_{fib}$ was set on 0.1 S/m. For small W_{p}, *V*_{myo}, and *V*_{fib} were greatly reduced.

### 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 Δ*V*_{myo} and Δ*V*_{fib} 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, ${\overline{\sigma}}_{fib}$ of 0.1 S/m and *R*_{myo,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 *V*_{myo} 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 *V*_{fib}. In this example, *R*_{myo, 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 *V*_{fib} at point 1 and depolarization at point 2, which is the opposite of what happens to *V*_{myo}. However, the coupling term between the two domains caused that the *V*_{myo} is pulling the *V*_{fib} in the same direction. At the very edge of the patch (Figure 9F) the two effects cancel out.

**Figure 9**. Spatial distribution of **(A)** ϕ_{e}, **(B)** *V*_{myo}, and **(C)** *V*_{fib} at the end of extracellular current application in the 2D model with the square-shaped average fibrosis patch, ${\overline{\sigma}}_{fib}$ of 0.1 S/m and *R*_{myo,fib} of 1 GΩ. Close-ups of the spatial distribution for **(D)** ϕ_{e}, **(E)** *V*_{myo}, and **(F)** *V*_{fib}. The extracellular current led to alterations of *V*_{myo} and *V*_{fib} 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 *R*_{myo,fib} = 1 TΩ, which represents negligible coupling between *V*_{fib} and *V*_{myo} (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 *V*_{myo}. 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 Δ*V*_{fib} are indeed opposite to the Δ*V*_{myo} at the patch boundary. Also, the spatial scale for the virtual electrode was much smaller than the virtual electrodes for the *V*_{myo} field, which agrees with the corresponding values for the space constants associated with the two fields.

**Figure 10**. Spatial distribution of **(A)** ϕ_{e}, **(B)** *V*_{myo}, and **(C)** *V*_{fib} at the end of extracellular current application in the 2D model with the disk-shaped average fibrosis patch, ${\overline{\sigma}}_{fib}$ of 0.1 S/m and *R*_{myo,fib} of 1 TΩ. Close-ups of the spatial distribution for **(D)** ϕ_{e}, **(E)** *V*_{myo}, and **(F)** *V*_{fib} at the end of application of extracellular currents. The extracellular current led to alterations of *V*_{myo} and *V*_{fib} 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 Δ*V*_{myo} and Δ*V*_{fib} at both locations of the patch boundaries as a function of *R*_{myo, 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 Δ*V*_{myo} 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 Δ*V*_{myo} for low values of the *R*_{myo,fib} and an approximate 5% difference for high values of the *R*_{myo,fib}. Also, Δ*V*_{myo} for the square-shaped were higher than for the disk-shaped patch, and close to the 1D patch geometry. Δ*V*_{fib} for the 1D and the two 2D geometries were similar (Figures 11C,D). For Δ*V*_{fib}, the influence of patch geometry was less pronounced.

**Figure 11**. Comparison of **(A,B)** *V*_{myo} and **(C,D)** *V*_{fib} 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):

with the injected current *f*_{e} due to extracellular current application, the myocyte membrane capacitance *C*_{myo} and the fibroblast membrane capacitance *C*_{fib}. 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):

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:

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 *V*_{myo} (hyperpolarization). If ∇σ_{myo} points in the opposite direction of *E*, the current application leads to an increase of *V*_{myo} (depolarization). We extend this description for *V*_{fib} using a decomposition of the second right hand term of Equation (14):

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 *V*_{myo}. On the contrary, ∇σ_{fib} at point 1 is aligned with E, which results in a hyperpolarization of *V*_{fib} at point 1.

However, this simple mechanistic explanation ignores the coupling term between the fibroblast and the myocyte domain. The current *I*_{myo, fib}, which is determined by the parameter *R*_{myo, fib} and the voltage between the domains, explains why in some cases, *V*_{fib} follows *V*_{myo}. In particular, at high coupling between myocyte and fibroblast domain (i.e., low values of *R*_{myo,fib}), *V*_{myo} is very similar to *V*_{fib}. Interestingly, at intermediate values of *R*_{myo,fib}, the simulations revealed a competition of the two antagonistic effects for *V*_{fib} (Figures 4, 6). These two effects can even cancel each other out in a sort of tug of war situation, leaving *V*_{fib} 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, *V*_{fib} and *V*_{myo} 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 *V*_{fib} 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, *V*_{fib} could pull *V*_{myo} 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 Δ*V*_{fib} during the extracellular current application. In Figures 4B,D, the reported Δ*V*_{fib} is the largest absolute value (keeping its sign) of the monitored difference. These indicators exhibited a discontinuous behavior for intermediate values for the *R*_{myo,fib} as exemplified in Figures 4B,D. The discontinuity is caused by the transition between dominant minima and maxima of *V*_{fib} as illustrated in Figure S2B. The transition occurred at decreased *R*_{myo,fib}, when increasing ${\overline{\sigma}}_{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 *V*_{myo}. In Figures 5A,C, discontinuities occurred at intermediate values of *R*_{myo, 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 *V*_{myo}. Though this effect is important for the dynamics of the system, in the context of our study it did affect only marginally Δ*V*_{fib} and Δ*V*_{myo}.

In our 2D simulations, the influence of patch geometry on Δ*V*_{fib} 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 *R*_{myo, 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 *R*_{myo, 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 *R*_{myo, 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 *R*_{myo, 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, *R*_{myo, 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 *R*_{myo, fib} limit (1 MΩ) in our studies appears small. We note that effects of *R*_{myo, 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 *R*_{myo, fib} and virtual electrodes.

We varied ${\overline{\sigma}}_{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 ${\overline{\sigma}}_{myo}$. Dependent on *Vol*_{fib} (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 *Vol*_{fib} 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 ${\overline{\sigma}}_{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 *Vol*_{fib} 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/m^{3}) of electrical connections between fibroblasts and myocytes is:

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, *R*_{myo, 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:

### 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 *Vol*_{fib} (Equation 9). Non-linear models might be more realistic. In particular, for small *Vol*_{fib}, 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 *V*_{myo} and *V*_{fib} 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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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, SpainReviewed by:

Bradley John Roth, Oakland University, United StatesEdward 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