Frontiers in Physiology

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:

$\begin{array}{rc}\nabla ·({\sigma }_{myo}\nabla {\varphi }_{myo})={-f}_{s,myo}+{\beta }_{myo,fib}{I}_{myo,fib}+{\beta }_{myo}{I}_{myo,e}& (1)\end{array}$

$\begin{array}{rc}\nabla ·({\sigma }_{fib}\nabla {\varphi }_{fib})={-f}_{s,fib}-{\beta }_{myo,fib}{I}_{myo,fib}+{\beta }_{fib}{I}_{fib,e}& (2)\end{array}$

$\begin{array}{rc}\nabla ·({\sigma }_e\nabla {\varphi }_e)={-f}_{s,e}-{\beta }_{myo}{I}_{myo,e}-{\beta }_{fib}{I}_{fib,e}& (3)\end{array}$

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³), 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³) were defined as:

$\begin{array}{rc}{\beta }_{myo}=\frac{Vo{l}_{myo}}{Vo{l}_{myo,single}}& (4)\end{array}$

$\begin{array}{rc}{\beta }_{fib}=\frac{Vo{l}_{fib}}{Vo{l}_{fib,single}}& (5)\end{array}$

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³). 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:

$\begin{array}{rc}{I}_{myo,fib}=\frac{{\varphi }_{myo}-{\varphi }_{fib}}{R_{myo,fib}}& (6)\end{array}$

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

$\begin{array}{rc}{\beta }_{myo,fib}=\frac{Vo{l}_{fib}}{Vo{l}_{fib,single}}& (7)\end{array}$

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

$\begin{array}{rc}{\sigma }_{myo}=Vo{l}_{myo}{\overline{\sigma }}_{myo}& (8)\end{array}$

$\begin{array}{rc}{\sigma }_{fib}=Vo{l}_{fib}{\overline{\sigma }}_{fib}& (9)\end{array}$

$\begin{array}{rc}{\sigma }_e=Vo{l}_e{\overline{\sigma }}_e& (10)\end{array}$

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⁻¹². 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² 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⁵ and −0.7 × 10⁵ A/m³, 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

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

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 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² 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⁵ and −0.7 × 10⁵ A/m³, 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:

$\begin{array}{rc}\Delta V_{myo}=a\left(1-e^{-\frac{W_p}{b}}\right)& (11)\end{array}$

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:

$\begin{array}{rc}\Delta V_{fib}=c-a{e}^{-\frac{W_p}{b}}& (12)\end{array}$

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

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⁴ 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

Figure 2. Representative plots of (A) ϕ_e, (B) V_myo, and (C) V_fibfor 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

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_fibwas 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

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

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

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

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

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.

TABLE 2

Table 2. Relationships of membrane voltage changes and W_p 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 Δ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

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

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

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):

$\begin{array}{l}\frac{\partial V_{myo}}{\partial t}=\frac{1}{C_{myo}}(\nabla ·\left({\sigma }_{myo}\nabla V_{myo}\right)+\nabla ·\left({\sigma }_{myo}\nabla {\varphi }_e\right)\\ +\cdots )& (13)\end{array}$

$\begin{array}{l}\frac{\partial V_{fib}}{\partial t}=\frac{1}{C_{fib}}(\nabla ·\left({\sigma }_{fib}\nabla V_{fib}\right)+\nabla ·\left({\sigma }_{fib}\nabla {\varphi }_e\right)\\ +\cdots )& (14)\end{array}$

$\begin{array}{l}\nabla ·\left\{({{\sigma }_e+{\sigma }_{myo}+\sigma }_{fib}\nabla {\varphi }_e)\right\}={-f}_e-\nabla ·({\sigma }_{myo}\nabla V_{myo})\\ -\nabla ·({\sigma }_{fib}\nabla V_{fib})& (15)\end{array}$

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):

$\begin{array}{rc}E=-\nabla {\varphi }_e& (16)\end{array}$

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:

$\begin{array}{rc}\frac{\partial V_{myo}}{\partial t}\propto -(\nabla {\sigma }_{myo})·E+{\sigma }_{myo}{\nabla }^2{\varphi }_e& (17)\end{array}$

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):

$\begin{array}{rc}\frac{\partial V_{fib}}{\partial t}\propto -(\nabla {\sigma }_{fib})·E+{\sigma }_{fib}{\nabla }^2{\varphi }_e& (18)\end{array}$

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³) of electrical connections between fibroblasts and myocytes is:

$\begin{array}{rc}{\beta }_{myo,fib}={\beta }^{(r)}\frac{Vo{l}_{fib}}{Vo{l}_{fib,single}}& (19)\end{array}$

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:

$\begin{array}{rc}{\beta }^{(fract)}=\frac{{\beta }^{(r)}}{R_{myo,fib}}& (20)\end{array}$

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

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