Critical and Supercritical Spatiotemporal Calcium Dynamics in Beta Cells

A coordinated functioning of beta cells within pancreatic islets is mediated by oscillatory membrane depolarization and subsequent changes in cytoplasmic calcium concentration. While gap junctions allow for intraislet information exchange, beta cells within islets form complex syncytia that are intrinsically nonlinear and highly heterogeneous. To study spatiotemporal calcium dynamics within these syncytia, we make use of computational modeling and confocal high-speed functional multicellular imaging. We show that model predictions are in good agreement with experimental data, especially if a high degree of heterogeneity in the intercellular coupling term is assumed. In particular, during the first few minutes after stimulation, the probability distribution of calcium wave sizes is characterized by a power law, thus indicating critical behavior. After this period, the dynamics changes qualitatively such that the number of global intercellular calcium events increases to the point where the behavior becomes supercritical. To better mimic normal in vivo conditions, we compare the described behavior during supraphysiological non-oscillatory stimulation with the behavior during exposure to a slightly lower and oscillatory glucose challenge. In the case of this protocol, we observe only critical behavior in both experiment and model. Our results indicate that the loss of oscillatory changes, along with the rise in plasma glucose observed in diabetes, could be associated with a switch to supercritical calcium dynamics and loss of beta cell functionality.

Introduction

Homeostasis of energy-rich nutrients in blood has to cope with behavioral and environmental extremes, such as ingestion of a large meal or prolonged fasting (Schmitz et al., 2008). The anabolic hormone insulin promotes postprandial storage of nutrients and tightly controls their consumption interprandially, thus playing a crucial homeostatic role, which becomes disrupted in obesity and diabetes (Kahn et al., 2014). Similarly to many other hormones, insulin concentration in blood oscillates, with a diurnal (meal-related) component, an ultradian component (period of 80–180 min), and a so called high-frequency component (period of 5–15 min), the last being evolutionary conserved in different mammals, such as humans and mice (Nunemaker, 2005; Satin et al., 2015). It has been established beyond doubt that the oscillations in blood are due to pulsatile release of insulin from the pancreas. In contrast, many questions remain to be answered with regard to how exactly this pulsatile release is brought about and regulated, since the beta cells which sense glucose (and other nutrients) and secrete insulin are scattered throughout the exocrine part of the gland in the form of small organs called islets of Langerhans, of which there are about a thousand in the mouse and about a million in the human pancreas, and each of which contains from a few to a couple of thousand beta cells, together with other endocrine and mesenchymal cells (Dolenšek et al., 2015).

The stimulus-secretion coupling process in beta cells involves entry of glucose (and other nutrients) into the cell and its metabolism to ATP, which in turn decreases the open probability of ATP-sensitive potassium channels, leading to depolarization of plasma membrane, opening of voltage-sensitive calcium channels, a rise in the cytosolic calcium concentration ([Ca²⁺]_c), and triggering exocytosis (Ashcroft and Rorsman, 2013). In addition to this canonical triggering pathway, there are probably an additional amplifying calcium-dependent (Henquin, 2011) and even a calcium-independent signaling pathway (Aizawa et al., 1998).

Individual, uncoupled beta cells display oscillations of membrane potential and [Ca²⁺]_c with a frequency close to the fastest of the abovementioned, but with a large degree of heterogeneity (Tengholm and Gylfe, 2009; Satin et al., 2015). Within islets of Langerhans, a strong intercellular coupling force in the form of intercellular gap junctions consisting of the connexin 36 protein and possibly other modes of cell-cell communication overcome the heterogeneity of individual beta cell oscillators (Bavamian et al., 2007; Konstantinova et al., 2007; Tengholm and Gylfe, 2009; Benninger et al., 2011; Rutter et al., 2017; Skelin Klemen et al., 2017). In whole islets, metabolism, membrane potential, [Ca²⁺]_c, and secretion of insulin oscillate at a frequency close to the abovementioned 5–15 min due to coupling, however islets also display so called fast oscillations (period of 1–15 s) in the form of bursts of membrane depolarizations (Dean and Matthews, 1970), with accompanying oscillations in [Ca²⁺]_c (Gilon et al., 1992; Dolenšek et al., 2013) and insulin secretion (Bergsten, 1995). This fast oscillatory component is responsive to changes in glucose concentration (see below), synchronized in different cells of an individual islet by means of membrane potential and [Ca²⁺] waves (Dolenšek et al., 2013), but does not seem to be entrained into a common rhythm among different islets in vivo (Valdeolmillos et al., 1996). In contrast, the slow oscillations of individual islets are not responsive to changes in glucose concentration, but probably are entrained into a common rhythm in vivo, yielding the 5–15 min oscillations observed in blood. It is believed that the 5–15 min oscillations are due to oscillations in metabolism and the 1–15 s oscillations are due to a feedback between calcium and potassium channels (Satin et al., 2015). This view is incorporated into a mathematical model called the dual-oscillator model, which in addition to the mixed pattern of carrying slow metabolic and superimposed fast electrical oscillations is also able to account for in vitro observations of only the fast or the slow component in absence of the other (Bertram et al., 2007).

Metabolic differences between individual islets (Nunemaker et al., 2005, 2009) can be overcome by weak coupling via an intrapancreatic neural network (Fendler et al., 2009), by negative feedback from the liver (Pedersen et al., 2005; Dhumpa et al., 2014), or both (Satin et al., 2015). According to the recent metronome model, glucose-responsive fast oscillations of individual islets determine the amplitude or pulse mass of the largely stable 5–15 min insulin oscillations (Satin et al., 2015).

Theoretically, an individual islet can respond to an increase in glucose concentration by recruiting more cells into a functional state, by enhancing the response of active cells, or both. Previous experiments have suggested that within a narrow range of glucose concentrations above the threshold concentration, recruitment rapidly saturates and that beyond that, all beta cells within an islet are active all the time, with synchronous membrane potential and [Ca²⁺]_c oscillations that increase in plateau fraction with increasing glucose concentrations (Henquin et al., 1982; Henquin, 1987; Valdeolmillos et al., 1989; Santos et al., 1991; Gilon and Henquin, 1995; Jonkers et al., 1999; Jonkers and Henquin, 2001). In sum, according to this view the pulse mass is more importantly determined by enhancing the responses of individual cells than by recruiting new cells (Jonkers et al., 1999; Jonkers and Henquin, 2001). The main shortcoming of previous studies aimed at quantitating the role of recruitment and enhancement is the fact that clusters of beta cells were used instead of islets to ensure spatial resolution at the level of individual cells, and that when whole islets were used, resolution at the level of individual cells was not achieved. Additionally, beta cells have traditionally been stimulated by elevating glucose to supraphysiological concentrations and in a constant, i.e., non-oscillatory manner, not to concentrations slightly above the threshold and in an oscillatory manner, as is probably the case in vivo.

