The Effect of Gap Junctional Coupling on the Spatiotemporal Patterns of Ca2+ Signals and the Harmonization of Ca2+-Related Cellular Responses

Calcium ions (Ca2+) are important mediators of a great variety of cellular activities e.g. in response to an agonist activation of a receptor. The magnitude of a cellular response is often encoded by frequency modulation of Ca2+ oscillations and correlated with the stimulation intensity. The stimulation intensity highly depends on the sensitivity of a cell to a certain agonist. In some cases, it is essential that neighboring cells produce a similar and synchronized response to an agonist despite their different sensitivity. In order to decipher the presumed function of Ca2+ waves spreading among connecting cells, a mathematical model was developed. This model allows to numerically modifying the connectivity probability between neighboring cells, the permeability of gap junctions and the individual sensitivity of cells to an agonist. Here, we show numerically that strong gap junctional coupling between neighbors ensures an equilibrated response to agonist stimulation via formation of Ca2+ phase waves, i.e. a less sensitive neighbor will produce the same or similar Ca2+ signal as its highly sensitive neighbor. The most sensitive cells within an ensemble are the wave initiator cells. The Ca2+ wave in the cytoplasm is driven by a sensitization wave front in the endoplasmic reticulum. The wave velocity is proportional to the cellular sensitivity and to the strength of the coupling. The waves can form different patterns including circular rings and spirals. The observed pattern depends on the strength of noise, gap junctional permeability and the connectivity probability between neighboring cells. Our simulations reveal that one highly sensitive region gradually takes the lead within the entire noisy system by generating directed circular phase waves originating from this region.


Introduction
Calcium ions (Ca 2+ ) play a crucial role for almost every aspect in the biology of organisms. Cells possess sophisticated machinery to precisely regulate the free Ca 2+ concentrations in the cytoplasm (c cyt ), the endoplasmic reticulum (c ER ) and the mitochondria (c mito ). Maintaining the low concentrations of Ca 2+ in the cytoplasm against a 10,000-fold higher extracellular Ca 2+ concentration, i.e. the strong trans-membrane electrochemical gradient of Ca 2+ ions needed for proper cell signaling [1] requires energy. Upon agonist stimulation, cytoplasmic Ca 2+ levels are elevated from two sources: (i) Ca 2+ influx from the extracellular space across the plasma membrane and (ii) Ca 2+ release from stores, mostly the endoplasmic reticulum (ER). Different types of Ca 2+ channels are responsible for the Ca 2+ influx across the plasma membrane including: voltage-dependent Ca 2+ channels, receptor-operated Ca 2+ channels including transient receptor potential channels (TRP), store-operated Ca 2+ channels, etc. [2]. The release of Ca 2+ from the ER results from activation of either the ryanodine receptors (RyR) or the inositol 1,4,5-trisphosphate (InsP 3 ) receptors (InsP 3 R). Previously, it was assumed that RyR are of primary importance for Ca 2+ release in excitable cells, while InsP 3 R were presumed to govern Ca 2+ release in non-excitable cells. However, both InsP 3 R and RyR are expressed in excitable as well as in non-excitable cells [3,4], indicating cooperation between the two types of receptors in all cell lines. RyR have structural and functional similarities with InsP 3 R, but show no sensitivity to InsP 3 [5]. One of the roles of RyR is to amplify the InsP 3 -mediated release of Ca 2+ [6]. Ca 2+ signals are often organized in specific temporal patterns. The rhythmic changes in c cyt are called Ca 2+ oscillations. Several ligand/receptor interactions together with the involvement of components of the intracellular Ca 2+ -signaling toolkit induce Ca 2+ oscillations [7,8] and as a result Ca 2+ oscillations can act as integrators of different stimuli [9]. The stimuli intensity was often found to be proportional to the oscillation frequency, which in turn was proportional to the evoked down-steam cellular response, e.g. histamine-dependent fluid secretion of blowfly salivary gland [10] or glucose-dependent insulin secretion of pancreatic islets of Langerhans [11].
An identical genotype, differentiation state and moreover a stable environment are not sufficient to guarantee the same phenotype for neighboring cells within a tissue or organ. Indeed, single cell analysis of genetically identical cells grown in vitro revealed rather large cell-to-cell variability [12]. The transcription of DNA to mRNA followed by translation to protein occurs stochastically, as a consequence of the low copy number of DNA and mRNA molecules involved [13][14][15]. Therefore, each cell is expected to have a stochastic number of receptors for a certain ligand or agonist. This leads to an individually different sensitivity to agonists and individually different cellular response. For instance, Ca 2+ responses in individual mesothelial cells show a wide range of different oscillatory patterns within a genetically homogenous cell population [16]. Nevertheless, in many organs, the neighboring cells have to overcome their individually different sensitivity and produce a synchronized response for instance, smooth muscle cells, in order to generate the contractile waves of the uterus or the peristaltic movement in the gastrointestinal tract.
Gap junctions are integral membrane structures that enable the direct exchange of cytoplasmic constituents (ions and low molecular weight metabolites) between neighboring cells. The core proteins of these channels are the connexins [17]. Gap junctions are permeable to small molecules including both Ca 2+ and InsP 3 . Thus, gap junctions are involved in intercellular Ca 2+ signaling. Besides forming gap junctions between the same cell types, numerous gap junctions are also known to exist between different cell types. For example, gap junction channels ensure heterocellular Ca 2+ waves between glia and neurons [18]. Intercellular Ca 2+ waves spreading via gap junctions occur in rat liver epithelial cells upon mechanical stimulation [19]. Besides of gap junctional transport of Ca 2+ and/or InsP 3 , ATP may serve as a coupling messenger. ATP is thought to be released into the extracellular space and subsequently activating adjacent cells through purinergic receptors [20]. The ATP-mediated Ca 2+ spreading within cell populations is slower than direct gap junctional coupling, but allows a connection between cells not connected by gap junctions [19].
Two types of Ca 2+ waves can be distinguished associated with gap junctional coupling: i) Ca 2+ "diffusion" or trigger waves and ii) Ca 2+ phase waves. Ca 2+ trigger waves arise, when a single cell is stimulated in the network. In this case the gap junctional transport of InsP 3 originating from the stimulated cell drives the Ca 2+ wave. In this case, local increases in c cyt may be considered as an indicator of intercellular diffusion of InsP 3 molecules and not of the bulk movement of Ca 2+ ions [21]. Ca 2+ trigger waves are slow and due to the dilution of InsP 3 , the intensities of the Ca 2+ signal decrease in distant cells and finally Ca 2+ waves fade out. In the case of Ca 2+ phase waves, all cells are stimulated in the network, yet to different extents. Ca 2+ phase waves are generated by a small shift in the phase between individual cells oscillating with the same or nearly the same frequencies. Evidently, a coupling agent is required that synchronizes the ensemble of cells. A Ca 2+ phase wave differs from a trigger wave, since the spreading of Ca 2+ phase waves can be much faster than diffusion and moreover can travel long distances without annihilation [22].
Many different models have been built to understand the versatile patterns of travelling Ca 2+ waves [23,24]. However, these models do not take into account the heterogeneous sensitivity of cells to stimulation by agonists. This heterogeneity is assumed to substantially contribute to different individual temporal patterns observed in physiological settings. In this article, we address this issue. In order to decipher the presumed function of Ca 2+ waves spreading among neighboring cells, a general, not cell-type specific mathematical model was developed. The concept of this general approach is based on the hypothesis of Fewtrell [25]: "Since a single cell type may exhibit most, if not all, of the different types of oscillatory patterns, it seems unlikely that each cell type has developed its own specific mechanism for the generation of Ca 2+ oscillations. Instead, each cell may either contain a number of different mechanisms or a single, rather complex mechanism that is capable of generating the full range of oscillatory patterns". Our model, although containing several simplifications related to the machinery implicated in Ca 2+ oscillations, allows to numerically modifying (i) the connectivity probability between neighboring cells, (ii) the permeability of gap junctions and (iii) the sensitivity of a single cell to a particular agonist. The intercellular heterogeneity in agonist sensitivity can manifest in the onset of a Ca 2+ response, as alterations in the steady-state levels of InsP 3 and/or changes in the Ca 2+ handling. The main goal of this work is thus to provide a coherent model for intercellular Ca 2+ waves propagation and analyze the effects of gap junctions on networks of cells with inhomogeneous properties. This model highlights possible situations leading to the formation of typical signaling wave patterns (such as rings or spirals). In order to conduct our investigations, we analyze them in a first step within a deterministic framework, where the parameter values on the network are fixed. In such a manner, prominent parameters can be isolated. In a second step, noise is added to the system to stick to more realistic situations and its influence is investigated.

