Multicellular automaticity of cardiac cell monolayers: effects of density and spatial distribution of pacemaker cells

Self-organization of pacemaker (PM) activity of interconnected elements is important to the general theory of reaction–diffusion systems as well as for applications such as PM activity in cardiac tissue to initiate beating of the heart. Monolayer cultures of neonatal rat ventricular myocytes (NRVMs) are often used as experimental models in studies on cardiac electrophysiology. These monolayers exhibit automaticity (spontaneous activation) of their electrical activity. At low plated density, cells usually show a heterogeneous population consisting of PM and quiescent excitable cells (QECs). It is therefore highly probable that monolayers of NRVMs consist of a heterogeneous network of the two cell types. However, the effects of density and spatial distribution of the PM cells on spontaneous activity of monolayers remain unknown. Thus, a simple stochastic pattern formation algorithm was implemented to distribute PM and QECs in a binary-like 2D network. A FitzHugh–Nagumo excitable medium was used to simulate electrical spontaneous and propagating activity. Simulations showed a clear nonlinear dependency of spontaneous activity (occurrence and amplitude of spontaneous period) on the spatial patterns of PM cells. In most simulations, the first initiation sites were found to be located near the substrate boundaries. Comparison with experimental data obtained from cardiomyocyte monolayers shows important similarities in the position of initiation site activity. However, limitations in the model that do not reflect the complex beat-to-beat variation found in experiments indicate the need for a more realistic cardiomyocyte representation.


Introduction
The study of collective behaviour of interconnected elements is important to the general theory of reaction-diffusion systems and is applicable to different areas of science from physics to physiology [1][2][3][4]. Diverse behaviours can be observed in networks of connected elements including a variety of synchronous regimes, pattern formation, spiral wave propagation, and spatio-temporal chaos. An example of a physiologically-relevant system is the monolayer cultures of neonatal rat ventricular myocytes (NRVMs) that is often used as experimental models in studies on multicellular cardiac electrophysiology [5][6][7][8]. These monolayers usually exhibit electrical automaticity [9][10][11], defined as the ability to generate action potentials without external (electrical, mechanical, or chemical) stimulations. Primary cultures of isolated NRVMs revealed the existence of pacemaker (PM) cells and quiescent excitable cells (QECs) [12]. Thus, it is highly probable that monolayers of NRVMs are heterogeneous networks of these two cell types. The spontaneous beating rate in these cultures is reportedly modulated by plating density [13]; however, the specific effects of PM cell density and spatial distribution on confluent monolayer spontaneous frequency are still unknown.
Understanding the impact of density and spatial distribution of PM cells on automaticity of the global cardiac network is highly relevant for several clinical situations, such as the reliability of heart activation by the sinoatrial node (SAN) and the initiation of arrhythmias by cells susceptible to delayed afterdepolarisations (DADs). When the heart is being activated, the SAN PM cells act as initial depolarizing current sources and the surrounding myocardial cells as related sinks. Considering the electrotonic depression of the SAN by electrical connection with the myocardium, the number of PM cells must be large enough to initiate this activation. As a matter of fact, the spontaneous rate and minimum diastolic potential of the SAN dramatically increase when it is isolated from atrial tissue [14]. This knowledge is exploited in the process of creating biological PMs, an alternative to electronic PMs in the treatment of bradycardia [15]. Indeed, to create a sustained PM function in canine ventricle, xenografted cell-containing fluid injected in the left ventricular myocardium must contain at least 700 000 adult human mesenchymal stem cells, and consist of ∼40% PM cells [16]. This is clear example of where the number of PM cells is important, and must be locally increased to activate the surrounding tissue. On the other hand, the presence of PM cells is also linked to arrhythmia initiation. Premature ventricular complexes, for example, can occur when a group of cells with DADs manage to initiate ectopic electrical propagation. These pathological cells act as secondary PM cells competing with the leading SAN. Since a ventricular cell is on average connected to 11 other cardiomyocytes [17], there must be a sufficient number of connected cells simultaneously generating DADs to overcome the electrotonic source-sink mismatch and initiate propagation [18].
Monolayer cultures of NRVMs may provide a useful biological tool to investigate the role of PM-cell-cluster patterns on automaticity and therefore gain better understanding of clinical issues, including creation of biopacemakers [19]. However, the seeding process is usually random and the final spatial distribution of the two cell types unknown. Experimental assessment of density and spatial distribution of PM cells within a monolayer remains an unsolved problem.
Mathematical modelling may help gain insight into the complex source-sink mechanism behind the effects of PM cell clusters on the spontaneous behaviour of the 2D network. Rather than concentrating on PM cell aggregate spatial characteristics, simulation studies on monolayers of excitable cells have predominantly focussed on ion channel properties [20,21], effects of coexistence of non-excitable cells [11,22,23], intercellular electrical connectivity [24][25][26], or sink (cells clamped at steady state potential) and break (regions with no conductivity) densities [27]. Theoretical studies on the spatial-temporal behaviour of random mixtures of excitable and oscillating cells have been performed. One notable study was performed by Kanako et al [28], who demonstrated the importance of diffusion constant on global activation pattern on a lattice. However, the specific issue of the impact of spatial distribution structure of PM cells on monolayer activity remains to be evaluated. To fill this gap, the present study aims to theoretically assess the effects of spatial distribution and density of PM cells on global automaticity of the cardiac monolayer.
Here, we hypothesized that: (a) The monolayer is a binary network of interconnected PMs and QECs (passive cells and cleft spaces were not included). (b) The diffusion is isotropic, since the gap junctions of NRVMs are evenly distributed [29]. (c) The spatial distribution and density of PM cells and QECs (independent variables) influence automaticity via modulation of spontaneous frequency and electrical propagation properties (dependent variables).
Based on these hypotheses, we created different 2D networks of PM cells and QECs. PM cells were spatially distributed using a stochastic model allowing the variation of aggregation (PM cells connecting to other PM cells) and nucleation (PM cells forming new clusters). The modulation of spontaneous frequency and electrical activation by PM cell spatial distribution and density was then studied with a FitzHugh-Nagumo (FHN)-type model of excitable media [30]. Simulations showed a clear nonlinear dependency of spontaneous activity (occurrence and amplitude of spontaneous period) on the spatial patterns of PM cells. In most simulations, the first initiation sites were found to be located near the substrate boundaries. Data were compared to fluorescence recordings of calcium waves in NRVM monolayers to highlight similarities and limitations of the model.

