Fibroblast mediated dynamics in diffusively uncoupled myocytes: a simulation study using 2-cell motifs

In healthy hearts myocytes are typically coupled to nearest neighbours through gap junctions. Under pathological conditions such as fibrosis, or in scar tissue, or across ablation lines myocytes can uncouple from their neighbours. Electrical conduction may still occur via fibroblasts that not only couple proximal myocytes but can also couple otherwise unconnected regions. We hypothesise that such coupling can alter conduction between myocytes via introduction of delays or by initiation of premature stimuli that can potentially result in reentry or conduction blocks. To test this hypothesis we have developed several 2-cell motifs and investigated the effect of fibroblast mediated electrical coupling between uncoupled myocytes. We have identified various regimes of myocyte behaviour that depend on the strength of gap-junctional conductance, connection topology, and parameters of the myocyte and fibroblast models. These motifs are useful in developing a mechanistic understanding of long-distance coupling on myocyte dynamics and enable the characterisation of interaction between different features such as myocyte and fibroblast properties, coupling strengths and pacing period. They are computationally inexpensive and allow for incorporation of spatial effects such as conduction velocity. They provide a framework for constructing scar tissue boundaries and enable linking of cellular level interactions with scar induced arrhythmia.


Introduction
The heart is an electro-mechanical pump, whose coordinated mechanical contraction is enabled by the propagation of synchronized waves of electrical excitation.In the mammalian heart, myocytes and fibroblasts constitute two of the most important type of cells, with the larger myocytes being responsible for cardiac electrical activity and the smaller but more numerous fibroblasts maintaining the electro-mechanical integrity of the heart [1].Fibroblasts play a crucial role in the repair of heart muscles especially post injury or disease [1,2].In aged or diseased hearts the number of fibroblasts may increase substantially (up to 40% [3]) resulting in increased collagen deposition causing fibrosis.
Furthermore fibroblasts themselves differentiate into much larger myofibroblasts in injured and diseased hearts.In this paper we have used the terms fibroblasts and myofibroblasts interchangeably even though we are mostly referring to myofibroblasts.For a long time fibroblasts were not believed to influence the electrical conduction in heart muscles and even today the exact nature of myocyte-fibroblast coupling in vivo is debated [2,4].However several in vitro experiments have reported the existence of gap-junctional coupling between myocytes and fibroblasts under both physiological and pathological conditions [5,6,7,8].Experiments have shown that coupling between myocytes and fibroblasts can significantly alter the conduction properties of the tissue [5,9], excitability of myocytes [10] and their resting membrane potential [6].At the level of the tissue and organ, fibrosis is known to substantially affect wave propagation in the heart and it is well understood that fibrosis can create a substrate for cardiac arrhythmia [11,12,13,14,15,16].
Several in silico studies have investigated how fibroblasts alter the electrical activity of individual myocytes and tissue under both normal and pathological conditions [17,18,19,20,21,22,23].These studies have usually simulated fibroblasts as being either attached to a myocyte or inserted in a tissue of coupled myocytes thereby coupling nearest neighbours.However heterocellular cell culture experiments have shown that fibroblasts can enable conduction up to 300µm [5].In vivo fibroblasts have been observed to form large sheet-like extensions having additional folds and elongated cytoplasmic processes [2].In the sinoatrial node it has been observed that an individual fibroblast could form membrane juxtaposition with a nearby myocytes covering up to 720 µm 2 [2,24].Fibroblasts that have such long extensions can potentially couple with multiple myocytes that are spatially distant.
Such long range interactions between distant myocytes mediated via fibroblasts have the potential to modify tissue electrophysiology and dynamics via conduction delays and by exciting resting and partially recovered regions.The possibility of such long range coupling increases in diseased or injured hearts where there are a larger number of myofibroblasts.Ablation lines are another scenario where the fibroblasts that aid the repair of the ablation scars could produce conduction pathways between regions that have been electrically isolated [25].Another possible scenario occurs when islands of myocytes are trapped in a scar within a sea of fibroblasts resulting in conduction pathways that connect distant uncoupled regions [26,27].
Most computational studies have approached the problem of fibroblast induced dynamics in terms of the effect of the fibroblast distribution, their density [22,23,28] and texture [29,30].On the other hand in this paper we describe a bottom-up approach by considering several simple motifs consisting of 1 or 2 fibroblast units coupled to a pair of mutually uncoupled myocytes.The idea of motifs has previously been used successfully to understand networks in complex biological domains including gene transcription, biochemical systems, neuronal networks [31,32,33].In this paper we have extended the idea of motif to model complex interactions between myocytes and fibroblasts at the scar boundaries.
We hypothesise that fibroblast mediated coupling that connects diffusively uncoupled myocytes can potentially alter conduction between myocytes via the introduction of delays or by the initiation of premature stimuli and can potentially result in reentry or conduction blocks in tissue.We have used the different motifs to investigate the effect of different configurations of fibroblast mediated electrical coupling between mutually uncoupled myocytes.We have also determined the effect of the different motif features such as connection topology, coupling strength, pacing period and cell parameters on the dynamics of the individual myocytes.
The 2-cell motifs allow for easy comparison of the influence of the intrinsic myocyte dynamics (characterized by its restitution) with the local coupling topology (characterized by the individual motifs) and pacing period.The advantage of using motifs is that they allow for a computationally inexpensive approach to study the effect of individual variation of the many features of the myocyte-fibroblast coupled system.By incorporating delay in the stimulation times of individual myocytes of the motif, we study the combined effect of conduction velocity in the tissue and the conduction delay arising from the coupling of disconnected myocytes via fibroblasts.Furthermore we have characterized the regimes of myocyte dynamics that can arise from this kind of fibroblast mediated long-range conduction.