Cell network and random graph
We generalize the model of Pecze and Schwaller [16,26] by embedding it into a network of N = n m cells on a two-dimensional graph of size n × m. Each cell v ij , 1 i n, 1 j m is composed of the endoplasmic reticulum (ER) and the cytosol and whose dynamics is similar as the one depicted in [16]. We consider that cytosolic Ca 2+ ions can pass from one cell to another via gap junction [17], but that ER lumen from neighboring cells are not connected, thus not allowing direct transfer of ER luminal Ca 2+ . The parameter d denotes the strength of gap junctional coupling; the stronger the gap junctional coupling, the faster the diffusion of Ca 2+ ions between linked cells. For simplicity, only the diffusion of Ca 2+ ions is considered in the models reported in the Main Text, in line with the study of Hofer [27] or Harris and Timofeeva [28]. In the supplemental S1 Text, we extend these models by also considering InsP 3 diffusion through gap junction and show that InsP 3 can also function as the molecule involved in synchronization of Ca 2+ oscillation in cell types, where Ca 2+ spikes are connected to InsP 3 fluctuations [29].
Synchronization means here that two randomly selected neighboring cells tend to adjust the times at which they produce Ca 2+ peaks (not the amplitude of these peaks). This is known as phase-synchronization [30,31]. First we select two neighboring cells in the network, u~v and estimate the phases of the time series X u and X v (Ca 2+ concentrations within cytosol) using the interpolation technique described in [32]. Denote by τ 1 , τ 2 , . . . the times at which X u attain its maxima. Between two maxima, the phase increases by 2π and in between a linear interpolation is used, so that the phase of X u for τ n t < τ n+1 is defined by We consider now the differences of the phase of X u (t) and X v (t) and rescale these angles in the interval [0, 2π] In a phase-locked situation, this difference would always be the same over time and would result in a constant. If this constant is 0, it means that X u (t) and X v (t) are perfectly synchronized, that is, all their peaks occur simultaneously. When the phases vary over time, so do their differences. Therefore, in order to quantify synchronization, one has to consider the distribution over time of the phase differences. The more the distribution is concentrated around 0, the more synchronized X u (t) and X v (t) are. A uniform distribution corresponds to the null case of no synchronization.
To quantify the concentration of the phase differences distribution Tass et al. [30] and Cazelles and Stone [32] use the following synchronization index, based on the Shannon entropy, where the Shannon entropy S u, v is estimated by À S n b ðu;vÞ k¼1 p k ðu; vÞlogp k ðu; vÞ and the maximal Shannon entropy is S max,u,v = log(n b (u, v)). The number of class intervals is n b (u, v) and p k (u, v) denotes the relative frequency that Δϕ u,v lies in the k th interval. The index Q u,v lies between 0 (no phase-synchronization) and 1 (perfect synchronization). Finally, we define a synchronization measure of the network by taking the mean of Q u,v over all neighboring cells u~v where |ε| is the number of edges, i.e. pairs of neighbors in the network. When m sync = 1, the synchronization between all cells is perfect and in particular no wave nor any pattern would appear. For medium ranges of m sync , synchronization between neighboring cell is efficient enough to enable waves to emerge. When m sync is close to zero, all cells tend to behave independently from each other. Our synchronization measure m sync indicates the degree of phase synchrony focusing on time periods at which peaks are produced within the whole network.
Amplitudes of the peaks are not taken into account for the parameter m sync . Two graph topologies were used: homogeneous, fully dense graph and heterogeneous sparse graph topologies (Fig 1).
I) The homogeneous, fully dense graph topology consists of a grid graph to which diagonals were added (G 0 , Fig 1A). In this setting, all neighboring cells are connected to each other. The set of nodes is V 0 = {v ij , 1 i n, 1 j m} and let ε 0 denote the set of links. Let v i 1 j 1 and v i 2 j 2 be two members of set V 0 . These nodes are connected in G 0 , if and only if they are neighbors, that is, II) Heterogeneous sparse graph (G) is a random sparse version generated from the homogeneous fully dense graph (G 0 ). In this situation, some links and nodes were randomly removed, exemplified in Fig 1B. Decreased gap junction coupling can be achieved by decreasing the permeability or numbers of the specific channels involved [21,33]. The latter situation can be modeled by link removal. Changes in the densities of gap junctions (strong decrease and recovery) was observed after partial hepatectomy [34]. Random sparse graphs were obtained by the following procedure. The link-removal probability p(v) is built with a colored noise in space, in order to create spatial correlations among the network. Let F be a Gaussian random field on the graph with standard deviation set to one and with a correlation function of the Gaussian form, that is is a correlation parameter. When randomized, the probabilities p(v) are as follows pðvÞ ¼ p maxðFÞ À FðvÞ maxðFÞ À minðFÞ ð7Þ for some p 2 [0, 1] in a way that in regions with nodes having low values of F, the probability to remove links is high, creating sparse areas, whereas regions with large values of F tend to be more strongly linked. Choosing a small quantile q of F, it is possible to create holes in the network by excluding nodes v from V 0 whenever F(v) < q. The resulting graph is denoted by G ¼ ðV; EÞ, ( Fig 1B).

