{'fraction_non_alphanumeric': 0.04503415903704048, 'fraction_numerical': 0.01589440907189664, 'mean_word_length': 4.4455879871773245, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Shaping bursting by electrical coupling and noise

Gap-junctional coupling is an important way of communication between neurons and other excitable cells. Strong electrical coupling synchronizes activity across cell ensembles. Surprisingly, in the presence of noise synchronous oscillations generated by an electrically coupled network may differ qualitatively from the oscillations produced by uncoupled individual cells forming the network. A prominent example of such behavior is the synchronized bursting in islets of Langerhans formed by pancreatic \beta-cells, which in isolation are known to exhibit irregular spiking. At the heart of this intriguing phenomenon lies denoising, a remarkable ability of electrical coupling to diminish the effects of noise acting on individual cells. In this paper, we derive quantitative estimates characterizing denoising in electrically coupled networks of conductance-based models of square wave bursting cells. Our analysis reveals the interplay of the intrinsic properties of the individual cells and network topology and their respective contributions to this important effect. In particular, we show that networks on graphs with large algebraic connectivity or small total effective resistance are better equipped for implementing denoising. As a by-product of the analysis of denoising, we analytically estimate the rate with which trajectories converge to the synchronization subspace and the stability of the latter to random perturbations. These estimates reveal the role of the network topology in synchronization. The analysis is complemented by numerical simulations of electrically coupled conductance-based networks. Taken together, these results explain the mechanisms underlying synchronization and denoising in an important class of biological models.

Square wave bursting, one of the most common firing patterns, is characterized by alternating periods of fast oscillations of the membrane potential and quiescence [58,21]. Typically, cells that generate bursting, under different conditions exhibit other firing patterns such as periodic or aperiodic spiking or reside in the excitable regime [5,51,52,70]. Transitions between different dynamical regimes in excitable cells often signal important physiological or cognitive events such as changes in the rate of hormone secretion or neurotransmitter release as in the cases of pancreatic β−cells [65] and midbrain dopamine neurons [27]; or changes in respiratory rhythm or attentional state as in the cases of neurons in the Pre-Botzinger complex [4] or Locus Coeruleus [72] respectively. Not surprisingly, mathematical models have been used extensively to explain the origins of different firing patterns. For single cell models, mechanisms underlying various modes of electrical activity have been thoroughly studied [41,21,51,58,59,69]. Some of the techniques developed for single cell models extend to cover small networks [46,60]. However, mathematical analysis of large conductance-based networks without special assumptions on network topology is an outstanding problem. Furthermore, there is a growing body of experimental and theoretical studies indicating the importance of noise in shaping neuronal dynamics [10,15,26,31,36,39,42,45,55,61,67,71]. In the presence of noise in the network dynamics, its analysis becomes even a more challenging problem. The goal of this paper is to elucidate principal factors shaping synchronous activity in large electrically coupled networks of bursting capable cells forced by small white noise.
Many ingredients contribute to the output of neuronal networks. Among them, intrinsic properties of the individual cells (local dynamics) and the type and the structure of connections between cells (network topology) are probably the most important ones. Direct electrical coupling through gap-junctions is a common way of communication between neurons as well as between cells of the heart, pancreas, and other physiological systems [11]. The role of electrical coupling in shaping firing patterns generated by neuronal networks has been studied using many different techniques: the theory for weakly connected networks [32,35,62], Poincare maps [29,6,38,53], and Lyapunov functions [54], to name a few. In the present study, we consider a relatively less studied case of strong electrical coupling [12,54], for which we develop two complementary approaches based on center-manifold reduction [37] and fast-slow decomposition [1,56]. Importantly, our method covers networks with arbitrary topology, which allows us to study a large class of models and to reveal the role of the network topology in shaping network dynamics.
Under fairly general conditions, strong electrical coupling synchronizes activity across the network [48,49]. Therefore, one might expect that dynamics of electrically coupled networks of bursting cells will closely resemble that of a single cell provided the coupling is strong enough. This is true in general for deterministic models. However, in the presence of noise network dynamics, while still synchronous, can be qualitatively different from that of a single cell uncoupled from its neighbors. For instance, single cell models, which exhibit irregular spiking in isolation (Fig. 1a) can generate very regular synchronous bursting when they are coupled electrically (Fig. 1b). Likewise, a coupled network exhibiting synchronous spiking for extremely long period of time (Fig. 1c) may be formed from bursting cells (Fig. 1d). The first scenario illustrated by Fig. 1a,b was proposed in [65,64] to explain why pancreatic β−cells burst within electrically coupled islets of Langerhans, but in isolation exhibit irregular spiking. Numerical experiments and formal analysis in [65,64] show that noise shaping the dynamics of individual β−cells becomes less effective when these cells are coupled electrically. It is certainly intuitive (albeit not obvious mathematically) that in a coupled network the effects of uncorrelated stochastic processes acting on individual cells may become weaker due to averaging. However, a remarkable property of electrical coupling, which as the analysis below shows should not be taken for granted, is that the variability of the coupled system can be fully controlled by varying two parameters: the coupling strength and the size of the network. We identify analytical conditions, which guarantee that the variability of the coupled systems (the meaning of variability will be explained below) goes to zero as the network size and the coupling strength tend to infinity.
Besides bursting in electrically coupled islets of Langerhans, other phenomena where denoising plays a central role, include episodes of phasic firing in the Locus Coeruleus network [72] and enhanced reliability of neural responses in gap-junctionally coupled networks [50,68]. In this paper, building on an earlier analysis of denoising in networks of integrate-and-fire neurons [50] and our recent study of spontaneous activity in a closely related model of the Locus Coeruleus network [55], we derive quantitative estimates characterizing denoising in electrically coupled networks of conductance-based models. We find that the results obtained for integrate-and-fire models for individual cells do not extend automatically to conductance-based models with higher-dimensional state phase. We identify additional features of the local dynamics and coupling architecture that are needed to guarantee denoising. In particular, our analysis highlights the role of the bifurcation structure of the bursting cell model for denoising. It also elucidates the contribution of the network topology to this important effect. We show that networks on the graphs with large algebraic connectivity [20] or small total effective resistance [3] are better equipped for implementing denoising. As a by-product of the analysis of denoising, we analytically estimate the rate with which trajectories converge to the synchronization subspace and the stability of the latter to random perturbations. Taken together, these results explain the mechanisms underlying synchronization and denoising in an important class of biological models.
The organization of the paper is as follows. In the next section, we formulate our assumptions on the single cell model ( §2.1) and explain how the network is formed ( §2.2). In § §2.3,2.4, we collect necessary information from the algebraic graph theory [2], which will be used for describing the role of the network topology in dynamical phenomena analyzed in this paper. In Section 3, we introduce two scenarios (A and B) leading to distinct firing patterns produced by the single cell models and by synchronized networks of these models. In the remainder of this section, we analyze the first of these scenarios illustrated in Fig. 1  (a,b). First, in Lemma 3.1, we show that in the single cell model, bursting can be destroyed with small noise (Figure 1a). This counter-intuitive result relies on the presence of a slow timescale in the dynamics of bursting. Even very small noise can have significant qualitative effects on the dynamics of a slow-fast system, because during the (long) periods of slow evolution there is sufficient time for large deviations (which are extremely unlikely on time intervals of order O(1)) to develop. The proof of the lemma uses large deviation estimates and is adapted from [17]. Thus, given a deterministic model of a bursting cell, one can switch its dynamics into irregular spiking by adding noise. We then show that when many such cells are coupled together, the effects of noise weaken and the bursting of the underlying deterministic model reemerges ( Figure 1b).
The analysis of the coupled system proceeds in two steps. First, we use the center-manifold reduction [37] to approximate the coupled system near an excitable equilibrium by a simpler lower-dimensional system of ordinary differential equations. Second, thanks to the gradient structure of the reduced problem, we can accurately estimate expected time that a trajectory of the fast subsystem of the coupled system spends near the excitable equilibrium. We show that this time is much longer for the network model than for the single cell one. This analysis (based on large deviation estimates [18]) yields one way of quantitative description of denoising.
In Section 4, we present a complementary method based on a slow-fast decomposition. We show that when the coupling is strong, network dynamics near the excitable equilibrium splits into two modes: fast synchronization and ultra-slow noise-driven escape from the basin of attraction of the equilibrium along a low-dimensional synchronization subspace. The results of this section yield valuable insights into synchronization properties of the coupled system. In particular, we estimate the rate of convergence of trajectories to the synchronization subspace and the stability of the latter against random perturbations. The estimates show explicitly the contribution of the structural properties of the network to stability properties of the synchronization subspace.
In Section 5, we take a look at denoising from a slightly different angle. Specifically, we study the linearization of the coupled system near an excitable equilibrium directly without invoking center manifold reduction. The analysis shows that when the dynamics of the individual cells lives in multidimensional phase space (unlike integrate-and-fire models or those of excitable cells near a saddle-node bifurcation), unless the coupling is full rank (see [48] for the definition of full versus partial rank coupling) denoising should not be expected. These results show that the proximity to a saddle-node bifurcation in the fast subsystem of a square wave bursting neuron model, which makes the dynamics near the excitable equilibrium essentially one-dimensional, is critical for observing distinct dynamics generated by the single cell and network models. Therefore, the bifurcation structure of the bursting cell models like those used in [65,64] and in this paper is important for the interplay between electrical coupling and noise in shaping the dynamics of the coupled network. Finally, the results of this study are summarized in Section 6.

The model
In this section, we introduce a model that will be studied in the remainder of this paper. First, we formulate the assumptions on the single cell model and describe its main dynamical regimes. Next, we introduce the coupled system. At the end of this section, we review some notions and facts from the algebraic graph theory [2], which will be useful for characterizing the contribution of the network topology to dynamics.

The single cell model
In numerical experiments throughout this paper, we are going to use a conductance based model of a pancreatic β−cell due to Chay [5]. The analysis in the following sections does not depend on any specific details of this model and more general assumptions will be formulated below. However, we believe that it is instructive to start from a concrete model to make the biological interpretation of the analysis that follows transparent. With this in mind, following [5], we introduce the a system of differential equations modeling electrochemical dynamics in the β−cells: Here v, n, and u stand for the cell membrane potential, gating variable, and the concentration of calcium respectively. I ion (v, n, u), I Ca (v), n ∞ (v) and τ (v) denote nonlinear functions, which are used for modeling ionic currents. C m denotes membrane capacitance. The small parameter > 0 multiplying the right hand side of the third equation reflects the separation of the timescales of the calcium dynamics and the fast variables v and n. The right hand sides of the first two equations also contain independent copies of Gaussian white noiseẇ 1,2 , which account for the deviations from the deterministic dynamics due to various fluctuations [66,71]. For further details of (2.1)-(2.2) including the values of parameters, we refer the reader to Appendix A and [5].
To describe the structure of (2.1)-(2.3), it is convenient to rewrite it in a more general forṁ Under the variation of y, the fast subsystem has the bifurcation structure as shown schematically in Fig.  2a. Specifically, (PO) There exists y hc ∈ R such that for each y < y hc , Equation (2.6) has an exponentially stable limit cycle of period T (y): The family of the limit cycles, L = y<y hc L y , forms a cylinder in R 3 , which terminates at a homoclinic loop at y = y hc [37] (Fig. 2a).
(EQ) There is a branch of asymptotically stable equilibria of (2.6), E = y>ysn E y , E y = {x = ψ(y)}, which terminates at a saddle-node bifurcation at y = y sn < y hc (Figure 2a).
(LS) For each y ∈ R, the ω−limit set of almost all trajectories of (2.6) belongs to L y E y .
Deterministic models of bursting are well understood (see, e.g., [21,41,51,58]). For small > 0, (2.4) 0 and (2.5) features three main regimes: bursting, spiking, and quiescence (or excitable). In the former, the trajectory alternates between drifting along the cylinder L foliated by periodic orbits of the fast subsystem and the curve of equilibria (see Fig. 2a,d). Alternatively, (2.4) 0 and (2.5) may have a stable limit cycle near L (spiking) or a stable fixed point near E (excitable). The latter two regimes are illustrated in Fig. 2 (b,e) and (c) respectively. The analytical conditions for bursting, spiking, and excitable regimes are given in Appendix B.

The electrically coupled network
Next, we consider a gap-junctionally coupled ensemble of n cells, whose dynamics is generated by (2.1)-(2.3). In the coupled network, Cell i, i ∈ [n] := {1, 2, . . . , n} receives current from other cells in the network. Conductance a ij > 0 if Cell i and Cell j are connected and a ij = 0, otherwise. Without loss of generality, we set a ii = 0, i ∈ [n], and denote A = (a ij ). Nondimensional parameter g > 0 is used to control the coupling strength.
By including the coupling current (2.8) into the models of individual cells (2.1)-(2.3), we obtain a differential equation model of the electrically coupled network where W (i) = (w (i,1) , w (i,2) ) T are independent copies of 2D Brownian motion. Using the notation, which we adopted for the single cell model in (2.4) and (2.5), we rewrite the coupled system in the following more general form:Ẋ (2) , y (2) ), . . . , g(x (n) , y (n) ) , and ⊗ stands for the Kronecker product.