Methods
The electrical activity of myocytes was described using the T N N P − T P 06 model of human ventricular cells [34,35], while the electrophysiological properties of the fibroblasts were described using the M acCannell "active" fibroblast model [36].The time variation of the transmembrane voltage V for myocytes coupled to n fibroblasts was given as, Here I ion describes the total of all ionic currents: ) where I N a is the sodium current, I to is the transient outward current, I K1 , I Kr and I Ks are the inward rectifier, delayed rectifier and slow delayed rectifier potassium currents, I CaL is the L-type Ca 2+ current, I N aK is the N a + /K + pump current, I N aCa is the N a + /Ca 2+ exchanger current, I pCa and I pK plateau calcium and potassium currents and I bCa and I bN a are the background N a + and Ca 2+ currents.V f is the fibroblast transmembrane potential while G gap is the strength of the gap junctional coupling between myocyte and fibroblast.
The M acCannell fibroblast model equations are used to describe the time evolution of the fibroblast membrane potential V f (similar to Equation. 1 and Equation.2) with the ionic currents comprised of inward rectifying potassium current I f K1 , the time-and voltage -dependent potassium currents I f Kv , I f N aK a sodium-potassium pump current and a background sodium current I bN a + .C m and C f are the cell capacitance per unit surface area of myocyte and fibroblast set to 150pF and 50pF (corresponding to the larger myofibroblast [37]) respectively.
For the myocytes, we used two parameter sets corresponding to Shallow and Steep restitution slopes (see Table 2, slope = 0.7 for Shallow and slope = 1.8 for Steep in ten Tusscher et al [35]).The uncoupled fibroblast resting membrane potential V F R were set to either −24.5 mV or −49.0 mV [12].Most of the results described in the paper were obtained with V F R set to −24.5 mV.A subset of results with V F R = −49.0mV is described in the supplementary material.The different resting membrane potentials were obtained by shifting the gating variable voltage dependence of the time dependent potassium current [19].
The action potential of the uncoupled myocyte stimulated at a period of 600 ms for both Shallow and Steep parameters is shown in Fig. 1(a) while Fig. 1(b,d) show the effect of weak fibroblast coupling on the myocyte and fibroblast trans-membrane potentials respectively.In Fig. 1(c) we have plotted the restitution curve for the Shallow and Steep uncoupled myocytes.
In our study we considered two identical myocytes coupled either to 1 or 2 fibroblast units via different connection motifs (Fig. 2).Each fibroblast unit represents a fixed number of fibroblasts (see below) connected in parallel and coupled to a neighbouring myocyte (similar to Figure .1 in [36] and Figure .2 in [38]).The different topology considered are (i) M otif − 1 where one fibroblast unit is coupled to two myocytes, (ii) M otif − 2 where there is an asymmetry in the number of passive cells connected to each myocyte and (iii) M otif − 3 where two fibroblast units are coupled in a symmetric fashion with the two myocytes.In the following text we have referred to each fibroblast unit as a single fibroblast (viz., Fib 1 and Fib 2 in Fig. 2).The myocyte-fibroblast interaction is modelled as occurring via a gap-junction.We considered two kinds of gap-junction conductance viz., G Loc and G Long corresponding to electrical conduction between fibroblast and the proximal and distal myocytes respectively.Note that for the results reported here the coupling strength of a fibroblast to a proximal myocyte G Loc , if a coupling exists between them, is the same for both fibroblast units.Similarly if coupling exists the coupling strength of a fibroblast to a distal myocyte G Long is the same for both myocyte fibroblast units.For the gap junctional coupling we have chosen values falling in the range (0−4 nS) that is considered to be representative of the effect of fibroblasts in cell-cultures [19].
We captured the effect of conduction velocity in tissue by the introduction of a delay (τ D ) in the stimulation of the distal myocyte (i.e., the myocyte coupled to a fibroblast with a coupling strength of G Long ).While propagation delays of 11 − 68 ms have been observed in heterocellular culture [5], in our study we performed simulations with two representative time delays of τ D = 10 ms and τ D = 25 ms.The resulting equation for the coupled system M otif − 3 (Fig. 2 (c)) is: where i = 1, 2, j = 1, 2 and i = j.Here V i and V fi represent the transmembrane voltage of the i th myocyte and i th fibroblast respectively, while N f represents the fibroblast density (number of fibroblasts in a fibroblast unit).The equations for M otif −1 are obtained by setting G Loc1 and G Long2 to zero.Similarly M otif −2 equations are obtained by setting G Long2 = 0 nS.The coupled ordinary differential equations were solved using an explicit adaptive time step method with a maximum time-step of 0.001 ms [39].In order to obtain an action potential, either or both the myocytes were individually stimulated every T ms by applying an external stimulus of strength −52 µA/mm 2 for a duration of 1 ms.While many of the results described in the paper were obtained for T = 600 ms, the effect of pacing period was also investigated by setting T = 500, 400 and 300 ms respectively.The fibroblast density N f was set to 4. Each simulation was performed by stimulating the myocytes 20 times, the first 10 action potentials generated were ignored and the mean action potential duration and fraction of action potentials in the Pacing and Response cells were calculated over the last 10 action potentials alone.The initial conditions of the variables describing both myocytes and fibroblasts were set to their uncoupled resting membrane values.

