In silico study of multicellular automaticity of heterogeneous cardiac cell monolayers: Effects of automaticity strength and structural linear anisotropy

The biological pacemaker approach is an alternative to cardiac electronic pacemakers. Its main objective is to create pacemaking activity from added or modified distribution of spontaneous cells in the myocardium. This paper aims to assess how automaticity strength of pacemaker cells (i.e. their ability to maintain robust spontaneous activity with fast rate and to drive neighboring quiescent cells) and structural linear anisotropy, combined with density and spatial distribution of pacemaker cells, may affect the macroscopic behavior of the biological pacemaker. A stochastic algorithm was used to randomly distribute pacemaker cells, with various densities and spatial distributions, in a semi-continuous mathematical model. Simulations of the model showed that stronger automaticity allows onset of spontaneous activity for lower densities and more homogeneous spatial distributions, displayed more central foci, less variability in cycle lengths and synchronization of electrical activation for similar spatial patterns, but more variability in those same variables for dissimilar spatial patterns. Compared to their isotropic counterparts, in silico anisotropic monolayers had less central foci and displayed more variability in cycle lengths and synchronization of electrical activation for both similar and dissimilar spatial patterns. The present study established a link between microscopic structure and macroscopic behavior of the biological pacemaker, and may provide crucial information for optimized biological pacemaker therapies.


Introduction
Oscillating, autonomous or spontaneous electrical activity is the basis of normal heart physiology [1], as well as some impaired rhythms triggered by ectopic activity [2]. Two oscillating mechanisms or clocks, the membrane and calcium clocks, are hypothesized to control the sinoatrial node (SAN) isolated cellular rate [3][4][5]. Membrane clock refers to the synergy of transmembrane ionic currents [6,7], and calcium clock to the oscillations of intracellular calcium concentration [8]. Developmental variations may change magnitudes of the respective clock components [9]. Interplay between these two strongly coupled mechanisms may be responsible for spontaneous activity and temporal fluctuation in heart rate [10]. At the cellular level, the clocks basically create an ionic imbalance during the diastolic period, leading to a net inward flux of ionic current that slowly increases membrane potential until the threshold (− 40 mV) to fire an action potential is reached. Inducing this net inward flux of ionic current during the diastole can actually generate automaticity in otherwise quiescent cardiomyocytes (CMs). This principle has been exploited in the design of biological pacemakers (BPs), a therapeutic alternative to overcome the shortcomings of cardiac electronic pacemakers [11] in the treatment of bradycardia. Different procedures have been proposed, including injection-based gene [12] and cell therapy [13], that locally modify cardiomyocyte phenotype or bring differentiated cells in the myocardium.
These concepts are limited by the lack of control on the spatial distribution and phenotype of pacemaker (PM) cells within the resting but excitable cellular network of the myocardium. We have shown that density and spatial distribution of PM cells can alter significantly the emergence and characteristics of multicellular spontaneous activity [14]. In fact, density and spatial distribution of PM cells, a priori unknown in BPs, may lead to a non-negligible intrinsic variability in the spontaneous activity of the overall network. Intrinsic variability is defined as behavioral discrepancies among BP samples that had undergone the exact same protocol. This phenomenon could eventually compromise the success of BP implantation in patients, and is observed even in in vitro BP models like monolayer cultures of neonatal rat ventricular myocytes (NRVMs), which are also heterogeneous network of autonomous and quiescent cardiomyocytes [15].
In the present simulation study, besides density and spatial distribution of PM cells, we introduce two additional variables: (a) automaticity strength and (b) anisotropy. Automaticity strength is defined as the ability of a pacemaker cell to maintain robust spontaneous activity with fast rate and to drive neighboring quiescent cells. It is strongly related to the amplitude of the net inward ionic current into the PM cell during the late diastolic period and the rising phase of the action potential (AP). For example, adding fetal bovine serum to monolayer cultures of NRVMs "strengthens" automaticity, i.e. favors higher firing rate, by upregulating inward long-lasting activation calcium current I CaL [16]. Early versions of engineered BPs have been created from quiescent monolayer cultures or quiescent in vivo CMs via the use of different techniques to upregulate inward pacemaking current I f [17]. The second newly introduced variable, anisotropy, can be created in cultures of NRVMs via several methods, notably by patterning the culture substrate [18,19] or directly seeding the cell into a thin slice of decellularized cardiac tissue [20]. These methods usually lead to functional cardiac network with elongated cell and faster propagation in the longitudinal direction [20]. It has been proposed that linear anisotropy could facilitate BP function [21]. However, the underlying mechanism remains unclear since most studies do not assess specific effects of anisotropy on spontaneous activity but instead focus on contractile function [20], electrical activation [22], or orientation-related response to stretch [23].
This study aims to assess modulation effects of automaticity strength and anisotropy on the spontaneous activity of cardiac monolayers with various densities and spatial distributions. The non-linear relationship between those two variables and automaticity will be characterized with simulation methods and discussed in details.