Finally, considering the formidable complexity of the mechanism supporting pulsatile insulin release, the teleological question seems appropriate, as to what evolutionary advantage is conferred by pulsatile insulin release. It has been suggested by modeling that pulsatile insulin secretion may be beneficial for the secretory capacity of beta cells in the sense that it allows the readily releasable pool of insulin granules enough time to refill during the resting periods between periods of activity (Pedersen and Sherman, 2009). Additionally, it has been proposed that pulsatile insulin may be important for physiological autocrine effects in islets and expansion of beta cell mass (Tengholm and Gylfe, 2009). On the other hand, experimental evidence points to a greater efficiency of pulsatile insulin on target tissues and continuous insulin leads to internalization, down-regulation of insulin receptors, and post-receptor signaling defects (Pørksen, 2002; Pørksen et al., 2002; Satin et al., 2015). Importantly, disrupted insulin pulsatility has been observed in prediabetes, type 2 diabetes mellitus (T2DM), and even in normoglycaemic relatives of patients with T2DM (Lang et al., 1981; O'Rahilly et al., 1988; Bingley et al., 1992). In T2DM beta cells are also unable to respond to entrainment by imposed oscillations in glucose, indicating a fundamental loss of entrainability (Hollingdal et al., 2000). Furthermore, it has been demonstrated that beta cell coupling is a target of diabetogenic insults, such as lipotoxicity (Hodson et al., 2013), glucotoxicity (Haefliger et al., 2013), and cytokines (Farnsworth et al., 2016; Johnston et al., 2016). Likewise, islets exposed to such insults and islets of Cx36 knockout mice display disrupted calcium waves and synchronization of oscillations (Benninger et al., 2008, 2014). Finally, in Cx36 knockout mice the in vivo observed insulin response and pattern of insulin oscillations are disrupted in a manner resembling typical changes in T2DM (Head et al., 2012).

Even though some mathematical models have addressed the collective behavior of coupled beta cells (Smolen et al., 1993; Zimliki et al., 2004) and the propagation of [Ca²⁺] waves (Nittala et al., 2007; Benninger et al., 2008), the development of multicellular computational models is just recently increasingly gaining attention. For the most part, this has been initiated by the growing evidence showing that cell–cell interactions via gap junctions are a prerequisite for proper hormone secretion (Charollais et al., 2000) and that impaired intercellular interactions disrupt normal oscillatory patterns of insulin secretion in a way similar to what occurs in diabetes (Head et al., 2012). Moreover, recent advances in imaging techniques along with the integration of complex system approaches in islet research revealed a complex functional organization of beta cells (Hodson et al., 2013; Stožer et al., 2013b; Cherubini et al., 2015; Gosak et al., 2015, in press; Markovič et al., 2015) that might be predominantly a consequence of cell-to-cell variability and a heterogeneous nature of intercellular interactions (Goel and Mehta, 2013; Barua and Goel, 2016; Cappon and Pedersen, 2016). Detailed analyses of pancreatic [Ca²⁺] waves have shown that the waves originate from specific yet rather randomly distributed sub-regions with elevated excitability (Benninger et al., 2014) or lower metabolic rates (Westacott et al., 2017), thereby giving emphasis also to the spatial aspect of beta cell heterogeneity and the sub-compartmental organization of islets (Markovič et al., 2015). Furthermore, percolating network models have been utilized with the aim to provide a phenomenological insight into the interplay between the coupling architecture and the beta cell activity in health and disease (Benninger et al., 2008; Hraha et al., 2014a,b; Stamper et al., 2014). Most importantly, by these means the existence of a critical behavior that reflects a phase transition between globally active and inactive states was predicted and confirmed both computationally and experimentally (Hraha et al., 2014b). Such a critical transition was suggested to be a general regulatory mechanism that islets utilize in order to leverage cellular heterogeneity, and is characteristic also for a variety of other complex real-life systems (Trefois et al., 2015).

In this vein, many natural systems were found to operate naturally near a critical point (Bak, 1996; Marković and Gros, 2014). The emergent dynamics in such systems is usually associated with the concept of self-organized criticality (SOC), which embraces a power-law distribution of systems' observables. SOC arises in complex systems that are far from equilibrium as a result of interactions of components and is increasingly gaining attention in the context of organizing principles of biological systems. For example, durations of brief awakenings during sleep exhibit a power-law distribution, indicating a scale-invariant dynamic that is typical for systems undergoing SOC (Lo et al., 2004, 2013; Allegrini et al., 2015). The power-law distribution is also characteristic for other human body dynamics, e.g., human gait and human heartbeat; however, the scale invariant feature characterized by 1/f scaling is only a hallmark of SOC. Ivanov et al. (Ivanov et al., 2009) showed that while both gait interstride interval and cardiac interbeat interval time series have comparable 1/f scaling, they are governed by different mechanisms and lead to different levels of complexity, characterized by monofractal and multifractal properties, respectively. Very recently, as a further example of SOC in human body, clinical evidence has been provided for self-organization of blood pressure regulation (Fortrat and Gharib, 2016). At the single cell level, hallmarks of SOC have been observed in the spatiotemporal organization of [Ca²⁺] waves in individual cardiac myocytes (Nivala et al., 2012) and oocytes (Lopez et al., 2012). However, the greatest progress in this framework has been done in the field of neuroscience where fingerprints of SOC have been identified at different levels of organization, ranging from interacting arrays of neurons or astrocytes to the entire brain (Jung et al., 1998; Beggs and Plenz, 2003; Plenz and Thiagarajan, 2007; Hesse and Gross, 2014). On different scales of observation, patterns of neuronal electrical activity show a high degree of diversity, characterized by a scale free distribution of event (i.e., neuronal avalanches) sizes. Theoretical and experimental work indicates that this form of activity reflects a critical state between a random and an ordered dynamical regime, which is believed to lead to optimal operational abilities (Beggs and Plenz, 2003; Kinouchi and Copelli, 2006). A balanced interplay between network dynamics and topology as well as an activity-dependent adaptability of neuronal networks were identified as leading neurobiological determinants that ensure a robust critical behavior (Levina et al., 2007; Rubinov et al., 2011; Hutt et al., 2014). However, in vitro experiments on cortical assemblies have shown that the neuronal activity does not necessarily fall into a critical regime, but can also show subcritical or supercritical states, e.g., during development (Tetzlaff et al., 2010), which has latter been associated with changes in the neuronal network connectivity patterns (Tetzlaff et al., 2010; Massobrio et al., 2015b). Moreover, an excess of large neuronal avalanches, as is characteristic for supercritical dynamical states, has been associated with pathological activity (Meisel et al., 2012; Massobrio et al., 2015a). It has been speculated that deviations from critical behavior reflect the abnormal synchronized firing of neurons involved in the epileptic process (Lehnertz et al., 2009), thereby substantiating the hypothesis that the normal, healthy brain resides in a critical or even slightly subcritical state (Priesemann et al., 2014; Massobrio et al., 2015a). However, in contrast to inanimate matter, for which the SOC principles are well understood, in living systems these mechanisms are still under debate and are conjectured to be a result of different evolutionary or adaptive processes, with structural disorder, order-parameter feedback, and extended criticality as possible explanations (Lovecchio et al., 2012; Moretti and Muñoz, 2013).

Despite growing evidence that critical-like dynamics and power-law scaling imply an efficient design of biological systems, only a few studies have investigated the presence of these features in other multicellular physiological systems. At least in part this lack can be explained by demanding experimental techniques that are required to precisely and noninvasively assess the function of a large number of cells simultaneously and over long periods of time. Previously, it has been demonstrated that the behavior of islets of Langerhans is not as deterministic as once thought (Dolenšek et al., 2013; Hodson et al., 2013; Stožer et al., 2013a,b; Benninger et al., 2014; Hraha et al., 2014a,b; Markovič et al., 2015; Rutter and Hodson, 2015). Thus, in the present study we examined whether fingerprints of SOC can be found in the spatiotemporal pattern of fast [Ca²⁺]_c dynamics in interconnected beta cells from islets of Langerhans. We first constructed a computational model of a network of heterogeneous and heterogeneously coupled beta cells and analyzed the statistical organization of intercellular [Ca²⁺]_c wave sizes. We simulated beta cell behavior after applying a constant stimulatory concentration of glucose, as well as under an oscillatory stimulation with a slightly lower average glucose concentration, in order to more closely mimic the physiological conditions with oscillating blood glucose and insulin levels. Finally, we compared the model predictions with experimental data obtained by means of confocal functional multicellular calcium imaging of [Ca²⁺]_c changes evoked in beta cells in acute tissue slices subjected to the same stimulation protocols as in simulations.

Materials and Methods

Single Cell Model

The dynamics of each beta cell is governed by the mathematical model proposed by Bertram et al. (2007). The model combines mitochondrial metabolism, glycolysis and plasma membrane electrical activity and Ca²⁺ activity in the cytosol and in the endoplasmic reticulum (Pedersen and Sherman, 2009). Interconnecting these three compartments represents a general and widely used theoretical framework that provides a firm description of the complex oscillatory patterns, such as compound bursting patterns, observed in experiments. A detailed description of the model equations along with parameter values is given in the Text S1.

Simulation of Glucose Stimulation

To simulate the stimulation with glucose, we increased the glucokinase reaction rate parameter J_GK (Bertram et al., 2007). In case of constant stimulation, i.e., switch from 6 to 8 mM glucose, the parameter was increased from J_{GK, L} = 0.04 μMms⁻¹ to J_{GK, H} = 0.38 μMms⁻¹. When we simulated an oscillatory stimulation protocol with glucose, glucokinase reaction rate was smoothly varied between these two values as follows:

$\begin{array}{cc}\begin{array}{l}{J}_{\text{GK}}(t)=J_{\text{GK},\text{L}}+\frac{1}{\pi }(J_{\text{GK},\text{H}}-J_{\text{GK},\text{L}})\\ (\frac{\pi }{2}-\frac{{\text{Tan}}^{-1}\left(k_w\left(\left(\left(t+T_{\text{lag}}\right)-T_{P1}\right)\text{ mod }{T}_L -T_{P2}\right)\right)}{2}\\ +\frac{{\text{Tan}}^{-1}\left(k_w\left(\left(\left(t+T_{\text{lag}}\right)\right)\text{ mod }{T}_L -T_{P1}\right)\right)}{2})\end{array}& (1)\end{array}$

where k_w is a pulse-smoothing parameter set to 0.0005. The period of the wave is given by T_L = 600, 000 ms. Durations of the low and high glucose concentration phases are given with the parameters T_P1 and T_P2. Both parameters were set to 300, 000 ms. Lastly, T_lag is used to adjust the onset of the first stimulation. In this manner an oscillatory stimulation protocol was stimulated with a period of 10 min with 5 min intervals of elevated glucose conditions. The course of the periodic variations in glucokinase reaction rate is shown in Figure S1.

Heterogeneity of Beta Cells

The pancreatic beta cells exhibit a high degree of heterogeneity, which manifests itself in cell-to-cell variability in their sizes, membrane capacitance, channel densities and conductances, in the rates in glucose-induced insulin synthesis and release, cellular thresholds for glucose utilization and oxidation, etc. (Pipeleers et al., 1994; Benninger and Piston, 2014). In our model we introduce heterogeneity of beta cells by means of a random distribution of some model parameters. Albeit cell-to-cell variability would result in diversity of various model parameters, in the present study we consider only heterogeneity of some crucial parameters in different parts of the machinery that governs the beta cells behavior. In particular, heterogeneity is introduced in the glyceraldehyde 3-P dehydrogenase (GPDH) reaction rate parameter k_{GPDH, i}, glycolytic flux affecting PDH activity J_{GPDHbas, i}, maximal PFK reaction rate V_{max, i}, the membrane capacitance C_i, and in g_{k(ATP), i} denoting the conductance of ATP-sensitive K⁺ channels. All these parameters were assumed to follow a normal distribution with a relative standard deviation of 30 % with a cut-off of 60 %.

Network of Beta Cells

We model the intercellular coupling between beta cells as random geometric graph. Initially N=150 cells are arranged randomly in a unit square with a prescribed minimal possible distance (0.04) to ensure a more homogeneous spatial distribution. Connections among the cells represent intercellular communication by means of electrical coupllig, defined as:

$\begin{array}{cc}{I}_{\text{cpl},i}=g_i{\displaystyle \sum _{j\ne i}^N}{d}_{ij}(V_j-V_i).& (2)\end{array}$

This coupling term I_{cpl, i} is added to the equation describing the dynamics of the membrane potential of the i-th cell (see Equation s21) in Text S1). The network structure is stored in the coupling matrix d. Its ij-th element d_ij is set to 1 if the i-th cells is connected with the j-th cell, whilst otherwise d_ij = 0. In particular, an identical radius range for all cells was chosen as $r_n=\sqrt{\langle k\rangle /(N\pi )}$, where 〈k〉 = 6 signifies the average number of connections per cell. Two cells are then connected if they fall within each other's range.

In Equation (M2) g_i stands for the electrical coupling coefficient. Previous studies have indicated a high degree of heterogeneity in the gap junctional conductances between beta cells that exceeds a normal distribution. To implement this feature into our simulations the values of g_i were distributed in accordance to an exponential distribution as:

$\begin{array}{cc}{g}_i=-\ln \left(\text{rand}\left(0,1\right)\right)\langle g_i\rangle ,& (3)\end{array}$

where 〈g_i〉 symbolizes the mean value of the distribution (〈g_i〉 = 200 pS) and rand(0, 1) is a uniformly distributed random number in the interval (0,1). The resulting network models are quite homogeneous and without unconnected components. The intercellular coupling, on the other hand, is rather heterogeneous (see Figure S2).

Confocal Calcium Imaging in Pancreatic Tissue

The preparation of the pancreatic tissue slices and the experimental protocols to monitor changes in intracellular calcium concentration were described in detail before (Speier and Rupnik, 2003; Dolenšek et al., 2013; Stožer et al., 2013b). Briefly, NMRI mice (age 10–20 weeks, both males and females) were sacrificed in order to have their abdomens exposed. 1.9 % low-melting point agarose was injected through the proximal bile duct into the ductal system of pancreas. After transferring the agarose blocks containing isolated pancreas tissue onto the wet ice, the hardened agarose allowed for cutting the soft pancreas tissue into 140 μm thick slices. Slices were incubated in the dye-loading solution containing 6 μM Oregon Green 488 BAPTA-1 AM (OGB-1, Invitrogen, Eugene, Oregon, USA), 0.03% Pluronic F-127 (w/v), and 0.12% dimethylsulphoxide (DMSO, v/v) at room temperature. Recordings of the calcium concentration changes were performed on a Leica TCS SP5 AOBS Tandem II upright confocal system using a 20x Leica HCX APO L water immersion objective (NA 1.0). The excitation wavelength was set to 488 nm and the emission collected in the range 500–650 nm. The sampling rate was 10 Hz. Regions of interest were selected based on cell morphology and exported as time series for off-line analysis. This study was carried out in accordance with the national recommendations. The protocol was approved by the Slovenian Ministry of Agriculture, Forestry and Food U34401-30/2016/U94-01.

Processing of the Recorded Time Series

Experimentally and numerically obtained [Ca²⁺] traces were initially leveled and smoothed to ensure a consistent and accurate binarization procedure in a two-step signal processing protocol. To extract the fast component of [Ca²⁺] oscillations and to exclude artifacts (i.e., photobleaching) a Butterworth band pass filter (Yao et al., 2012) has been applied to the data sets. The order of the filter has been set to 2 to achieve a steeper frequency cut-off and not to destabilize the filter (resonance disaster). The upper F_HIGH and lower F_LOW cut-off frequencies have been determined by visually inspecting all the filtered time traces of a given sample. What followed was the smoothing of the time traces, by means of a sliding window algorithm (Yaroslavsky et al., 2001). The temporal width of the sliding window Δt_SW has been adjusted to the sampling rate in order to avoid over smoothing of the data (Δt_SW = ±2 frames). The solely prepared data has then been used for the binarization of individual beta cell Ca²⁺ activity. The onset and ending of an activation pulse have been determined by computing the first time derivate of individual time traces $(\stackrel{.}{C{a}_i^{2+}})$ and the corresponding standard deviations of the time traces std $(C{a}_i^{2+})$ as well as time derivates std $(\stackrel{.}{C{a}_i^{2+}})$. The combined three information's were then used for the binarization. Whenever a local maxima in $(\stackrel{.}{C{a}_i^{2+}})(t)$ is >1.5$*\text{std }(C{a}_i^{2+})$ and the corresponding local maxima in $C{a}_i^{2+}(t)$ is >1.5$*\text{std }(C{a}_i^{2+})$, we treat the time t as the onset of the n-th activation time t_{START, i}(n). Between the n-th local maxima t_{START, i}(n) and its first successor t_{START, i}(n+1), we search for local minima in $(\stackrel{.}{C{a}_i^{2+}})$. The discrete time, corresponding to the minima is then set as the end time of the n-th activation pulse t_{END, i}(n). Values of the binary time trace matrix $C{a}_{BIN,i}^{2+}$ which lie within the interval t_{START, i}(n) and t_{END, i}(n) are set to 1 and other values are set to 0. The process is schematically shown in Figures 1C,D.