Dynamics of Ca 2+ concentrations on the network
Denote by (X v (t), Y v (t)) the free Ca 2+ concentrations in the cytosol (c cyt ) and in the ER (c ER ) of the cell v at time t. The corresponding dynamical system on the graph G is driven by the following system of differential equations: for all v 2 V with initial conditions X v (0) = c cyt,ini = 110nM and Y v (0) = c ER,ini = 260μM according to [16]. At the boundaries or in the case of holes, fluxes follow only existing directions, and there are no fluxes via the missing links. The function J INF (v,t) comprises the following movements of Ca 2+ ions that increase the free Ca 2+ concentration in the cytoplasmic compartment: (i) Ca 2+ ions entering the cytoplasm from the extracellular space or from organelles (mitochondria, Golgi apparatus, lysosomes) and (ii) Ca 2+ ions released from cytoplasmic Ca 2+ -binding proteins, which are considered as Ca 2+ buffers in a broad sense. In a stricter sense, Ca 2+ buffers comprise specific subsets of mobile Ca 2+ -binding proteins such as calretinin, calbindin D-28k and parvalbumin. The particular case of calretinin is presented in the S2 Text). The diffusion of Ca 2+ through gap junctions from a cell's cytosol to the one of an adjacent cell is given where the constant γ is the ratio between the changes in X v and Y v caused by the same amount of Ca 2+ ions transported through the ER membrane. This value is derived from the difference in the effective volume of the ER lumen and the cytosolic volume and from the different fractions of free and protein-bound Ca 2+ in these compartments. The function J EFF represents a combination of Ca 2+ fluxes: (i) from the cytosol to the extracellular space, (ii) from cytosol to organelles and (iii) Ca 2+ ions temporarily sequestered by cytosolic Ca 2+ buffers. The common features of these Ca 2+ -ion movements are their dependence on c cyt . The higher c cyt is, the more Ca 2+ ions are removed by mitochondria or plasma membrane extrusion systems, or are bound to cytosolic Ca 2+ buffers. We simulated J EFF by a linear equation, in line with the work of Fink et al [35] and based on the experimental results of Herrington et al. [36]. The different dissociation constant (K d ) values of individual components (plasma membrane Ca 2+ ATPase's, exchangers, mitochondria) ensure that the extrusion fluxes will never reach their saturation points in the range of biologically relevant values of c cyt .
where the indicator function 1 {x>0} returns 1 if x > 0 and 0 otherwise. The function J SERCA denotes the Ca 2+ flux from the cytosol to the ER. For simplicity, we assume that it depends only on c cyt , in line with our previous study [16]. Here, we also assumed that the different SERCA variants with different K d values [37] ensure that this flux will never reach its saturation point in the range of biologically relevant values of c cyt . It is given by Modeling of J SERCA and J EFF fluxes with saturating Hill equations, did not modify qualitatively the behavior of our model. The simulations are presented in S3 Text. The parameter J ERLEAK represents a small flux of Ca 2+ ions diffusing from the ER to the cytosol. The origin of this flux is unknown as well as its main characteristics [38]. For simplicity, we considered this flux to be constant, J EREFF describes the flux of Ca 2+ passing from the ER to the cytosol through InsP 3 R. In our model, InsP 3 R are influenced both by c cyt and by c ER , however without an allosteric regulation between the two. We consider J EREFF as the sum of two individual functions. The first one depends on c cyt and has a bell-shaped form when considering concentrations in logarithmic scale [39], where M is the concentration mean on a logarithmic scale, R i,max the maximum and where the variance σ 2 is a positive constant. The second function expresses the experimental fact that an increase in [InsP 3 ] causes significant Ca 2+ release from the ER even if c cyt = 0 [40]. The shape of J EREFF depends on [InsP 3 ] and is encapsulated in the functions M, R i,max and R i2 . Indeed, [InsP 3 ] modifies the sensitivity of InsP 3 R to changes in c cyt and in c ER . Elevating [InsP 3 ] mainly changes the mean (M on a logarithmic scale) and the maximum (R i,max ) of the bell-shaped curve of c cyt dependence [41]. In our model, we did not consider that the open probability curve at high [InsP 3 ] is not precisely bell-shaped [41]. For simplicity, we didn't consider that the Ca 2+ flux from the ER, either through leak channels or InsP 3 R depends on the driving force originating from an electrochemical gradient across the ER membrane. We used stationary open probabilities, not considering the binding/unbinding kinetics of Ca 2+ to activating or inhibitory sites of InsP 3 R. Taking the abovementioned factors into account did not modify critically the behavior of the system. The simulations are presented in S3 Text. Based on the experimental data presented in [42,43], elevating [InsP 3 ] also has an effect on the loading of the ER. Increased [InsP 3 ] reduces the amount of Ca 2+ ions stored in the ER. We simulated this effect by changing the parameter R i2 . The functions M, R i,max and R i,2 have similar forms, given by We assume that within cell populations, each cell has a different sensitivity to agonists. A stimulation induces two processes that are important for Ca 2+ oscillations: (i) it increases the levels of InsP 3 by G-protein-regulated phospholipase C and (ii) it increases the Ca 2+ flux characterized by J INF , mainly due to the opening of plasma membrane Ca 2+ channels. We assume that both processes are positively related to the stimulus intensity, but two cells with the same [InsP 3 ] can have different J INF values. Hence, we used two input parameters: and In accordance with experiments [44], we assumed that stronger activation results in smaller t 1 value. Thus, the values of i JINF,max (v) and i IP3 , max (v) are positively correlated among cells, but the values of i JINF,max (v) and t 1 (v) are negatively correlated, as exemplified in Fig 1D and 1E. In our model, we did not consider that J INF is partially sensitive to changes in c cyt . All parameter values are presented in Table 1. All these values were set to reproduce as closely as possible the Ca 2+ concentration changes in cytosol and ER measured experimentally by fluorescent Ca 2+ indicators.