The graph of the network
Dynamics of the coupled system depends on the spectrum of matrix L appearing in the coupling operator in (2.12). The eigenvalues of L in turn depend on the structure of the graph of the network. Throughout this paper, we will repeatedly use the relation between the spectrum of L and the structural properties of the network to study how the network topology affects its dynamics. To this end, we will need certain constructions and results from the algebraic graph theory, which we review below following [2].
Let G = (V (G), E(G)) denote the graph of interactions between the cells in the network. Here, V (G) = {v 1 , v 2 , . . . , v n } and E(G) = {e 1 , e 2 , . . . , e m } denote the sets of vertices (cells) and edges (pairs of connected cells), respectively. Throughout this paper, we assume that G is an undirected connected graph.
It is instructive to consider first the case when all nonzero conductances are equal to 1: 1 Cell i and Cell j are connected, 0 otherwise. (2.15) As before, we set the diagonal elements of A to zero, a ii = 0. Matrix A = (a ij ) in (2.15) is called the adjacency matrix of G and is called the Laplacian of G. By deg(v i ) we denote the degree of v i , i.e., the number of edges incident to v i . Alternatively, the graph Laplacian can be defined by where H ∈ R m×n is the coboundary matrix of G [2]. The definition of H uses an orientation of the edges of G. For each edge e j = (v j 1 , v j 2 ) ∈ V (G) × V (G), we specify the positive and negative ends; e.g., let v j 1 be the negative end of e j 2 . Then the coboundary matrix of G is defined as follows (cf. [2]) Definitions (2.16) and (2.17) are equivalent (cf. [2]). From either of these definitions, it is easy to see that λ 1 (L) = 0 is an eigenvalue of L. If G is a connected graph then the zero eigenvalue is simple and all other eigenvalues are positive (cf. [20]) Remark 2.1. Following a common in the algebraic graph theory convention, we will refer to the eigenvalues of L as the eigenvalues of G.
The second eigenvalue a(G) = λ 2 (L) is called the algebraic connectivity of G, because it yields a lower bound for the edge and the vertex connectivity of G [20]. The algebraic connectivity is important for a variety of combinatorial, probabilistic, and dynamical aspects of the network analysis. In particular, it is used in the studies of the graph expansion [33], random walks [3], and synchronization of dynamical networks [23,47,55]. Below, we show that a(G) determines the rate of convergence to synchrony in the coupled model (2.12) and (2.13).
Another spectral function of G, which will be useful in the analysis of the coupled system is the total effective resistance of G (cf. [40,73]) For electrical and graph-theoretic interpretations of R(G) as well as for numerous applications, we refer the reader to [3,25].
All definitions and constructions, which we reviewed above, naturally extend to cover nonhomogeneous networks with different conductances a ij . To this end, we define conductance matrix The algebraic connectivity and the total effective resistance of G are defined through the EVs of L(G) as before.