Generation of spatial patterns for PM cells
Patterns were created assuming the existence of only two types of NRVMs (PM cells and QECs) and this led to a binary 2D network. Each PM cell made connections with all its neighbours (the eight nearest cells in a square lattice). The spatial pattern formation algorithm randomly determined the positions of PM cells, and QECs filled the remaining available space.
The algorithm is presented in figure 1(a) and is summarized as follows. The starting point was an empty square lattice of dimension N × N. For the first iteration, an initial PM cell was randomly seeded in one of the N 2 available sites. The status of a place containing an autonomous cell was referred to as the autonomous site ('aut' site). During the second iteration, a random number ⩽ ⩽ p p (0 1) was drawn from a uniform distribution and compared with a user-defined iteration-independent parameter P thr ⩽ ⩽ ( ) p 0 1 thr that decided whether the new PM cell would or would not be in contact an existing cluster (a group of one or more interconnected PM cells).
If ⩽ p p , thr then the new cell would not contact any existing cluster and would be randomly deposited in a free site from the M 2 set (blue available sites in figure 1(b)). Inversely, if > p p thr , then the new cell would be part of the M 1 set (red available sites in figure 1(b) corresponding to the 8-connected free neighbourhood). Repeating this procedure for n iterations increased the density of autonomous cells (D aut = number of PM cells/N 2 ) to the desired level. There was a D aut limit (D aut > D aut,max (p thr )) beyond which no new sites could be added to the lattice without converging on an existing cluster.
In summary, spatial distribution was modulated by parameter p thr while the density was modulated by the number n of iterations.