FIGURE 1

Figure 1. The procedure for assessing spatiotemporal [Ca²⁺]_c dynamics in beta cells. (A,B): We used computational simulation (A, red) and functional multicellular calcium imaging (B, green) to determine beta cell spatial coordinates within networks of beta cells and the individual [Ca²⁺]_c activity either in experiments or in simulations (C,D). Oscillations were binarized (D) and further processed to extract individual spatio-temporal clusters of [Ca²⁺]_c activity (E,F).

Space-Time Cluster Analysis

In order to identify clusters of active cells we implemented the space-time cluster (STC) analysis, similar to the one proposed by Jung (Jung, 1997; Jung et al., 1998). Individual time frames were stacked together to obtain a large space-time cube. Each frame embeds the systems spatial organization, i.e., positions of cells, and the corresponding binary states of the included cells (active or inactive). The time interval between two frames is given by Δt_F. We start the algorithm by creating a cube around every active cell with a spatial side length and a temporal side length Δt_TSL. The spatial side length Δs for a given dataset was computed as the average distance between the 6 nearest cells. Typical values of Δs for experimental data were between 20 μm and 25 μm and for the computational model around 0.11. The temporal side length Δt_TSL was set to 0.5 s. A STC was then defined as a group of cells, for which the created cubes overlap in a time forward direction. In other words, a cluster results for an isolated intercellular [Ca²⁺] wave that propagates between nearby cells, and its size p reflects the number of involved cells and the durations of their [Ca²⁺] pulses in this given event. If two such clusters, i.e., [Ca²⁺] waves, collide, the incoming cluster is joined to the cluster it collides with. The whole procedure from time series preparation to space-time cluster analysis is schematically presented in Figure 1. It should be noted that the size p of an individual wave (indicated by the color) actually reflects the volume of a given space-time cluster.

To characterize the spatiotemporal characteristics of intercellular [Ca²⁺] waves, we measured the distribution of cluster sizes. More specifically, we divided the range of sizes' values in a series of non-overlapping intervals with size p and counted how many values N(p) fell into each interval. For both simulated and measured dynamics, we calculated the distribution separately for the activation and for the plateau phase in case of constant stimulation, and over the whole interval (simulated or measured) in case of periodic stimulation. For the representation of the computational results we pooled the data from three (constant stimulation) and four (periodic stimulation) independent simulations runs and for the experimental results we merged the data from three (constant stimulation) and four (periodic stimulation) different tissue slices originating from three different mice. In this manner, the number N(p) reflects the average number of detected waves with size p out of three/four settings. Since the number of cells in the pancreatic slices in different experimental recordings was different, we normalized all slices with respect to the largest detected cluster, i.e., the largest cluster has a size p = 1. Finally, the data (experimental and computational) were fitted with a power-law function to qualitatively evaluate if the activity patterns can be treated as critical or supercritical.

Results

Previous experimental investigations and modeling endeavors have shown both non-trivial temporal activity patterns of beta cell populations and a very complex spatiotemporal organization at the multicellular level (Benninger and Piston, 2014). More specifically, despite gap junctional coupling that fosters the intrinsically heterogenoeus cells to operate in synchrony, a very heterogeneous activity was observed especially in the activation phase when the cells respond to stimulation (Stožer et al., 2013a). Moreover, recent studies also indicate that the beta cell syncytium is functionally organized in a quite complex manner (Stožer et al., 2013b; Markovič et al., 2015; Rutter and Hodson, 2015; Cappon and Pedersen, 2016; Johnston et al., 2016), which provides a basis for non-trivial activity patterns that are less synchronized than once thought. With the aim to explore the spatiotemporal behavior of beta cells, we first built a computational model of interconnected beta cells. The dynamics of individual beta cells is driven by the comprehensive mathematical model proposed by Bertram et al. (2007) that combines glycolysis, and mitochondrial metabolism with plasma membrane electrical activity and [Ca²⁺]_i activity. In accordance with previous findings that suggested an extensive heterogeneity among beta cells (Pérez-Armendariz et al., 1991), we introduced variability of some model parameters that govern the beta cell behavior. The cell-to-cell electrical coupling between beta cells was modeleld by means of interactions formed by the random geometric network model, whereby the electrical coupling was assumed to be highly heterogeneous (see Materials and Methods for details). Namely, previous reports and our own observations have shown a high degree of heterogeneity in the gap junctional conductances between beta cells, which do not follow a normal distribution and the values span over a broad interval with a mean value around 200 nS (Pérez-Armendariz et al., 1991). In the continuation, we compare the simulated behavior of a network of interconnected beta cells with experimentally monitored beta cell activity assessed by means of confocal functional multicellular calcium imaging (fMCI).

As already mentioned above, in both simulations and experiments we used two different stimulation protocols: a constant and an oscillatory stimulation. In both cases a typical response to stimulation consisted of a delayed elevation of [Ca²⁺]_i with superimposed oscillatory changes of [Ca²⁺]_i. These oscillations were detected in the form of coordinated [Ca²⁺]_i waves spreading across cells. [Ca²⁺]_i traces either simulations or experiments were then processed and binarized in order to provide ground for characterization of [Ca²⁺]_i signal propagation. To this purpose, we detected and labeled individual intercellular [Ca²⁺]_i waves by means of space-time cluster analysis. The whole procedure is schematically presented in Figure 1. Finally, we quantified the intercellular [Ca²⁺]_i activity by measuring the spatiotemporal size of waves and looked for the distribution of the number of [Ca²⁺]_i waves of a given size [N(p)] as a function of size (p).

Constant Stimulation: Computational Results

Stimulation with glucose was modeled by simultaneously increasing the glucokinase rate parameter J_GK from 0.04 to 0.38 at 100 s after the beginning of simulation in all cells in the network. Stimulation initiated an activation phase with a characteristic progressive recruitment of active beta cells with delays spanning 2–10 min, which is in accordance with previous experimental findings (Stožer et al., 2013a). In this stage, the spatiotemporal [Ca²⁺]_i activity encompassed smaller clusters of active cells. The activation phase was followed by a plateau phase with characteristic [Ca²⁺]_i waves that often encompassed the majority of cells, and, in contrast to the activation phase, displayed a temporally quite ordered and deterministic-like sequence. These observations are illustrated in Figure 2A, whereas a more precise insight into the temporal evolution of the spatiotemporal behavior of the model after stimulation can be observed in Video S1. Next, we quantitatively characterized the [Ca²⁺]_i waves, separately for the activation and the plateau phase. To this purpose we first extracted and labeled individual waves, as presented in Figures 2B,C. We then calculated relative wave sizes p during the activation (Figure 2D) and the plateau phase (Figure 2E). To ensure a better statistical accuracy, we calculated the distributions as the average of three independent simulation runs. The distribution N(p) was found to obey a power law in the activation phase, whereas an excess of larger waves was detected in the plateau phase. In sum, this indicates that the persistent stimulation initiates a critical state during the activation phase, which then switches to a supercritical state during the plateau phase.

FIGURE 2

Figure 2. Simulation of constant stimulation: spatiotemporal organization of intercellular [Ca²⁺]_i waves in islets. (A) Typical computed [Ca²⁺]_c responses of four different beta cells after switching to stimulatory conditions (upper panel, gray area indicates stimulatory conditions) and binarization of the computed oscillations of all cells (lower panel). Extracted individual waves as denoted by different colors in space-time graphs during the activation phase (B) and in the plateau phase (C). The purple dots in the x-y plane denote the coordinates of individual cells. [Ca²⁺]_i wave size distribution of the computed data during the activation (D) and plateau phase (E). The gray line indicates a power-law fit with slopes −1.41 and −1.58 for the activation and plateau phase, respectively.

Oscillatory Stimulation: Computational Results

Using the same computational model, we simulated the beta cell activity during an oscillatory stimulation regime by smoothly changing the glucokinase reaction rate with a period of 10 min, as described in Materials and Methods section. By this means we mimicked oscillations of glucose concentration observed in humans and mice in vivo (Head et al., 2012; Satin et al., 2015) in a manner that is compatible with physical limitations of our experimental perifusion setup and is also comparable with other recent reports employing more specialized microperifusion chambers (Pedersen et al., 2013; Dhumpa et al., 2015; Sun et al., 2015). The results are shown in Figure 3. The activity of cells displayed a phase shift of a few minutes with regard to stimulation. Except for the first stimulation pulse, the majority of cells were active at least once in each pulse. Notably, even in the low-stimulatory phases between the pulses (white areas in Figure 3A), cellular activity was observed, but was typically limited to [Ca²⁺]_c increases in smaller subgroups of cells. A more detailed insight into the spatiotemporal beta cell dynamics can be seen in Video S2. To evaluate the observed behavior that appeared to be qualitatively different from the one observed in the plateau phase during a constant stimulation, we quantified the intercellular [Ca²⁺]_i activity. Figure 3B depicts [Ca²⁺]_i waves during the 4th stimulation pulse and the colors denote individual clusters. The average distribution N(p) of the wave sizes p was calculated on the basis of four independent simulation runs and was found to obey a power law, indicating critical behavior of the system during oscillatory stimulation.

FIGURE 3

Figure 3. Simulation of oscillatory stimulation: spatiotemporal organization of intercellular [Ca²⁺]_c waves in islets. (A) Computed [Ca²⁺]_c responses of four typical cells during oscillatory glucose stimulation (upper panel) and binarized dynamics of all cells in the network (lower panel). The gray areas denote the stimulatory pulses realized by periodic increases of glucokinase reaction rates. (B) Space-time clusters of [Ca²⁺]_i activity during the four 5-min glucose stimulations, the colors denote different [Ca²⁺]_i waves. (C) The distributions N(p) of spatiotemporal [Ca²⁺]_i wave sizes p. The gray line indicates a power-law fit with a slope of −1.48.

Constant Stimulation: Experimental Results

We stimulated the beta cells within an islet of Langerhans with 8 mM glucose. As in simulations, the cells responded to stimulation with an initial activation phase that was followed by a plateau phase. Characteristic for the activation phase was a transient increase in [Ca²⁺]_c and occurrence of fast oscillations, but beta cells were being recruited only gradually during this phase (Figure 4A, 300 s < t < 800 s). During the plateau phase that followed the activation phase, all cells were active and displayed repeated and more regular oscillations (Figure 4A, t > 800 s). Video S3 shows the temporal evolution of the measured and binarized [Ca²⁺]_c activity from the onset of stimulation.

FIGURE 4