Examples of the network connectivity
In this subsection, we present several examples of network connectivity, which will be useful for illuminating the role of the network topology in the dynamical phenomena analyzed in this paper.
Example 2.2. One of the most common examples of local connectivity is a 1D nearest neighbor coupling. It has been used in many studies of the coupled ensembles of β−cells (see, e.g., [65,64]). We use the 1D nearest-neighbor coupling in most numerical experiments throughout this paper.
In this configuration, each cell in the interior of the array is coupled to two nearest neighbors. This leads to the following expression for the coupling current: The coupling currents for the cells on the boundary are given by The corresponding graph Laplacian is (2.23) The undirected graph corresponding to the 1D nearest neighbor coupling scheme is called a path and is denoted by P n [2].
Clearly, the total number of edges in the graph (connections in the network) is one of the important factors shaping the dynamics of the coupled system. However, the way how these connections are distributed, i.e. the connectivity of the network, is also important as the analysis of the following example will show. Example 2.3. Consider two graphs on n nodes of degree d = 4 (i.e., each node in each of these graphs has precisely four edges incident to it). Such graphs are referred to as (n, d)−graphs. Thus, each of these graphs has 2n edges. To illuminate the role of connectivity, we assign to these graphs two different connectivity patterns as shown schematically in Fig. 3. The graph in Fig. 3a has symmetric connections. The edges of the graph in Fig. 3b were selected randomly.
Specifically, the second graph was generated using so-called permutation model of a (n, d) random graph of even degree d = 2d In this model, one chooses at randomd permutations π 1 , π 2 , . . . , πd in the symmetric group S n . Then the edges between n vertices v 1 , v 2 , . . . , v n of G are generated as follows Spectral properties of the random graphs similar to the one constructed in the previous example have important implications for analyzing dynamical phenomena like synchronization and denoising in large networks. In the remainder of this subsection, we will discuss several facts about the spectra of random graphs that are particularly relevant to the analysis that follows. Our review is very brief and we refer an interested reader for more information and extensive bibliography to an excellent survey by Hoory et al [33].
First, we note that the random graph constructed in Example 2.3 is a (spectral) expander, which means that for some positive α λ 2 (G n ) ≥ α, uniformly in n ∈ N, (2.24) where n stands for the power of the set of vertices [33,63]. For those readers who have not seen this condition before, we note that this property fails to hold for any family of lattices (see Section 5 in [47] for a related discussion). In particular, for a path on n vertices (cf. Example 2.2), the second eigenvalue tends to zero as n → ∞. Thus, the existence of a uniform bound for the second eigenvalue postulated by 2.24 is a special, if not counter-intuitive, property. Nonetheless, it holds for random graphs. There are also known explicit (nonrandom) algorithms generating expanders, including the celebrated Ramanujan graphs [44,43].
In many applications, it is desirable to have a large bound on the expansion constant α. However, for [33]. A remarkable property of the family of (n, d) random graphs is that for n 1 they posses nearly optimal expansion constant with overwhelming probability. In particular, it is known that for any > 0 where G n stands for the family of (n, d) random graphs of degree d ≥ 3 and n 1 [19].

Transitions to bursting in the coupled model
Having reviewed bursting in deterministic systems and the definition of the coupled network, we now turn to the main theme of this work -the roles of noise and electrical coupling in shaping firing patterns of the coupled system (2.12) and (2.13). Fig. 1 shows transitions in the coupled system's dynamics observed over long intervals of time upon increasing the coupling strength. In the first case, which in the sequel we will refer to as Scenario A, irregular spiking patterns for weak coupling are transformed into fairly regular bursting patterns for sufficiently strong coupling (see Fig. 1a,b). In the second case, Scenario B, robust very irregular bursting becomes synchronous spiking when coupling strength is increased.
We show that at the heart of both transitions lies denoising, the mechanism responsible for greatly diminishing the effects of noise on network dynamics. Below, we discuss both scenarios illustrating them with numerical results. For Scenario A, we also develop analytical estimates characterizing denoising. Our goal is to show how statistical features of the firing patterns depend on the coupling strength, excitability, and network topology. The analysis of Scenario B can be developed along the same lines, but it requires certain additional techniques for dealing with the analysis of trajectories near a limit cycle. These extensions will be considered in the future work.