Cardiac network model
Semi-discrete microstructure models are more suitable than continuous models when individual cell sizes, shapes and orientations are variables under investigation [24,25]. For this reason, a previously described semi-discrete microstructure model [26] was used to simulate two 2D network geometries corresponding to isotropic and anisotropic monolayers. The two network geometries were identical in all aspects, except: (a) aspect ratio of cells (AR, length divided by width of the cell), and (b) distribution of gap junctions. As summarized in Table 1, a grid of 920 x 920 nodes was created with 6 μm resolution, and assigned to 42,642 CMs to create a 5.5 mm x 5.5 mm monolayer.
Each cell included~20 nodes. CMs for anisotropic geometry had an average AR of 3 compared to 1 for isotropic geometry. Longitudinal and transverse intercellular conductivities were adjusted to fit experimental conduction velocities found in NRVMs monolayer cultures [18]. The experimental isotropic conduction velocity was reported to be 16.8 ± 2.1 cm/s in all directions; and for anisotropic monolayer cultures, the longitudinal and transverse conduction Table 1. Characterization of the monolayers: isotropic vs. anisotropic. Isotropic and anisotropic monolayers are identical for most of the features. The cells differ only in aspect ratio and intercellular conductivities.

Isotropic Anisotropic
Total number of nodes (#) 920 x 920 920 x 920 Resolution (μm) 6 6 Length of monolayer (mm) 5.5 5.5 Width of monolayer (mm) 5.  [27] was used to represent the CMs, with the application of a constant inward bias current (I bias ) to generate spontaneous activity [28,29]. The LR1 model is simulated at every node of each cell, and examples of APs and total ionic currents obtained in a single cell with I bias = 2.6 μA/cm 2 and I bias = 3.5 μA/ cm 2 are illustrated in Fig 1A and 1B  This cell model was chosen because its bifurcation structure related to oscillatory behavior has been fully characterized. Indeed, bifurcation analysis undertaken with AUTO continuation software [30] is displayed in Fig 1C. The S-shape curve of fixed points has a lower and upper branch connected by an intermediate branch of unstable fixed points. Both the lower and upper branches change stability through subcritical Hopf Bifurcations (H1 and H2, magenta square). Stable cycle exist in between the two Cycle Saddle nodes bifurcation (blue lines from SNC1 and SNC2). At high I bias , the branch of unstable cycles created at SNC2 (green line) connect the Hopf bifurcation H2. The branch of unstable cycles created at SNC2 exist only on a small interval of I bias and ends through a Homoclinic bifurcation with the intermediate branch of fixed point. Similarly, the branch of unstable cycles created at H1 exist on a tiny interval of I bias and also disappears through a Homoclinic bifurcation with the intermediate branch of fixed point. The cycle length of stable spontaneous activity decreased with I bias and ranged from 1989 ms to 516 ms ( Fig 1D).

Stochastic distribution of pacemaker cells
The 2D cardiac network was assumed to contain two populations of cells: PM (I bias = 2.6 μA/ cm 2 or I bias = 3.5 μA/cm 2 ) and quiescent (I bias = 0 μA/cm 2 ) excitable cells. The density of pacemaker cells (D aut 2 [0,1]) was defined as the percentage of PM cells within the network. The spatial distribution was dependent on a variable (p thr 2 [0,1]) determining how homogeneous PM cells were spread in the network. Fig 2 provides intuitive disambiguation between density and spatial distribution. Density described the total number of PM cells in the network, regardless of their scattering within the network.
Spatial distribution described the scattering of PM cells, and hence completed the spatial information from the overall density. Both variables were required to fully specify a spatial pattern. Two networks might have the same density of PM cells but different spatial distributions, or conversely they might share similar spatial distribution of PM cells but had different densities. Higher values of D aut and p thr correlate with higher percentage and more homogeneous spatial distribution of PM cells respectively. A previously described stochastic algorithm [14] was modified and implemented to randomly attribute positions to PM cells in the isotropic and anisotropic networks. Aggregation occurred when a PM cell was placed in the immediate neighborhood of another PM cell. Otherwise nucleation was said to occur, i.e. the PM cell occupied a position where only quiescent CMs were in its immediate neighborhood. Any random process of the algorithm followed continuous uniform probability distribution. For each pair (D aut , p thr ), the algorithm worked as stated below and as illustrated in Fig 3: Step 1: randomly place the 1 st PM cell in the network (equal probability for each cell of the network); Step 2: determine M 1 (aggregation available sites) and M 2 (nucleation available sites); Step 3: randomly generate decision number p 2 [0,1]; Step 4: perform aggregation in M 1 (if p > p thr ) or nucleation in M 2 (if p p thr ); Step 5: perform aggregation in M 1 (if M 2 is empty); Step 6: repeat steps 2 to 5 until the required D aut is reached. Automaticity strength, anisotropy, and cardiac spontaneous activity Default I bias value for all cells was 0, i.e. all cells are quiescent unless set otherwise. Placing a PM cell consisted in setting I bias for all nodes of a cell to a single non-zero value inside the interval [2.6-4.7] μA/cm 2 . The same procedure was used to distribute spontaneous cells in the isotropic and anisotropic cell layouts.
In Fig 4 are presented the first 50 x 50 nodes of two geometries. Cells in network with isotropic geometry demonstrated no preferential orientation compared to cells in network with anisotropic geometry which are clearly oriented along the longitudinal axis.
Within the cardiac 2D network, a PM cluster was a subgroup of interconnected PM cells. Two PM cells were considered interconnected if they shared gap junctions. The size of a PM cluster was the number of PM cells in that cluster, and the maximum PM cluster size S cluster was the size of the network's biggest PM cluster. Porosity was the fraction of quiescent cells in a PM cluster. In fact, any subgroup of interconnected PM cells was a PM cluster, but quiescent cells might also be enclosed within the cluster, i.e. totally surrounded by the cluster's PM cells. Given S Tcluster , the average of maximum PM cluster size including both PM and quiescent cells, porosity was defined as follows: Porosity: average fraction of quiescent cells in the largest PM cluster S cluster : average of maximum PM cluster size, counting only PM cells S Tcluster : average of maximum PM cluster size, counting both PM and quiescent cells

Simulation protocols and data analysis
As previously described [26], a 2D monodomain approximation with fine discretization was used to formulate the microstructure model. No-flux boundary conditions were applied to the four sides of the network. Initial conditions for all cells corresponded to the resting state of quiescent cells. The total simulation duration was 10 s, and the steady-state behaviors were reached rapidly within two to three autonomous period for most spontaneous cases, the longest transient behaviors found at the transition between non-autonomous to autonomous multicellular activity. Analysis was done on simulations after removing the first action potential thus including the time from the 2 nd action potential up until 10 s. Simulations were performed to study the effect of D aut and p thr on the spontaneous activity of 4 groups: Post-simulation analysis was performed in Matlab (The Mathworks, Natick, MA). The network was said to have spontaneous activity if two complete activations or more were detected during 10 s of simulation. Conversely, the simulation was labeled as non-automatic if a single AP or no activation was detected. The activation time of an action potential (AP) was defined as the time when the transmembrane voltage depolarizes beyond -40 mV. For the i th AP of a given simulation, the activation map (M tact,i in ms) was a matrix constructed from detected activation times for all nodes.
A first set of measures was computed for each 10 s simulation with spontaneous activity.
Ã Normalized activation map (M ntact,i in ms) for the i th AP was obtained as follows: Ã Spontaneous cycle length map (ΔM tact,i in ms) for the i th AP was defined as: The spontaneous cycle length value (Δt act,i in ms) for that map was: And the average spontaneous cycle length Δt act,i in ms) for the series of N APs was: Ã The map of synchronization times (M τsync,i in s/cm) for the i th AP was the inverse of conduction velocities, calculated from the spatial gradient of M tact,i using a previously described method [31]. It quantified activation delays between cells. If all cells were synchronized, meaning they activated at the same time, M τsync,i would be zero for each position on the lattice. The average synchronization time for N APs was: A second set of measures were defined for each pair (D aut ,p thr ) to assess how spontaneous activity behaves for similar spatial patterns and hence the variability between realizations of the same random process of pattern generation. In fact, as previously stated, eight monolayers have been produced for every pairs (D aut ,p thr ). The monolayers produced with the same (D aut , p thr ) had similar spatial patterns, compared to monolayers produced with different (D aut ,p thr ) which had dissimilar spatial patterns. Ã For each pair (D aut ,p thr ), given n simulations (up to 8) with spontaneous activity, average and standard deviation of spontaneous cycle length for similar spatial patterns were defined as: Dt act;j ð6Þ DT act : average cycle length over n simulations for a pair (D aut , p thr ) σ ΔTact : standard deviation of cycle length over n simulations for a pair (D aut , p thr ) Dt act;j : cycle length for the j th simulation, as stated in Eq (4) Ã For each pair (D aut ,p thr ), average and standard deviation of synchronization time for similar spatial patterns were defined as: T sync : average synchronization time over n simulations for a pair (D aut , p thr ) σ Tsync : standard deviation of synchronization time over n simulations for a pair (D aut , p thr ) τ sync,j : synchronization time for the j th simulation, as stated in Eq (5) The same process is used to calculate T sync;x , T sync;y , longitudinal and transverse components of T sync .
The third set of measures were defined for dissimilar spatial patterns, i.e. monolayers with different values of (D aut ,p thr ).
Ã The cycle length range for dissimilar patterns was defined as the percentage variation between the minimum and maximum values of T sync over all pairs of (D aut ,p thr ): DT act : average cycle length over n simulations for a pair (D aut , p thr ), as stated in Eq (6) Ã The synchronization time range for dissimilar patterns was defined as the percentage variation between the minimum and maximum values of T sync over all pairs of (D aut ,p thr ): T sync : average synchronization time over n simulations for a pair (D aut , p thr ), as stated in Eq (8) Furthermore, each monolayer was divided in two areas separated by a square of 2.75 mm side (half of the side of the monolayer) to distinguish between border and central foci (first initiation site of AP). Border foci behavior was assessed by anisotropy ratio. The border foci anisotropy ratio (r) was defined as: η L : number of border foci in longitudinal x-direction (anisotropy direction) η T : number of border foci in transverse y-direction