Figure 4. Experimental constant stimulation: spatiotemporal organization of intercellular [Ca²⁺]_c waves in islets. (A) [Ca²⁺]_i responses of four typical cells within an islet to stimulation with 8 mM glucose (upper panel). Binarized oscillations for all cells within the same islet (lower panel). The gray shaded area indicates the switch from 6 to 8 mM glucose. Space-time clusters of [Ca²⁺]_i activity in the activation phase (B; 300 s < t < 800 s) and in the plateau phase (C; t > 800), the colors denote different [Ca²⁺]_i waves. The distributions N(p) of relative spatiotemporal [Ca²⁺]_i wave sizes p for the activation phase (D) and in the plateau regime (E). The gray line indicates a power-law fit with slopes −1.82 and −1.92 for the activation and plateau phase, respectively.

Next, we separately determined [Ca²⁺]_c waves during the activation and the plateau phase. Figures 4B,C show such waves during a part of the activation and the plateau phase, respectively, individual waves and their respective sizes are denoted with different colors. Distribution N(p) of relative wave sizes p was determined and plotted on a log-log scale for both the activation (Figure 4D) and the plateau phase (Figure 4E). For the former, N(p) follows a linear relationship on the log-log scale, whereas for the latter a larger proportion of larger and more global waves could be detected. The distribution signifies the average of three different slices subjected to the same protocol. It should be noted that all of them showed a conceptually very similar behavior, i.e., an activation phase that was followed by a plateau phase with dominating global events, even though the average firing rate and durations of oscillations might be different in different islets. Since the number of cells in each slice was different (111, 125, and 98), we normalized the [Ca²⁺]_c wave sizes with respect to the largest one detected in a given slice. Apparently, under constant stimulation, the nature of the spatiotemporal organization changed from a critical regime in the activation phase to a supercritical regime in the plateau phase. These results overlap well with model predictions and, most importantly, substantiate our choice of the model particularities, i.e., a high degree of heterogeneity in properties of individual cells as well as in intercellular coupling.

Oscillatory Stimulation: Experimental Results

Finally, we checked if the model correctly predicted the critical spatiotemporal organization of [Ca²⁺]_i waves under the oscillatory stimulation protocol. To this purpose, we varied the glucose concentration periodically from 6 to 8 mM with a period of 10 min, whereby each concentration was applied for 5 min. The period of the resulting oscillations lies at the upper end and the amplitude is roughly an order of magnitude larger than what has been described in vivo. However, shorter-period and smaller-amplitude oscillatory stimulations are practically impossible to achieve with our perifusion system and this seems to be a challenge also for dedicated microperifusion systems (Roper Industries, US). Results are presented in Figure 5 and Video S4. As in the simulation, a smaller proportion of cells were active already during the first stimulation pulse, however most cells were active during subsequent stimulations. Moreover, different sizes of [Ca²⁺]_i waves can be observed and again, a phase shift between glucose stimulation and maximal beta cell activity was observed. We used the same methodological paradigm as with the persistent stimulation to determine [Ca²⁺]_i waves. Figure 5B shows spreading of the [Ca²⁺]_i waves during the second stimulus, each wave being denoted with a different color. For this stimulatory regime, the distribution N(p) of the relative sizes p of the [Ca²⁺]_i waves follows a power law (Figure 5C), as it was predicted by the model. The distribution was calculated on the basis of four different slices subjected to the same protocol and qualitatively identical behavior was observed in all recordings. Since the number of cells in each slice was different (33, 68, 76, and 35), we normalized the [Ca²⁺]_i wave sizes with respect to the largest one detected in a given slice. Apparently, oscillations in glucose concentration that mimicked in vivo conditions, trapped the beta cell population in a critical regime.

FIGURE 5

Figure 5. Experimental oscillatory stimulation: spatiotemporal organization of intercellular [Ca²⁺]_i waves in islets. (A) Typical [Ca²⁺]_i responses in four typical cells to an oscillatory stimulation with 8 mM glucose (upper panel). Oscillations are binarized and depicted for all cells within the same islet (lower panel). The gray and white areas denote glucose concentrations of 6 and 8 mM, respectively. (B) Space-time clusters of [Ca²⁺]_i activity during the second 5 min glucose stimulations, the colors denote different [Ca²⁺]_i waves. (C) The distributions N(p) of relative spatiotemporal [Ca²⁺]_i wave sizes p. The gray line indicates the power-law fit with a slope of −1.78. Note that during glucose nadirs, beta cells remain active.

Discussion

Beta cells respond to stimulation by glucose with an oscillatory activity pattern (Farnsworth and Benninger, 2014). Electrical coupling between beta cells provides the necessary but probably not the only substrate for their coordinated activity and regulated hormone release (Rutter and Hodson, 2013; Benninger and Piston, 2014; Skelin Klemen et al., 2017). However, beta cells within an islet are not completely synchronized. First, measuring a number of different parameters of beta cell function, it has been shown that not all beta cells within an islet respond to stimulation simultaneously, but are progressively recruited into an active state (Schuit et al., 1988; Hiriart and Ramirez-Medeles, 1991; Kiekens et al., 1992; Jonkers and Henquin, 2001; Zarkovic, 2004; Stožer et al., 2013a). Second, the membrane potential and [Ca²⁺] changes spread over the islet in a wave-like manner with a finite speed and do not necessarily always encompass all beta cells in a given islet (Benninger et al., 2008; Dolenšek et al., 2013; Stožer et al., 2013a; Cappon and Pedersen, 2016). These two features reflect the intrinsically heterogeneous nature of beta cells (Benninger et al., 2014; Cappon and Pedersen, 2016), which is retained despite intercellular coupling, most probably for the sake of an at least partly selective and gradual regulation of their function.

To get a more detailed insight into the spatiotemporal organization and the subsequent physiological function of this complex multicellular system, we made use of computational modeling approaches in combination with advanced high spatial and temporal resolution confocal imaging. We developed a multicellular computational model of interconnected beta cells that was based on the theoretical framework of Bertram et al. (2007). By incorporating known particularities, such as cell-to-cell variability in glucose metabolism and conductance of ATP-sensitive K⁺ channels, and a high degree of heterogeneity in the intercellular coupling, we obtained a nice agreement of theoretical results with experimental measurements. A high degree of heterogeneity in coupling was necessary to omit a nearly complete synchronization of cells (see Figure S3). Most importantly, our results showed that [Ca²⁺]_c responses and spreading of [Ca²⁺]_c events between beta cells after stimulation with non-oscillatory glucose are characterized by a two-phased dynamics: the activation phase and the plateau phase. During the activation phase that lasts around 5 min in case of stimulation with 8 mM glucose, beta cells are recruited and the distribution of [Ca²⁺]_c event sizes is characterized by a power-law distribution, suggesting critical behavior. After this initial period, the dynamical nature changes qualitatively and becomes more stable and organized. In the plateau regime that follows, a high number of bigger, more global [Ca²⁺]_c events are observed, indicating supercritical behavior of [Ca²⁺]_c dynamics.

A constant, i.e., a non-oscillatory elevation in plasma glucose, may not reflect physiological stimulation of beta cells, since blood glucose levels in fasting man and other mammals are oscillating with periods of around 5–15 min (Goodner et al., 1977; Lang et al., 1979). Additionally, the average glucose concentration measured in vivo in mice (Kjems et al., 2002; Matyšková et al., 2008) lies just slightly above the in vitro determined threshold for glucose-induced metabolic, membrane potential, and [Ca²⁺]_c changes (Jonkers and Henquin, 2001; Zarkovic, 2004; Stožer et al., 2013a; Skelin Klemen et al., 2017), and stimulation with 8 mM glucose, albeit being low compared with concentrations that are usually used (e.g., 11.1 or 16.7 mM), is probably still supraphysiologically high. Thus, we attempted to more closely mimic physiological conditions and applied in our simulations and experiments a pulsatile stimulation with 5 min glucose pulses, with an average concentration of 7 mM glucose. This way, the total glucose load that the beta cells received and the maximum concentration reached in the perifusion chamber was the same in the two stimulation protocols (i.e., 20 min × 8 mM vs. 4 × 5 min × 8 mM, maximum = 8 mM). Moreover, during constant stimulation, the glucose concentration was just slightly above the average during oscillatory stimulation and we deliberately chose this value to mimic the continuous rise to hyperglycemia.

Noteworthy, although this was by no means the main aim of our study, our oscillatory protocol also addressed entrainability of islets. We wish to point out that switching from a non-stimulatory to a stimulatory concentration differs significantly from switching between two stimulatory concentrations when it comes to studying entrainability of nonlinear oscillators entrainment (Pedersen et al., 2013). In particular, the former scenario switches the system on and off and always results in entrainment, whereas the latter enables real islet entrainment (Pedersen et al., 2013). Since 6 mM glucose is usually regarded as non-stimulatory, at first our oscillatory protocol seems to correspond to the case of switching the system on and off. However, in our case, even during nadirs of glucose concentration, beta cell activity did not cease and it seems that our protocol therefore enables a true assessment of entrainability and our findings further substantiate recent reports that islets are in fact entrainable by stimuli of periods and amplitudes comparable to ours (Pedersen et al., 2013; Dhumpa et al., 2015; Sun et al., 2015).

During oscillatory stimulation, spatiotemporal organization of [Ca²⁺] waves remained in the critical regime without an excess of global [Ca²⁺] events. Noteworthy, the computational approach is, in contrast to experimental measurements, not restricted by technical limitations and can therefore account for very long observation times. The theoretical results have shown that the periodic stimulation does not lengthen the activation phase and delays the onset of supercritical behavior, but keeps the system essentially trapped in the critical regime. In general, the computational and experimental results are in very good agreement for both constant and periodic stimulation protocols. The only discrepancy occurs in the power-law exponents reflecting critical behavior, whose values are in simulations ranging between −1.41 and −1.58, whereas in experiments the values are a bit more negative and span between −1.78 and −1.94. This difference indicates that despite high levels of physiological complexity, the theoretical model probably still does not cover all the details about the intra- and inter-cellular aspects of beta cell signaling. We believe that especially a more detailed description of beta cell heterogeneity, both metabolic and electrophysiological, could provide a further step toward reality, but future experimental and theoretical work will be needed to address this issue. Moreover, in both computational and experimental results the exponents are found to be a bit more negative in the plateau phase than in the activation phase, which most probably reflects a sharper decay of smaller to intermediate events on account of global calcium waves.

Apparently, periodic entraining of the glucokinase reaction rates confines the supply of energy in the form of ATP, which in turn adjusts the rhythmic activity of ATP-dependent ion channels. This provides in combination with the intercellular coupling, which has either a suppressive effect that prevents uncoordinated and spontaneous activations, or a regulatory role in terms of mediating the synchronizing signal across the islet, the necessary substrate for self-organized activity. To be more precise, as a result of heterogeneities, dynamic and confined regions with elevated excitability emerge, from which [Ca²⁺]_c waves are triggered, which also goes in hand with previous findings, where the initiation of waves was associated with regions with a higher glucose metabolism (Benninger and Piston, 2014). The outreach of these waves was found to depend on the metabolic state of surrounding cells and on the local intercellular connectivity. On the other hand, if the glucokinase activity and the subsequent supply of metabolic energy is high, as in prolonged stimulatory conditions, all cells become on average more excitable and the level of their cell-to-cell variability decreases. Consequently, the regulatory role of gap junctional coupling fades and hence an activation of a few cells frequently leads to [Ca²⁺] waves that propagate throughout the whole islet, as is characteristic for supercritical behavior. However, due to the interplay between heterogeneities in cell metabolism, in the level of excitability, and in the intercellular coupling, the transition between the inactive and active beta cell network after switching to stimulatory conditions is not abrupt. Instead, the beta cell recruitment is a time dependent process that evokes an emergent spatiotemporal dynamics characterized by power-law scaling. In this vein, the concept of a functionally heterogeneous beta cell population can be regarded as a key determinant of critical behavior in islets.