Mathematical model of autonomous cells and QECs
Modelling could be performed with various types of excitable representations ranging from the simplest FHN model to cardiac realistic ionic models [31,32]. Here, we selected a simple modified FHN model of PM cells integrated into a 2D excitable network similar to cultured monolayers. The model was described previously by Borek et al [30]: where V and w represent the activation (similar to the transmembrane voltage) and inhibition variables, respectively. The parameter values were fixed to obtain an excitable but nonautonomous medium: diffusion coefficient D x = D y = 0.2, κ = 0.6, β = 0.7, w h = w l = 0.1, I P = 0. A PM current I P = 1 was applied to achieve autonomous activity. Spatial distribution of the autonomous cell was thus implemented by setting I P = 1 in 'aut' sites (all the locations occupied by an autonomous cell) of the lattice built with the stochastic model described in section 2.1.
The model was solved on the N 2 lattice (N = 200, unless stated otherwise in the text) using an Euler forward finite difference in time (Δt = 0.001) coupled with a central difference in space of Δx = Δy = 0.05 and no-flux boundary conditions applied to the four sides of the lattice. The estimated space constant of the system with the above parameters was greater than 10Δx. Initial conditions were the resting state for QECs and the state of minimum V for PM cells. All PM cells were set to the same phase in this study for simplicity.

Experiments with NRVMs monolayer cultures
All animal handling procedures were concordant with the Canadian Council on Animal Care guidelines and were approved by the Montreal Heart Institute Animal Research Ethics Committee. Sprague-Dawley rats aged 1-3 d were sacrificed by decapitation. Beating hearts were removed immediately and kept in cold Ca 2+ -and Mg 2+ -free Hank's balanced salt solution. Ventricular muscle was excised and tissue was minced on ice into 1-3 mm 3 pieces. The mixture was subjected to purified enzymatic digestion using a Neonatal Cardiomyocytes Isolation System (Worthington Biochemical, Lakewood, NJ). Isolated cells (enriched cardiomyocytes) were counted and seeded at a density of 1 × 10 6  washed out twice with DMEM and replaced with 1.5 mL fresh DMEM (without FBS and P/S). A 15 min equilibration period was allowed for the de-esterification process of Fluo-4 AM to occur before initiating acquisitions.
Mapping experiments were performed with a setup developed in-house using a high-speed CCD camera system (80 × 80 pixels, RedShirtImaging, LLC, Decatur, GA). The dye was excited with a quartz tungsten halogen lamp (Oriel Instruments, Stratford, CT). The filters used for excitation and emission were λ excitation ≈ 480 ± 20 nm (Chroma Technology, Bellows Falls, VT) and λ emission ≈ 535 ± 25 nm (Semrock, Rochester, NY), respectively. The system was set to image a field of view of 1.38 cm 2 and the acquisition frame rate was 125 Hz for all experiments. Signals were filtered and analysed using a program developed in-house using Matlab software (R2008, MathWorks, Natick, MA). In summary, raw acquisitions were normalized by the minimum fluorescence for each pixel and the dF/dt (the first time derivative of fluorescence) was approximated using a finite difference approach. The resulting signals were then spatially filtered using a Gaussian kernel (5 × 5 pixels with σ = 1.5) followed by a temporal moving average with a 2-sample window.