Particular models
We present and analyze four particular types of models and networks, each type allowing to highlight different phenomena and to isolate the key features leading to the various patterns arising. Model G 0 : Individual cell. Oscillations in a single cell are investigated, shedding light on the parameter values leading to the oscillatory response to an input signal. The parameters are fixed according to Table 1 unless specified.
Model G D : Noise-free graph with distinct regions. We consider the graph G 0 . All parameters are deterministic and fixed for all nodes v 2 V 0 . Onto this network, we introduce 1, 2 or 3 regions with different sensitivities and inputs. More formally, we define s 2 N subgraphs ðG j Þ 1;...;s of sizesñ j Âm j , within the graph G 0 . The subgraphs are disjoint (i.e. with no common nodes). Thus i JINF,max (v) takes the form: with i JINF,max,j > 0 for all j = 0, . . ., s. In order to maintain the correlations between  (v), i IP3,max (v) and t 1 (v), we randomize them with a white noise (identical and independent randomizations), so that  (v) and t 1 (v), the parameters t 1 (v) and i JINF,max (v) were randomized as follows. Let (ψ v ) be a centered white noise with standard deviation s t 1 and (φ v ) a centered white noise with standard deviation s i JINF;max . Let denote the order statistics of the noises (ε v ), (φ v ) and (ψ v ), such that ε [1] is the minimum of (ε v ) and ψ [1] the maximum of (ψ v ). For a given The resulting noises (ψ v ) and (φ v ) are white, but (φ v ) is negatively correlated with (ε v ) and (ψ v ) is positively correlated with (ε v ). Finally, set All analyses were performed in MATLAB version 2013b (Mathworks Inc., Natick, MA, USA) and the code for the scripts is available in the Supplementary Material (S1 Code).