Characterization of the cardiac 2D network
S cluster , average of S cluster over 8 monolayers for each pair (D aut ,p thr ), increased with D aut , independently of p 1=4 thr , with a transition from below 10,000 PM cells to over 30,000 PM cells around D aut = 0.5 (Fig 5A and 5D).
The extent of the transition phase, i.e. the number of pairs (D aut ,p thr ) with S cluster between 10,000 and 30,000 PM cells, correlated with increased standard deviation of S cluster (Std.S cluster ). In fact, Std.S cluster (Fig 5B and 5E) was below 1,000 PM cells for all pairs (D aut ,p thr ) except for the transition phase where Std.S cluster > 2,000 PM cells. As shown in Fig 5C and 5F, the number of clusters (N clusters ) increased with D aut as long as D aut D aut,max (white solid line), and then decreased for D aut > D aut,max . For each p 1=4 thr , D aut,max was the maximum D aut beyond which M 2 became empty, i.e. there was no more available site to perform cluster nucleation during the creation of the spatial distribution of PM cells. Thus, once D aut,max has been reached, only aggregation was possible for increasing D aut . In Fig 6, a relationship was established between the transition in S cluster and D aut,max .
The maximum of the derivative of S cluster as a function of D aut , i.e. max[S 0 cluster ] = max[ΔS cluster / ΔD aut ], representing the sharpness of the transition of S cluster from below 10,000 PM cells to over 30,000 PM cells, was calculated for all p 1=4 thr and plotted against D aut,max . The sharpness of the transition was inversely proportional to D aut,max (Fig 6D). Differences between isotropic and anisotropic geometries were unsurprisingly negligible for S cluster , N clusters , and D aut,max , since the number of neighbors was approximately the same for both geometries. As such, spatial distribution of spontaneous cells created by our aggregation and nucleation process are not affected by cell preferential orientation.