Characterization of the spatial distribution of autonomous cells based on the stochastic model
Patterns of PM cells obtained by simulating the stochastic model are shown in figure 2 for D aut = 0.35 corresponding to 35% occupation of the lattice by PM cells. When p thr = 0, only a single cluster (normalized area of 0.35) is created by aggregation of the new PM cells around the initial random seeding. All patterns obtained with p thr = 0 have a single Eden-like cluster and are thus similar to the example presented in panel (a). For this special case, the only difference between created patterns stems predominantly from the initial random seeding around which the aggregation process occurs. Patterns in figure 2 were obtained with > p 0 thr and resulted in a higher number of clusters (B, 72; C, 1020; D, 5400; E, 14015; F, 14602); thus, decreasing the mean normalized area per cluster (B, 486.1 × 10 -5 ; C, 34.3 × 10 -5 ; D, 6.5 × 10 -5 ; E, 2.5 × 10 -5 ; F, 2.4 × 10 -5 ).
As D aut increases when building patterns, the number of positions in set M 2 , i.e. sites not touching an existing clusters, decreases such that there exists a limit for each pattern when M 2 is empty and < D 1 aut . Except for = p 0 thr for which the limit is D aut = 1, the limit is represented by a distribution of D aut,max (p thr ) corresponding to the distribution of maximum D occ for individual patterns. This limit is shown in figure 3(a) as three curves representing the 5th percentile line (bottom curve), median (middle curve), and 95th percentile (top curves) obtained for a set of 1000 patterns per (p thr ,D aut ) pair. The result is shown with scaling of p thr ( ) p thr 1/4 to facilitate viewing since most of variations occurred at low p thr . Median D aut,max (p thr ) decreases with increasing p thr . The interval between the 5th and 95th percentiles of D occ,max (p thr ) is not constant as a function of p thr 1/4 . It is narrow for p thr 1/4 close to zero, widest around The average number of clusters (N clusters ) obtained for dimension N = 100 and 1000 samples per (p thr , D aut ) pair is shown in figure 3(b). Two general results can be stated. Firstly, for a given > p 0 thr , the number of clusters increases and then decreases with D aut in a biphasic manner, as long as D aut < D aut,max (p thr ) (the median shown by a white line in the figure).
Secondly, for a given D aut , the number of clusters is higher, again with greater p thr , as long as D aut < D aut,max (p thr ). Thus, the maximum N clusters is found for p thr = 1 and D aut immediately below D aut,max (1). The average maximum cluster area (mean S cluster ) was calculated from the set of maximum cluster area for each simulated patterns and is shown in figure 3(c). p thr 1/4 = 1 is a special case where all clusters consist of a single site. Otherwise, increasing D aut resulted in higher mean S cluster independent of p thr 1/4 . However, a difference was observed in the slope for increasing D aut .
Indeed, for low p thr 1/4 the increase is more linear and starts at low D aut while the transition to high-area cluster occurs at high D aut (around 0.6) with a rapid augmentation. Interestingly, for a given D aut , increasing p thr 1/4 results in a decrease of mean S cluster . However, biphasic variation occurs for D aut between 0.6 and 0.8. This is also the region of greatest variability between samples as highlighted by the higher standard deviation in the maximum cluster size (std. S cluster ) depicted in figure 3(d).