Oscillatory regimes depending on the sensitivity of individual cell in model G O
An isolated cell (model G 0 ) oscillates in response to a given input encapsulated into the functions J INF and InsP 3 for well-chosen parameters (see [16]). Unless specified, all parameters are set according to Table 1 and the ranges of parameters leading to oscillations are depicted in Fig  2. The period of oscillations (inverse of the frequency) is determined by i JINF,max and i IP3,max ; three zones can be distinguished. In zone I, the Ca 2+ concentration oscillates with a constant frequency as exemplified in Fig 2B. In zone II (the boundary of zone I) any small oscillation rapidly vanishes after the initiating stimulation signal (Fig 2C). Finally, in zone III, no oscillations occur and Ca 2+ concentrations saturate (Fig 2D). The model for each cell has a Hopf bifurcation that separates region III from I and a homoclinic one that separates region I from region II. The system oscillates in a certain range of i JINF,max and i IP3,max values. We thus differentiate two types of cells: self-oscillating ones and non-self-oscillating ones.
Effect of the gap junctional coupling (d) on Ca 2+ waves in the noise-free models G D Consider a noise-free graph (model G D ) with gap junctional coupling d > 0 and one particular region at the center of the homogenous, fully dense graph G 0 , where the sensitivity is increased (i JINF,max,1 >i JINF,max,0 , i IP3,max,1 > i IP3,max,0 and t 1,1 < t 1,0 ). In this setting, a Ca 2+ wave emerges from the center of the graph, indicating that the most sensitive cells are the wave initiators (Fig 3 and S1 Movie). The wave front is initiated in the ER as illustrated in Fig 4A. Hence, intercellular Ca 2+ waves are driven by ER Ca 2+ release as reported in [29]. This is also in accordance with the findings in [16], where in single cells, oscillations are initiated from within the ER. With a strictly positive gap junctional coupling, single oscillations thus propagate through all connected cells with a velocity controlled by the strength of the gap junctional coupling d, as illustrated in Fig 4C. In Fig 4B, we report that i JINF,max,0 and i IP3,max,0 additionally contribute to the wave speed. Large values of the parameter d diminish the period of the wave and thus increase the oscillation frequency. Note however that if parameters i JINF,max,0 and i IP3,max,0 are too large, the network homogenizes rapidly, so that c cyt reaches saturating values in every cell, thus dampening any oscillatory behavior.
Indeed, parameters i JINF,max,0 and i IP3,max,0 negatively correlate with the period of oscillation in all cells of the graph as illustrated in Fig 4B. In this figure, the period of oscillation has been calculated for various parameter values of i JINF,max,0 , i IP3,max,0 and a moderate d. Both parameters participate to the homogenization of the patterning through the graph. If both parameters i JINF,max,0 and i IP3,max,0 are sufficiently large, no oscillation is observed. All cells fill up with Ca 2+ and remain on a stationary state. In Fig 4, which is seen as an added percentage of sensitivity in the sub-graphG. This parameter k has no remarkable global effect except that for sufficiently large values of i JINF,max,0 , it slightly negatively contributes to the wave speed (Fig 4D). In the S2 Text, an extension of this model, which explicitly takes into account the presence of calretinin, shows that this specific buffer has an effects on wave propagation: an increase in calretinin tends to increase the oscillation frequency and the wave speed, as illustrated in Fig. A in S2 Text. Interestingly, the synchronization of neighboring cells occurs even if the network contains non-self-oscillating cells. To verify this, we add to the previous network a single region Gap Junctions and Ca 2+ Signals Z of cells with parameters (i JINF,max,2 , i IP3,max,2 ) outside of the range that would allow them to oscillate independently (zone III in Fig 2D). Any oscillations in those cells would rapidly die out (Fig 2C). But since they are coupled to neighboring oscillating cells, they synchronize and oscillate in response to the behavior of their neighbors (Fig 5 and S2 Movie). Ca 2+ ions transported to non-self-oscillating cells thus evoke changes in their oscillations-properties. This phenomenon is accelerated by higher gap junctional coupling values. In such a configuration, we observe that the wave initiators are shifted from the center to cells adjacent to the region Z.   Effect of randomization and holes in the random model G R When considering noise in the system, the same local effects depicted previously (Figs 4 and 5) occur. Neighboring cells, although having different sensitivities, try to synchronize their operation. This is illustrated in Fig 6 and S3-S5 Movies, where we used the random model G R described in the previous section for a system with strong noise and with p = 0, i.e. no links were removed (graph G 0 ). The extreme case of no coupling (d = 0) leads to independent Ca 2+ oscillations in each cell (Fig 6A and S3 Movie). For low and moderate couplings, small bursting phenomena are visible before synchronization of two neighboring cells (Fig 6B and  S4 Movie). By "burst" we mean that a cell changes its phase during the development of a Ca 2+ spike resulting in a prolonged irregular Ca 2+ spike. Moreover, no coherent behavior or typical wave patterns emerge from such situations. Incorporating calretinin into the model has the very interesting effect of promoting the formation of coherent wave patterns. This phenomenon is exemplified in S18 Movie, where the same framework as in Fig 6B is used, with the addition of calretinin (see S2 Text). More generally in Fig. B in S2 Text, we illustrate that a deficiency in coupling can be generally compensated by sufficiently high levels of calretinin from the point of view of synchronization. Notice that strong coupling leads to well-synchronized oscillations, even in the presence of strong noise (Fig 6C and S5 Movie). We observed that when phase synchronization occurs, the initial differences in the amplitudes of individual cells decrease, i.e. the harmonization of the network is also visible as a harmonization of the amplitudes. Highly sensitive cells still initiate Ca 2+ waves travelling through the network. In such noisy systems, waves can arise from different random places. When they collide, they aggregate or split depending on their velocity and on local properties of the graph.
In the S1 Text, we also consider the effect of InsP 3 coupling within this random model. In these particular settings, Fig. A in S1 Text illustrates how InsP 3 and Ca 2+ oscillations synchronize in an arbitrary cell. S14-S16 Movies show the evolution of Ca 2+ concentrations when (i) Ca 2+ coupling is active, but there is no InsP 3 coupling (S14 Movie), (ii) InsP 3 coupling is active but there is no Ca 2+ coupling (S15 Movie) and (iii) Both couplings are active (S16 Movie). There are two important phases in order to compare the three situations. In the beginning, around t 1 (v), InsP 3 coupling is very efficiently involved in the synchronization process, but afterwards acts poorly as synchronizing agent (see Fig. B in S1 Text). The reverse is true for Ca 2+ coupling (Fig. B in S1 Text). All results are reported in Table B in S1 Text.
Next we examined the effect of link removal on Ca 2+ wave propagation i.e. the heterogeneous sparse graph was used ðGÞ. Besides the gap junctional coupling (d), the density of the connection p(v) influences the wave propagation. Remember, the parameter p(v) indicates, whether two neighboring cells are connected or not, while the parameter d represents the strength of the coupling if such a connection exists. Hence it is not surprising that increasing the link removal probability p(v) has a similar effect on the wave propagation as decreasing the strength of the gap junctional coupling d. Highly connected cells (regions) enable waves to spread easily, while poorly linked regions will act as a barrier deviating the wave front to another direction. In the extreme case, holes were added to the network (Fig 7 and S6-S8 Movies).