Scenario A: The single-cell model
In this and in the following subsections, we discuss Scenario A in more detail. We start with the single-cell model (2.4) and (2.5).
Assume that the deterministic model (2.4) 0 and (2.5) is in the bursting regime, i.e., a typical trajectory alternates between drifting along a curve of stable equilibria of the fast subsystem, E, and a cylinder of limit cycles, L (see Fig. 2a). Let us focus now on a slow evolution along E. Sinceẏ = O( ) (cf. (2.5)), on time intervals o( −1 ) long, the dynamics of the slow-fast system (2.4)-(2.5) is approximated by the fast subsystemẋ = f (x, y) + Σẇ, where y ∈ (y sn , y hc ) is fixed, Note that the transitions take place along the unstable manifold of the saddle. This shows that the center manifold reduction used in Sections 3 and 4 adequately describes the dynamics observed numerically.
The frozen system (3.1) 0 is bistable. It has two co-existing attractors: the stable equilibrium E y and the limit cycle L y (see Fig. 4 a). A trajectory of the randomly perturbed system (3.1) alternates between the basins of E y and L y , B(E y ) and B(L y ), separated by the stable manifold of the saddle point O (see Fig. 4b). In fact, most of the time it spends in small neighborhoods of E y and L y . The plot in Fig. 4b shows that a typical trajectory leaves the neighborhood of E y along the weak stable manifold of the sink, i.e., the dynamics near E y is effectively one-dimensional. We will use this observation to simplify the analysis of the coupled system by reducing it (via the center manifold theorem [8]) to a simpler system.
The following exit problem is instrumental for understanding the effects of noise on the dynamics of the bistable system (3.1). For a trajectory of (3.1), which starts from a deterministic initial condition Below we show that one can choose σ( ) > 0 such that with overwhelming probability for small > 0, and at the same time σ( ) → 0 as → 0 (see Lemma 3.1). Therefore, with arbitrarily small noise one can make the trajectory of the frozen system (3.1) leave B(E y ) in time insufficient for y to change by O(1) amount thus destroying bursting.
In Fig. 5a, we show a typical trajectory of (2.1)-(2.3) superimposed on the bifurcation diagram of the fast subsystem. Note that the trajectory can not advance far enough along the curve of stable equilibria E, because it is forced to jump to the vicinity of L by noise. Thus, it lands on L near y hc every time and does not have enough room to generate many spikes before it reaches the right end of L and jumps back to E. This results in very irregular bursting patterns with very few spikes in one burst (see Fig. 5b,c). In theory, for sufficiently small > 0 and small σ > 0 one can make the probability of clusters of 2 and more spikes in the timeseries of (2.1)-(2.3) arbitrarily low. This means, that one can transform bursting into irregular spiking by adding small noise. Showing this numerically, however, requires integrating stiff stochastic differential equations over very long time. Thus, we did not strive to achieve irregular spiking in our experiments resorting to irregular bursting patterns with very few spikes in Fig. 5b, which already illustrate this effect.
The noise-induced irregular firing pattern shown in Fig. 5b is characterized by the approximetelly geometric distribution of spikes in one cluster (see Fig. 5c). The geometric distribution has its origins in the exponential distribution of the first exit time τ (E y , x 0 ) [13], which implies roughly that the distance from the landing point on L, y l , to the right hand end of L, y hc , y l − y hc is distributed approximately exponentially. The exponential distribution of y l − y hc translates into the geometric distribution of the number of spikes in one burst. Later we will see that this distribution is qualitatively different for the coupled system (see Fig. 5f). The distinct probability distributions in Fig. 5c and Fig. 5f corresponding to different values of the coupling strength show that the transition in the network dynamics has taken place.
In conclusion, we note that while the fact that we were able to destroy deterministic bursting with noise may not seem very surprising, the possibility of doing this with small noise is far from obvious. Note that the equilibria of the fast subsystem, E y , for y near y hc are extremely stable. The slow-fast structure of the vector field is the key to this important effect. Extremely slow evolution along E gives the trajectory of the frozen system enough time to develop large deviations (necessary to leave B(E y )), which are highly unlikely on time intervals O(1) long. This is a general mechanism by which adding noise to slow-fast systems may create qualitatively new dynamical regimes [17].
Before we move on to discuss the coupled system, we prove the following lemma, which provides an estimate of the noise intensity sufficient for destroying bursting.
Lemma 3.1. Let k > 0 and Σ = diag (σ, kσ) be a matrix defining the noise intensities in (3.1). Then there exists C 1 > 0 such that for every C 2 ≥ C 1 and a trajectory of (3.1) with initial condition x(0) ∈ B (E y ) leaves B (E y ) in time o( −1 ) with probability converging to 1 as → 0.
Proof. By rescaling x 2 in (3.1), one can always achieve k = 1, which we assume without loss of generality.
Thanks to the large deviation theory, we have the following estimate for the first exit time τ (E y , x 0 ) (cf. Theorem 4.4.2, [18]): for any h > 0, where positive constant V y is the minimum of the quasipotential associated with the randomly perturbed system (3.1). The definition and the properties of the quasipotential can be found in [18].

Scenario A: The coupled system
Next, we turn to the analysis of the coupled system. We want to understand how bursting is recovered for larger values of the coupling strength (see Fig. 5 d-f). Thus, we consider the coupled system (2.12) and (2.13). Below, we use two simplifying asuumptions, which let us avoid certain technical details, which are peripheral to the mechanism analyzed below, First, we focus on the fast subsystem ignoring O( ) changes in the slow variables: where X = (x (1) , x (2) , . . . , x (n) ) T ∈ R d × · · · × R d ∼ = R nd . Furthermore, it can be shown that with the diffusive coupling like in (3.6), y i s synchronize and remain close after some initial transients, provided the coupling is sufficiently strong (see [48,49]). Thus, we set the frozen slow variables for all subsystems to the same value y ∈ (y sn , y hc ) so that For the coupled system (3.6), we consider the problem of exit from the basin of the rest state E y = E y × · · · × E y ∈ R 2n . To this end, we define The key difference between the coupled system (3.6) and the individual subsystems (3.1) is that the former is much less susceptible to the effects of noise. The robustness of the coupled system to noise manifests itself in the disparity of the exit times: More precisely, the asymptotic relation between the exit times corresponding to the single cell and coupled models is given in the following theorem.
Theorem 3.2. Suppose that (3.1) is close to a nondegenerate saddle-node bifurcation 1 Then for some σ 0 > 0 and g 0 > 0, and for all 0 < σ < σ 0 and g > g 0 the following asymptotic relations hold:

7)
where denotes logarithmic asymptotics (cf. (3.4)), n is the number of the cells in the network, and C 3 is a positive constant independent of σ and n.  .7) show that for the level of noise chosen in Scenario A, by taking sufficiently strong coupling with sufficiently many cells in the network, one can make τ (B(E y ), x 0 ) longer than the time necessary for y to reach the vicinity of y sn . This means that the level of noise, which prevents a single cell model from bursting, does not affect bursting in the coupled system (compare Fig. 5 a and d).
Proof. The proof of the theorem follows from the analysis of a closely related model in [55]. For completeness, we outline the main steps of the proof and refer the interested reader to [55] for further details. The proof consists of several steps, the main of which are the center manifold reduction of the single cell and coupled models near the saddle-node bifurcation (cf. (3.12)) and the variational interpretation of the reduced systems (cf. (3.15)).
1. By our assumptions on the fast subsystem, the Jacobian Df (0, 0) has a 1D kernel. Denote e ∈ ker Df (0, 0)/{0} and p ∈ ker(Df (0, 0)) T such that p T e = 1. (3.8) In addition to standard nondegeneracy and transversality conditions for a saddle-node bifurcation we assume that a 3 = p T Je = 0, (3.11) where J 1 is the matrix involved in the coupling operator (see (3.6)). Condition (3.11) guarantees that the projection of the coupling onto the center subspace is not trivial. Without loss of generality, we assume that nonzero coefficients a 1,2,3 are positive.
2. Under the conditions in 1., near the bifurcating equilibrium the coupled model can be reduced to an n−dimensional slow manifold. The reduced system (after appropriate rescaling of the dependent and independent variables and dropping higher order terms) has the following form: where γ is the rescaled coupling strength, 1 n = (1, 1, . . . , 1) T ∈ R n , and W is the standard Brownian motion in R n 1 .
Similarly, the single cell model is reduced to the following 1D equation: 3. We recast (3.13) and (3.12) as randomly perturbed gradient systems. The former is rewritten as (3.14) For the latter, we use the structure of the coupling matrix (2.22), to reduce it to the following forṁ is a square root of the nonnegative definite conductance matrix C (cf. 2.21). 4. The large deviation estimates yield the logarithmic asymptotics for the exit times associated with the stable fixed point of (3.14) (cf. Theorems 2.1 & 3.1 in [18]): where Φ(1) = 4/3 is the value of the potential at the barrier at ζ = 1. Thus, which shows the first relation in (3.7).
5. The (deterministic) coupled system (3.12) 0 has a stable fixed point at z = −1 n . The basin of attraction of this fixed point is bounded by The estimate of the exit times from D follows from the analysis of the minima of the potential function U γ (·) on ∂D. In [55], we prove that for γ ≥ 2λ 1 (L 1 ) −1 , U γ (z) achieves its minimal value on ∂D at z = 1 n : (Here, λ 1 (L 1 ) denotes the smallest eigenvalue of matrix L 1 , obtained from L by deleting the first row and the first column.) Knowing the minimum of the potential function on the boundary of the basin of attraction of the stable fixed point, the large deviation estimates for the randomly perturbed gradient systems yield (3.20)

Scenario B
In this scenario, the deterministic single cell model is tuned to be in the spiking regime, i.e., it has a stable limit cycle L yc , y c < y hc in the vicinity of L (see Fig. 2b,e). When noise is added to the modeling equations (2.1)-(2.3), it forces the trajectory to leave the basin of attraction of the limit cycle B(L yc ) once in a while (see Fig. 6a), producing irregular bursting (Fig. 6). This mechanism of irregular bursting was studied in detail in [31]. In particular, it is shown in [31] that the number of spikes in one burst has asymptotically geometric distribution with parameter p ≈ exp{C 4 σ −2 } for some positive constant C 4 .
The normalized histogram in Fig. 6c shows that with the level of noise chosen for this experiment, the system generates very long irregular bursts. In the second row in Fig 6d-f, we show the results for two coupled cells and the coupling strength g = 50. Taking just two cells already changes the statistics of bursting significantly (compare Fig. 6 c and f). For three coupled cells, the bursts become even longer preserving the geometric distribution (see Fig. 6 i). For ten coupled cells and the coupling strength g = 500, we were unable to detect a single termination of spiking activity. While the coupled system is technically still in the bursting regime, the bursts become so long that for all practical purposes, one can speak about effective transition to spiking (see Fig. 6 g,h).
The changes in statistical properties of the firing patterns generated by the coupled system for increasing values of the coupling strength are due to denoising. Here, electrical coupling diminishes the effects of noise a b c d e f g h i Figure 6: Scenario B. The underlying deterministic system is in the spiking regime (see Fig. 2b). Adding small noise transforms spiking into irregular bursting shown in plots (a,b). Coupling two cells together already increases the duration of bursts significantly (d,e). For ten coupled cells the duration of typical bursts is extremely long so that for all practical purposes (g,h), it can be considered as a transition to spiking. The normalized histograms for the number of spikes in one burst are plotted for a single cell, two, and three coupled cells in c, f, and i respectively. We were unable to numerically detect any bursts for the 10 coupled cells due to their extremely long duration. All histograms show approximately geometric distribution confirming the large deviations nature of the bursting patterns. The parameters of the geometric distribution (i.e., the expected value of the number spikes per burst) increases significantly as we go from one cell to three cell network (c,f,i).
on synchronous limit cycle oscillations and results in longer times that the trajectory of the coupled system spends in B(L yc ) for larger values of g. Qualitatively, this is the situation that we have already analyzed for Scenario A. However, the study of the Scenario B involves certain additional technical details needed for extending the analysis to systems of coupled limit cycle oscillators (see [49,48]). We will consider the analysis of Scenario B in the future work.

Synchronization
In this section, we study synchronization in electrically coupled network (2.12) and (2.13). Our motivation is twofold. On the one hand, synchronization is an important aspect of the network dynamics, because for sufficiently strong coupling the network activity becomes synchronized. In particular, all numerical results shown in this paper for the coupled system feature synchronous activity. On the other hand, we show that the mechanism of synchronization is closely related to that of denoising. In particular, the estimates of stochastic stability of synchrony, which we derive below, reveal the contribution of the network topology and the strength of connections distribution to synchronization and denoising. We show that networks with larger algebraic connectivity and smaller total effective resistance are more effective for implementing denoising.
To illustrate this point, we use random and symmetric degree-4 graphs defined in Example 2.3. Our analysis and numerics show that random graphs and, more generally, expanders have good synchronization and denoising properties. The analytical approach, which we develop in this section, complements that of the previous section.

