Large variability in dynamical transitions in biological systems with quenched disorder

Coherent oscillatory activity can arise spontaneously as a result of increased coupling in a system of excitable and passive cells, each being quiescent in isolation. This can potentially explain the appearance of spontaneous rhythmic contractions in the pregnant uterus close to term. We investigate the transition to periodic activity using a model system comprising a lattice of excitable cells, each being coupled to a variable number of passive cells whose distribution defines a quenched realization (replica) of spatial disorder. Close to the transition between quiescent state and sustained oscillations in the system we observe large fluctuations between different replicas induced by variations in the local density of passive cells around an excitable cell. We demonstrate that the disorder-induced fluctuations can be described in terms of a simple scaling relation which involves the strength of coupling between excitable cells, the mean passive cell density, as well as the logarithm of the system size. Our results can be interpreted as suggesting that larger organs will have greater variability in the onset of persistent activity.


I. INTRODUCTION
Biological phenomena at different length scales often exhibit periodic activity whose frequency can span many time scales, and where the rhythmic behavior is crucial to the functioning of the relevant systems [1]. For instance, in the physiological context, oscillations are observed in variations of intracellular molecular concentrations [2,3], activation of pancreatic β-cell islets that control the release of insulin [4], changes in the activity levels of different brain areas [5], circadian rhythms that control the daily sleep-wake cycle [6] and patterns of locomotion [7]. In a system comprising many different oscillating elements, it is crucial for their activity to be synchronized in order for the generation of rhythmic activity at the macro-or systems-level. One of the most significant biological contexts in which such synchronization is observed is in the regular electromechanical contraction of the heart. A specialized group of pacemaker cells located in the sino-atrial node of the heart coordinates this periodic activity in a robust manner [8,9]. Disruption of the coherent collective behavior can result in arrhythmia that may be fatal if untreated in many cases [10]. However, such clock-like centralized coordinating agencies, although seen in several contexts [11], cannot be identified for many biological processes, which raises the possibility that coherence appears through self-organization among the many interacting oscillating elements.
A biological system whose function is critically dependent on the occurrence of rhythmic activity is the pregnant uterus, where coherent contractions are observed immediately preceding delivery. However, the available evidence does not suggest that this is centrally coordinated by a specialized cluster of pacemaker-like cells, as in the heart. For most of pregnancy, the uterus does not exhibit any sustained spontaneously generated excitation and can be well described as an excitable medium [6]. However, on approaching term, episodes of periodic activity are observed more and more frequently, and they increase in magnitude as well as in duration until system-wide powerful contractions eventually expel the fetus [12]. As the mechanical activity of the uterus is governed by electrical excitation of the tissue, the rhythmic contractions are related to the synchronization of oscillations in the transmembrane potentials of the cells in the myometrium that forms the bulk of uterine wall. In pathological cases (occurring in more than 10% of all pregnancies), rhythmic contractions may be initiated much earlier than normal leading to preterm birth [13] which account for a third of infant mortality in the USA [14]. It is not yet known why in many cases rhythmic activity is initiated significantly earlier than normal. As there is currently no effective treatment for events leading to preterm labor, understanding the dynamical processes underlying the onset of self-organized coherent rhythmic activity in excitable medium may have potential clinical benefits.
The present work is concerned with the emergence of synchronized oscillations in a system of coupled excitable and passive cells [15][16][17][18][19] that is motivated by the cellular organization of the mammalian uterus, where none of the constitutive cells can spontaneously oscillate in isolation [12,13,20]. Muscle tissue in the uterus, like that of many other biological organs, is heterogeneous in composition with electrically excitable myocytes (smooth muscle cells) forming the bulk but also containing small fraction of electrically passive cells like fibroblasts and telocytes [13,14,21]. The role of inter-cellular communication, including that between dissimilar cell types, in generating spontaneous periodic activity in tissue has recently come into focus [15,16]. This is of particular relevance for the uterus as the electrical coupling between arXiv:1212.3466v1 [cond-mat.dis-nn] 14 Dec 2012 cells through specialized proteins known as gap junctions is seen to increase remarkably during pregnancy [22][23][24].
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 [14], we consider the number of passive cells connected with an excitable cell to be randomly distributed about a given mean density. 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 [25], reminiscent of what has been observed in other systems of coupled electrically excitable and passive cells [16,26,27]. More generally, our results are relevant for a broad class of systems that comprise coupled heterogeneous elements, e.g., in the cardiac context [28,29]. In qualitative agreement with biological observations of uterine tissue, it is seen from the model that strong cellular coupling promotes synchronized oscillations [25].
While our earlier work provided a global description of the overall dynamics of the model system [25], using a more detailed approach to understand the genesis of the various dynamical regimes we now point out the key role played by disorder. In this paper, we show that the local fluctuation in passive cell density can result in 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, whereas individual members of the self-organized oscillatory groups referred to above are not inherently different from the other cells in the system. Thus, in the heart, pacemaker function is an intrinsic property of certain specialized cells, whereas, in the uterus it is an outcome of interactions between heterogeneous cell types.
A given distribution of passive cells attached to each myocyte effectively represents a particular realization of a quenched disordered system. The important role played by disorder in various phase transitions has been extensively investigated in physical systems [30,31]. Here, we investigate how the disorder described above influences the transition to sustained coherent activity in a vital biological organ. Specifically, we consider fluctuations in the macroscopic behavior of the system, viz., emergence of synchronized oscillations in different replicas, where a replica refers to a particular realization of the disorder. If excitable and passive cells are coupled with strength C r , the combination is capable of oscillating when 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 this transition between quiescent state and oscillatory activity occurs through a subcritical bifurcation [32,33], so 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 just the mean value of f over the entire system but also the spatially fluctuating values of local passive cell density that are important. The definition of the local density involves a coarse-graining length, d, which is expected to be related to the 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 II we describe the model system followed in Section III by a discussion of the properties of an excitable cell coupled to a variable number of passive cells. Section IV 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 V 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 VI, we conclude with a summary of our results.