Spiral waves
The most interesting feature occurring in random systems is the waves-patterning transitions. As explained above, waves initially emerge from highly sensitive cells and spread radially away from them. Such travelling waves commonly appear in many oscillating systems [45]. In inhomogeneous systems built with our model, such patterns do not necessarily stabilize with time. Circle centers move, some waves die out as other gain in strength, collide and even spiraling phenomena are observed (S8 Movie). In this section the different frameworks resulting to spiraling phenomena are explored.
In our settings, surprisingly, even homogenous noise-free networks can exhibit transitions from circles to spiral waves. The cause of this has to be looked for in the inhomogeneity of cell sensitivities. Indeed, if we consider the same framework as in [46], by considering model G D , with three additional small regions where cells are more sensitive, spiral waves appear (see Fig  8A and S9 Movie). These spirals are of transitory nature and the system finally ends up with stable circular waves emanating from the three sensitive zones.
Notice that defining two particular zones is sufficient to create spirals as demonstrated in Fig 8B and S2 Movie. In this figure, the central region is very sensitive and a second area (upper-left) is low sensitive (with parameters in zone III of Fig 2A). After a sufficiently long time, a spiral develops around this second area and moves along the network. This sheds light on the effect of gap junctional coupling. Here its value (d = 0.0045) is moderate. An identical framework, but with a high gap junctional coupling, would prevent spiraling behaviors. Considering two highly sensitive areas leads to the same conclusion, as illustrated in S10 Movie.
Beyond collisions of different travelling circular waves and topological considerations, gap junctional coupling and inhomogeneous activation (as defined by e.g. noise strength) play key roles in spiral generation. Independently of its strength, noise acts as a catalyst for producing spirals. For a low level of noise, moderate gap junctional coupling enables the appearance of transitions between circles and spirals. In highly disordered systems, the level of gap junctional coupling needed to smooth the spirals and obtain stable circles increases. Indeed, a high gap junctional coupling allows waves to spread rapidly among the cells and thus enables a fast synchronization of their behavior ( Fig 6C). As such, any inhomogeneity that can be generated from different traveling waves breaks up.
To explore the particular effect of inhomogeneous activation, the noise is set to a low value. In Fig 9A (and S11 Movie) the random model G R is used with a relatively high sensitivity. Due to the noise, the system is inhomogeneous enough to produce spirals, as locally explained in Fig 8A-8B. Turning to model G R,C , by adding a small central zone with higher sensitivity regulates the behavior of the system and concentric circles develop (Fig 9B and S12 Movie). In this case the sensitive center enables the homogenization of the whole network. This occurs in a similar way as what is shown in the noise-free framework of Fig 3. Increasing simultaneously the two main frequency-determining parameters, m i IP3;max and m i JINF;max , results in similar concentric circles, but with higher frequencies (see S19 Movie, where the same parameters as in Fig 9C are used with m i JINF;max;0 ¼ 0:9 and m i JINF;max;1 ¼ 1). However, very high values for one parameter (e.g. m i IP3;max as in Fig 9C and S13 Movie) shifts the model to non-organized wave propagation and creates spirals again; the sensitive center is no more sufficient to homogenize the whole network. Our explanation for this phenomenon is the following: increasing one parameter excessively increases the number of non-self-oscillating cells (Fig 2A Zone III), but when increasing simultaneously both parameters, most cells remain in the oscillating Zone I and the network is able to produce rhythmic concentric circles The effect of specific buffers (calretinin) on m i IP3;max À overstimulated systems is investigated in S2 Text. Our general finding is that Ca 2+ buffers promote coherent behavior in these situations.