Trapping beta cells in a self-organized critical state might be of crucial physiological importance. Namely, criticality seems to be crucial for enabling living tissues reasonable handling of energy and providing efficient functioning and optimized responses to external stimuli. For example, in the field of neuroscience it has been hypothesized that a normal, healthy brain resides in a critical state, which provides the fastest and most flexible adaptation to different challenges from the environment (Chialvo, 2010; Hesse and Gross, 2014; Massobrio et al., 2015a). Tomen et al. (2014) have also shown that cortical networks exhibit optimal neuronal information processing at a near-critical state, i.e., in a narrow region in the phase space at the transition from subcritical to supercritical dynamics. Some studies have even put forward the idea that an optimal physiological neural activity does not perfectly reflect the SOC state but can be characterized as a subcritical regime slightly below the SOC without a separation of time scales, exhibiting the so-called ≫mélange≪ of avalanches (Priesemann et al., 2013, 2014). Potential advantages of this marginally subcritical regime slightly below the SOC may be a more efficient information processing, and a safety margin from supercriticality, which has been linked to some pathophysiological disorders. Interestingly, very recently it has been shown that also social systems are poised in the proximity of critical points and by tuning the distance to this point facilitates the system to favor either stability or flexibility (Daniels et al., 2017).

Even within nadirs of our oscillatory stimulation protocol, we were unable to detect any convincing evidence of subcriticality. However, this does not exclude the possibility that in vivo, subcriticality in the form of localized individual [Ca²⁺]_c might be present. Glucose is a signal that rapidly reaches all beta cells also in vivo and it is hard to believe that differences in supply of glucose to different parts of islets or differences in sensitivity of beta cells to glucose could be much higher in vivo than in our in situ preparation and bring about subcritical behavior (Michau et al., 2016). However, in case of incretins which are also able to evoke [Ca²⁺]_c changes in beta cells, evidence exists that especially in human islets, not all beta cells respond to GLP-1 equally well (Hodson et al., 2013, 2014). If differences in sensitivity are great enough, one can imagine that a local [Ca²⁺]_c response of a highly sensitive cell might remain localized, especially if the momentary glucose concentration is not great enough to support propagation to neighboring cells (Eddlestone et al., 1984). Likewise, in vivo, subcritical behavior of spatiotemporal [Ca²⁺]_c dynamics could be brought about by signals able to evoke beta cell [Ca²⁺]_c responses that stem from local sources and conceivably stimulate some beta cells more than others, for instance ATP, acetylcholine, or glucagon released from parasympathetic nerve endings or alpha cells (Rodriguez-Diaz et al., 2011a,b; Gylfe et al., 2012). These possibilities remain to be investigated experimentally.