II. 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 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 article we focus on phenomena that occur independent of the finer-scale details of specific models and therefore use the generic FitzHugh-Nagumo (FHN) model [9]. The functional form for the ionic current in the FHN system 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/ [9]. For the chosen set of parameter values, the dynamics of the FHN system involves either relaxation towards a fixed point for a subthreshold perturbation or a finite amplitude excursion before decaying to 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 [34]: so that 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. Similar to Γ e , Γ p 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 2-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 Eq. (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 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 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 do not result in qualitatively different phenomena [25]. 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.

III. GLOBAL OSCILLATIONS AND LARGE FLUCTUATIONS
A. 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 Eqs. (1)(2)(3)(4) [15,16,25]. 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 Eq. (4) can be replaced by a real number f , which corresponds to an exact representation of the system dynamics in the mean-field limit. Fig. 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 [15]. 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 Eqs. (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 Eq. (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 [ Fig. 1  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 [ Fig. 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 Fig. 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 Cr→∞ 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 [ Fig. 1 (a)], while the minimal value is f = 0.484.

B. Two-dimensional 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 system of coupled excitable and passive cells with periodic boundaries (Fig. 2). Fig. 3 shows snapshots of activity in two different replicas of such a system as the diffusion constant D is increased. 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 the medium does not show any oscillation; however, for finite values of D, it is possible to observe oscillatory behavior either in local clusters or globally as traveling waves (Fig. 3). We define a nondimensionalized 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 (Fig. 3). This is manifested as large fluctuations in the system behavior, measured by various order parameters defined later. For example, while 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 (Fig. 4). While for low values of C r , increasing D results in 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 were obtained by using suitable order parameters (defined in the next subsection), and choosing appropriate threshold values. As can be seen by comparing Fig. 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 signifi-cant variability from one replica to the other. By averaging over many such replicas, we can obtain a "mean" phase diagram. The diagram obtained here, Fig. 4 (c), for f = 0.5 is qualitatively similar to the one obtained in Ref. [25] for a higher value of f (= 0.7). The region corresponding to Coherence or "COH" phase also exists in the present case, appearing at large value of D which is outside the region of interest of the present paper.

C. 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 amplitude higher than a chosen threshold, to the total number of cells N . Note that the oscillating cells have approximately the same amplitude which is a consequence of the subcritical nature of the transition across f l c (Fig.1). This enables a clear distinction between oscillating and non-oscillating cells, so that the value of f osc does not depend sensitively on the choice of the threshold.
The different phases shown in Fig. 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 distribution of f osc over an ensemble of replicas (Fig. 5 a). At very low values of D, the distribution of f osc is peaked about a value close to f osc = f . Upon increasing D, the distribution broadens with a bias towards 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 Fig. 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 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.
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 V.

IV. 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.

A. Averaging procedure
We define local passive cell densityn d coarse-grained over a length scale d as the convolutionn 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., K d (i, j) = k d (i/d) × k d (j/d). 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 ST H d in section V 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 .

B. Emergence of pacemaker-like regions through diffusion
As can be seen from Fig. 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 Fig. 3 in agreement with the mean-field analysis. Indeed, we observe that at larger values of D, the number of oscillating cells is reduced on average. The coarse-graining procedure reflects this phenomenon as with increased kernel width d the variations inn d are reduced significantly. 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 coarsegraining 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. Fig. 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 V.

V. STATISTICAL DESCRIPTION
In this section, we investigate the effect of coarsegraining 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, Eq. (6), for coarse-graining. domly 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 coarse-graining 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.
B. 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 which is given by the cumulative distribution corresponding to Eq. (7): where I x (a, b) is the regularized (incomplete) beta function [35] 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 [36], .
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 Fig. 7.
In the limit of large system size N , P ∞ N d can be rewritten by expressing the incomplete beta function in Eq. (8) as a continuous fraction [36] 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 Eq. (9) by expanding the regularized Relation between the largest values of coarsegraining length d * and diffusion coefficient D * for which a given 2D system possesses regions with oscillatory activity. Each point corresponds to a different replica whereas different symbols represent different system sizes N = L 2 . Error bars centered around averages of replicas for the same system size express the standard deviation. The linear fit between d * 2 and D * , shown as a broken line, has a slope T 36.6 that lies in the middle of the range of observed oscillation periods (lying between 31 and 40) for different values of the diffusive coupling strength in a 2D system with L = 50.
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 Figs. 7 and 9). In addition, for finite N , corrections of order N d /N are expected in the expression involving the reduced variable X , Eq. (11). These effects are responsible for deviations between the predictions of our analysis and the numerical results. C. Probability of the occurrence of a pacemaker-like region in the 2D system As mentioned earlier, Fig. 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 coarse-graining is done over a neighborhood of N d cells. Fig. 8 shows the sigmoid form of this probability as a function of N d , computed by 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 higher probability of having a pacemaker-region for a fixed value of d (or N d ) which is explicitly shown in Fig. 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 neighbouring blocks introduces significant correlations between them, leading to an overestimation of Π N N d in Eq. (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 the following argument. Consider the covariance of the two variables, f d (i, j) = N p (i, j)/N d and f d (i , j ), i.e., of the average number of passive cells at the two sites (i, j) and (i , j ) separated by a distance l = (i − i ) 2 + (j − j ) 2 . It has the exact expression σρ(l)/N d where σ = (1 − ζ)f and ρ(l) = exp[−(l/2d) 2 ] is the correlation between the two sites. Thus, how fast the two values of f d decorrelate is simply given by ρ(l). 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 Fig. 8 with the theoretical expression of Eq. (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π. We now obtain the cutoff value N * d for the coarsegraining 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 , Eq. (10). This implies that the only system-size dependence of Π N N d is from the exponent mN/N d . We thus obtain a linear relation between log(N ) and N * d (Fig. 9): which effectively defines N * d . For large systems, the behavior of N * d is well described by the log(N ) term in Eq. (14) shown by a straight line in Fig. 9. It is worth noting that the log(N ) scaling appears very naturally in the general context of extreme value statistics [37]. The constant term in the R.H.S. of Eq. (14) is responsible for the deviations from log(N ) scaling seen for N < ∼ 10 4 .

VI. SCALING IN THE TRANSITION TO ACTIVITY IN SPATIALLY EXTENDED SYSTEMS
The analysis presented in the previous section allows us to identify a characteristic coarse-graining size, sug- gesting 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 Eq. (14). In addition, Fig. 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 (subsection IV B). Based on this, we expect that the fraction of oscillating cells will depend on the reduced variable 4πDT /N * d . Fig. 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, Eq. (10), is not valid when N d is small; see in particular the discussion at the end of subsection V B), 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 Fig. 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 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 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.

VII. CONCLUSION
In this paper we explore the role of heterogeneity in the spatial distribution of passive cells, which introduces quenched disorder in a system of coupled excitable and passive cells, on the transition between quiescent state and oscillatory activity. We begin our analysis with the mean-field limit of the spatially extended system, which is equivalent to a single excitable cell coupled to a number of passive cells. We observe here that the transition from quiescent to oscillatory activity is subcritical, i.e., oscillations appear at onset with a finite amplitude. We next investigate the system at finite values of the coupling between excitable cells and characterize the coarsegraining 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. We observe large variability in the dynamical behavior of different replicas, i.e., different realizations of the spatial distribution of passive cells. The strong effect of this quenched disorder on the dynamics is reminiscent of glassy systems. To characterize the properties aver-aged 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 has been motivated by the biological phenomenon of onset of coherent oscillations 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 occurring in several physical systems, including spin glasses.