Analysis of the stochastic model
The parameter p thr in the model represents the threshold in the interval [0, 1] selected between creating a new cluster (nucleation) and adding an autonomous site to an existing cluster (aggregation). Thus, the probability for a new cluster to be created should be equal to p thr as long as D aut < D aut,max (p thr ). The average creation rate estimated by the mean change in the number of clusters per iteration is depicted for three increasing D aut intervals in figure 4(a). As expected, the average creation rate correlates almost perfectly with p thr for small D aut (black . Thus, on average, the total number of clusters are less than the value predicted by p thr . The decrease in the average creation rate can be explained from the following model. Starting from the initial single autonomous site in an empty square lattice, the number of clusters (n) at time zero is equal to 1 (n 0 = 1). The change in average number of clusters at iteration i + 1 (n i+1 ), can then be approximated by  As expected, increasing D aut yields a higher area. However, the increase is more linear for low p thr 1/4 , while high p thr 1/4 shows minimal area until D aut ≈ 0.6 where a rapid transition occurs due to aggregation of individual clusters. Note that there is a slight shift to higher D aut for mean S cluster around 60 for p thr 1/4 between 0.2 and 0.4. (d) Standard deviation of the maximum cluster area (std. S cluster ) showing that intrinsic variation between generated patterns occurs over a different range of D aut as a function of p thr 1/4 . The maximum range with higher variation is found for p thr 1/4 between 0.2 and 0.4.
Here ε p D ( , ) thr aut represents the probability that adding the new site results in the merger or aggregation of two or more clusters. The average of ε(p thr , D aut ) obtained for 1000 simulations per pair of parameters are plotted in figure 4(b). At low D aut , ε(p thr , D aut )→0 adding a new autonomous site to an existing cluster should not connect to existing clusters because the number of clusters is small compared to the distance between clusters. As D aut increases, so does ε(p thr , D aut ) since the average distance between neighbouring clusters decreases. Finally, ε → 0 with increasing D aut as the number of clusters decreases to a single cluster because of aggregation. The variation in ε(p thr , D aut ) with p thr is shown in figure 4(b) with values increasing as p thr increases. For the special case of p thr = 1, ε(p thr , D aut ) is not defined for D aut < D aut,max (p thr ) because of the impossibility of adding an autonomous site to an already existing cluster.
It is still possible to estimate the average total number of clusters for this model. For example, for p thr = 0, ε = 0 since there is a single cluster and no new cluster is created on every iteration (see figure 2 with D aut = 0.35). In contrast, for p thr = 1, the added new site will never come into contact with an existing cluster as long as D aut < D aut,max (1). We can deduce from [0,1] thr exhibiting biphasic variation with an initial increase followed by a decrease for D aut > D aut , max (p thr ). The black line represents the median of D aut , max (p thr ). equation (2) that the average number of clusters should be: where the first term corresponds to the average number of clusters created for the density D aut in a medium with N 2 sites and the second term corresponds to the total loss (over the addition process from 0 to D aut ) of clusters caused by cluster aggregation or fusion.

Modulation of the autonomous period by PM cell density in simulations
The period of activity (normalized to the single-cell period) is increased in a heterogeneous spatially extended system compared to single-cell autonomous activity of the FHN for equation (1) with I P = 1. The change in period when simulating a circular patch with increasing radius is depicted in figure 5(a). Two main characteristics can be found: (1) there is a minimal radius for automaticity occurrence, and (2) the normalized period is increased by ∼32% when automaticity starts to occur and decreases asymptotically to the period of the single-cell model.
We then asked what effects stochastic patterns would have on the minimum density of autonomous sites (D aut,min ) for autonomous activity to occur. Simulation results are shown in figure 5(b) where a clear transition between low D aut,min and high D aut,min is found with ∼10 times larger density needed for high p thr 1/4 (D aut,min = 0.015 for p thr 1/4 close to zero and D aut,min = 0.155 forp thr 1/4 close to 1). The transition occurs around p thr 1/4 = 0.5 as estimated by fitting a sigmoidal function (black curve).
Low density of autonomous cells also always results in longer periods of activity. Simulation results for In comparison, for p thr 1/4 = 0.6 (middle panel), the decrease in period is slower. While some variability in period is detected between samples, the level seems to be less and restricted to lower D aut . A high p thr 1/4 of 0.9 results in a slow decrease in period T, but also negligible variability between samples, a hallmark of the low variability observed in the maximum cluster area (as seen in figure 3(d)). In all cases of figure 5(c), a greater minimum density is required to generate spontaneous activity. Global mean of the period (mean T) of activity is presented as a function of p thr 1/4 and D aut in figure 6(a) and the corresponding standard deviation of the period (std. T) in figure 6(b). A closer look at the results in figure 6(a) shows that a given delay occurs at a small D aut for small p thr 1/4 while more PM cells are needed to obtain a similar rate for high p thr 1/4 . For example for T = 1.15 (black curve), D aut is around 0.03 for p thr 1/4 = 0.1 compared to 0.39 for p thr 1/4 = 0.9; thus, giving rise to a 36% difference in cell density. As T approaches unity, the difference increases such that for T = 1.10 (white curve), the difference in cell density is 49%. In most of the parameters space, the transition from low to high D aut occurs within the interval < < p 0. where patterns consist of a combination of small and large clusters (see figure 2(c)). As expected, the variability between results is generally higher at low D aut followed by a monotone decrease with increasing density. However, the change in variability as a function of p thr 1/4 for a given D aut has biphasic variation within the same interval . This indicates an instability in the generated patterns, and a correlation between variability and increased sensitivity of the period to changes in p thr 1/4 . 2). It clearly shows that a minimum cluster area is needed to drive the monolayer (highlighted by the grey area labelled no automaticity). While a small cluster yielded ∼32% slowing, increasing the cluster area decreased the period towards the isolated cell period (T = 1). (b) Minimum D aut (D aut,min ) as a function of p thr 1/4 is presented. A clear transition between low D aut,min close to 0.015 is found for p thr 1/4 < 0.4 and rapid transition to values around 0.15 for p thr 1/4 > 0.6. A sigmoidal fit with p thr 1/4 at half value of 0.5 is superimposed (black line). (c) Normalized period of activity (T) for three p thr 1/4 values is shown where a faster decrease in T is found for p thr 1/4 = 0.05 (top) compared with p thr 1/4 = 0.6 (middle) and p thr 1/4 = 0.9 (bottom). In addition, variability between samples occurs over a larger D aut interval for smaller p thr 1/4 .