The new coordinates
To study synchronization, we introduce a special system of coordinates, which to leading order decouples the two principal modes of the system's dynamics: fast synchronization and ultra slow excitation due to noise. Remarkably, the latter is captured by a scalar stochastic ordinary differential equation. The slowfast decomposition is used to show that synchronization takes place for sufficiently strong coupling and to quantify various aspects of the synchronized dynamics.
The new coordinate system takes into account the structure of the network. To this end, we will need a spanning tree of G,G = (V,Ẽ), i.e., a subgraph of G with n = |V (G)| vertices, n − 1 edges, and containing no loops [2]. Having chosen the spanning treeG, letH ∈ R (n−1)×n denote the coboundary matrix corresponding toG. Matrix S, which we introduce in the next lemma will be useful for constructing the new coordinates.  Then SS T = I n−1 and S1 n = 0. Proof. Properties (4.2) follow from the definition (4.1).
The synchronization subspace spanned by 1 n coincides with the kernel of S, while the columns of S T form an orthonormal basis of 1 n ⊥ . These two subspaces are important for studying synchronization, which motivates the following coordinate transformation where ξ = Sz and η = n −1 1 n T z. Note that |ξ| = |P 1n ⊥ z| measures the distance of the solution of the coupled system (3.12) to the synchronization subspace corresponding to ξ = 0. Here, P 1n ⊥ stands for the orthogonal projector onto 1 n ⊥ .

The slow-fast system
Throughout this section, we work with the reduced (rescaled) system (3.12), which we derived using the center manifold approximation of (3.6) near the excitable equilibrium. In Section 3, we analyzed (3.12) using its gradient structure and the large deviations estimates. This time, we use a complementary approach by first identifying and then exploiting the slow-fast structure of (3.12).
By projecting (3.12) onto 1 n ⊥ and 1 n , we obtaiṅ ξ = (−γL + 2ηI n−1 )ξ + σSẆ + O(|ξ| 2 ), (4.5) where matrixL = SLS T is the unique solution of the matrix equation Here,Ẇ andẇ denote the Gaussian white noise processes in R n and R 1 respectively. In the derivation of (4.6), we use the fact that 1 n TẆ and n −1/2ẇ are identically distributed. For details of the derivation of (4.5) and (4.6), we refer the interested reader to Lemma 5.3 of [55]. About the reduced matrixL, the following is known:L is a positive definite matrix, whose spectrum consists of the nonzero eigenvalues of L (cf. Lemma 2.6, [47]). Moreover, S maps the generalized eigenspaces of the nonzero eigenvalues of L bijectively onto those ofL [47]. Equation (4.5) captures the dynamics along the orthogonal complement of the synchronization subspace 1 n . Thus, it describes the process of synchronization. On the other hand, Equation (4.6) tracks the motion along the synchronization subspace. SinceL is positive definite, in the strong coupling regime (γ 1), it follows from (4.5) that the trajectory of the full system (4.5) and (4.6) relaxes to an O(σ) neighborhood of the synchronization subspace, 1 n , at O(γ) rate. Equation (4.6) to leading order (with ξ ≈ 0) is a standard model of a particle in a potential well forced by noise. On the time intervals not exceeding the Kramer's time O(exp{C 5 σ −2 }) for some C 5 > 0, η remains close to −1 for the most of the time. This sums up the qualitative dynamics of (4.5) and (4.6). In the remainder of this section, we study it in more detail.

The fast subsystem: synchronization
The analysis of the fast subsystem elucidates several important aspects of synchronization in electrically coupled networks. In particular, we estimate the rate of synchronization in terms of the network connectivity. We then study robustness of synchrony to noise.
The stability of the synchronization subspace is determined by the linear part of (4.5): (4.7) Assuming for simplicity the deterministic initial condition ξ(0) = ξ 0 ∈ R n−1 , we compute the mean vector and the covariance matrix of ξ(t): It follows from (4.8) and the geometric interpretation of ξ, that the rate of convergence to the synchronization subspace is set by the product of the strength of coupling γ and the algebraic connectivity a(G): where P 1n ⊥ · stands for the orthogonal projection onto 1 n ⊥ and C 6 is a positive constant independent from γ and a. Therefore, networks with larger algebraic connectivity synchronize faster. There are many fine results in the spectral graph theory relating algebraic connectivity to various structural properties of the network (see [33] and references therein). These results can be used via (4.10) to elucidate the contribution of the network topology to synchronization properties of the coupled system. In particular, (4.10) shows that networks on random graphs and on expanders in general (cf. Example 2.3) admit a lower bound on the synchronization rate, which is uniform in the size of the network n.
In the presence of noise, it is important to know how well the synchrony can withstand stochastic perturbations. This leads to the question of stochastic stability. There are many ways for measuring stability of synchrony to random perturbations [30]. In this paper, we use the mean square stability, which provides a natural metric for the problem at hand. Specifically, we are interested in transverse stability of the synchronization subspace. From (4.10), we know that Thus, we next look at the second moments, i.e., we estimate the dispersion of the trajectories around 1 n : var P 1n ⊥ z(t) = n−1 i var ξ i (t) =: var ξ(t). (4.11) From (4.9), we find that var ξ(t) = Tr cov ξ(t) = σ 2 2γ Tr L −1 I n−1 − exp{−2tγL} .
SinceL is positive definite, var ξ(t) has a finite asymptotic value where R(G) is the total effective resistance of the weighted graph of the network G (cf. 2.20). Thus, var ξ provides a convenient measure of stochastic stability of the synchronization subspace. The smaller the value of var ξ, the more stable the synchrony is. In the next section, we will use the asymptotic value of the variance to estimate the effectiveness of denoising.
Estimate (4.12) shows that the stochastic stability of the synchronization subspace is fully determined by the strength of coupling and the total effective resistance of the network. Similarly to the algebraic connectivity of the graph, the value of the total effective resistance R(G) can be related to the structure of the graph and the weight distribution [25]. However, while the rate of convergence to the synchronization subspace depends only on the leading nonzero eigenvalue of the graph Laplacian, the description of stochastic stability requires the entire spectrum of G. The information about higher eigenvalues of G in general is hard to obtain. However, as one can see from the following crude estimate of R(G) graphs with larger algebraic connectivity enjoy better bounds on the total effective resistance. In fact, for expanders, from (4.13) one gets a bound on R(G) = O(n 2 ), which can not be improved (up to a multiplicative constant). Therefore, we expect that graphs with good bounds on algebraic connectivity, like a random graph in Example 2.3, are robust against random perturbations. Below, we will illustrate this point with numerical results.
Above we have used the analysis of the fast subsystem to gain useful information about synchronization properties of the coupled network. Further, with these results at hand one can easily understand the 1D slow equation (4.6) and get a complete description of the slow-fast system (4.5) and (4.6). This is our next step.
Since on time intervals not exceeding the Kramer's time, with high probability ξ(t) = O(σ), we approximate (4.6) byη = f (η) + σ √ nẇ . (4.14) From (4.14), we estimate the time that a typical trajectory spends in the basin of the rest state η = −1 Comparison of (3.20) and (4.15) shows that the results of the fast-slow analysis are consistent with those of the previous section. Furthermore, through the estimates of stability of the synchronization subspace (4.10) and (4.12) the analysis of this section elucidates the contribution of the structural properties of the network to its synchronization properties.