We might further hypothesize that islet tissue operates as a so-called driven subcritical system that is switched to SOC by an entrainment with glucose dynamics. It is well established that glucose oscillations are a hallmark of normal glucose tolerance and islet function (Lang et al., 1979, 1981; Mao et al., 1999; Ritzel et al., 2005). Blood glucose concentration oscillates in both monkeys and humans with a period in the range of the slow insulin oscillations and an amplitude of approximately 1-10 % with respect to the average concentration (Goodner et al., 1977; Lang et al., 1979; Mao et al., 1999). These glucose oscillations seem to be strongly correlated with pulsatile secretion of insulin (O'Meara et al., 1993; Mao et al., 1999; Ritzel et al., 2005; Pedersen et al., 2013; Nunemaker and Satin, 2014; Satin et al., 2015). In fact, it was shown that oscillations of insulin secretion can be entrained by imposed small changes in glucose concentration in vivo in normal subjects and that this ability is lost in T2DM (Mao et al., 1999; Hollingdal et al., 2000). It was further suggested that oscillatory glucose actually amplifies the mass of insulin secretory pulses that coincide with imposed glucose stimuli and that the intrinsic frequency of insulin oscillations does not change (Ritzel et al., 2005). This seems in contrast with recent findings that slow oscillations in metabolism and [Ca²⁺]_c that are believed to underlie insulin pulses (Pedersen et al., 2013; Sun et al., 2015), as well as insulin secretion can indeed be entrained to glucose (Dhumpa et al., 2015). However, all these contradictions might only be a consequence of still unknown underlying mechanisms for self-organization in human body and the relation of these micro-mechanisms with the emergent macro-phenomena. Recently, Lo et al. (2013) have demonstrated that in contrast to the macroscopic homeostatic equilibrium that describes sleep at the circadian time scale of several hours, the sleep micro-architecture at scales from seconds to minutes exhibits a non-equilibrium behavior of SOC type. They argue that the asymmetry in the transitions between quiet states and avalanches is important as the energy can slowly build up during quiet states, i.e., slowly approaching the critical point, and dissipates rapidly when avalanches occur. Nevertheless, methodologies and novel, more holistic approaches to studying such interconnected physiological processes, occurring in a broad range of space and time scales, are only beginning to emerge. The recently emerging fields of network physiology and network medicine show great potential to provide new insights into how global behavior at the organism level can arise out of micro-mechanisms on the cellular and tissue level, along with networked interactions among different organ systems, to generate health or disease (Bashan et al., 2012; Ivanov et al., 2016).

Our results show that a change in glucose stimulus similar to the one during development of T2DM, i.e., exposure of islets to non-oscillatory and slightly elevated glucose concentration, leads to an excess of large events in the spatiotemporal pattern of fast [Ca²⁺]_c oscillations that are believed to determine the pulse mass of insulin oscillations, suggesting that a switch to supercriticality might be an important pathophysiological mechanism in T2DM. This finding corresponds with recent investigations in the brain, showing that epilepsy, for example, is characterized by a supercritical behavior of neurons (Meisel et al., 2012; Priesemann et al., 2014), as well as in a number of other tissues and disease states, such as in obesity and irritable bowel syndrome (caused by a transition in microbial composition), asthma and other pulmonary diseases, depression, inflammation, cancer, and cardiovascular events (for review see Trefois et al., 2015).

In general, an improved understanding of such transitions to supercriticality from normal healthy states during onset and progression of diseases could have important applications in health care. From the viewpoint of preventive, the identification and characterization of early warning signals could predict upcoming critical transitions; and from the view point of curative, the understanding of critical transitions might help us develop more effective therapeutic applications.

Author Contributions

MG, AS, RM, JD, MP, MR, and MM designed and performed the research as well as wrote the paper.

Funding

The authors acknowledge the financial support from the Slovenian Research Agency (research core funding, Nos. P3-0396 and I0-0029-0552, as well as research projects, Nos. N3-0048, J3-7177, J1-7009, and J7-7226). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2017.01106/full#supplementary-material

Text S1. Mathematical model and parameter values for beta cells.

Figure S1. The course of the simulated oscillatory stimulation protocol. Variations in the glucokinase reaction rate reflect the oscillatory changes in glucose.

Figure S2. Coupling in the multicellular beta cell model. (A) A typical structure of the intercellular network of beta cells. (B) The corresponding degree distribution. The beta cell network is quite homogeneous with a mean degree around 6. (C) The distribution of the electrical coupling coefficient. The coupling strength is rather heterogeneous and follows an exponential distribution with a mean of 200 pS.

Figure S3. Simulated spatio-temporal activity under constant stimulation with homogeneous coupling. Typical computed [Ca²⁺]_c responses of four different beta cells after switching to stimulatory conditions (upper panel, gray area indicates stimulatory conditions) and binarization of the computed oscillations of all cells (lower panel). The electrical coupling coefficient was distributed normally with mean 200 pS and relative SD of 30%. In this case very synchronized behavior is obtained without progressive and heterogeneous activations of cells, as observed in experiments.

Video S1. Representative animation of computed and binarized spatiotemporal [Ca²⁺]_i activity under constant stimulation with glucose.

Video S2. Representative animation of computed and binarized spatiotemporal [Ca²⁺]_i activity under periodic stimulation with glucose.

Video S3. Movie of experimentally measured and binarized [Ca²⁺]_i activity under constant stimulation with 8 mM glucose from the onset of glucose increase.

Video S4. Movie of experimentally measured and binarized [Ca²⁺]_i activity under periodic stimulation with 6-8-6-8-6-8-6 mM glucose.

References

Aizawa, T., Komatsu, M., Asanuma, N., Sato, Y., and Sharp, G. W. (1998). Glucose action ‘beyond ionic events’ in the pancreatic β cell. Trends Pharmacol. Sci. 19, 496–499.

PubMed Abstract | Google Scholar

Allegrini, P., Paradisi, P., Menicucci, D., Laurino, M., Piarulli, A., and Gemignani, A. (2015). Self-organized dynamical complexity in human wakefulness and sleep: different critical brain-activity feedback for conscious and unconscious states. Phys. Rev. E 92, 32808–32820. doi: 10.1103/PhysRevE.92.032808

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashcroft, F. M., and Rorsman, P. (2013). KATP channels and islet hormone secretion: new insights and controversies. Nat. Rev. Endocrinol. 9, 660–669. doi: 10.1038/nrendo.2013.166

CrossRef Full Text | Google Scholar

Bak, P. (1996). How Nature Works. New York, NY: Springer New York.

Google Scholar

Barua, A. K., and Goel, P. (2016). Isles within islets: the lattice origin of small-world networks in pancreatic tissues. Phys. D 315, 49–57. doi: 10.1016/j.physd.2015.07.009

CrossRef Full Text | Google Scholar

Bashan, A., Bartsch, R. P., Kantelhardt, J. W., Havlin, S., and Ivanov, P. C. (2012). Network physiology reveals relations between network topology and physiological function. Nat. Commun. 3, 702–722. doi: 10.1038/ncomms1705

PubMed Abstract | CrossRef Full Text | Google Scholar

Bavamian, S., Klee, P., Britan, A., Populaire, C., Caille, D., Cancela, J., et al. (2007). Islet-cell-to-cell communication as basis for normal insulin secretion. Diabetes, Obes. Metab. 9, 118–132. doi: 10.1111/j.1463-1326.2007.00780.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Beggs, J. M., and Plenz, D. (2003). Neuronal avalanches in neocortical circuits. J. Neurosci. 23, 11167–11177.

PubMed Abstract | Google Scholar

Benninger, R. K. P., Head, W. S., Zhang, M., Satin, L. S., and Piston, D. W. (2011). Gap junctions and other mechanisms of cell–cell communication regulate basal insulin secretion in the pancreatic islet. J. Physiol. 58922, 5453–5466. doi: 10.1113/jphysiol.2011.218909

CrossRef Full Text | Google Scholar

Benninger, R. K. P. P., Hutchens, T., Head, W. S., McCaughey, M. J., Zhang, M., Le Marchand, S. J., et al. (2014). Intrinsic islet heterogeneity and gap junction coupling determine spatiotemporal Ca²⁺ wave dynamics. Biophys. J. 107, 2723–2733. doi: 10.1016/j.bpj.2014.10.048

CrossRef Full Text | Google Scholar

Benninger, R. K. P., and Piston, D. W. (2014). Cellular communication and heterogeneity in pancreatic islet insulin secretion dynamics. Trends Endocrinol. Metab. 25, 399–406. doi: 10.1016/j.tem.2014.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Benninger, R. K. P., Zhang, M., Head, W. S., Satin, L. S., and Piston, D. W. (2008). Gap junction coupling and calcium waves in the pancreatic islet. Biophys. J. 95, 5048–5061. doi: 10.1529/biophysj.108.140863

PubMed Abstract | CrossRef Full Text | Google Scholar

Bergsten, P. (1995). Slow and fast oscillations of cytoplasmic Ca²⁺ in pancreatic islets correspond to pulsatile insulin release. Am. J. Physiol. 268, E282–E287.

PubMed Abstract | Google Scholar

Bertram, R., Satin, L. S., Pedersen, M. G., Luciani, D. S., and Sherman, A. (2007). Interaction of glycolysis and mitochondrial respiration in metabolic oscillations of pancreatic islets. Biophys. J. 92, 1544–1555. doi: 10.1529/biophysj.106.097154

PubMed Abstract | CrossRef Full Text | Google Scholar

Bingley, P. J., Matthews, D. R., Williams, A. J. K., Bottazzo, G. F., and Gale, E. A. M. (1992). Loss of regular oscillatory insulin secretion in islet cell antibody positive non-diabetic subjects. Diabetologia 35, 32–38. doi: 10.1007/BF00400849

PubMed Abstract | CrossRef Full Text | Google Scholar

Cappon, G., and Pedersen, M. G. (2016). Heterogeneity and nearest-neighbor coupling can explain small-worldness and wave properties in pancreatic islets. Chaos 26, 53103–53107. doi: 10.1063/1.4949020

PubMed Abstract | CrossRef Full Text | Google Scholar

Charollais, A., Gjinovci, A., Huarte, J., Bauquis, J., Nadal, A., Martín, F., et al. (2000). Junctional communication of pancreatic β cells contributes to the control of insulin secretion and glucose tolerance. J. Clin. Invest. 106, 235–243. doi: 10.1172/JCI9398

PubMed Abstract | CrossRef Full Text | Google Scholar

Cherubini, C., Filippi, S., Gizzi, A., and Loppini, A. (2015). Role of topology in complex functional networks of beta cells. Phys. Rev. E 92, 42702–42712. doi: 10.1103/PhysRevE.92.042702

PubMed Abstract | CrossRef Full Text | Google Scholar

Chialvo, D. R. (2010). Emergent complex neural dynamics. Nat. Phys. 6, 744–750. doi: 10.1038/nphys1803

CrossRef Full Text | Google Scholar

Daniels, B. C., Krakauer, D. C., and Flack, J. C. (2017). Control of finite critical behaviour in a small-scale social system. Nat. Commun. 8, 14301–14308. doi: 10.1038/ncomms14301

PubMed Abstract | CrossRef Full Text | Google Scholar

Dean, P. M., and Matthews, E. K. (1970). Glucose-induced electrical activity in pancreatic islet cells. J. Physiol. 210, 255–264. doi: 10.1113/jphysiol.1970.sp009207

PubMed Abstract | CrossRef Full Text | Google Scholar

Dhumpa, R., Truong, T. M., Wang, X., Bertram, R., and Roper, M. G. (2014). Negative feedback synchronizes islets of langerhans. Biophys. J. 106, 2275–2282. doi: 10.1016/j.bpj.2014.04.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Dhumpa, R., Truong, T. M., Wang, X., and Roper, M. G. (2015). Measurement of the entrainment window of islets of Langerhans by microfluidic delivery of a chirped glucose waveform. Integr. Biol. 7, 1061–1067. doi: 10.1039/C5IB00156K

PubMed Abstract | CrossRef Full Text | Google Scholar

Dolenšek, J., Rupnik, M. S., and Stožer, A. (2015). Structural similarities and differences between the human and the mouse pancreas. Islets 7, e1024405–e1024416. doi: 10.1080/19382014.2015.1024405

PubMed Abstract | CrossRef Full Text | Google Scholar

Dolenšek, J., Stožer, A., Skelin Klemen, M., Miller, E. W., and Slak Rupnik, M. (2013). The relationship between membrane potential and calcium dynamics in glucose-stimulated beta cell syncytium in acute mouse pancreas tissue slices. PLoS ONE 8:e82374. doi: 10.1371/journal.pone.0082374

PubMed Abstract | CrossRef Full Text | Google Scholar

Eddlestone, G. T., Gonçalves, A., Bangham, J. A., and Rojas, E. (1984). Electrical coupling between cells in islets of Langerhans from mouse. J. Membr. Biol. 77, 1–14. doi: 10.1007/BF01871095

PubMed Abstract | CrossRef Full Text | Google Scholar

Farnsworth, N. L., and Benninger, R. K. P. (2014). New insights into the role of connexins in pancreatic islet function and diabetes. FEBS Lett. 588, 1278–1287. doi: 10.1016/j.febslet.2014.02.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Farnsworth, N. L., Walter, R. L., Hemmati, A., Westacott, M. J., and Benninger, R. K. P. (2016). low level pro-inflammatory cytokines decrease connexin36 gap junction coupling in mouse and human islets through nitric oxide-mediated protein kinase Cδ. J. Biol. Chem. 291, 3184–3196. doi: 10.1074/jbc.M115.679506

PubMed Abstract | CrossRef Full Text | Google Scholar

Fendler, B., Zhang, M., Satin, L., and Bertram, R. (2009). Synchronization of pancreatic islet oscillations by intrapancreatic ganglia: a modeling study. Biophys. J. 97, 722–729. doi: 10.1016/j.bpj.2009.05.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Fortrat, J.-O., and Gharib, C. (2016). Self-Organization of blood pressure regulation: clinical evidence. Front. Physiol. 7:113. doi: 10.3389/fphys.2016.00113

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilon, P., and Henquin, J. C. (1995). Distinct effects of glucose on the synchronous oscillations of insulin release and cytoplasmic Ca²⁺ concentration measured simultaneously in single mouse islets. Endocrinology 136, 5725–5730. doi: 10.1210/endo.136.12.7588329

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilon, P., Henquin, J., Gilontj, P., and Henquintll, J.-C. (1992). Influence of membrane potential changes on cytoplasmic Ca²⁺ concentration in an electrically excitable cell, the insulin-secreting pancreatic B-cell. J. Biol. Chem. 267, 20713–20720.

PubMed Abstract | Google Scholar

Goel, P., and Mehta, A. (2013). Learning theories reveal loss of pancreatic electrical connectivity in diabetes as an adaptive response. PLoS ONE 8:e70366. doi: 10.1371/journal.pone.0070366

PubMed Abstract | CrossRef Full Text | Google Scholar

Goodner, C., Walike, B., Koerker, D., Ensinck, J., Brown, A., Chideckel, E., et al. (1977). Insulin, glucagon, and glucose exhibit synchronous, sustained oscillations in fasting monkeys. Science 195, 177–179. doi: 10.1126/science.401543

PubMed Abstract | CrossRef Full Text | Google Scholar

Gosak, M., Dolenšek, J., Markovič, R., Slak Rupnik, M., Marhl, M., and Stožer, A. (2015). Multilayer network representation of membrane potential and cytosolic calcium concentration dynamics in beta cells. Chaos Solitons Fractals 80, 76–82. doi: 10.1016/j.chaos.2015.06.009

CrossRef Full Text | Google Scholar

Gosak, M., Markovič, R., Dolenšek, J., Slak Rupnik, M., Marhl, M., and Stožer, A. (in press). Network science of biological systems at different scales: a review. Phys. Life Rev. doi: 10.1016/j.plrev.2017.11.003.

PubMed Abstract | CrossRef Full Text | Google Scholar

Gylfe, E., Grapengiesser, E., Dansk, H., and Hellman, B. (2012). The neurotransmitter ATP triggers Ca²⁺ Responses promoting coordination of pancreatic islet oscillations. Pancreas 41, 258–263. doi: 10.1097/MPA.0b013e3182240586

PubMed Abstract | CrossRef Full Text | Google Scholar

Haefliger, J. A. J.-A., Rohner-Jeanrenaud, F., Caille, D., Charollais, A., Meda, P., and Allagnat, F. (2013). Hyperglycemia downregulates Connexin36 in pancreatic islets via the upregulation of ICER-1/ICER-1γ. J. Mol. Endocrinol. 51, 49–58. doi: 10.1530/JME-13-0054

PubMed Abstract | CrossRef Full Text | Google Scholar

Head, W. S., Orseth, M. L., Nunemaker, C. S., Satin, L. S., Piston, D. W., and Benninger, R. K. P. (2012). Connexin-36 gap junctions regulate in vivo first- and second-phase insulin secretion dynamics and glucose tolerance in the conscious mouse. Diabetes 61, 1700–1707. doi: 10.2337/db11-1312

PubMed Abstract | CrossRef Full Text | Google Scholar

Henquin, J. C. (1987). Regulation of insulin release by ionic and electrical events in B Cells. Horm. Res. 27, 168–178. doi: 10.1159/000180806

PubMed Abstract | CrossRef Full Text | Google Scholar

Henquin, J.-C. (2011). The dual control of insulin secretion by glucose involves triggering and amplifying pathways in β-cells. Diabetes Res. Clin. Pract. 93, S27–S31. doi: 10.1016/S0168-8227(11)70010-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Henquin, J. C., Meissner, H. P., and Schmeer, W. (1982). Cyclic variations of glucose-induced electrical activity in pancreatic B cells. Pflugers Arch. 393, 322–327. doi: 10.1007/BF00581418

PubMed Abstract | CrossRef Full Text | Google Scholar

Hesse, J., and Gross, T. (2014). Self-organized criticality as a fundamental property of neural systems. Front. Syst. Neurosci. 8:166. doi: 10.3389/fnsys.2014.00166

PubMed Abstract | CrossRef Full Text | Google Scholar

Hiriart, M., and Ramirez-Medeles, M. C. (1991). Functional subpopulations of individual pancreatic B-cells in culture. Endocrinology 128, 3193–3198. doi: 10.1210/endo-128-6-3193

PubMed Abstract | CrossRef Full Text | Google Scholar

Hodson, D. J., Mitchell, R. K., Bellomo, E. A., Sun, G., Vinet, L., Meda, P., et al. (2013). Lipotoxicity disrupts incretin-regulated human β cell connectivity. J. Clin. Invest. 123, 4182–4194. doi: 10.1172/JCI68459

PubMed Abstract | CrossRef Full Text | Google Scholar

Hodson, D. J., Tarasov, A. I., Gimeno Brias, S., Mitchell, R. K., Johnston, N. R., Haghollahi, S., et al. (2014). Incretin-modulated beta cell energetics in intact islets of langerhans. Mol. Endocrinol. 28, 860–871. doi: 10.1210/me.2014-1038

PubMed Abstract | CrossRef Full Text | Google Scholar

Hollingdal, M., Juhl, C. B., Pincus, S. M., Sturis, J., Veldhuis, J. D., Polonsky, K. S., et al. (2000). Failure of physiological plasma glucose excursions to entrain high-frequency pulsatile insulin secretion in type 2 diabetes. Diabetes 49, 1334–1340. doi: 10.2337/diabetes.49.8.1334

PubMed Abstract | CrossRef Full Text | Google Scholar

Hraha, T. H., Bernard, A. B., Nguyen, L. M., Anseth, K. S., and Benninger, R. K. P. (2014a). Dimensionality and size scaling of coordinated Ca²⁺ dynamics in MIN6 β-cell clusters. Biophys. J. 106, 299–309. doi: 10.1016/j.bpj.2013.11.026

CrossRef Full Text | Google Scholar

Hraha, T. H., Westacott, M. J., Pozzoli, M., Notary, A. M., McClatchey, P. M., and Benninger, R. K. P. (2014b). Phase transitions in the multi-cellular regulatory behavior of pancreatic islet excitability. PLoS Comput. Biol. 10:e1003819. doi: 10.1371/journal.pcbi.1003819

PubMed Abstract | CrossRef Full Text | Google Scholar

Hutt, M.-T., Kaiser, M., and Hilgetag, C. C. (2014). Perspective: network-guided pattern formation of neural dynamics. Philos. Trans. R. Soc. B Biol. Sci. 369, 20130522–20130522. doi: 10.1098/rstb.2013.0522

PubMed Abstract | CrossRef Full Text

Ivanov, P. C., Liu, K. K. L., and Bartsch, R. P. (2016). Focus on the emerging new fields of network physiology and network medicine. New J. Phys. 18, 100201–100209. doi: 10.1088/1367-2630/18/10/100201

CrossRef Full Text | Google Scholar

Ivanov, P. C., Ma, Q. D. Y., Bartsch, R. P., Hausdorff, J. M., Nunes Amaral, L. A., Schulte-Frohlinde, V., et al. (2009). Levels of complexity in scale-invariant neural signals. Phys. Rev. E 79, 41920–41912. doi: 10.1103/PhysRevE.79.041920

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnston, N. R., Mitchell, R. K., Haythorne, E., Pessoa, M. P., Semplici, F., Ferrer, J., et al. (2016). Beta cell hubs dictate pancreatic islet responses to glucose. Cell Metab. 24, 389–401. doi: 10.1016/j.cmet.2016.06.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Jonkers, F. C., and Henquin, J.-C. (2001). Measurements of cytoplasmic Ca²⁺ in islet cell clusters show that glucose rapidly recruits β-cells and gradually increases the individual cell response. Diabetes 50, 540–550. doi: 10.2337/diabetes.50.3.540

PubMed Abstract | CrossRef Full Text | Google Scholar

Jonkers, F. C., Jonas, J.-C., Gilon, P., and Henquin, J.-C. (1999). Influence of cell number on the characteristics and synchrony of Ca²⁺ oscillations in clusters of mouse pancreatic islet cells. J. Physiol. 520, 839–849. doi: 10.1111/j.1469-7793.1999.00839.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jung, P. (1997). Thermal waves, criticality, and self-organization in excitable media. Phys. Rev. Lett. 78, 1723–1726. doi: 10.1103/PhysRevLett.78.1723

CrossRef Full Text | Google Scholar

Jung, P., Cornell-Bell, A., Madden, K. S., and Moss, F. (1998). Noise-induced spiral waves in astrocyte syncytia show evidence of self-organized criticality. J. Neurophysiol. 79, 1098–1101.

PubMed Abstract | Google Scholar

Kahn, S. E., Cooper, M. E., and Del Prato, S. (2014). Pathophysiology and treatment of type 2 diabetes: perspectives on the past, present, and future. Lancet 383, 1068–1083. doi: 10.1016/S0140-6736(13)62154-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Kiekens, R., In 't Veld, P., Mahler, T., Schuit, F., Van De Winkel, M., and Pipeleers, D. (1992). Differences in glucose recognition by individual rat pancreatic B cells are associated with intercellular differences in glucose-induced biosynthetic activity. J. Clin. Invest. 89, 117–125.

PubMed Abstract | Google Scholar

Kinouchi, O., and Copelli, M. (2006). Optimal dynamical range of excitable networks at criticality. Nat. Phys. 2, 348–351. doi: 10.1038/nphys289

CrossRef Full Text | Google Scholar

Kjems, L. L., Ravier, M. A., Jonas, J.-C., and Henquin, J.-C. (2002). Do oscillations of insulin secretion occur in the absence of cytoplasmic Ca²⁺ oscillations in β-Cells? Diabetes 51, S177–S182. doi: 10.2337/diabetes.51.2007.S177

PubMed Abstract | CrossRef Full Text | Google Scholar

Konstantinova, I., Nikolova, G., Ohara-Imaizumi, M., Meda, P., Kucera, T., Zarbalis, K., et al. (2007). EphA-Ephrin-A-Mediated β cell communication regulates insulin secretion from pancreatic islets. Cell 129, 359–370. doi: 10.1016/j.cell.2007.02.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Lang, D. A., Matthews, D. R., Burnett, M., and Turner, R. C. (1981). Brief, irregular oscillations of basal plasma insulin and glucose concentrations in diabetic man. Diabetes 30, 435–439. doi: 10.2337/diab.30.5.435

PubMed Abstract | CrossRef Full Text | Google Scholar

Lang, D. A., Matthews, D. R., Peto, J., and Turner, R. C. (1979). Cyclic oscillations of basal plasma glucose and insulin concentrations in human beings. N. Engl. J. Med. 301, 1023–1027. doi: 10.1056/NEJM197911083011903

PubMed Abstract | CrossRef Full Text | Google Scholar

Lehnertz, K., Bialonski, S., Horstmann, M.-T., Krug, D., Rothkegel, A., Staniek, M., et al. (2009). Synchronization phenomena in human epileptic brain networks. J. Neurosci. Methods 183, 42–48. doi: 10.1016/j.jneumeth.2009.05.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Levina, A., Herrmann, J. M., and Geisel, T. (2007). Dynamical synapses causing self-organized criticality in neural networks. Nat. Phys. 3, 857–860. doi: 10.1038/nphys758

CrossRef Full Text | Google Scholar

Lo, C.-C., Bartsch, R. P., and Ivanov, P. C. (2013). Asymmetry and basic pathways in sleep-stage transitions. Europhys. Lett. 102:10008. doi: 10.1209/0295-5075/102/10008

PubMed Abstract | CrossRef Full Text | Google Scholar

Lo, C.-C., Chou, T., Penzel, T., Scammell, T. E., Strecker, R. E., Stanley, H. E., et al. (2004). Common scale-invariant patterns of sleep-wake transitions across mammalian species. Proc. Natl. Acad. Sci. U.S.A. 101, 17545–17548. doi: 10.1073/pnas.0408242101

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopez, L., Piegari, E., Sigaut, L., and Ponce Dawson, S. (2012). Intracellular calcium signals display an avalanche-like behavior over multiple lengthscales. Front. Physiol. 3:350. doi: 10.3389/fphys.2012.00350

PubMed Abstract | CrossRef Full Text | Google Scholar

Lovecchio, E., Allegrini, P., Geneston, E., West, B. J., and Grigolini, P. (2012). From self-organized to extended criticality. Front. Physiol. 3:98. doi: 10.3389/fphys.2012.00098

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, C. S., Berman, N., Roberts, K., and Ipp, E. (1999). Glucose entrainment of high-frequency plasma insulin oscillations in control and type 2 diabetic subjects. Diabetes 48, 714–721. doi: 10.2337/diabetes.48.4.714

PubMed Abstract | CrossRef Full Text | Google Scholar

Marković, D., and Gros, C. (2014). Power laws and self-organized criticality in theory and nature. Phys. Rep. 536, 41–74. doi: 10.1016/j.physrep.2013.11.002

CrossRef Full Text | Google Scholar

Markovič, R., Stožer, A., Gosak, M., Dolenšek, J., Marhl, M., and Rupnik, M. S. (2015). Progressive glucose stimulation of islet beta cells reveals a transition from segregated to integrated modular functional connectivity patterns. Sci. Rep. 5, 7845–7810. doi: 10.1038/srep07845

PubMed Abstract | CrossRef Full Text | Google Scholar

Massobrio, P., de Arcangelis, L., Pasquale, V., Jensen, H. J., and Plenz, D. (2015a). Criticality as a signature of healthy neural systems. Front. Syst. Neurosci. 9:22. doi: 10.3389/fnsys.2015.00022

PubMed Abstract | CrossRef Full Text | Google Scholar

Massobrio, P., Pasquale, V., and Martinoia, S. (2015b). Self-organized criticality in cortical assemblies occurs in concurrent scale-free and small-world networks. Sci. Rep. 5, 10578–10516. doi: 10.1038/srep10578

PubMed Abstract | CrossRef Full Text | Google Scholar

Matyšková, R., Maletínská, L., Maixnerová, J., Pirník, Z., Kiss, A., Železná, B., et al. (2008). Comparison of the obesity phenotypes related to monosodium glutamate effect on arcuate nucleus and/or the high fat diet feeding in C57BL/6 and NMRI mice. Physiol. Res 57, 727–734.

PubMed Abstract | Google Scholar

Meisel, C., Storch, A., Hallmeyer-Elgner, S., Bullmore, E., and Gross, T. (2012). Failure of adaptive self-organized criticality during epileptic seizure attacks. PLoS Comput. Biol. 8:e1002312. doi: 10.1371/journal.pcbi.1002312

PubMed Abstract | CrossRef Full Text | Google Scholar

Michau, A., Hodson, D. J., Fontanaud, P., Guillou, A., Espinosa-Carrasco, G., Molino, F., et al. (2016). Metabolism regulates exposure of pancreatic islets to circulating molecules in vivo. Diabetes 65, 463–475. doi: 10.2337/db15-1168

PubMed Abstract | CrossRef Full Text | Google Scholar

Moretti, P., and Muñoz, M. A. (2013). Griffiths phases and the stretching of criticality in brain networks. Nat. Commun. 4, 2521–2510. doi: 10.1038/ncomms3521

PubMed Abstract | CrossRef Full Text | Google Scholar

Nittala, A., Ghosh, S., Wang, X., Kurths, J., and Tass, P. (2007). Investigating the role of islet cytoarchitecture in its oscillation using a new β-cell cluster model. PLoS ONE 2:e983. doi: 10.1371/journal.pone.0000983

PubMed Abstract | CrossRef Full Text | Google Scholar

Nivala, M., Ko, C. Y., Nivala, M., Weiss, J. N., and Qu, Z. (2012). Criticality in intracellular calcium signaling in cardiac myocytes. Biophys. J. 102, 2433–2442. doi: 10.1016/j.bpj.2012.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Nunemaker, C. S. (2005). Insulin secretion in the conscious mouse is biphasic and pulsatile. Am. J. Physiol. 290, E523–E529. doi: 10.1152/ajpendo.00392.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Nunemaker, C. S., Dishinger, J. F., Dula, S. B., Wu, R., Merrins, M. J., Reid, K. R., et al. (2009). Glucose metabolism, islet architecture, and genetic homogeneity in imprinting of [Ca²⁺]_i and insulin rhythms in mouse islets. PLoS ONE 4:e8428. doi: 10.1371/journal.pone.0008428

CrossRef Full Text | Google Scholar

Nunemaker, C. S., and Satin, L. S. (2014). Episodic hormone secretion: a comparison of the basis of pulsatile secretion of insulin and GnRH. Endocrine 47, 49–63. doi: 10.1007/s12020-014-0212-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Nunemaker, C. S., Zhang, M., Wasserman, D. H., McGuinness, O. P., Powers, A. C., Bertram, R., et al. (2005). Individual mice can be distinguished by the period of their islet calcium oscillations: is there an intrinsic islet period that is imprinted in vivo? Diabetes 54, 3517–3522. doi: 10.2337/diabetes.54.12.3517

PubMed Abstract | CrossRef Full Text | Google Scholar

O'Meara, N. M., Sturis, J., Blackman, J. D., Byrne, M. M., Jaspan, J. B., Roland, D. C., et al. (1993). Oscillatory insulin secretion after pancreas transplant. Diabetes 42, 855–861.

PubMed Abstract | Google Scholar

O'Rahilly, S., Turner, R. C., and Matthews, D. R. (1988). Impaired pulsatile secretion of insulin in relatives of patients with non-insulin-dependent diabetes. N. Engl. J. Med. 318, 1225–1230.

PubMed Abstract | Google Scholar

Pedersen, M. G., Bertram, R., and Sherman, A. (2005). Intra- and inter-islet synchronization of metabolically driven insulin secretion. Biophys. J. 89, 107–119. doi: 10.1529/biophysj.104.055681

PubMed Abstract | CrossRef Full Text | Google Scholar

Pedersen, M. G., Mosekilde, E., Polonsky, K. S., and Luciani, D. S. (2013). Complex patterns of metabolic and Ca²⁺ entrainment in pancreatic islets by oscillatory glucose. Biophys. J. 105, 29–39. doi: 10.1016/j.bpj.2013.05.036

CrossRef Full Text | Google Scholar

Pedersen, M. G., and Sherman, A. (2009). Newcomer insulin secretory granules as a highly calcium-sensitive pool. Proc. Natl. Acad. Sci. U.S.A. 106, 7432–7436. doi: 10.1073/pnas.0901202106

PubMed Abstract | CrossRef Full Text | Google Scholar

Pérez-Armendariz, M., Roy, C., Spray, D. C., and Bennett, M. V. (1991). Biophysical properties of gap junctions between freshly dispersed pairs of mouse pancreatic beta cells. Biophys. J. 59, 76–92. doi: 10.1016/S0006-3495(91)82200-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Pipeleers, D., Kiekens, R., Ling, Z., Wilikens, A., and Schuit, F. (1994). Physiologic relevance of heterogeneity in the pancreatic beta-cell population. Diabetologia 37(Suppl. 2), S57–S64. doi: 10.1007/BF00400827

PubMed Abstract | CrossRef Full Text | Google Scholar

Plenz, D., and Thiagarajan, T. C. (2007). The organizing principles of neuronal avalanches: cell assemblies in the cortex? Trends Neurosci. 30, 101–110. doi: 10.1016/j.tins.2007.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Pørksen, N. (2002). The in vivo regulation of pulsatile insulin secretion. Diabetologia 45, 3–20. doi: 10.1007/s125-002-8240-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Pørksen, N., Hollingdal, M., Juhl, C., Butler, P., Veldhuis, J. D., and Schmitz, O. (2002). Pulsatile insulin secretion: detection, regulation, and role in diabetes. Diabetes 51, S245–S254. doi: 10.2337/diabetes.51.2007.S245

PubMed Abstract | CrossRef Full Text | Google Scholar

Priesemann, V., Valderrama, M., Wibral, M., and Le Van Quyen, M. (2013). Neuronal avalanches differ from wakefulness to deep sleep – evidence from intracranial depth recordings in humans. PLoS Comput. Biol. 9:e1002985. doi: 10.1371/journal.pcbi.1002985

PubMed Abstract | CrossRef Full Text | Google Scholar

Priesemann, V., Wibral, M., Valderrama, M., Pröpper, R., Le Van Quyen, M., Geisel, T., et al. (2014). Spike avalanches in vivo suggest a driven, slightly subcritical brain state. Front. Syst. Neurosci. 8:108. doi: 10.3389/fnsys.2014.00108

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritzel, R. A., Veldhuis, J. D., and Butler, P. C. (2005). The mass, but not the frequency, of insulin secretory bursts in isolated human islets is entrained by oscillatory glucose exposure. Am. J. Physiol. Endocrinol. Metab. 290, E750–E756. doi: 10.1152/ajpendo.00381.2005

CrossRef Full Text | Google Scholar

Rodriguez-Diaz, R., Abdulreda, M., Formoso, A., Gans, I., Ricordi, C., Berggren, P.-O., et al. (2011a). Innervation patterns of autonomic axons in the human endocrine pancreas. Cell Metab. 14, 45–54. doi: 10.1016/j.cmet.2011.05.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Rodriguez-Diaz, R., Dando, R., Jacques-Silva, M. C., Fachado, A., Molina, J., Abdulreda, M. H., et al. (2011b). Alpha cells secrete acetylcholine as a non-neuronal paracrine signal priming beta cell function in humans. Nat. Med. 17, 888–892. doi: 10.1038/nm.2371

PubMed Abstract | CrossRef Full Text | Google Scholar

Rubinov, M., Sporns, O., Thivierge, J.-P., Breakspear, M., and Weinberger, D. (2011). Neurobiologically realistic determinants of self-organized criticality in networks of spiking neurons. PLoS Comput. Biol. 7:e1002038. doi: 10.1371/journal.pcbi.1002038

PubMed Abstract | CrossRef Full Text | Google Scholar

Rutter, G. A., and Hodson, D. J. (2013). Minireview: intraislet regulation of insulin secretion in humans. Mol. Endocrinol. 27, 1984–1995. doi: 10.1210/me.2013-1278

PubMed Abstract | CrossRef Full Text | Google Scholar

Rutter, G. A., and Hodson, D. J. (2015). Beta cell connectivity in pancreatic islets: a type 2 diabetes target? Cell. Mol. Life Sci. 72, 453–467. doi: 10.1007/s00018-014-1755-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Rutter, G. A., Hodson, D. J., Chabosseau, P., Haythorne, E., Pullen, T. J., and Leclerc, I. (2017). Local and regional control of calcium dynamics in the pancreatic islet. Diabetes, Obes. Metab. 19, 30–41. doi: 10.1111/dom.12990

PubMed Abstract | CrossRef Full Text | Google Scholar

Santos, R. M., Rosario, L. M., Nadal, A., Garcia-Sancho, J., Soria, B., and Valdeolmillos, M. (1991). Widespread synchronous [Ca²⁺]_i oscillations due to bursting electrical activity in single pancreatic islets. Pflugers Arch. 418, 417–422. doi: 10.1007/BF00550880

PubMed Abstract | CrossRef Full Text | Google Scholar

Satin, L. S., Butler, P. C., Ha, J., and Sherman, A. S. (2015). Pulsatile insulin secretion, impaired glucose tolerance and type 2 diabetes. Mol. Aspects Med. 42, 61–77. doi: 10.1016/j.mam.2015.01.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmitz, O., Rungby, J., Edge, L., and Juhl, C. B. (2008). On high-frequency insulin oscillations. Ageing Res. Rev. 7, 301–305. doi: 10.1016/j.arr.2008.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Schuit, F. C., In ', P. A., Veldt, T., and Pipeleers, D. G. (1988). Glucose stimulates proinsulin biosynthesis by a dose-dependent recruitment of pancreatic beta cells (protein synthesis/beta cell heterogeneity/insulin). Proc. Natl. Acad. Sci. U.S.A. 85, 3865–3869.

Google Scholar

Skelin Klemen, M., Dolenšek, J., Slak Rupnik, M., and Stožer, A. (2017). The triggering pathway to insulin secretion: functional similarities and differences between the human and the mouse β cells and their translational relevance. Islets 9, 109–139. doi: 10.1080/19382014.2017.1342022

PubMed Abstract | CrossRef Full Text | Google Scholar

Smolen, P., Rinzel, J., and Sherman, A. (1993). Why pancreatic islets burst but single beta cells do not. The heterogeneity hypothesis. Biophys. J. 64, 1668–1680. doi: 10.1016/S0006-3495(93)81539-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Speier, S., and Rupnik, M. (2003). A novel approach to in situ characterization of pancreatic beta-cells. Pflugers Arch. 446, 553–558. doi: 10.1007/s00424-003-1097-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Stamper, I. J., Jackson, E., and Wang, X. (2014). Phase transitions in pancreatic islet cellular networks and implications for type-1 diabetes. Phys. Rev. E 89, 12719–12715. doi: 10.1103/PhysRevE.89.012719

PubMed Abstract | CrossRef Full Text | Google Scholar

Stožer, A., Dolenšek, J., Rupnik, M. S., Arkhammar, P., and Thastrup, O. (2013a). Glucose-stimulated calcium dynamics in islets of langerhans in acute mouse pancreas tissue slices. PLoS ONE 8:e54638. doi: 10.1371/journal.pone.0054638

PubMed Abstract | CrossRef Full Text | Google Scholar

Stožer, A., Gosak, M., Dolenšek, J., Perc, M., Marhl, M., Rupnik, M. S., et al. (2013b). Functional connectivity in islets of langerhans from mouse pancreas tissue slices. PLoS Comput. Biol. 9:e1002923. doi: 10.1371/journal.pcbi.1002923

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, Y., Quyen, T. L., Hung, T. Q., Chin, W. H., Wolff, A., and Bang, D. D. (2015). A lab-on-a-chip system with integrated sample preparation and loop-mediated isothermal amplification for rapid and quantitative detection of Salmonella spp. in food samples. Lab. Chip 15, 1898–1904. doi: 10.1039/C4LC01459F

PubMed Abstract | CrossRef Full Text | Google Scholar

Tengholm, A., and Gylfe, E. (2009). Oscillatory control of insulin secretion. Mol. Cell. Endocrinol. 297, 58–72. doi: 10.1016/j.mce.2008.07.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Tetzlaff, C., Okujeni, S., Egert, U., Wörgötter, F., and Butz, M. (2010). Self-organized criticality in developing neuronal networks. PLoS Comput. Biol. 6:e1001013. doi: 10.1371/journal.pcbi.1001013

PubMed Abstract | CrossRef Full Text | Google Scholar

Tomen, N., Rotermund, D., and Ernst, U. (2014). Marginally subcritical dynamics explain enhanced stimulus discriminability under attention. Front. Syst. Neurosci. 8:151. doi: 10.3389/fnsys.2014.00151

PubMed Abstract | CrossRef Full Text | Google Scholar

Trefois, C., Antony, P. M., Goncalves, J., Skupin, A., and Balling, R. (2015). Critical transitions in chronic disease: transferring concepts from ecology to systems medicine. Curr. Opin. Biotechnol. 34, 48–55. doi: 10.1016/j.copbio.2014.11.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Valdeolmillos, M., Gomis, A., and Sainchez-Andres, J. V. (1996). in vivo synchronous membrane potential oscillations in mouse pancreatic beta-cells: lack of co-ordination between islets M. J. Physiol. 493, 9–18. doi: 10.1113/jphysiol.1996.sp021361

CrossRef Full Text | Google Scholar

Valdeolmillos, M., Santos, R. M., Contreras, D., Soria, B., and Rosario, L. M. (1989). Glucose-induced oscillations of intracellular Ca²⁺ concentration resembling bursting electrical activity in single mouse islets of Langerhans. FEBS Lett. 259, 19–23. doi: 10.1016/0014-5793(89)81484-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Westacott, M. J., Ludin, N. W. F., and Benninger, R. K. P. (2017). Spatially organized β-Cell subpopulations control electrical dynamics across islets of langerhans. Biophys. J. 113, 1093–1108. doi: 10.1016/j.bpj.2017.07.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, S. N., Collins, T., and Jančovič, P. (2012). Hybrid method for designing digital butterworth filters. Comput. Electr. Eng. 38, 811–818. doi: 10.1016/j.compeleceng.2012.03.015

CrossRef Full Text | Google Scholar

Yaroslavsky, L. P., Egiazarian, K. O., and Astola, J. T. (2001). “Transform domain image restoration methods: review, comparison, and interpretation,” in Nonlinear Image Processing and Pattern Analysis XII, eds E. R. Dougherty and J. T. Astola (San Jose, CA, United States), 155–169.

Google Scholar

Zarkovic, M. (2004). Synchronization and entrainment of cytoplasmic Ca²⁺ oscillations in cell clusters prepared from single or multiple mouse pancreatic islets. Am. J. Physiol. Endocrinol. Metab. 287, E340–E347. doi: 10.1152/ajpendo.00069.2004

CrossRef Full Text | Google Scholar

Zimliki, C. L., Mears, D., and Sherman, A. (2004). Three roads to islet bursting: emergent oscillations in coupled phantom bursters. Biophys. J. 87, 193–206. doi: 10.1529/biophysj.103.038471

PubMed Abstract | CrossRef Full Text | Google Scholar