Impact of the spatial distribution of autonomous cells on synchronous activation
The simulations were analysed and the total activation time (t act,tot ) for the last simulated beat calculated. For a specific beat, the total activation time is the activation delay between the first and last nodes and represents the time required for the initial activation to fully propagate in the monolayer. It is calculated as time of the last activation-time of the first activation for the same beat. The mean t act,tot (figure 7(a)) and the standard deviation of t act,tot (std. t act,tot , figure 7(b)) show a variation as a function of p thr 1/4 and D aut . Small p thr 1/4 yields the largest delay (example for D aut = 0.25 with t act,tot = 0.37 is shown in the top panel of figure 7(c)). The delay decreases with increasing p thr 1/4 (examples with t act,tot = 0.26 for p thr 1/4 = 0.5 and 0.31 for p thr 1/4 = 0.9). Thus, while large clusters of autonomous cells yield the smaller period of activity for a given density, the activity is less synchronous and more dependent on propagation. On the other hand, a more uniform distribution of autonomous cells (as created with large p thr 1/4 ) is slower, being more affected by electrotonic effect; however, it yields a more synchronized activation. Importantly, as depicted in figure 7(d), the normalized period of activity T shows no interbeat variation regardless of pattern type (e.g. p thr 1/4 values).

Localization of the first activation sites
Another interesting result can be seen in figure 7(c), which shows that the sites where activation starts (site of first activation, dark blue from the colour scale) are located close to the boundaries. Indeed, the example for p thr 1/4 = 0.2 (top panel of figure 7(c)) and p thr 1/4 = 0.5 (middle panel) show a single first activation site respectively located at the bottom and top of the substrate. Although only a single first activation site is observed for the example with p thr 1/4 = 0.9 (bottom panel), a second region starts to activate rapidly at the bottom left corner.
Pooling and plotting the first activation sites of the last three activations of each simulation (all < D 0.35 aut and p thr 1/4 ) yields figure 8(a). Clearly, the vast majority of the first activation sites is found close to the boundaries. Indeed only 6.3% of sites are located within the red square. Simulations with the spatial patterns of autonomous cells were repeated after adding a border of excitable cells (with a layer of 0.65 cm). The distribution of the first activation sites is displayed in figure 8(b), which shows a more homogeneous distribution of sites. By adding the excitable layer around the pattern of autonomous cells (width of 13 sites all around the previously simulated distribution patterns), the proportion of sites within the red square increased to 9.2% (corresponding to an increase of 46% compared with control conditions in panel (a)). A comparison of the density of first initiation sites in the red square was made between these two sets of simulations at cumulative D aut . The result is shown in figure 8(c) and indicates that the presence of no-flux boundary conditions near autonomous cells favour their activation by constraining local diffusion; thus, decreasing the electrotonic effect on spontaneous activity.

Complex dynamic behaviour: insights from experimental data
A typical post-culture recording of calcium waves is shown in figure 9. In this example, three different initiation sites were identified and located at the bottom right (panel (a)), upper right (panel (b)), and bottom left (panel (c)) corners of calculated activation maps, similar to what is shown in our simulations. However, contrary to the simulations, the pacemaking loci compete with each other. During transitions between first activation sites, some beats show multiloci quasi-simultaneous activation (panel (b)). In panel (d), the initiation site of each beat is specified with a coloured star on a single-pixel signal. Large beat-to-beat variations are observed in the dynamics, such as the interbeat interval and the position of the initiation site. The complex dynamical behaviour in autonomous activity is clear when looking at the interbeat sequence (panel (e)), which displays a large variation in intervals, ranging from 664 to 1392 ms.