The mechanism of denoising
In this section, we continue to study denoising, the mechanism responsible for variability of activity patterns in electrically coupled network (2.9)-(2.10). The analysis of synchronization and denoising in the previous sections relied on the proximity of the models of individual cells to the saddle-node bifurcation. To clarify to what extent the mechanism of denoising depends on this assumption and to elucidate the scope of its applicability, in this section, we analyze the coupled system around the stable fixed point without assuming anything about the bifurcation structure of the problem.
Because the dependence of the single cell models on the slow variable is not essential for the mechanism of denoising per se, we omit to indicate it explicitly to simplify notation. Thus, in this section, we consideṙ where f : R d → R d is a smooth function, w is a standard Brownian motion in R d , and Σ ∈ R d×d is a nonsingular matrix. We assume that (5.1) 0 has a stable equilibrium at the origin, i.e., where A s := 0.5(A + A T ) > 0 is a positive definite matrix. The coupled system has the following forṁ where as before . . , f (x (n) )) T , L = H T CH, and J is d × d matrix, such that J s is positive definite. Finally, W stands for the Brownian motion in R nd The coupled system has an equilibrium at the origin in R nd . The linearization of (5.3) about this equilibrium is given bẏ In the previous sections, to quantify the effect of denoising we used the times of the first exit from the basins of stable equilibria of the single cell and coupled models (cf. (3.20) and (4.15)). Since in this section we do not assume that the local systems are located near a saddle-node bifurcation, we no longer can rely on explicit estimation of the first exit times. Thus, we seek other means for measuring the effectiveness of denoising. To this end, we recall that in the previous section we found that the asymptotic value of the variance of the trajectories near the synchronization subspace to provide a convenient measure for estimating stochastic stability. We thus adapt it to our current purpose. Specifically, we use to measure the variability of the trajectories of the coupled system X(t) = (x (1) (t), x (2) (t), . . . , x (n) (t)). We estimate the effect of denoising by comparing max i∈[n] var x (i) to var x, where x(t) is the solution of the local subsystem (5.1). We identify conditions which guarantee that and show how the former quantity depends on the coupling strength, network size and topology, as well as on the intrinsic properties of the local system (5.1).
The following lemma is the key for understanding the mechanism of denoising.
In the process of proving Lemma 5.1, we derive the following estimates, which are of independent interest.
Proof. The proof of Lemma 5.1 consists of the following steps.
We want to understand how the variability of the coupled system can be smaller than that of a local subsystem. Estimate (5.5) suggests a possible scenario. Note that the first term on the right hand side of (5.5) can be made arbitrarily small by increasing n. The second term can be controlled by g provided that λ 1 (B s ) remains O(1) as g → ∞. Thus, we need to understand the behavior of λ 1 (B s ) for increasing values of g. In the following lemma, we show that the latter depends on the coupling architecture and identify two cases of full and partial coupling, which are important in the context of denoising.
where J s ≥ 0, A s > 0, and 0 < δ 1 (cf. (5.4)). If J s is a full rank matrix then Otherwise, let k = dim ker J s > 0, choose {q 1 , q 2 , . . . , q k } an orthonormal basis for ker J s , and define an Then Proof. For small δ > 0, the EVs of B perturb smoothly from the EVs ofL ⊗ J s (cf. [24]). If J s > 0 then λ 1 (L ⊗ J s ) = λ 1 (L)λ 1 (J s ) and (5.21) follows. Otherwise, 0 is an EV ofL ⊗ J s of multiplicity k(n − 1). We construct an orthonormal basis for the 0−eigenspace ofL ⊗ J s where e i = (δ i 1 , δ i 2 , . . . , δ i n−1 ) and δ i j stands for the Kronecker delta. By the perturbation results for the multiple eigenvalues of symmetric matrices (cf. Appendix, [24]), we have It is easy to see that G = I n−1 ⊗ G with G defined in (5.22). Therefore, the EVs of G and G coincide. This shows (5.23).
In conclusion, we discuss the implications of the analysis of this section. There are two effects that contribute to denoising: averaging and synchronization. The former can be interpreted as a manifestation of the law of large numbers: the combined effect of independent stochastic processes acting on individual cells vanishes as the size of the network goes to infinity. It is captured by the first term on the right hand side of (5.5). From (5.5) it follows that averaging depends on the network size, n, and the dissipativity of the local dynamics through λ 1 (A s ) but is independent from the properties of the coupling operator. The second term on the right hand side of (5.5), which captures synchronization, depends on the dissipativity of the coupling operator through λ 1 (B s ), the coupling strength g, and the size of the network. Denoising takes place when the variance of the trajectories of the coupled system can be controlled by n and g. In fact, we can always make the first term in (5.5) arbitrarily small by taking n large enough. To control the second term on the right hand side of (5.5), we can use g provided lim g→∞ λ 1 (B s ) > 0. Thus, one can make the right hand side of (5.5) arbitrarily small by increasing n and g provided λ 1 (B s ) is bounded away from 0 as g → ∞.
Lemma 5.3 identifies two cases important in the context of denoising and synchronization. Depending on the rank of J s , the leading eigenvalue of B s , λ 1 (B s ), is either O(1) or O(δ) for g 1. If the coupling is full rank then λ 1 (B s ) = O(1) and (5.5) shows that the variability of large coupled systems can be effectively controlled by the coupling strength and the network size. However, if the coupling is partial (rank J s < d), as in the model of gap-junctionally coupled network (2.9)-(2.11), (5.5) becomes This shows that in the partial coupling case the mechanism of denoising may fail. Our numerical experiments and the analysis of (2.9)-(2.11) show that denoising is at work despite the fact that the coupling is partial in this model. We are able to observe denosing in (2.9)-(2.11) because of the additional structure of this problem -the proximity of the models of individual cells to the saddle-node bifurcation. This additional feature of the problem affords center manifold reduction, and the coupling of the reduced system is already full rank. This reconciles the results of this section with the previous analysis and our numerical results. This also highlights the importance of the bifurcation structure of the fast subsystem of the bursting cell model (2.1)-(2.3) for the mechanism of denoising. More generally, this discussion implies that while the mechanism of denoising can be already studied at the level of the linearized system, the structure of the nonlinear system is nonetheless important for the realization of denosing in concrete (biophysical) models.

