The effect of quenched disorder on dynamical transitions in systems of coupled cells

Non-equilibrium systems are characterized by a rich variety of dynamical behaviors, which may sensitively depend on control parameters. Here, we investigate and provide a quantitative analysis of the role of disorder on the transitions between different dynamical regimes in extended heterogeneous systems of excitable and passive cells, induced by varying the strength of the coupling between cells. The random distribution of passive cells provides a quenched disorder in important biological contexts, such as the appearance of contractions in the pregnant uterus. We observe a large variability between different realizations of the disorder (replicas) in a lattice of excitable cells, each cell being coupled to a random number of passive cells. The statistics of these large disorder-induced fluctuations are related to the properties of the coarse-grained distribution of passive cells, in particular, to its extreme values. We show that these fluctuations can be characterized by a simple scaling relation, involving the strength of the coupling between excitable cells, the mean passive cell density and the logarithm of the system size. Our results provide a quantitative understanding of the important effect of a quenched disorder in the transition between dynamical regimes in extended dynamical systems.


Introduction
Varying control parameters in a wide class of non-equilibrium systems, in particular in physical contexts, can lead to very different dynamical regimes [1,2]. In biological systems, the transition to a regime of rhythmic activity is important to determine the functioning of crucial physiological properties, with a wide range of observed periods [3][4][5][6][7][8][9][10][11]. The case of the pregnant uterus close to delivery leads to a particularly intriguing problem. None of the cells of the tissue spontaneously oscillate [12][13][14]; the origin of the transition to a regime with strong contractions is not fully understood. This is a serious impediment in any attempt to control pre-term contractions, which are a major cause of infant mortality [15]. It has been noticed that the electric coupling between cells very significantly increases shortly before delivery [16][17][18]; many results provide evidence that the coupling between cells plays a crucial role in the triggering of rhythmic contractions [19][20][21].
The presence of disorder in the structural composition in many realistic natural systems could affect the transitions between dynamical regimes, as it does affect the equilibrium phase transitions [22,23]. Such disorder is frequently seen in biological systems, including the uterine muscle which has a heterogeneous composition comprising electrically excitable myocytes (smooth muscle cells) that form the bulk of the tissue and a small fraction of electrically passive cells, like fibroblasts and telocytes [13,15,24]. The passive cells are much fewer compared to the myocytes and experimental evidence suggests that the (frozen) pattern of their connections to excitable cells is highly irregular [15], and can be treated as a quenched disorder in the system. This work is concerned with the emergence of oscillations from a quiescent state in such a disordered system of coupled excitable and passive cells [19,20,[25][26][27][28][29][30]. The main contribution of this paper is to provide a quantitative understanding and analysis of the role played by disorder in the transition between dynamical regimes in a simple model, proposed to describe the pregnant uterus [21]. Our approach should also lead to important insight to understand other models of biophysical significance [30].
Our model, inspired by the heterogeneous architecture described above, comprises a lattice of diffusively coupled excitable cells, each of which may be connected to a number of passive cells. As observations made on biological tissue do not reveal any regular organization in the connection density between excitable and passive cells [15], we consider the number of passive cells connected with an excitable cell to be randomly distributed about a given mean density. A specific distribution of passive cells attached to the myocytes effectively represents a particular realization of quenched disorder in the system. Numerical exploration of the model dynamics exhibits a variety of spatio-temporal regimes at different values of inter-cellular coupling and average density of passive cells [21], reminiscent of what has been observed in other systems of coupled electrically excitable and passive cells [20,31,32]. More generally, our results are relevant for a broad class of systems that comprise coupled heterogeneous elements, e.g. in the cardiac context [33,34]. In qualitative agreement with biological observations of uterine tissue, it is seen from the model that strong cellular coupling promotes synchronized oscillations [21].
While, in earlier work [21], we have provided a global description of the overall dynamics of the model system, this paper provides a detailed understanding of the genesis of the different dynamical regimes by revealing the key role played by disorder in the transition between them. In particular, we show that local fluctuations in passive cell density can result in the self-organized emergence of 'pacemaker-like regions', i.e. clusters of coupled cells that spontaneously oscillate, eventually activating the entire system. We stress that these are very different from the specialized pacemaker cells of the heart whose role cannot be adopted by any of the other cells in cardiac tissue under normal circumstances. In contrast, individual members of the self-organized groups of oscillating cells referred to above are not inherently distinct from other cells in the system. Thus, unlike the heart, where pacemaker function can be thought of as an intrinsic property of certain specialized cells, in the uterus this function is an emergent property of interactions between heterogeneous cell types.
Variability in dynamical transitions has been observed in spatially extended excitable media, e.g. in certain chemical systems [35][36][37][38], even in the absence of any spatial disorder. In such problems, the regimes observed over long times depend on the initial state of the system, as well as on the boundary conditions. While such effects are also observed in our model, this study focuses on how a replica, i.e. a particular realization of the disorder, affects the transition between the dynamical regimes observed in the system. Specifically, we consider fluctuations between different replicas in the macroscopic behavior of the system, e.g. in the emergence of oscillations. We use the earlier observation [19,21] that when excitable and passive cells are coupled with strength C r , the combination is capable of oscillating provided the number density of passive cells f coupled to an excitable cell is within a critical range: A key aspect of the model is that the transition between quiescent state and oscillatory activity occurs through a subcritical bifurcation [39,40], such that the oscillations always have a finite amplitude. For f f l c , the system exhibits very strong fluctuations from replica to replica as the coupling between excitable cells D is varied. This is because close to f l c , it is not only the mean value of f for the entire system, but also, the spatially fluctuating values of the local passive cell densities that are important. The definition of the local density involves a coarse-graining length, d, which is expected to be related to the diffusive coupling between excited cells as d ∝ √ D on general physical grounds. Therefore, whether a system oscillates for given values of C r and D will depend critically on the details of the passive cell distribution in a specific realization. More specifically, oscillating groups will emerge from regions where the local passive cell density is greater than f l c . The probability of these extreme events, i.e. the occurrence of such regions through statistical fluctuations, increases with system size. Based on these arguments, we derive a scaling relation of the system dynamics as a function of system size and cellular coupling strengths that we have verified through numerical simulations.
The paper is organized as follows. In section 2 we describe the model system followed in section 3 by a discussion of the properties of an excitable cell coupled to a variable number of passive cells. Section 4 reports our numerical observations of fluctuations in the system behavior across different replicas and explains how regions with high local passive cell densities can act as organizing centers for oscillatory activity (i.e. 'pacemakers') in the system. Section 5 contains the key theoretical result of our paper where we derive the probability of occurrence of such pacemaker regions as a function of different system parameters. We show that it has a logarithmic dependence on system size and obtain a scaling relation. Finally, in section 6, we conclude with a summary of our results. We have included a list of the variables and symbols used throughout the paper in the appendix for ease of reference.

The model
As mentioned above, the system under investigation comprises electrically excitable as well as passive cells. The dynamics of excitable cells can be described by equations having the general form where V e (mV) represents the transmembrane potential, C m = 1 µF cm −2 is the membrane capacitance density, I ion (µA cm −2 ) represents the total density of currents transported across the cell membrane by different ion channels, g i are variables describing the gating dynamics of the ion channels and e corresponds to external stimulation that could be due to coupling with neighboring cells. While the explicit form for the ionic current density I ion depends on the type of cell being considered and the level of biological realism desired, in this paper we focus on phenomena that occur independently of the finer-scale details of specific models and therefore use the generic FitzHugh-Nagumo (FHN) model [11]. The functional form for the ionic current in the FHN system, expressed in dimensionless form, is is a parameter governing the fast activation kinetics, α (= 0.2) is the activation threshold and g represents an effective repolarization current. The time-evolution of g is described by with (= 0.08) corresponding to a relatively slow rate of recovery of the system from the excited state. If the system is capable of oscillation, its period is governed by the time scale ∼1/ [11]. For the chosen set of parameter values, the dynamics of the FHN system involves either relaxation toward a fixed point for a subthreshold perturbation or a finite amplitude excursion before decaying to a fixed point after an interval (action potential) for a strong enough or suprathreshold stimulus. Thus, the active cells are in an excitable regime in which they are not capable of exhibiting oscillations spontaneously (i.e. in the absence of external stimulation). In addition to excitable cells, the system also contains passive cells described by the dynamics of the transmembrane potential V p [41]: so that the activity of the cell tends to decay at an exponential rate to its resting value V R p (= 1.5) with K (= 0.25) specifying the characteristic time scale of relaxation. The term e represents the external stimulus to the passive cell. When excitable and passive cells are coupled to each other via an electrical conductance (representing gap junctions in biological tissue), the external stimulus of each cell corresponds to the input received from the other cell so that the respective external stimuli terms are The strength of coupling is represented by C r while n p corresponds to the number of passive cells coupled to an excitable cell.
To investigate spatially extended patterns that are seen in biological tissue such as that of the uterus, we consider a two-dimensional (2D) system of coupled excitable and passive cells. The system comprises a square lattice of L × L (= N ) excitable cells coupled diffusively to their nearest neighbors on the lattice. The coupling results in an additional term D∇ 2 V e in the term e in equation (1) where D is the effective diffusion constant (the lattice spacing being set to 1) and the Laplacian being numerically implemented by a five-point stencil.
Each excitable cell is coupled to a variable number n p of passive cells, with each passive cell being connected to a single excitable cell. The total number of passive cells in the system is M = N × f , f being the passive cell number density. Each passive cell is assigned to a randomly chosen excitable cell drawn uniformly from the population of N cells. We use a binomial distribution for n p , which at large N converges to a Poisson distribution with mean f . Each realization of M passive cells randomly distributed over N excitable cells, which is an instance of quenched disorder in the system, is referred to as a replica.
We have integrated the system numerically using a fourth-order Runge-Kutta scheme, using the time step dt = 0.05. Periodic boundary conditions have been used to obtain all the results reported here; our earlier results have shown that using no-flux boundary conditions does not result in qualitatively different phenomena [21]. After starting from random initial conditions, the system is first allowed to evolve for 2 × 10 4 time steps (corresponding to transient behaviors) before recording the data.

Single excitable cell coupled to passive cells: the mean-field limit
Before proceeding to investigate the spatially extended system, we consider the dynamical system comprising an excitable cell coupled to one or more passive cells, given by equations (1)-(4) [19][20][21]. This problem, which formally corresponds to zero spatial dimension, can also be viewed as the mean-field limit of extremely strong diffusive coupling between excitable cells, such that spatial fluctuations can be ignored and the lattice effectively behaves as a single excitable cell coupled to the average number of passive cells, f . Even though an excitable cell couples to only an integer number n p of passive cells, the mean f will in general not be an integer. Thus, the integer value of n p in equation (4) can be replaced by a real number f , which corresponds to an exact representation of the system dynamics in the mean-field limit. Figure 1(a) shows that the system is capable of exhibiting spontaneous oscillation over a range of values of C r and f , as has been reported earlier [19]. For a given value of coupling C r , the fixed point solution of the system is unstable when f l c < f < f u c , where the lower and upper bounds can be obtained analytically by linear stability analysis of equations (1)-(4) around the steady state. A simple relation between f l,u c and C r is obtained in the limit of large K , i.e. when the passive cell relaxes rapidly. Under this condition, equation (3) is replaced by an algebraic equation such that e = f K (V R p − V e )/(1 + K /C r ). Thus, as f and C r appear only via e , the critical value of f depends on C r as f l,u c = c l,u + m l,u /C r , where c l,u , m l,u are independent of C r . Note that, although the above argument is strictly valid only in the large K limit, we observe from our numerical results that the dependence of the critical values of f on C r follows the above relation even for small K ( figure 1(a)). Thus, f l c (C r ) decreases monotonically with increasing C r and approaches a finite value as C r → ∞.
While the Hopf bifurcation resulting in loss of stability of the fixed point solution at f = f u c (C r ) is supercritical, the bifurcation at f = f l c (C r ) is subcritical so that oscillating solutions of finite amplitude exist for f < f l c (C r ). For example, for C r = 1.9 we observe oscillations when f > f l c = 0.545 (figure 1(b), inset), while they are never observed for f < f 0 0.529. Thus, for mean passive cell density in the interval f 0 < f < f l c , the system is multistable so that both fixed point and oscillatory solutions can be observed numerically (see the narrow interval enclosed by the dash-dotted line in figure 1(b)).
In this paper we focus on the situation when f < f l c (C r ). As discussed in detail later for this condition we observe large fluctuations between replicas when the system is made to undergo dynamical transitions by increasing the coupling between excitable cells, D. In contrast, when C r is increased so that f > f l c (C r ), increasing the coupling D results in a much simpler transition that approaches the mean-field behavior. The replica variations are most pronounced when f just exceeds lim C r →∞ f l c (C r ), the minimal value of f l c (C r ). Most of the numerical results reported here are obtained with f = 0.5 (figure 1(a)), while the minimal value is f = 0.484.

Two-dimensional (2D) media
Biological organs exhibiting electrical activity that are often functionally important, such as the heart or the uterus, are spatially extended objects, being made of tissue comprising a large number of excitable and passive cells. As mentioned earlier, we have modeled these by a 2D . The corresponding time-averaged spatial correlation, which provides a characteristic size of the active regions, is shown in the second row. The frequencies of individual oscillators in the medium are shown in the pseudocolor plots in the third row (white corresponds to absence of oscillation). The last row shows the local density of passive cells averaged over a length scale d indicated below each frame. For the replica shown in (a), we observe global synchronization at D = 2 followed by progressive cessation of spontaneous periodic activity in the system indicated as a shrinking region of oscillating cells when coupling is increased further. However, for the replica in (b), coherent oscillation is not observed as D is increased, with the existing localized oscillating clusters having distinct frequencies gradually decreasing in size. system of coupled excitable and passive cells with periodic boundaries (figure 2). Figure 3 shows snapshots of activity in two different replicas of such a system as the diffusion constant D is increased. The replicas are obtained by taking two different distributions of passive cells attached to the myocytes. Following the schematic representation in figure 2, in the two replicas each myocyte (large circle) is coupled to a different number of passive cells (small circles) while keeping the total number of passive cells constant. We have investigated the system for a value of the average passive cell density f = 0.5 and coupling strength between excitable and passive cells C r = 1.9. For these values, f < f l c ≈ 0.545, so that in the mean-field limit We define a non-dimensionalized distance µ = ( f l c − f )/ f l c between the system under study and the critical system where the fixed point becomes unstable in the mean-field limit (for f = 0.5, µ ≈ 5.7 × 10 −2 ). For small µ, the behavior of the system is sensitively dependent on the particular replica chosen (figure 3). This is manifested as large fluctuations in the system behavior, measured by various order parameters defined later. For example, while global synchronization where all cells oscillate with the same frequency is seen in one replica (figure 3(a), D = 2), another replica for the same set of parameters exhibits only localized clusters of cells oscillating with different frequencies ( figure 3 This sensitive dependence of the system dynamics on the distribution of n p can also be seen from the (D, C r )-parameter space diagrams of the 2D system for different replicas (figure 4). The determination of such phase diagrams has been used in the general context of non-equilibrium phase transitions [1,2,42,43], including problems of extended excitable media [35][36][37][38], to classify the various dynamical possible regimes. Here, we fix all the parameters that define the cells, and vary the coupling strengths between excitable cells, D and between excitable and passive cells, C r . While for low values of C r , increasing D results in a complete absence of oscillations in the system ('no oscillation' or NO phase), for larger values of C r we observe clusters of oscillating cells ('cluster synchronization' or CS phase) at low D, with each cluster having a characteristic frequency that may differ from other clusters. As D is increased, the clusters gradually synchronize with each other although isolated regions of non-oscillating cells can exist ('local synchronization' or LS phase), until all cells eventually oscillate at the same frequency at a sufficiently high D ('global synchronization' or GS phase). In numerical simulations, the phases are obtained by using suitable order parameters (defined in section 3.3), and choosing appropriate threshold values. As can be seen by comparing figures 4(a) and (b), the diagrams differ quite significantly in terms of the actual dynamical regimes that are observed for the same set of values of D and C r , illustrating the significant variability from one replica to the other.
In obtaining the two phase diagrams shown in figures 4(a) and (b), we have used for each value of the parameters D and C r an initial condition chosen randomly, such that each variable is close to its resting state value. It is known that the dynamical regimes observed in spatially extended excitable systems can depend sensitively on initial conditions and this feature is indeed seen in our system. This is manifested as noisy patches visible in figures 4(a) and (b). It reflects the fact that for some values of the parameters (D, C r ) the dynamics can belong either to a LS or to a CS regime, thereby defining a region of parameter space where the system is multi-stable. We have indeed verified that for several values of C r and D, the nature of the asymptotic regime depends on the precise value of the initial condition. Thus, if other initial conditions had been used, the phase diagrams shown in figures 4(a) and (b) would have looked slightly different in the region where multistability is observed, as revealed by the seemingly random appearance of black (CS) and white (LS) patches (regimes), although the principal features of the diagrams would have been unchanged. Another manifestation of this phenomenon is shown in [21], where it was observed that in the CS regime as one increases D very slowly for a fixed value of C r , traveling waves are obtained in a system with periodic boundary conditions, the wavenumber depending on the initial conditions used.
It may be of interest to note here that the behavior of systems having relatively small size (L = 50) illustrated in figures 4(a) and (b) do not exhibit as rich a variety as seen, for example, in related models of diffuse fibrosis in cardiac tissue [28][29][30]. In these studies, spiral waves are observed to play an important role, and transitions between different regimes involve several intervening states. While we cannot rule out the possibility that spiral waves may be observed in our model for the larger system size, they do not play any role in the phase diagrams shown in figures 4(a) and (b).
By averaging over many replicas, we can obtain a 'mean' phase diagram. The diagram obtained here, figure 4(c), for f = 0.5 is qualitatively similar to the one obtained in [21] for a higher value of f (=0.7). The region corresponding to coherence or 'COH' phase, defined by the condition that all cells in the system oscillate with the same frequency, and with a very small phase difference between them, also exists in the present case, appearing at large value of D which is outside the region of interest of this paper.

Order parameters
For a detailed quantitative analysis of the spatiotemporal dynamics of the 2D system, we have used two order parameters: (i) the fraction of oscillating cells in the medium f osc and (ii) the width of the frequency distribution as measured by the standard deviation σ ν , where ν denotes the frequencies of the oscillating cells. The amplitude of an oscillating cell is obtained from the square root of the integral of power spectral density (PSD) of the corresponding V e time-series, ν being the frequency (>0) at which the PSD is maximum. The fraction of oscillating cells in the system f osc = N osc /N is the ratio of the number of oscillating cells N osc , i.e. those having an The different phases shown in figure 4 are defined in terms of the two order parameters as follows. The LS and GS phases both have σ ν → 0; however, for the former f osc < 1 while for the latter f osc 1. The CS phase has a finite σ ν while the NO phase is characterized by f osc 0. The order parameters for a given system depend on the exact form of the quenched spatial disorder and we can investigate the statistical properties of the distribution of f osc over an ensemble of replicas ( figure 5(a)). At very low values of D, the distribution of f osc peaks about a value close to f osc = f . Upon increasing D, the distribution broadens with a bias toward large values of f osc . A peak at f osc = 1 develops around D 1. Increasing D further, one observes a higher probability for very small values of f osc . At D 2, the distribution has an almost bimodal form, corresponding to realizations with either very few oscillating cells, f osc 1, or almost all cells oscillating f osc ≈ 1. At yet larger values of D, the peak at f osc ≈ 1 disappears, and the distribution concentrates around f osc = 0, implying a strong reduction in spontaneous activity of the system. This is reflected in figure 5(b), which shows the ensemble average f osc (where denotes an averaging over replicas) and the fluctuations about the mean (inset) as a function of the coupling D for 2D systems with L = 50 and 100 at C r = 1.9. We observe that at large values of D, the fraction of oscillating cells in the system decreases on average with D and hence is consistent with the mean-field result. The standard deviation of f osc becomes extremely high in the range 1 < D < 2 which is a signature of the large fluctuations between different replicas. Comparing between the two curves for L = 50 and 100 shows that the larger system has a higher probability of showing oscillations than the smaller one for the same parameter values. This dependence on the size of the medium is in fact systematic and is discussed in section 5.

Fluctuations of passive cell density and local activity
The behavior described in the previous section can be qualitatively understood by noticing that diffusion effectively couples excitable cells that are within a distance ∝ √ D of each other. As a consequence, one can expect that the behavior of a given excitable cell depends on the density of passive cells in its neighborhood, characterized by the distance √ D. With this motivation, we begin by describing the procedure for computing the local density of passive cells surrounding a given excitable cell.

Averaging procedure
We define local passive cell densityn d coarse-grained over a length scale d as the convolution n d = K d * n p of the spatial distribution of passive cells n p (i, j) (i, j = 1, . . . , L) over the lattice with a 2D averaging kernel K d (i, j). For simplicity, we have considered separable kernels, i.e.
As the coupling term is effectively described by a Laplacian, a natural choice for the coarse-graining kernel K d is the Gaussian kernel K G d given by which has been used for the numerical investigation of the model system presented in this section. For analytical convenience, we have also used the simpler 'square top hat' kernel, K STH d in section 5, which is based on the classical 1D top hat filter: where H (.) is the Heaviside step function, i.e. H (x) = 0 when x 0 and H (x) = 1 otherwise. For the square top hat kernel,n d (i, j) is the number of passive cells, averaged over a square subregion of size N d = d 2 centered at the point (i, j). In the case of the Gaussian kernel K G d , the contribution from a site (i , j ) to n p at site (i, j) depends on the distance between the two. The variance ofn d can be written as σ 2 d = σ 2 /N d , where σ 2 is the variance of n p and N d (=4π d 2 ) is the effective number of sites contributing ton d .

Emergence of pacemaker-like regions through diffusion
As can be seen from figure 3, large-amplitude oscillations occur in regions that have the highest local density of passive cells. In addition, the features observed on increasing the diffusion coefficient D resemble the patterns of the local passive cell density seen upon increasing the averaging size d. As f < f l c , increasing D is expected to ultimately suppress oscillation in This decreases the probability that a cell will have sufficiently high local passive cell density to generate oscillations. We define a 'pacemaker-like region' to be a group of adjacent cells with a local passive cell densityn d larger than the threshold f l c (C r ). The oscillatory activity arising in these regions may propagate in the form of waves to the rest of the system, so that they effectively act as pacemakers that are seen in systems having specialized coordination centers.
The number of cells comprising the pacemaker-like region varies from replica to replica. When f < f l c , we expect the size of the region to shrink with increasing coarse-graining length d, eventually disappearing at large d (as predicted by mean-field analysis). Thus, for a given replica, we can measure the largest value of the coarse-graining length, d * , for which a pacemaker-like region still exists. This is the smallest value of d for whichn d < f l c everywhere in the 2D system.
As mentioned above, the coarse-graining length d is related to the diffusion coefficient D, expected to be of the form d ∼ √ D. Therefore, in analogy with d * , we can define a value D * of the diffusion coefficient above which the number of oscillating cells in the system goes to zero. Figure 6 demonstrates that the relation between d * and D * (shown for different replicas and system sizes) follows d * 2 ≈ T D * , where T is the slope of the linear fit and defines a characteristic time of the order of the slow time scale in the system, ∼ 1/ . We note that as system size increases, both D * and d * increase. Thus, for systems just below the critical threshold for oscillatory activity, i.e. f f l c , oscillations are seen for a wider range of D in larger systems. This size dependence is investigated quantitatively in section 5.

Statistical description
In this section, we investigate the effect of coarse-graining length d on the presence and spatial extent of pacemaker-like regions in the 2D system. For analytical convenience, we use the 'square top hat' kernel, equation (6), for coarse-graining.

Distribution of passive cells
Consider a 2D lattice with N = L 2 excitable cells in which a total number of M = f N passive cells are randomly distributed. The probability that there are N p passive cells in a region containing N d = ζ N (0 < ζ < 1) excitable cells is given by the binomial distribution: When N → ∞, this reduces to a Poisson distribution with mean f N d For a square top hat kernel, the averaging occurs over N d = d 2 sites where d is the coarsegraining length and the local passive cell densityn d = N p /N d . For a Gaussian kernel, the corresponding coarse-grained region comprises N d = 4π d 2 sites.

Probability of the occurrence of a spontaneously oscillating cell in a neighborhood of N d cells
Given the probability distribution of passive cells, we now determine the probability that a given cell is capable of spontaneous oscillations when it is effectively coupled through diffusion to a neighborhood consisting of N d cells. This is obtained by considering the probability that the local passive cell densityn d = N p /N d f l c , which is given by the cumulative distribution corresponding to equation (7) where I x (a, b) is the regularized (incomplete) beta function [44] and P N N d depends on system size N . In the limit N → ∞, P N N d , tends to a regularized incomplete gamma function [45], .
The probabilities obtained for different values of N d from the above expression agrees well with the corresponding values numerically obtained from different replicas of the 2D system in figure 7. kernel. The broken curves indicate the corresponding probabilities obtained from the analytical expression of P N N d , equation (8). The solid curve shows P ∞ N d , the probability values obtained for an infinitely large system, equation (9).
In the limit of large system size N , P ∞ N d can be rewritten by expressing the incomplete beta function in equation (8) as a continuous fraction [45] and keeping only lower order terms in N −1 d and N −1 . In addition, neglecting ζ = N d /N we obtain for any fixed value of µ where λ = −µ − log(1 − µ). This expression can also be obtained from equation (9) by expanding the regularized incomplete gamma function and using the Stirling approximation. For small µ, we have λ µ 2 /2 > 0 and on introducing a reduced variable X = µ 2 f l c N d one obtains Note that in order to ensure that all replicas for a given value of f have exactly the same mean number of passive cells f =n p connected to an excitable cell, we have used a binomial distribution of n p for generating the different replicas in our numerical study. In comparison, the Poisson distribution used in our analysis leads to a mean fraction of passive cells attached to an excitable cell that varies from one replica to another for a given value of f . In the limit of large system size, the fluctuations in the number of passive cells become very small, so that both binomial and Poisson distributions lead to the same behavior. However, when N is not very large, these fluctuations are appreciable when using the Poisson distribution, which introduces biases in the results (see figures 7 and 9). In addition, for finite N , corrections of order N d /N are expected in the expression involving the reduced variable X , equation (11). These effects are responsible for deviations between the predictions of our analysis and the numerical results.

Probability of the occurrence of a pacemaker-like region in the 2D system
As mentioned earlier, figure 6 shows d * , the largest value of the coarse-graining length d for which a pacemaker-like region exists in the 2D system, as a function of the diffusion coefficient. We define the probability N N d of having a pacemaker-like region in the system as the probability of having at least one cell withn d f l c in a replica of size N = L 2 cells when the coarsegraining is done over a neighborhood of N d cells. Figure 8 shows the sigmoid form of this probability as a function of N d , computed using 400 replicas. To obtain an effective value of d * (or N * d ) for an ensemble of many replicas, we define the quantityd * (orÑ * d ) by the condition that Ñ N * d = 1/2. As mentioned earlier, larger systems have a higher probability of having a pacemaker-region for a fixed value of d (or N d ), which is explicitly shown in figure 8.
If we assume that there are η uncorrelated, i.e. effectively independent, subsets of size N d in the 2D system, we have One possible estimate of η is to assume that the independent subsets are obtained by tiling the system using non-overlapping blocks of size N d , i.e. η = N /N d . However, a single pacemaker region may be split among several neighboring tiles, which effectively leads to an underestimation of the probability N N d . Another possible estimate of η is to assume that all blocks of size N d are independent, so that η = N . However, the high degree of overlap between neighboring blocks introduces significant correlations between them, leading to an overestimation of N N d in equation (12). In general, we expect that the value of η will be related to N and N d by η = m N /N d . The existence of a number m independent of the system size can be qualitatively understood from To estimate the number of independent subsets containing N d points, we introduce κ, defined as the maximum possible correlation between two independent subsets. The correlation length l 0 is then defined by ρ(l 0 ) = κ, so that two subsets separated by a distance l > l 0 are independent. This implies the existence of N /l 2 0 independent blocks, suggesting in turn that m = N d /l 2 0 . Intuitively, one expects l 0 to be of the order of the width of the Gaussian kernel, d, yielding m ≈ 4π. We have numerically obtained the effective number m by least-square fitting of the numerical data shown in figure 8 with the theoretical expression of equation (12). This gives values of m lying in the range 9.5 m 16.5, which is in fact consistent with the heuristic estimate m = 4π.

System size dependence of the probability of the occurrence of a pacemaker-like region
We now obtain the cutoff value N * d for the coarse-graining size below which oscillations are present in the system by solving To solve this equation, we assume that the system is large enough, so that P N N * d can be replaced by P ∞ N * d , equation (10). This implies that the only system-size dependence of N N d is from the  (14). The symbols indicated in the legend correspond to different values of L and f , i.e. N and µ. The different curves collapse to the same form (within statistical accuracy) when N d is large, the regime for which the analysis described in the text is valid. exponent m N /N d . We thus obtain a linear relation between log(N ) and N * d ( figure 9) which effectively defines N * d . For large systems, the behavior of N * d is well described by the log(N ) term in equation (14) shown by a straight line in figure 9. It is worth noting that the log(N ) scaling appears very naturally in the general context of extreme value statistics [46]. The constant term in the rhs of equation (14) is responsible for the deviations from log(N ) scaling seen for N 10 4 .

Scaling in the transition to activity in spatially extended systems
The analysis presented in the previous section allows us to identify a characteristic coarsegraining size, suggesting that the number of pacemaker-like regions in a given system is determined by the ratio N d /N * d , where N * d is given by equation (14). In addition, figure 6 shows that diffusive coupling D results in coarse-graining of the passive cell density over a region of size N d ≈ 4π DT , with T being the typical oscillation period (section 4.2). Based on this, we expect that the fraction of oscillating cells will depend on the reduced variable 4π DT /N * d . Figure 10 shows the average (a) and the standard deviation (b) of the fraction of oscillating cells, f osc , averaged over many (∼100) replicas as a function of 4π DT /N * d , for several values of f (or equivalently, of µ) and L. The curves for the mean value of f osc seem to collapse to a common form at large N d where our analysis applies, equation (10), is not valid when N d is small; see in particular the discussion at the end of section 5.2), supporting the conclusion that the properties of the system indeed vary with the scaling parameter D/N * d . The deviations seen in the different curves can be explained as an outcome of the large standard deviation in f osc , as seen from figure 10(b). Remarkably, apart from the curves corresponding to the average values, we also observe relatively good collapse in the curves corresponding to the standard deviation of f osc , at least when D/N * d is not too small. This suggests that our analysis may apply not only to the averaged value but also to the second moment of the distribution of the fraction of oscillating cells, f osc .
As mentioned in the previous section, in the limit of very large system size, N * d is a function of the reduced variable µ 2 D/log(N ). This scaling implies that for a fixed value of the coupling between excitable cells, as the system size is increased and/or the mean value f of the passive cell distributions comes closer to the critical value f l c , there is a higher probability of observing sustained oscillatory activity in the system. Alternatively, if f is kept fixed at a value smaller than f l c , then the larger the system size, the larger the range of values of D over which oscillations can be observed.

Conclusion
In this paper we explore the role of quenched disorder on the transition between quiescent state and oscillatory dynamics in extended systems, composed of passive and excitable cells. Such a situation is particularly important in a biological context; the model we have considered is relevant to the understanding of spontaneous activity in the heart [30,31,33,34], or in the uterus [21]. While it is clear that in these systems, the intrinsic randomness must be taken into account, how the specific realization of disorder (replica) affects their behavior has rarely been studied quantitatively.
Our study centers around the key observation that different realizations of the disorder lead to large variability in the corresponding phase diagrams. While these diagrams give an indication of the rich variety of possible dynamical regimes seen for different parameter values as observed in a broad class of non-equilibrium systems [1,2,42], our study highlights the remarkable variability of the phase diagram as a function of the particular replica being realized. We emphasize that this variability is not merely an effect of multistability previously reported in the context of homogeneous excitable systems [35][36][37][38]. The results reported here demonstrate that the precise role of disorder in non-equilibrium transitions deserves more attention.
In this paper, we first analyzed the mean-field limit of the model system, which is equivalent to a single excitable cell coupled to a number of passive cells. We observe that the transition from quiescent to oscillatory activity is subcritical, i.e. oscillations appear at onset with a finite amplitude. We next investigated the spatially extended system at finite values of the coupling between excitable cells and characterize the coarse-graining effect of diffusion. We show that the dynamical behavior of the system is related to the local passive cell density obtained using a suitable coarse-graining length scale. The local density of passive cells can be different for each replica and the extreme values of this density distribution controls the observed variability in the phase diagram between the replicas. The strong effect of this quenched disorder on the dynamics is reminiscent of glassy systems. To characterize the properties averaged over many realizations of the disorder, we assume that the effect of diffusion of strength D is to coarsegrain the local passive cell distribution over N d sites. We indeed empirically establish that D ∝ N d . From this perspective, we analyze the transition between quiescent state and sustained oscillations. We obtain a scaling relation that describes the transition as a function of (i) system size, (ii) coupling between excitable cells and (iii) average passive cell density. One of the important implications of this relation is that the occurrence of oscillatory behavior depends on the logarithm of the system size N so that increasing N enhances the probability of observing oscillations.
As mentioned in the introduction, the model system that we have analyzed was motivated by the biological phenomenon of the onset of coherent oscillations in the pregnant uterus close to term. A possible limitation of our work is the use of a simplified model of an excitable cell. While the FHN model studied here qualitatively describes many natural phenomena reasonably well (e.g. see [47]), in the context of the uterine system, two important issues have to be addressed, viz., (i) whether the results concerning the transitions, discussed in [21] and in this work (see figure 4), continue to hold when using more realistic models of uterine myocytes and (ii) whether the transitions observed are in the biologically relevant range of coupling strengths. Our preliminary results obtained using a biologically detailed model of uterine myocytes developed recently [48], suggest that the phase diagram observed in this work is indeed seen in more realistic models. We also observe that the transition to global activity occurs when the coupling strength approaches the high values measured experimentally in the pregnant uterus close to term.
One of the implications of our work is that larger organs may show greater variability in their dynamical behavior for a given set of parameters describing the state of the system. In particular, they may be more likely to exhibit oscillations even prior to the transition point as a result of spatial fluctuations, potentially implying that mammals having bigger uteri will be at higher risk of having pre-term rhythmic activity. As it is not yet well-understood why in some cases periodic dynamical behavior is initiated in uterine tissue significantly earlier than normal, our study of the role of disorder in creating an effective pacemaker-like region giving rise to rhythmic activity in such systems may be of potential significance for possible clinical applications. Our work also connects the dynamical phenomena seen in such biological systems with the study of the role of disorder in phase transitions that occur in several physical systems, e.g. spin glasses.