Occurrence of spontaneous activity
In general, 2D cardiac networks with isotropic geometry demonstrated circular-shaped electrical activation, as illustrated in Fig 7B for the pattern shown in panel a. Networks with the anisotropic geometry typically had ellipse-shape electrical activation ( Fig 7D for the pattern in  panel c).  For example, 51.5% of pairs (D aut ,p thr ) demonstrated [n = 8] in ISO-3.5 versus 13.3% in ISO-2.6. Interestingly, proportions of pairs (D aut ,p thr ) with [0<n<8] were very similar for all four groups (~9%). Automaticity was thus more likely to be observed for higher values of D aut and lower values of p 1=4 thr , and no difference in occurrence of autonomous activity was found between isotropic and anisotropic geometries.
A S cluster combined with porosity offered crucial insights on the morphology of the transition curves. Typically, automaticity did not appear when S cluster was below 5,000 PM cells and when porosity was over 0. 35  For all 4 groups, transition curves to [0<n<8] were subtracted from transition curves to [n = 8], and the differences were displayed in Fig 9E and 9F. Peak differences were higher for ISO-2.6 (0.50) and ANISO-2.6 (0.45) compared to ISO-3.5 (0.35) and ANISO-3.5 (0.40) respectively. Furthermore peak differences occurred for lower p 1=4 thr for ISO-2.6 and ANISO-2.6 (~0.2) compared to ISO-3.5 and ANISO-3.5 (~0.40).