Discussion
Over the last 25 years, many different models have been developed to describe Ca 2+ oscillations in cells [23,24,47]; new methods including high-resolution spatiotemporal recordings of Ca 2+ signals enabled the reconsideration and further development of the different models. As an example, the simultaneous monitoring of c cyt and InsP 3 production has revealed that Ca 2+ oscillations are not the direct consequence of fluctuations in [InsP 3 ] [29]. Simultaneous monitoring of c cyt and c ER showed a time shift between the maximum of c cyt and the minimum in c ER during a Ca 2+ spike [48][49][50][51], findings that have allowed for a better understanding of the mechanisms implicated in oscillations. The continuous loading of the ER with Ca 2+ during the interspike periods followed by a rapid release during the Ca 2+ spike in c cyt results in sawtooth- like Ca 2+ oscillations in c ER [16]. This indicates that the loading state of the ER is an essential parameter to understand Ca 2+ oscillations in the cytosol. A next Ca 2+ spike can be generated only, if c ER reaches a certain threshold value [52]. This threshold is determined by the prevailing InsP 3 concentration [42]. The replenishment of the ER Ca 2+ store is modulated by a constant Ca 2+ influx across the plasma membrane. An increase in the Ca 2+ influx rate leads to a higher frequency of Ca 2+ oscillations, while decreasing the Ca 2+ influx reduces the frequency [16,26,53]. In some conditions, mitochondrial Ca 2+ transport (uptake and release) was found to substitute for the plasmalemmal Ca 2+ exchange function, thus rendering the oscillations independent of extracellular Ca 2+ [26]. The magnitude of the Ca 2+ transport into mitochondria was also found to influence the Ca 2+ oscillation frequencies [54]. Often it is the refilling of the ER that sets the oscillation period (frequency), not the InsP 3 R dynamics. The model proposed in our previous works [16,26] has demonstrated to be coherent with the above-mentioned experimental findings at the single-cell level (zero dimension). In this study, we showed that this single-cell model in a 2D framework is a useful tool for the prediction and understanding of several phenomena in naturally-occurring multicellular noisy systems. Among existing models for intercellular Ca 2+ waves, the key limitations are the size of the system (two [27] or three cells [55]), restriction of cell activation to a single cell [56] or the fact that cells with altered sensitivities are restricted to few areas in an otherwise noise-free system [46]. To the authors' best knowledge, the model that we propose in the present study is the first to address simultaneously those limitations.
Our model revealed the ER to be the initiation site of the Ca 2+ phase wave front (Figs 3A and 4A). Thus, intercellular Ca 2+ phase waves are driven by an initiative wave front starting in the ER; even if the neighboring cells might communicate with each other by exchanging their cytosolic but not luminal Ca 2+ ions via gap junctions. This model prediction has already been proven experimentally. Keller and coworkers have found in guinea heart myocytes that Ca 2+ waves in c cyt are driven by "sensitization" wave fronts in c ER [57]. Differently from Ca 2+ phase waves, the initial wave front of Ca 2+ trigger waves starts from within the cytosol. Another difference between the two types of waves is that two Ca 2+ phase waves annihilate each other when they collide, while two Ca 2+ trigger waves add up to generate a new wave of greater amplitude.
It has been known for a long time that gap junctions are permeable both to InsP 3 and Ca 2+ ions [58]. However, because Ca 2+ but not InsP 3 is strongly buffered by cytoplasmic proteins and/or unidentified immobile buffers, Ca 2+ movement within a cell is very restricted and slower than that of InsP 3 [59]. Thus InsP 3 is more likely to diffuse to greater distances and subsequently to mediate intercellular Ca 2+ waves [21]. Although this concept is quite attractive, one has to take into account that (i) Albritton et al. [59] measured the Ca 2+ diffusion in cell extracts, where the Ca 2+ pumps and InsP 3 -metabolizing enzymes were blocked, i.e. not in physiological conditions. (ii) To generate a long distance Ca 2+ wave, Ca 2+ ions do not need to travel for long distances; it is enough to diffuse to the next InsP 3 R or RyR. More precisely, since these receptors are organized in clusters, Ca 2+ ions only have to diffuse from one cluster to the neighboring cluster. Since the ER consists of a network in the cytoplasm filling almost the entire cell, the cluster-to-cluster distances between adjacent cells should be of similar magnitude than that of intracellular cluster-to-cluster distances, i.e. based on the images presented in Chalmers et al. [60], approximately 1 μm. (iii) Ca 2+ buffer proteins not only take up Ca 2+ ions, but to the same extent, also release Ca 2+ ions. Depending on their Ca 2+ -binding parameters (fast or slow kinetics) they can promote or inhibit intracellular Ca 2+ wave formation [61]. Most probably, they have similar effects on intercellular wave formation, yet no experimental data are available. As a preliminary result in the S2 Text, we have already simulated the effect of a specific Ca 2+ buffer, calretinin, on wave propagation. Our model predicts that calretinin promotes intercellular synchronization. (iv) As presented in the introduction, there are two types of Ca 2+ waves: Ca 2+ trigger waves and Ca 2+ phase waves. If we complete our model with the gap junctional transport of InsP 3, this results in Ca 2+ trigger waves, since only few cells get activated before activation of the large majority of the other cells. In this case the gap junctional transport of InsP 3 is the main agent harmonizing the initial Ca 2+ signal, in agreement with previous experimental results [21]. However, the existence of individual Ca 2+ oscillatory machinery and the individual stimulation of each cell, as in our model, allows the generation of phase waves in which the emergence of the waveform is due to Ca 2+ release from adjacent oscillating cells slightly differing in their oscillation phase.
In summary, our model shows that both Ca 2+ ions and InsP 3 are likely the synchronizing agents, possibly in a synergistic way. InsP 3 is more involved in the formation of the initial Ca 2+ trigger waves, while Ca 2+ ions serving as a coupling agent are more implicated in the later ones indicating that Ca 2+ phase waves are dominant at the later stages of a model experiment (See S1 Text). Of note, in other studies, the authors have concluded that only InsP 3 may serve as the coupling agent [55,62], reporting that the pertaining actual [InsP 3 ] is the only frequency-determining factor. In our model, we identified two factors determining the Ca 2+ oscillation frequencies, i.e. i JINF,max and i IP3,max .
Stimulation of primary mesothelial cells induces Ca 2+ responses showing a wide range of different oscillatory patterns within a single, supposedly homogenous cell population. Since a single cell type may exhibit most, if not all, of the different types of oscillatory patterns, most likely each cell contains all components of the Ca 2+ signaling toolkit (possibly to different extents) required to generate the full range of oscillatory patterns and spreading Ca 2+ waves [25]. Our model indicates that different spatiotemporal patterns of intercellular Ca 2+ signals are mostly the consequence of different strengths of coupling via gap junctions. Non-synchronous oscillations in individual cells is favored when gap junctional coupling is weak as was observed in cultured primary mesothelial cells [16]. Moderate levels of gap junctional coupling led to temporally and spatially restricted Ca 2+ bursts (Fig 6). Nevertheless, other mechanisms can also be involved in the formation of Ca 2+ bursts [63][64][65]. This type of oscillations and waves was found in pancreatic beta cells upon glucose stimulation [66]. Strong gap junctional coupling resulted in Ca 2+ waves, even if individual cells had different sensitivities with respect to stimulation. Synchronous smooth muscle cell-mediated contractions of the uterus driven by Ca 2+ waves [67,68] are a typical example of strong coupling. The influence of gap junctional coupling on the generation and spreading of intercellular Ca 2+ waves was experimentally revealed by analysis of GT-1 cells, a cell line derived from immortalized Luteinizing Hormone-Releasing Hormone neurons. GT-1 cells were then further subcloned to result in lines GT1-1, GT1-3 and GT1-7 cell lines differing in expression levels of connexin 26 (Cx26) [69,70]. Low Cx26-expressing GT1-7 cells mostly displayed frequent spontaneous asynchronous Ca 2+ oscillations, while high Cx26-expressing GT1-1 cells showed spontaneous intercellular Ca 2+ waves, completely in line with our model. The occurrence of Ca 2+ waves strongly depends on cell-cell contact probability and moreover the strength of gap junctional coupling. Our model also predicts that in a noisy system, in which each cell has an individual sensitivity to stimulation, the most sensitive cells act as the wave initiator cells.
Physiological cellular responses to evoked Ca 2+ signaling are cell-type dependent: Ca 2+ signals elicit contraction in muscle cells [71], neurotransmitter release in neurons [72] and e.g. insulin secretion in pancreatic beta cells [73]. Ca 2+ signals in immune cells participate in the regulation of cell differentiation, gene transcription and effector functions [74]. They are also involved in the regulation of cell proliferation of cancer cells [75]. Intracellular Ca 2+ signaling, frequently in the form of Ca 2+ oscillations, activates specific enzymes and transcription factors in a cell, which are often involved in cell proliferation [76]. Ca 2+ signals are usually not restricted to individual cancer cells, but are propagated to neighboring cells in the form of intercellular Ca 2+ waves, usually Ca2+ trigger waves are observed [77]. Assuming that intraand intercellular Ca 2+ signaling is the major way by which cells encode and transmit information, it is likely that during the passageway of a Ca 2+ wave, several Ca 2+ -dependent targets in affected cells would be activated and/or deactivated. Our model predicts that Ca 2+ waves play an important role in the harmonization of evoked responses i.e. Ca 2+ wave initiator cells are capable of activating neighboring cells, even when those cells would not oscillate by themselves; however they start to oscillate and hence support waves when coupled to other cells. Thus, a highly sensitive cell may trigger cell division even if the neighboring cells are not sensitive enough to mitogenic stimuli. The exploration of the details on Ca 2+ wave-dependent harmonization of Ca 2+ -related cellular responses remains an interesting topic to be investigated experimentally in the future.
Our model is capable of producing both circular waves and spirals within a specific range of parameters. Systems producing spirals, such as the famous Fitzhugh-Nagumo equations have been mathematically analyzed [45,[78][79][80] and several ways to create spirals in oscillating systems was proposed in the work of Mckenzie [81]. In biological systems, Ca 2+ spirals are observed in "overstimulated" conditions. Frog oocyte overexpressing muscarinic acetylcholine receptor produce spiral waves in c cyt upon stimulation [82].
During normal physiological function of many organs, directed rhythmic Ca 2+ waves are required, i.e. Ca 2+ phase waves propagating along a certain direction. For instance, rhythmic Ca 2+ phase waves are required for the contraction of the gastric pylorus [83], uterus [67,68], intestine [84] or urinary bladder [85]. Heart contraction is even more orchestrated with very fast directed Ca 2+ phase waves, but in this case the gap junctional transmission of action potentials is thought to be the relevant synchronizing process [22]. Nevertheless, it is not excluded that gap junctional transport of Ca 2+ and InsP 3 also plays role in the synchronization. This might explain the phenomena of delayed after-depolarization [86] or arrhythmias associated with altered operation of InsP 3 R [87]. In our model, directed waves observed as circular rings can be generated even in noisy systems, if there is one highly sensitive, wave initiator region. Incorporating one highly sensitive region is an effective way to harmonize networks and to determine the wave direction within a broad frequency range. However, the regulated network collapses, if the number of non-oscillating "signal-plateau" cells increase. Cells usually show signal-plateau response at high stimulation intensity, e.g. by strongly elevated levels of InsP 3 [29]. This corroborates the experimental observations in biological systems, i.e. spirals appear as a consequence of very extensive stimulation in spatially connected systems [82,88].
From a general viewpoint, noise can arise in any other excitable oscillating network, in which the individual units are connected by different synchronizing agents. Examples for coupled excitable units exist in the field of Biology (synchronization of cellular clocks in the suprachiasmatic nucleus [89]), Physics (Josephon junction circuits [90]), Chemistry (e.g. Belousov-Zhabotinsky reaction [91] or catalytic oxidation of carbon monoxide on a platinum surface [92]) or Social Sciences (social interaction such as waving of a human crowd during a football match [93], an audience clapping in synchrony [94], flashing of fireflies [95] or cricket chirps [96]). We propose that phenomena that we observed in our model might be extrapolated to other systems. For example the Kuramoto model [97] predicts a transition with increasing global coupling strength, at which the oscillators with originally different frequencies become coherent. Also, the global coupling strength was investigated in the limit of large number of oscillating sites [98,99]. Critical coupling was found to differentiate coherent and non-coherent macroscopic behaviors. Their general methods allowed to reducing the involved dynamical description of coupled equations to a finite number of differential equations for the macroscopic state of a system. We demonstrated how local coupling is involved as a globally synchronizing agent in our model, quantified by our synchronization index. Locally, this transition would occur around the "highly sensitive region" and this region would be the initiator (pacemaker) of the subsequent phase waves in originally non-synchronized noisy systems. The "highly sensitive region" might represent the area of higher frequencies or the area with a higher density of links depending on the system. Based on our observation that in all cases, a highly sensitive region gradually determined the behavior of the entire system, as is also observed in many model structures (S1-S3 Texts), one can deduce some general conclusions. This phenomenon may occur if I) each unit has its own machinery for the generation of oscillations, II) there is a certain level of noise within the system, III) there is one or more coupling agents modifying the oscillation frequencies of the coupled units and IV) there is a highly sensitive region allowing for a global synchronization.