Effect of fibroblast coupling on APD
In order to characterize the effect of connection topology and strength of myocyte-fibroblast coupling, we plotted 2-parameter portraits for the mean APD of the myocytes.The two parameters G Loc and G Long were varied over the range 0 − 4 nS in steps of 0.5 nS. Figure .3 shows the effect of strength of coupling on mean APD for both myocyte − 1 (Fig. 3 (a,c,e) and myocyte − 2 (Fig. 3 (b,d,f)) for the different motifs with Steep restitution parameters.(See Supplementary Fig. 10 for the equivalent figure for Shallow parameters).
While the general effect of fibroblast coupling was to reduce the myocyte AP D, there was significant variation across motifs.M otif − 1 was the simplest motif with one fibroblast coupled to two myocytes, and showed a reduction in APD of nearly 70 ms for the case of the strongest coupling (G Loc = G Long = 4.0 nS) as compared to the case of no coupling (Fig. 3 (a,b)).However the largest reduction in APD in myocyte − 1 occurred for the case of (G Loc , G Long ) = (0.0, 4.0) nS (For myocyte − 2, (G Loc , G Long ) = (4.0,0.0) nS).The reduction in APD for Shallow parameters was smaller than that of the Steep parameters, with a maximum reduction of APD compared to that of no coupling being around 44 ms (see Supplementary Fig. 10).In the case of Shallow parameters the maximum reduction occurred for the case of (G Loc = G Long = 4.0 nS).
We next considered the effect of APD for M otif − 2 which has a structural asymmetry with each myocyte exposed to a different number of fibroblasts.Figure . 3(c,d) shows the effect of M otif −2 coupling on APD for both myocyte− 1 and myocyte − 2 for the Steep parameters.Due to the asymmetry in their coupling the magnitude of APD was not the same for both myocytes for a given (G Loc , G Long ) pair.However APD of both myocytes became equal as G Long decreased.The reduction of APD in myocyte−1 was greater than in myocyte−2.This motif too showed a greater reduction in APD for the Steep parameters as compared to the Shallow parameters (Supplementary Fig. 10(c,d)).M otif − 3 is symmetric in terms of coupling and the reduction in APD for a given coupling strength was the same for both myocytes, i.e. the difference in APD between the myocytes was zero.However the decrease in APD with coupling strength was larger for the Steep parameter set (Fig. 3 (e,f)) compared to the Shallow parameter case(Supplementary Fig. 10 (e,f)).Comparing across the different motifs we find that the reduction in APD with coupling increased with complexity of the topology, with the largest reduction in APD occurring for the case of M otif − 3.
For a representative set of coupling strengths we have compared the action potential profiles of both myocyte−1 and myocyte−2 in M otif −2 for both Shallow and Steep parameters in Fig. 4. For each of the motifs, there was significant change in features of the myocyte action potential shape including peak, dome and duration of recovery depending on the coupling strengths (G Loc , G Long ).The corresponding action potential profiles for M otif − 1 and M otif − 3 are plotted in the Supplementary section (Fig. 11 and Fig. 12).

Effect of delay τ D in stimulation
In the result described above, both myocytes were stimulated simultaneously.To capture the effect of conduction velocity across tissue we stimulated myocyte−2 with a delay after the stimulation of myocyte − 1. Supplementary Fig. 14  We next considered the effect of delay for the case of M otif − 2. In the case of zero delay in stimulation for Shallow parameters myocyte − 2 showed a greater APD than myocyte − 1 only for G Loc ≤ 1nS.For stronger local coupling (G Loc > 1 nS), ∆AP D was positive indicating that the reduction in APD of myocyte − 2 due to coupling was greater than the reduction in that of myocyte − 1 (Supplementary figure Fig. 13 (d)).However for the case of Steep parameters with no delay, ∆AP D ≥ 0 only for G Long = 0 or when G Loc was large (Fig. 5(d)).For small G Loc the reduction in APD was greater for myocyte − 1. Furthermore the magnitude of ∆AP D was more negative for Steep compared to the Shallow parameter set.While for stronger coupling, the magnitude of ∆AP D became more positive for longer delays for both Shallow and Steep parameters, with a greater increase for Shallow compared to Steep parameters (Supplementary Fig. 13 (e,f)).On the other hand for the Steep case (Fig. 5 (e,f), the reduction in APD of myocyte − 1 compared to myocyte − 2 for small G Loc increased with the magnitude of delay.However for all values of delay, ∆AP D < 0 ms occurred only when G Loc ≤ 1.5 nS.
In M otif − 3, for a delay of 0 ms both the myocytes had the same AP D as seen in Fig. 5 (g).For a delay of 10 ms, both Steep (Fig. 5(h)) and Shallow parameters (Supplementary figure Fig. 13(h)) show a reduction in APD for myocyte − 2 (corresponding to a positive value of ∆AP D) for all pairs of (G Loc , G Long ).However for a delay of 25 ms alone the Shallow parameter set showed a reduction in AP D in myocyte − 2 for all conductance pairs (Supplementary Fig. 13(i)).On the other hand for the Steep parameter set (Fig. 5(i

Initiating action potentials in a resting cell via fibroblasts
We next considered the case of infinite delay in stimulating one of the cells.In other words only one myocyte in a motif was stimulated (Pacing Cell) and the effect of this stimulus on the other myocyte (Response Cell) via the fibroblast was determined.The goal was to identify regions in the 2-parameter conductance space that would result in excitation of the quiescent Response Cell.In this experiment, the Pacing Cell was stimulated at a fixed period and the effect on the Response Cell was characterized in terms of the fraction of stimuli in the Pacing Cell that elicited a response in the Response Cell.
In Figure .6 (a-b), the black region (N R) corresponds to conductance value (small G Loc or G Long ) that did not produce any excitation in the Response Cell, even though the Pacing Cell produced an action potential for every applied stimulus.However for sufficiently strong coupling (larger values of G Loc and G Long ) every excitation in the Pacing Cell, was followed by an excitation in the Response Cell.This parameter region was marked 1 : 1. Between N R and 1 : 1 regions, there was an intermediate region of parameter values where the Response Cell did not respond to every excitation of Pacing Cell.Rather there was an intermediate response (IR), with only some of the action potentials in Pacing Cell resulting in an action potential in the Response Cell.
Figure .6(c-f) shows the action potentials for the set of parameter points marked c, d, e, f on Fig. 6(a).The initiation of excitation at a distal (not directly coupled) myocyte was critically dependent on the (G Loc , G Long ) conductance values in the motif.For G Long = 0.5 nS there was no response (Fig. 6)(c).When the value of G Long was further increased to 1.0, a shorter action potential was elicited in Response Cell for every stimulation of the Pacing Cell (Fig. 6)(d).However further increase of G Long to 1.5 nS gave rise to more complex intermittent dynamics in the Response Cell (Fig. 6)(e).Finally for G Long ≥ 2.5 nS, the Pacing Cell elicited a 1 : 1 response in the resting cell (Fig. 6)(f).
For the cases where there was a depolarization in the Response Cell, there was always a time delay with respect to the depolarization of Pacing Cell.However this delay was a function of the conductance strength (as seen in (Fig. 6)(d,f).The action potential profiles in (Fig. 6)(d,f) suggest a reduction in APD of the Response Cell compared to the Pacing Cell.Furthermore for large conductance there was an overlap in the repolarization time-series of both myocytes suggesting a synchronization of the recovery process in both the myocytes.Similar 2-parameter plots for the Pacing and Response cells are plotted for both M otif − 1 and M otif − 2 (Supplementary figures Fig. 15 and Fig. 16).
Since M otif −2 had an asymmetry, we considered two scenarios, viz., stimulating either myocyte − 2 (Supplementary figures Fig. 15(c,d) and Fig. 16 (c,d)) or myocyte − 1 (Supplementary figures Fig. 15(e-f) and Fig. 16(e-f)).Stimulating only myocyte − 2 could initiate action potentials in myocyte − 1 for a large range of conductance for both Steep and Shallow parameter sets.While for both Steep and Shallow parameters, small G Loc or G Long values did not produce a response in myocyte − 1, at larger conductance 1 : 1 response was obtained for both the parameter sets.The region of intermediate response (IR) was larger for Steep compared to Shallow.On the other hand a very small number of conductance pairs could initiate action potential in myocyte − 2 when myocyte − 1 alone was stimulated for the Shallow parameter set.For the Steep parameter set, 1 : 1 response was not observed for any of the conductance values.
In order to identify the effect of the pacing period on the initiation of action potentials via fibroblasts in a resting cell, we performed the above simulations for pacing periods T = 500 ms, T = 400 ms and T = 300 ms.The parameter space diagrams for T = 500 ms and T = 400 ms (not shown) were found to be qualitatively similar to those obtained for T = 600 ms, with only the boundaries between the different regimes and size of regimes varying marginally depending on the pacing period.On the other hand pacing at T = 300 ms produced regimes that were spatially more patchy suggesting sensitive dependence to coupling strength at very rapid pacing.Furthermore at T = 300 ms pacing IR regimes were also observed in the Pacing Cell implying that at very rapid pacing 1 : 1 response was not always guaranteed especially when coupled to other cells that can act as a current sink.For T = 300, supplementary figures Fig. 17 and Fig. 18 describe the different dynamical regimes obtained by stimulating either myocyte − 1 or myocyte − 2 for the case of M otif − 2. Supplementary Fig. 19 and Fig. 20 describe the different regimes obtained by pacing one cell with T = 300 ms for M otif − 1 and M otif − 3.
The 2-parameter plots are useful in identifying the different dynamical regimes and their relation to myocyte-fibroblast coupling strength.However in order to determine the influence of each of the features (viz., connection topology, myocyte parameters, pacing period) on the myocyte dynamics, we determined the fraction of instances for every regime, keeping one feature fixed at a time.In Fig. 7 we plot the fraction of occurrence of each regime for Steep and Shallow parameters individually, summing across the other features for both Pacing Cell (P Cell) and Response Cell (RCell).Similarly Fig. 8 and Fig. 9 describe the fraction of occurrence of the different regimes for the different pacing periods and motifs respectively.
Based on Fig. 7, Fig. 8 and Fig. 9, we can summarise that (i) Obtaining 1 : 1 response was more difficult with Steep parameters compared to Shallow.On the other hand Steep showed IR for more parameters than Shallow.(ii) With increase in pacing period, there were fewer instances of 1 : 1 and more We also repeated our simulations with a more negative fibroblast resting membrane potential V F R = −49.0mV for M otif − 2 and M otif − 3. The results for a subset of these simulations are described in the Supplementary section.
A more negative fibroblast resting membrane potential resulted in a greater one myocyte was stimulated for Steep and Shallow regimes respectively.While the different dynamical regimes are qualitatively similar to those obtained in V F R = −24.5mV, it can be observed that due to the more negative fibroblast resting potential a stronger coupling was required to elicit an action potential in the Response Cell.

Discussion
In this paper we have described a simple way to capture complex cellular dynamics that can occur due to the interaction between 2 myocytes that are coupled only via fibroblasts.In the heart the absence of gap-junctional coupling between myocytes could be due to the presence of ablation lines, scars or fibrosis resulting in spatial and electrical separation of myocytes.However under such conditions there is some evidence that conduction can still occur via fibroblasts across novel pathways coupling otherwise disconnected myocytes [26,27,25,40].
As we have described in this paper, such connections can give rise to a range of dynamical behaviour including reduced APD, synchronization of repolarization in uncoupled myocytes, initiation of action potentials in resting cells and conduction delays.At the level of cardiac tissue such non-local coupling might possibly impact wave propagation and lead to reentry or conduction blocks.
We have considered 3 different topological arrangements of 2 myocytes cou-pled to either 1 or 2 fibroblast units and investigated their effect on myocyte dynamics.We have chosen myocyte parameters corresponding to Shallow and Steep restitution slopes [35] and for the fibroblast we have used the MacCannell "active" fibroblast model [36] with modifications made to obtain different resting membrane potentials [19].Since the setup we have considered here is more likely in a heart tissue undergoing repair from either surgery, injury or disease we have modified the MacCannell model to simulate myofibroblasts by setting Irrespective of the connection topology, the primary effect of the coupling was a decrease in the APD of the myocyte with increase in coupling strength (Fig. 3) with the myofibroblast acting as a leaky capacitor [12].However the magnitude of the change in APD depended on the type of motifs and the myocyte parameters (i.e., whether the restitution is Steep or Shallow).The decrease in APD is more pronounced for the case with Steep parameters (Fig. 3) compared to Shallow parameters (Supplementary Fig. 10) across all the motifs.Also the magnitude of the change in APD increased with the complexity of the connections with M otif − 3 (the most complex and tightly coupled of the motifs considered) showing the greatest decrease in APD at maximum coupling.Furthermore we observed that the decrease in APD also increases for a more negative resting membrane potential of the fibroblast (Supplementary figures Fig. 23 and Fig. 21), an observation consistent with earlier studies [12].
Since the two myocytes considered in the motifs are uncoupled and spatially separated, there is a time difference in the excitation of the two cells.This could be due to the difference in times that the cells are excited in the myocardium (i.e., the effect of conduction velocity across the heart tissue) and/or due to the effect of conduction delay in the propagation of current via fibroblasts.We first studied the case where the predominant delay is the time difference (τ D ) between the stimulation of both myocytes (Fig. 5 and Supplementary Fig. 13).The delay in the stimulation of myocyte − 2 breaks the symmetry of coupling and results in differential myocyte APD even for M otif − 1 and M otif − 3 (Supplementary Fig. 14) for the action potential profiles for cells stimulated at different times).We observed that the sign and magnitude of the difference in APD of the two myocytes (∆AP D) is sensitive to one or more of the features varied here viz., connection topology, coupling strengths, delay (τ D ) and the myocyte parameters (Shallow or Steep).For conditions of finite stimulus delay, τ D acts as a proxy to the effect of conduction velocity across the myocardium.
Another type of delay that can occur in these systems is the time taken by the excitation to propagate from one myocyte to another purely via a fibroblast.In order to investigate the effect of this conduction we considered the case of τ D → ∞.In practical terms this means that one of the cells in the motif is never stimulated externally for the duration of the simulation, while the other myocyte is stimulated periodically as before.Figure . 6 (a-b) shows the 2-parameter phase space characterizing the different dynamical regimes as a function of the gap-junctional coupling conductance.Figure . 6 (c-f) highlights the sensitive dependence of the myocyte dynamics on the coupling strength, with small changes in coupling parameters resulting in widely different myocyte responses.The idea of infinite delay in stimulation is especially useful to illustrate the dynamics of initiation of action potential in regions that are isolated from its neighbours except for conduction via fibroblasts.
The intermittent sequence of action potentials observed in the IR regime in the Response Cell (and in Pacing Cell at very rapid pacing) suggests a plausible dynamical mechanism that can result in conduction block and initiation of reentry in tissue.Fig. 7, Fig. 8 and Fig. 9 summarize the role of the individual features in determining the dynamical regimes.While steep restitution and rapid pacing favour complex or irregular dynamics in both the Response and Pacing Cells, the motif structure and the location of the stimulated myocyte (with respect to fibroblasts connected to it) also play a critical role in determining the dynamical regime.In particular M otif − 1 (where one fibroblast couples to two uncoupled myocytes) with its very low number of 1 : 1 response is a very plausible connection topology in tissue that can give rise to conduction blocks and reentry.
While heterocellular coupling between myofibroblasts and myocytes have been reported to initiate ectopic activity in vitro [7], to the best of our knowledge this is the first in silico study to systematically investigate the different features that can potentially influence the myocyte dynamics when coupled purely via fibroblasts.We hypothesised that non-local coupling could potentially result in the initiation of action potentials in the quiescent myocyte and identified the parameters that described the various dynamical regimes possible for different connection topology.
Motifs are a simple prototype to investigate the effect of long-distance connections between uncoupled myocytes connected only via fibroblasts.Many studies have looked at the effect of fibroblast distributions in simulated tissue on wave dynamics [22,23,28].The idea of fibrotic and functional clusters [41,42] have been developed based on percolation theory to investigate the interaction of wave-propagation with fibrosis.Interstitial fibrosis is associated with non-ischemic cardiomyopathy and has been modelled as infinitesimal splits in a finite element mesh [43].Machine learning algorithms have been employed to understand the effect of local fibrosis patterns especially at border zones [44].More recently homogenisation techniques have been applied to model fibrosis as spatially repeating structures [45], or by using graph theoretical [46] and volume averaging approaches [47] to incorporate microscopic structures into a macro-scale problem.
Our methodology in this study differs significantly from the methods used in the above papers.We have adopted a bottom-up approach where we have developed simple structural motifs that can in principle be scaled to build scar boundary zones.While motifs have been used extensively in other areas [31,32,33], to the best of our knowledge this is the first paper that uses the idea of motifs to describe the various dynamics arising from myocyte-fibroblast interaction.Our approach is especially suitable to study the effect of fibroblast connections that form across ablation lines.Fibroblasts that are involved in the repair of the surgically ablated zones could enable conduction across the separated regions [25].The motifs described here provide a possible structural mechanism to couple disconnected regions in tissue.
We conclude by stating the limitations of our study and the scope of future work.While the motifs developed here are a useful prototype to simulate non-local coupling in cells across ablation line or fibrotic regions, they do not account for the effect of electrotonic diffusion, which acts to smooth wavefronts and reduce the effect of local variations in topology.So in as much as motifs can be a useful tool to quickly explore dynamics while varying several factors their maximal utility would be realised when they are used to build scar tissue boundaries and simulate conduction across non-conducting ablation lines in heart tissue.While we have considered homogeneous myocytes in our motifs, myocyte properties are heterogeneous in diseased hearts.Expanding the motif prototype to 2D tissue and incorporating heterogeneous myocytes in the motifs are two areas of future work.In our simulations we fixed the number of fibroblasts in a unit (N f ) connected to a myocyte to be 4.But the number of fibroblasts coupled to a myocyte is an important factor that affects the action potential and can be systematically varied to investigate its effect on wave propagation [48].Lastly an important limitation of this paper is that we have only simulated the electrophysiology of myocytes and not the mechanics of their contraction and relaxation.Mechanical contractions and the resultant change in tissue geometry have significant effect on wave-propagation.Another factor that influences wave-propagation is the mechano-electric feedback.These are aspects that will be incorporated in future studies.

Figure 1 :
Figure 1: Time series of transmembrane potentials for myocytes and fibroblasts.(a) Transmembrane voltage for uncoupled myocyte (a) and the corresponding S1S2 restitution curves generated at T = 600 ms pacing for both Shallow (broken) and Steep (solid) parameters (c).The time series of the transmembrane potential of myocyte (b) and fibroblast (resting potential V F R = −49 mV) (d) for the case of weak coupling (G gap = 0.5 nS) for both Shallow (broken) and Steep (solid) parameters.

Figure 3 :
Figure 3: Effect of coupling on APD.Two-parameter plots for all the 3 motifs describing the effect of the coupling strengths G Loc , G Long on the APD of both myocyte − 1 (a, c, e) and myocyte − 2 (b,d,f) with Steep parameters.

Figure 4 :
Figure 4: Action potential time series for coupled myocytes.Action potential profiles for both myocytes in M otif − 2 while pacing the cells at T = 600 ms at G Loc = G Long = 0.5 nS (a), 2.0 nS (b) and 4.0 nS (c) respectively.The solid and broken blue lines correspond to myocyte − 1 and myocyte − 2 for the Shallow parameter set while solid and broken red traces correspond to the same for the Steep parameter set.

Figure 5 :
Figure 5: Effect of time delay in the stimulation of myocyte − 2. The difference in APD (∆AP D) between the two myocytes for M otif − 1 (a, b, c), M otif − 2 (d,e,f) and M otif − 3 (g, h, i) for Steep parameters.The values of the delay are τ D = 0 ms (a,d,g), = 10 ms (b, e, h) and = 25 ms (c, f, i).
)) ∆AP D was positive at large (G Loc , G Long ) values, while a large G Long and small G Loc conductance results in a greater reduction of APD for myocyte − 1 than for myocyte − 2 (as seen by the negative ∆AP D values).

Figure 6 :
Figure 6: Dynamical regimes described by two-parameter conductance maps.Dynamical regimes characterised as the fraction of stimulus that elicit an action potential in the Response Cell for both Shallow (a) and Steep (b) cell parameters for M otif − 3. (c-f) The action potential profiles for the parameters marked (c,d,e,f) in the case of Shallow parameter set (a) comparing the different dynamical behaviour in the Pacing and Response cells.

Figure 7 :Figure 8 :Figure 9 :
Figure 7: Effect of restitution.Fraction of occurrence of each dynamical regime, viz., N R, IR and 1 : 1 summed over all pacing cyles and motifs at both Pacing Cell (P Cell) and Response Cell (RCell) for Shallow (a) and Steep (b) parameter sets.

Figure 11 :
Figure 11: Action potential profiles for both myocytes in M otif −1 while pacing the cells at T = 600 ms at G Loc = G Long = 0.5 nS (a), 2.0 nS (b) and 4.0 nS (c) respectively.The blue and red traces correspond to Shallow and Steep parameter sets.

Figure 12 :
Figure 12: Action potential profiles for both myocytes in M otif −3 while pacing the cells at T = 600 ms at G Loc = G Long = 0.5 nS (a), 2.0 nS (b) and 4.0 nS (c) respectively.The blue and red traces correpsond to Shallow and Steep parameter sets.

Figure 15 : 2 -
Figure 15: 2-parameter conductance map for pacing period T = 600 ms describing the different dynamical regimes for the Steep parameter set, viz., N R, IR and 1 : 1, characterised as the fraction of stimuli that elicit an action potential in the Response Cell for M otif − 1 (a,b) and M otif − 2 (c-f).

Figure 21 :
Figure 21: Effect of coupling motifs on APD for Steep parameters with V F R = −49 mV.

Figure 22 :
Figure 22: Identifying regimes for Steep parameters at T = 600 ms that initiate action potentials in non-stimulated cells for V F R = −49 mV.

Figure 23 :
Figure 23: Effect of coupling motifs with Shallow cell parameters on APD for V F R = −49 mV.

Figure 24 :
Figure 24: Identifying conductance parameters that initiate action potentials in non-stimulated myocyte with Shallow restitution for V F R = −49 mV.