Rate of spontaneous activity
In Fig 10A-10D, DT act was calculated for each pair (D aut ,p thr ) with [n = 8] and displayed as color scale map vs. D aut and p 1=4 thr . Independently of groups and for any given p 1=4 thr , DT act decreased with increasing D aut . And ranges, i.e. intrinsic variability for dissimilar patterns, of DT act for ISO-2.6 (16.19%) and ANISO-2.6 (14.56%) were much smaller than ISO-3.5 (98.86%) and ANISO-3.5 (100.96%).
In Fig 10E is presented mean[DT act ], calculated as the average DT act over all pairs (D aut ,p thr ) with [n = 8] shown in Fig 10A-10D. Detailed values can be found in Table 2.
Mean[DT act ] for ISO-2.6 and ANISO-2.6 were respectively 110.26% and 110.96% higher compared to those of ISO-3.5 and ANISO-3.5. This behavior was not surprising, considering the difference observed in single cell simulations (1428 ms for I bias = 2.6 μA/cm 2 compared to 599 ms for I bias = 3.5 μA/cm 2 , in Fig 1D). Besides, it is important to notice that, for a given I bias value, single pacemaker cells always display lower cycle length than monolayers. No difference in mean[DT act ] was found between isotropic and anisotropic geometries for [n = 8], but focusing on the set with [0<n<8] yielded interesting results. Pooling together all pairs (D aut , p thr ) satisfying the condition [0<n<8] ( Table 2), there was a higher mean[DT act ] ANISO-2.6 versus ISO-2.6 (3.08%), and ANISO-3.5 versus ISO-3.5 (5.50%).
Similarly to DT act , σ ΔTact was calculated for each pair (D aut ,p thr ) with [n = 8]. Mean[σ ΔTact ], the average value of σ ΔTact over all pairs (D aut ,p thr ) with [n = 8] is also shown in Fig 10F and Table 2 for all groups. Mean[σ ΔTact ] was a measure of intrinsic variability between monolayers with similar spatial patterns of PM cells. Mean[σ ΔTact ] for ISO-2.6 and ANISO-2.6 were respectively 78.52% and 116.68% higher compared ISO-3.5 and ANISO-3.5. A difference between isotropic and anisotropic geometries was found where mean[σ ΔTact ] in ANISO-2.6 was 30% higher vs. ISO-2.6, while a more limited increase of 8% was found for ANISO-3.5 compared to ISO-3.5.

Spatial characteristics of spontaneous activity
Difference in spatial characteristics of spontaneous activity between isotropic and anisotropic geometry was evaluated. The position of foci (i.e. first initiation sites of electrical activation) was estimated to determine if differences existed between the two network geometries. Independently of geometry, focal activation was highly stable in time for a given spatial pattern of spontaneous cells, demonstrating no beat-to-beat variability. For direct comparison between geometries, the focal position of the last simulated spontaneous beat was selected for all simulations with automaticity (i.e. [n>0]). Pooled positions were plotted in Fig 11A-11D.
The proportion of central foci over all foci is shown in Fig 11E and 11F and Table 2  The border foci in the longitudinal x-direction were the foci located outside the red box, exclusively to the left and to the right. The border foci in the transverse y-direction were the foci located outside the red box, exclusively at the top and the bottom. Non-exclusive border foci at the corners, i.e. foci that are common to longitudinal and transverse directions (blue areas in Fig 11A-11D), were not considered in the calculations. Ratio with values greater than one suggested that there were more border foci in the longitudinal direction compared to the transverse direction. Values of r were calculated for [n = 8] and [0<n<8] and are presented in Table 2. Interestingly, r in anisotropic geometry are consistently greater than one and always higher compared to the value obtained in isotropic geometry. In fact, for [n = 8] case, r raised by 27% from ISO-2.6 to ANISO-2.6, and by 204% from ISO-3.5 to ANISO-3.5. For [0<n<8], r raised by 297% from ISO-2.6 to ANISO-2.6, and by 119% from ISO-3.5 to ANISO-3.5. No clear preference in border foci position was found for the isotropic network.
Synchronization of electrical activation is an important marker of spontaneously beating multicellular monolayer, either in silico or in vitro. T sync was displayed as color map vs. D aut and p 1=4 thr for each pair (D aut ,p thr ) with [n = 8] in Fig 12A-12D.   No particular tendencies were observed in T sync map for ISO-2.6 and ANISO-2.6. For ISO-3.5 and ANISO-3.5, and for D aut ! 0.8, T sync increased and then decreased as a function of p 1=4 thr . I bias = 3.5 μA/cm 2 led to much higher ranges compared to I bias = 2.6 μA/cm 2 . Indeed, range[T sync ] for ISO-2.6 (6.5%) and ANISO-2.6 (16.9%) were much smaller than ISO-3.5 (414%) and ANISO-3.5 (485%). Moreover, the increase in ranges could be noticed in anisotropic geometries versus isotropic geometries: 160% increase from ISO-2.6 to ANISO-2.6 and more moderate 17.15% increase from ISO-3.5 to ANISO-3.5.
Anisotropy consistently yielded higher synchronization times.