Discussion
In this work, we have analyzed synchronous regimes in electrically coupled networks of bursting capable excitable cells. The individual cells in the network can be tuned to one of the three main activity states: periodic spiking, bursting, or quiescence. Strong electrical coupling synchronizes activity across the network. However, in the presence of noise, the synchronous patterns generated by the network can differ qualitatively from the patterns that the cells comprising the network exhibit in isolation. We have identified two scenarios of such behavior: Scenario A and Scenario B. In the former, the network of irregularly spiking cells generates very regular bursting provided the coupling is sufficiently strong. In the latter, the network formed from irregularly bursting cells is effectively switched to spiking once the coupling strength and the size of the network become sufficiently large.
In constructing these scenarios featuring the disparity of the firing patterns generated by the single cell and coupled models, we used the large deviations type mechanisms to generate irregular patterns in the single cell models and denoising for shaping the patterns of collective behavior in the network. In particular, for irregular spiking in Scenario A, we used noise to perturb the system from slowly evolving stable equilibrium. A closely related mechanism was studied in the context of emerging regular dynamics in randomly perturbed slow-fast systems [17] and self-induced stochastic resonance [57,14]. The irregular bursting in Scenario B was organized by perturbing a system with a stable limit cycle via the mechanism proposed in [31]. In both cases, the slow-fast structure of the single cell models was essential. The emerging firing patterns in both scenarios are very sensitive to the variations in the intensity of noise. Utilizing the ability of electrical coupling to synchronize the activity and to reduce the effects of noise on network dynamics [50,48,49], for the coupled system we achieved effective control of the firing patterns by varying the strength of coupling and the size of the network.
In this paper, we analyzed in detail the mechanism of denoising for coupled systems with (slowly evolving) stable equilibria (Scenario A). We expect that the analysis of Scenario B can be done along the same lines by combining the techniques of Section 5 and local analysis near the synchronous limit cycle of the coupled system as in [48,49]. We will address this problem in the future work. For Scenario A, we considered two cases: one -when the individual subsystems are close the saddle-node bifurcation as in the network of model bursting neurons (2.12) and (2.13), our motivating example, and another without assuming the proximity to the saddle-node bifurcation. For systems near the bifurcation, we adapted two complementary approaches from [55]. The former of which utilizes the gradient structure of the reduced system , which suggests that denoising is more effective in a network on a random graph. The corresponding timeseries are shown in (c) and (d); and the histograms for the number of spikes in one burst -in (e) and (f). The latter shows that the number of spikes in one burst generated by the random network has approximately Gaussian distribution tightly localized around 15. This is in a stark contrast with a much broader distribution corresponding to the symmetric network in (e).
to yield accurate estimates of the first exit times from the basins of the stable equilibria of the single cell and coupled models; while the latter relies on the algebraic graph theory techniques [47], which elucidate the contribution of the network topology to synchronization and denoising. Furthermore, we extended the analysis of denoising in systems of coupled integrate and fire models in [50] to conductance based models with multidimensional phase spaces. We have shown that in this setting the coupling architecture is very important for implementing denoising. In particular, for the common in applications partial coupling case, we have shown that denosing generically does not take place, in contrast to what one might expect from the analysis of the coupled one-dimensional systems [50]. This result highlights the significance of the bifurcation structure of the biophysical models of square wave bursting neurons for the realization of denoising. The proximity to the saddle node bifurcation makes the dynamics of the individual subsystems effectively one-dimensional and the coupling effectively -full rank, thus circumventing the problems for implementing denosing in partially coupled systems identified in Section 5. Therefore, the results of this study extend and complement the existing results characterizing denoising in electrically coupled networks [50,68] and highlight the importance of the bifurcation structure of the local dynamical systems comprising the network for implementing denoising.
In analyzing the coupled system, we paid special attention to the role of the network topology in shaping the network dynamics. We have found that two spectral functions: the algebraic connectivity and the total effective resistance, feature prominently in the quantitative descriptions of synchronization and denoising in electrically coupled networks. The algebraic connectivity of the weighted graph of the network sets the rate of convergence to the synchronization subspace (cf. (4.10)), while the total effective resistance is involved in the estimates of stochastic stability of the synchronous regime (cf. (4.12)). These analytical estimates allow one to use many known results relating the spectra of the graphs and their structural properties (see [9,33] and references therein) to elucidate the contribution of the network structure to its dynamics. As an example of such application, we used the spectral properties of random graphs [33] to show that networks on random graphs feature fast synchronization and robustness of synchrony to noise. They are also more effective for implementing denoising than their symmetric counterparts (see Example 2.3), as shown in numerical experiments in Fig. 7.
Understanding dynamical mechanisms for different patterns of electrical activity in excitable cells and transitions between them has been long recognized as a fundamental problem in mathematical neuroscience. The results of this paper describe the transition from irregular spiking to nearly periodic bursting in networks of electrically coupled cells in the presence of noise. This transition was used in [65,64] to explain why pancreatic β−cells burst in electrically coupled islets of Langerhans but not in isolation. Our results support and extend the previous analysis in [65,64]. We believe that the results and techniques of this work will be useful for understanding dynamics in many other biophysical models of gap-junctionally coupled networks. (2.1)-(2.3).
The first term on the right hand side of (2.1) models the combined effect of sodium and calcium currents, I N a+Ca , the calcium-dependent potassium current, I KCa , delayed rectifier I K and a small leak current, I l I ion = I N a+Ca + I KCa + I K + I l , (A. 1) where The steady state functions used to model the ionic currents above are given by The values of the remaining parameters are summarized in the following table.

Appendix B. Dynamical regimes of the conductance-based model
The geometric theory for singularly perturbed systems implies the existence of the exponentially stable locally invariant manifolds E and L , which are O( ) close to E {(x, y) : y > y sn + δ} and L {(x, y) : y < y hc − δ}, respectively, for arbitrary fixed δ > 0 and sufficiently small > 0 [1,22]. Manifolds E and L are called slow manifolds. For small > 0, the dynamics of (2.4) and (2.5) on the slow manifolds is approximated by L :ẏ = G(y), y < y hc − δ, (B.1) E :ẏ = g(ψ(y), y), y > y sn + δ, and T (y) stands for the period of the limit cycle L y . Depending on the location of the null surface S = {(x, y) ∈ R 3 : g(x, y) = 0}, the slow-fast system (2.4) and (2.5) can be in one of the following states: bursting, spiking, and quiescent (see Fig. 2).
The following conditions on the slow subsystem yield bursting. For some c > 0 independent of , Under these assumptions, for sufficiently small > 0 a typical trajectory of (2.4) and (2.5) consists of the alternating segments closely following L and E and fast transitions between them (see Fig. 2a). The corresponding timeseries is shown in Fig. 2d.
Substituting (SB) with the condition that follows will switch (2.4) and (2.5) 0 to a spiking regime.
(SS) G(y) has a unique simple zero at y = y c ∈ (y sn , y hc ): G(y c ) = 0 and G (y c ) < 0. (B.6) In this case, the Pontryagin-Rodygin theorem [56] yields the existence of an exponentially stable limit cycle L (y c ) of period T (y c ) + O( ) lying in an O( ) neighborhood of L yc , provided (SS) holds and > 0 is sufficiently small (see Fig. 2b,e).
Finally, the slow-fast system (2.4) 0 and (2.5) is said to be in the excitable regime if it has a stable fixed point lying on E (see Fig. 2c): (Q) g(y) := g(ψ(y), y) has a unique simple zero at y = y q ∈ (y sn , y hc ): g(y q ) = 0 and g (y q ) < 0. (B.7)