Discussion and conclusion
In the present study we proposed, to the best of our knowledge, the first theoretical study assessing the effects of spatial distribution and density of PM cells on the global automaticity of the strongly-coupled cardiac monolayer.
The spatial distribution of PM cells on a × N N lattice was generated by a very simple stochastic model (depicted in figure 1), with a single parameter p thr determining the rate of nucleation (formation of new clusters) and aggregation (growth of existing clusters). Using this model, we were able to obtain various spatial distributions of PM cells. The remaining available sites of the lattice where filled with excitable but resting cells, to form a 100% confluent monolayer. The spatial distribution homogeneity is modulated by p thr and the density by the number n of iterations. It is important to note that the algorithm is based on an 8-neighbour rule. Using a 4-neighbour rule would affect the spatial characteristics of PM cell distribution. Furthermore, this rule may impact pacemaking activity with increasing D aut , which is a subject of a future study. The pattern formation model shows similarity to other models for cooperative sequential adsorption [33], invader-resident competition [34], and reaction-diffusion controlled cluster-cluster aggregation [35].
With this approach to build spatially distributed patterns, the number and size of PM cell aggregates can be varied while still mimicking the random nature of the spatial distribution and density of PM cells in experimental cardiomyocyte monolayers. Created patterns demonstrate complex and nonlinear changes in the number and spatial dispersion of PM cell clusters as a function of the parameter p thr (figure 2). The main interest of this pattern formation method is that it easily permits the effects of spatial patterns of PM cells to be studied within an excitable but resting lattice. Currently, this question cannot be easily addressed with available experimental techniques.
A previous study showed that a random distribution of PM cells with varying frequencies coupled to excitable cells exhibit from incoherent behaviour to global synchronization with increasing coupling between cells [28]. Here, we were interested at looking at the effects of having different spatial correlations in PM clusters in a condition of highly-coupled cells where global synchronisation is favoured. Thus, we found in our simulations (and similarly in experiments) all or none spontaneous activity of the monolayer but no asynchronous or reentrant activity as seen in monolayers with lower intercellular coupling [23,25,28,36,37]. Patterns obtained with our stochastic model confirmed the existence of a pattern-dependent minimum PM cluster size to generate automaticity ( figure 5(a)). This was observed even for the simplest FHN excitable tissue where the PM cluster must contain a sufficiently large number of PM cells to initiate electrical activation of the monolayer. Indeed, around 280 000 PM cells are necessary for a biological PM to drive a canine ventricle [16]. Early-after depolarisations (EADs), which underlie abnormal ectopic activity, have been found to be electrotonically suppressed in well-connected networks of NRVMs [38]. Propagation cannot be initiated when EAD-susceptible cells are connected with too many non-EAD-susceptible myocytes, because the transmembrane voltage gradient generating the diffusion current becomes too low to rapidly depolarize above threshold in non-EAD cardiomyocytes.
For our model, a smaller density is required to generate automaticity (figure 5(b)) with large aggregates of PM cells (low p thr ). This result emphasizes the importance of the granularity (local density) of PM spatial distribution on spontaneous activity, a crucial aspect to take into account in regenerative cardiology for in vivo cell injection procedures as that large aggregates of injected cells prone to spontaneous activity could induce tissue activity. A similar condition determines the capacity of DAD-susceptible foci to initiate an activation wave that can propagate throughout the tissue [18] and highlights the impact of non-excitable cells on silencing of spontaneous activity [23]. The sigmoid-like variation in the minimum density for spontaneous activity to occur (D aut,min ) and the spatial distribution correlation through the parameter p thr remains a key result emphasizing the nonlinear relation between cluster spatial characteristics and intercellular interaction to overcome the electrotonic effect. Unsurprisingly, the period of spontaneous activity always decreases for increasing PM cell density ( figure 5(c)). However, a clear p thr -dependent decay rate in period is accompanied by a change in variability between samples. This result requires further experimental validation. For a given plating density (same D aut ), the variation with patterns obtained for different p thr can have a strong effect on autonomous activity (notably the period) of the monolayer ( figure 6). This result may explain, at least in part, the variability in spontaneous behaviour observed experimentally among NRVM monolayers with identical plating densities and similar culture conditions. This knowledge could be important in the process of creating biological PMs, an alternative to electronic PMs in the treatment of bradycardia [15].
Interestingly, simulations revealed that the initiation sites of spontaneous activity are mostly located at the edge of the monolayer (figure 7(c)) compared with ∼20% in experiments from a previous study [21]. This finding may be explained by a source-sink mechanism [18,39]. At the edge, which has no-flux boundary conditions, PM aggregates are effectively connected with fewer QECs; as such, they are less electrotonically depressed and possess increased capacity to initiate activation. It is also particularly interesting to note that seemingly competing initiation sites appear when PM cells are distributed in a large number of small clusters (large values of p thr ). This result shows some similarities with experimental results (figure 9) that demonstrate the existence of localized PM foci and a tendency for the first initiation sites to be located at the boundary of the culture dish where sink effects would be decreased. As mechanical effects have shown to induce spontaneous activity at a certain distance of the boundaries [40], how spatially-distributed PM cells in an electrical-mechanical model would modify this property.
However, a major difference between experimental and modelling results exists in the temporal behaviour where, for the same acquisition, experimental data show complex competition between PM loci. No simulations with the FHN model, a simple excitable model, exhibited the complex sequence of activation depicted in figure 9. In simulations, there is no beat-to-beat variability in autonomous frequency and initiation site position for a given PM cell pattern as depicted in figure 7(d). However, experimentally, three initiation sites competed. Furthermore, the fastest locus seemed to experience fatigue, occasionally leaving the competition open to the two slower loci. Similar complex behaviour has also been documented by others [21]. The dynamical behaviour highlights limitations of the FHN model and the need to develop and study the behaviour in monolayers with a more realistic NRVM ionic model. However, the level of model and its characteristics needed to reproduce this results is unknown since comparison between the FHN and Luo-Rudy models (a more realistic cardiac model) showed qualitatively similar results [23]. In the present study, the focus has been placed on the effect of the spatial structure of PM cells on electrical activity. However, this activity also relies on the excitability of the membrane, which depends on the properties of ion channels; these are not well represented by FHN-like models and have important effects on propagation [41]. A quiescent version of the NRVM ionic model has been published [42] and could be modified to exhibit spontaneous activity more similar to real cardiomyocytes in culture with complex spontaneous activity [43].
In the present study, monolayers were modelled as a 100% confluent and isotropic network composed solely of two types of cardiomyocytes (quiescent and PM cells) identically connected with strong coupling between cells. It is clear that experimental frequencies of activity will vary between PM cells [12], a factor that was not taken into account in this study, but this has been carefully studied previously for random distributions [28]. Such variation in intrinsic frequencies may have a physiological role even for cardiomyocyte monolayers. For example, to perform normal PM function, the SAN is known to require a heterogeneous architecture that involves multiple cell types with a gradient of cellular and intercellular characteristics [44]. Also, more realistic spatial characteristics of monolayers of cardiomyocytes have been shown to affect global electrophysiological behaviour [45]. The spatial characteristics not included in our model include cell shape, cell-to-cell coupling, and degree of cleft space. The simplifications (single diffusion coefficient, identical number of neighbours, single-frequency PM cells with identical initial conditions) were performed to isolate the specific effects of spatial distribution of PM cells. It remains to be seen how inclusion of these more realistic characteristics of monolayers would influence monolayer spontaneous activity as a function of spatial distribution of PM cells. Such future studies will undoubtedly improve our understanding of the complex temporal behaviour found in experiments and will help in the development of biopacemakers as a cardiac therapy.