Discussion
To our knowledge, this paper presents the first in silico study in a microstructure model describing how automaticity strength and structural linear anisotropy may modulate the effects of density and spatial distribution of PM cells on the spontaneous activity of the BP. Our previous study [14] demonstrated that not only the density but also the spatial distribution of PM cells may induce important intrinsic variability in the BP dynamical behavior. The study focused on a simple continuous and isotropic 2D substrate with a single period of activity for the autonomous cells. These results quickly raised a very important question: can this intrinsic variability be limited despite the lack of control on density and spatial distribution of autonomous cells? Two straightforward solutions may be considered: (a) either the BP is tested ex vivo and then implanted into the myocardium, or (b) the current methods are combined with new processes intended to minimize the intrinsic variability. Investigating the second candidate solution is the main motivation behind the mathematical modeling optimization process that is underway. Here, we studied how different automaticity strengths and structural linear anisotropy can influence the spatial-temporal activity of the BP.
The main new contributions are described as follows: 1. Spatial patterns of PM cells were mathematically defined and characterized within a semidiscrete 2D model, describing both automaticity strength and structural linear anisotropy.
2. Automaticity strength enhanced occurrence of spontaneous activity, decreased variability in cycle length DT act and synchronization time T sync for similar spatial patterns of PM cells, but increased the number of central foci, and the variability of DT act and T sync for dissimilar spatial patterns PM cells.
3. Structural linear anisotropy had no important effect on occurrence of spontaneous activity, increased variability in DT act and T sync for both similar and dissimilar spatial patterns of PM cells, and decreased the proportion of central foci.
4. Intrinsic variability was modulated but not eliminated by neither automaticity strength nor structural linear anisotropy, since there was still a proportion of pairs (D aut ,p thr ) with no corners, i.e. foci that are common to longitudinal and transverse direction are in the blue areas and are not considered in the calculation of border foci anisotropy ratio (r) in Eq (12). spontaneous activity and important performance discrepancies notably for dissimilar spatial patterns of PM cells.
PM cells were randomly placed in a semi-discrete 2D microstructure via a stochastic algorithm with a parameter D aut controlling density and a parameter p 1=4 thr determining homogeneity of spatial distribution. Two levels of automaticity strength (weak and strong) were achieved by two different values of I bias (2.6 μA/cm 2 and 3.5 μA/cm 2 ). Linear anisotropy was structurally created by: (1) geometrically increasing aspect ratio of cells from one to three, and (2) fitting longitudinal and transverse gap junction distribution to published conduction velocities from anisotropic monolayer cultures of NRVMs [18].
Isolated PM cell dynamics were described by LR1 cardiac ventricle myocyte model [27], with a constant inward bias current to generate automaticity, as described elsewhere [28,29]. More complex ionic models (for example the published NRVM ionic model [32] or a mathematical models representing stem cell-derived human cardiomyocytes [33]) were not used: (1) to avoid important computational cost because of the associated long transient dynamics and (2) because increasing dimensions of mathematical ionic models does not necessarily lead to more predictive power. As a matter of fact, in many high dimension models very different set of parameters can lead to the same AP; in this case the model is said to have identifiability issue, a major flaw to predictive reliability [34]. The LR1 model obviously does not include the detailed ionic currents of cardiomyocytes (pacemaking and resting) nor can it reproduce the intricate interaction between the membrane/voltage clocks [6,28]. However it remains an appropriate trade-off option since it is identifiable [34] and can reproduce the basic physiological behaviors considered in the present study, namely automaticity and AP propagation.
To avoid confounding effects, some simplifications were made: 1. We considered only two populations of cells: PM and quiescent excitable cells. Other types of heterogeneities such as sinks (spots with voltage fixed at resting potential like fibroblasts) and breaks (spots with no conduction, like monolayer damage) were not considered although they may have had an effect on conduction [35].
2. I bias was the same for all PM cells of a specific monolayers, i.e. all PM cells were identical. All cells had the same initial conditions and no time delay. It is important to notice however that clinical BPs with all identical PM cells would lack robustness to external perturbations. For example they may completely stop firing after acetylcholine stimulation. Cellular diversity is indeed an important aspect to preserve robustness in the native sinus node [36].
3. The number of neighbors was approximately the same for each cell. Different number of neighbors could create electrotonic disparities between PM clusters and induce confusion with effects of porosity on occurrence of automaticity.
No reentry or asynchronous activation was detected in the simulations; this fact is an indication of strong intercellular coupling [28,29,[37][38][39] but also may be limited by the homogeneous initial conditions used in the simulations. But the strongest coupling is not necessarily the best for clinical in situ BPs, without no-flux boundary conditions. In fact, to maintain source-sink balance with the atrium, native SAN cells are coupled with low conductance connexins (Cx-45) instead of high conductance Cx-43 [40]. There is indeed an inverse proportional relationship between coupling strength and the safety factor of propagation [41].
S cluster and N cluster as a function of D aut and p 1=4 thr displayed remarkable similarities with our previous study [14], despite different model structure and different number of neighbors. In fact, average of the maximum PM cluster size S cluster monotonically increased with D aut , with a transition around D aut = 0.5 (Fig 5). However, in this study we brought deeper spatial characterizations compared to the previous one. For example, the sharpness of this transition was demonstrated to be inversely proportional to D aut,max (Fig 6). The sharpness of the transition relied on two related phenomena: aggregation of clusters at low p 1=4 thr and «fusion» [14] of clusters at high p 1=4 thr . For low p 1=4 thr (heterogeneous spatial distribution of PM cells) S cluster grew mainly because of aggregation and was proportional to D aut which increased linearly, hence showing a smooth transition. For higher p 1=4 thr (homogeneous spatial distribution) S cluster grew mainly because of «fusion» since there was more nucleation and less aggregation. In fact, as D aut increases, clusters start having common neighboring available sites. When a PM cell is placed on one of those sites, previously separated clusters «fuse» into one bigger cluster. Decrease of N cluster indicates a start in the «fusion» process, systematically above D aut,max . So S cluster started growing slowly for high p 1=4 thr and low D aut because there was almost no aggregation and no «fusion». Once D aut,max was reached, «fusion» process provoked a sharp increase of S cluster and hence a sharp transition.
Higher D aut and lower p 1=4 thr led to more chance of having automaticity ([n>0]) (Fig 8). Even higher D aut and lower p 1=4 thr were necessary for all simulations to show automaticity ([n = 8]). The lower p 1=4 thr was, the lower D aut was required to generate automaticity. The transition curves are also a major characterization improvement, compared to our previous study [14]. Superimposing transition curves over S cluster and porosity color scale maps (Fig 9) gave crucial information about relationship between the monolayer microstructure organization and occurrence of automaticity. It was really interesting to observe that small S cluster at low p 1=4 thr could generate spontaneous activity while large S cluster at high p 1=4 thr could not. A big cluster at high p 1=4 thr may contain too much inside spots with quiescent cells and not be able to fire an AP, while a small cluster at low p 1=4 thr containing no inside quiescent spot at all will be able to do it. This is where porosity remarkably came into play. Porosity established the importance of electrotonic source-sink balance inside the PM cluster itself.
Spontaneous activity occurs more often with strong automaticity (I bias = 3.5 μA/cm 2 ). Basically, lower S cluster was needed and higher porosity was tolerated for automaticity with stronger PM cells to be observed.
For each p 1=4 thr value, average autonomous cycle length DT act decreased with increasing D aut (Fig 10) since bigger clusters (with reasonably low porosity) have more electrotonic driving force. Unsurprisingly, strong automaticity brought smaller period of activity (lower DT act ), but also resulted in higher ranges of DT act , i.e. more intrinsic variability for dissimilar spatial patterns. Stronger automaticity allowed occurrence of spontaneous activity in more difficult conditions (lower S cluster , and higher porosity), but the resulting periods would be noticeably higher compared to more favorable conditions. On a clinical point of view, this result means that creating BPs with stronger automaticity would not solve the problem of not having control over the fate of the PM cells in situ (intrinsic variability). The generated BP may have better chance to display automaticity, but would have very disparate performances from one patient to another. The best way to bypass the problem is to create multiple BPs ex vivo, pick the cultured BP with the targeted performance and then graft it to the patient myocardium.
In accordance with the literature [14,42,43], foci were usually located at the border (Fig 11) since, with no-flux boundary condition, border clusters had less cells to depolarize. Strong PM cells led to more central foci because more clusters had enough electrotonic driving force to depolarize neighboring cells without interacting with the no-flux boundary condition. It was interesting to notice that anisotropy yielded more border foci in the longitudinal direction compared to the transverse direction. In fact, because of the aspect ratio, there were actually more cells at the two longitudinal borders and thus more chance of having a sufficiently bigger S cluster to fire an AP.
Electrical activations were less synchronous in anisotropic monolayers (Fig 12). Very low T sync in the longitudinal direction were hindered by very high T sync in the transverse direction, leading to an electrical activation that was globally less synchronous. The evolution of T sync as a function of p 1=4 thr for D aut ! 0.8 is very interesting and requires further investigations. Further investigations are also needed to explain how structural linear anisotropy consistently induced more intrinsic variability for both similar and dissimilar spatial patterns. Anisotropy had also limited effect on occurrence of automaticity and DT act , compared to strength of automaticity. If there was indeed an effect, it eventually happened below the 5% resolution of D aut and p 1=4 thr . Automaticity strength and anisotropy may not be as independent as displayed in this study. Anisotropy has been shown to induce changes in intracellular calcium transients ([Ca 2+ ] i ) dynamics, decreasing the diastolic [Ca 2+ ] i levels and increasing [Ca 2+ ] i influx per cardiac cycle [44]. Anisotropy may also alter the properties of voltage-gated ion channels, notably expression and regulatory properties of voltage-gated calcium channels [45]. So, on a strict structural point of view (i.e. aspect ratio and distribution of gap junctions), anisotropy may not have an impact, but it may indirectly affect automaticity strength, which then will have an effect on spontaneous activity. As such, even if linear structural anisotropy does not increase occurrence of automaticity and does not decrease intrinsic variability, it may still be clinically interesting.
This study offered crucial insights on the relationships between the microstructure of the BP and its macroscopic behavior. Automaticity strength and structural linear anisotropy may modulate effects of density and spatial distribution of PM cells on the spontaneous activity of the BP, but will probably not eliminate intrinsic variability among BPs. Therapeutic procedures that do not take characteristics of spatial distribution into account (eg. cell and gene therapies) may end up with a non-negligible intrinsic variability, despite standardized protocols. Increasing the number of pacemaker cells may not solve the issue. As a matter of fact, the native SAN tissue exhibits a specific architecture with gradual transition of intermediate cells and gradually decreasing density of pacemaker cells from the center to the periphery. That spatial arrangement promotes automaticity by maintaining the delicate balance in the source-sink relationship between the SAN and the surrounding atrial tissue [36]. With an unknown spatial arrangement, that balance may randomly be compromised, hence the intrinsic variability. Furthermore, unlike the SAN which is electrically isolated from the rest of the atrium with the exception of few exit pathways, BPs lack electrotonic barrier at the macroscopic level, another challenge to their performance. All those facts stress the importance of ex vivo design and performance assessments on BPs in bioreactors before implantation to patients.
A sheet-based BP as a replacement to the normal sinus node has important intrinsic differences with the physiological structure. The normal pacemaker is a compact node of heterogeneous spontaneous cells believed to be connected at specific exit points thus limiting contact to a restricted number of atrial cells [46]. The simulated BP in this study are very different where a population of autonomous cells are connected to resting ventricular myocytes, both having similar morphology and dimensions. Cell dimensions and morphologies vary between pacemaker cells (assuming spindle-like cell morphology as in the sinus node) and resting but excitable cardiomyocyte (with also differences between atrial and ventricular cells). Cardiomyocyte dimensions is well known to affect electrical propagation [47]. Thus, important work is still needed to study the differences in spontaneous activity of BP monolayer when considering morphological differences between cell types although the differences between cell types (sinus node cell, atrial and ventricular-like derived cardiomyocytes) for derived cardiomyocytes may not be as important as in adult hearts. More importantly, BP activity when electrically coupled to myocardium in order to drive the tissue will be depressed by the electrotonic effect [48] and biomimetism of the sinus node coupling structure to the atria may need to be considered for coupling the BP to the myocardium.

Conclusion
In summary, a pure change from isotropic to anisotropic substrate modelled by an elongated cell shape and anisotropic intercellular conductivity without modifications of ion channel expression nor spatial distribution has limited effects on spontaneous activity. However, increasing the intrinsic rate of autonomous cells has a much stronger effects. Although the two were studied together and independently, it is of importance to note that there is a strong possibility that both changes (anisotropy and autonomous strength) could be coupled [44]. Further work is thus needed to uncover this importance of the interaction (and by how much methods to induce cellular anisotropy can increase the cellular automaticity strength) and to elucidate how it could favor the BP development.