Self-organized critical dynamics is assumed to be an attractive mode of functioning for several real-life systems and entails an emergent activity in which the extent of observables follows a power-law distribution. The hallmarks of criticality have recently been observed in a plethora of biological systems, including beta cell populations within pancreatic islets of Langerhans. In the present study, we systematically explored the mechanisms that drive the critical and supercritical behavior in networks of coupled beta cells under different circumstances by means of experimental and computational approaches. Experimentally, we employed high-speed functional multicellular calcium imaging of fluorescently labeled acute mouse pancreas tissue slices to record calcium signals in a large number of beta cells simultaneously, and with a high spatiotemporal resolution. Our experimental results revealed that the cellular responses to stimulation with glucose are biphasic and glucose-dependent. Under physiological as well as under supraphysiological levels of stimulation, an initial activation phase was followed by a supercritical plateau phase with a high number of global intercellular calcium waves. However, the activation phase displayed fingerprints of critical behavior under lower stimulation levels, with a progressive recruitment of cells and a power-law distribution of calcium wave sizes. On the other hand, the activation phase provoked by pathophysiologically high glucose concentrations, differed considerably and was more rapid, less continuous, and supercritical. To gain a deeper insight into the experimentally observed complex dynamical patterns, we built up a phenomenological model of coupled excitable cells and explored empirically the model’s necessities that ensured a good overlap between computational and experimental results. It turned out that such a good agreement between experimental and computational findings was attained when both heterogeneous and stimulus-dependent time lags, variability in excitability levels, as well as a heterogeneous cell-cell coupling were included into the model. Most importantly, since our phenomenological approach involved only a few parameters, it naturally lends itself not only for determining key mechanisms of self-organized criticality at the tissue level, but also points out various features for comprehensive and realistic modeling of different excitable systems in nature.

Introduction

Self-organized collective dynamics is a remarkable phenomenon observed in various natural and man-made systems, in which collective behavior emerges from local interactions between individual elements (Bak, 1996; Marković and Gros, 2014; Ellis and Kopel, 2019). Regardless of the specific mechanisms responsible for self-organization, the resulting coherent global structures or dynamics are characterized by scale-invariant properties and a power-law distribution of systems’ observables (Khaluf et al., 2017; Muñoz, 2018). Such emergent behavior is often associated with critical dynamics and is assumed to be particularly beneficial for the functioning of several living systems, from the microscopic to the macroscopic (Tamayo et al., 1999; Nykter et al., 2008; Chialvo, 2010; Bialek et al., 2012; Furusawa and Kaneko, 2012; Sasai, 2013; Allegrini et al., 2015; Muñoz, 2018). Criticality has been argued to originate from the fact that many biological systems operate in the vicinity of a critical point of a phase transition between an ordered and disordered phase, which ensures a balance between robustness against perturbations and flexibility to adapt to a changing environment. However, the exact reasons why signatures of criticality can be conjectured to emerge in living systems are still under debate and the underlying principles are incompletely understood (Lovecchio et al., 2012; Moretti and Muñoz, 2013; Nonnenmacher et al., 2017). Most importantly, despite some skepticism and limitations, the evidently increasing amount of empirical evidence, fostered also by technological and computational advances, is nowadays inspiring more and more researchers to investigate the complexity of biological systems through the lens of phase transition behavior and criticality.

In the domain of biological networks, the concepts of self-organization and criticality have received the most attention in the field of neuroscience. On the smallest scales, patterns of activity in neuronal populations have been found to be very heterogeneous, with sizes of so-called neuronal avalanches following a power law distribution (Beggs and Plenz, 2003; Pasquale et al., 2008; Timme et al., 2016). Empirical evidence for criticality has been reported in both different in vivo preparations and on larges scales of whole-brain imaging (Plenz and Thiagarajan, 2007; Haimovici et al., 2013; Hesse and Gross, 2014). The presence of emergent critical dynamics in the nervous system is theoretically appealing and consequently computational models and tools from the realms of statistical physics have been utilized to unveil the mechanisms and correlations between phase transition behavior and the occurrence of scale-invariant neuronal avalanches (Plenz and Thiagarajan, 2007; Rubinov et al., 2011; Friedman et al., 2012; Zare and Grigolini, 2013; Tkačik et al., 2015; Brochini et al., 2016; di Santo et al., 2018). Importantly, several studies have underscored the emergence of critical dynamics in neuronal networks as one of the key pillars for their optimal operational abilities (Kinouchi and Copelli, 2006; Shew and Plenz, 2013; Shew et al., 2015; Stoop and Gomez, 2016). Moreover, complex and hierarchically organized network structures along with neuronal plasticity were identified as the main neurophysiological determinants that ensure robust critical behavior (Levina et al., 2007; Rubinov et al., 2011; Moretti and Muñoz, 2013; Hutt et al., 2014; Massobrio et al., 2015b). It should be noted that deviations from critical behavior occur in neuronal networks during development (Tetzlaff et al., 2010; Massobrio et al., 2015b) and under pathological conditions (Massobrio et al., 2015a; Tagliazucchi et al., 2016; Hahn et al., 2017). Especially during epileptic seizures (Hobbs et al., 2010; Meisel et al., 2012) or by pharmacological disruptions of the excitation-inhibition balance (Barral and D Reyes, 2016), an excess of large system-spanning avalanches occur, as is characteristic for supercritical dynamical states. Consequently, it has been hypothesized that the healthy brain resides near a critical or even slightly subcritical state, thereby ensuring a safety margin from supercriticality, which has been linked to some pathophysiological disorders (Priesemann et al., 2014; Tomen et al., 2014; Massobrio et al., 2015a).

Notably, recent research indicates that the concept of critical dynamics and power-law scaling in living beings applies well beyond the spatiotemporal activity patterns of neurons. At the (sub)cellular level, mitochondrial network of heart myocytes was reported to operate at the edge of dynamic instability characterized by a fractal scaling of depolarized mitochondrial clusters (Aon et al., 2004). In this regime, constancy in terms of a steady supply of ATP is provided in combination with flexibility, which ensures the adaptation of energy production in accordance with metabolic demands (Aon et al., 2006). Moreover, hallmarks of self-organized criticality have also been observed in the spatiotemporal organization of Ca²⁺ waves. Jung et al. (1998) reported a power law distribution of noise-induced spiral Ca²⁺ wave sizes in cultured networks of astrocytes. Heavy-tailed distributions and an avalanche-like behavior have also been observed in intracellular Ca²⁺ signalization in cardiac myocytes (Nivala et al., 2012) and in immature oocytes (Lopez et al., 2012). In both studies Ca²⁺ waves in individual cells resulted from random local Ca²⁺ events, reflecting small Ca²⁺ release events from individual channels or a cluster of channels, which can occasionally integrate to global events reflecting a whole-cell Ca²⁺ signal (Berridge et al., 2000). As the localized subcellular Ca²⁺ events interact, e.g., via diffusion, they can self-organize and lead to avalanches of activity that propagate through the cell. The concept does not only assert that the extent of such events is characterized by scale invariance, but also makes the global Ca²⁺ signals appear rather deterministic in spite of their stochastic origin (Skupin et al., 2008).

Even though information processing in living organisms is often performed by large networks of interacting cells, little attention has been devoted to the principles underlying critical dynamics on multicellular and tissue levels of organization. Recently, we have empirically shown that fingerprints of criticality are also found in the spatiotemporal dynamics of interconnected beta cells from islets of Langerhans (Gosak et al., 2017). These endocrine cells synthesize and release insulin, the anabolic hormone which promotes postprandial storage of nutrients, and thus serves a crucial role in homeostasis of energy that becomes disrupted in diabesity (Kahn et al., 2014). Insulin concentration in the blood displays inherent multimodal oscillations (Satin et al., 2015), and several studies have attempted to reveal the underlying mechanism by providing links to oscillations in beta cells metabolism and to a feedback between ion channels and electrical activity. This was further corroborated by modeling the interplay of the two signaling aspects (Bertram et al., 2007). Moreover, recent theoretical studies emphasized the role of biphasic feedback circuits in controlling functional beta cell mass (Karin et al., 2016) and progression of diabetes mellitus (Karin and Alon, 2017). On the organizational level of a single islet, beta cells respond to nutrient stimulation with an initial transient depolarization, followed by fast oscillations in membrane potential that are superimposed on a plateau phase (Gilon and Henquin, 1992; Rorsman and Braun, 2013; Skelin Klemen et al., 2017). Intracellular Ca²⁺ concentration ([Ca²⁺]_IC) closely follows changes in membrane potential due to tight coupling between electrical and calcium dynamics in beta cells (Gilon and Henquin, 1992; Dolenšek et al., 2013). However, all beta cells within an islet do not show identical electrical or [Ca²⁺]_IC activity and can therefore not be regarded as uniformly and strongly coupled identical units or even as a single supercell. Rather, the collective activity of beta cells is characterized by a phase shift between individual cellular oscillations, ultimately resulting in heterogeneous waves of membrane potential and [Ca²⁺]_IC changes. These waves spread repetitively over an islet, but not always from the same source and not always throughout the whole syncytium (Benninger et al., 2008; Dolenšek et al., 2013; Stožer et al., 2013a). They are thought to originate in specific sub-regions with elevated excitability (Benninger et al., 2014) or higher intrinsic oscillation frequency (Westacott et al., 2017). A plethora of evidence demonstrates that an essential prerequisite for the coordinated beta cell activity and formation of waves is intact intercellular connectivity mediated via gap junctions (Calabrese et al., 2003; Ravier et al., 2005; Bavamian et al., 2007; Bosco et al., 2011; Benninger et al., 2014) and probably other modes of communication, such as paracrine, contact-dependent, and ciliary signaling (Squires et al., 2002; Konstantinova et al., 2007; Yang et al., 2011; Gerdes et al., 2014). Most importantly, intercellular connectivity is not only necessary for normal islet function, its perturbations were also linked to metabolic diseases and impaired insulin secretion (Hamelin et al., 2009; Carvalho et al., 2012; Head et al., 2012; Hodson et al., 2013; Benninger and Piston, 2014; Benninger et al., 2018; Nasteska and Hodson, 2018).

Moreover, individual beta cells are intrinsically highly heterogeneous (Gutierrez et al., 2017). Several different approaches have demonstrated relatively large differences in the extent of coupling between beta cells (Pérez-Armendariz et al., 1991; Farnsworth et al., 2014), as well as different levels of excitability (Jonkers and Henquin, 2001; Benninger et al., 2011) and rates of glucose metabolism (Benninger et al., 2014; Benninger and Hodson, 2018; Nasteska and Hodson, 2018). Because of these features, the spatiotemporal responses of beta cells are very complex and to understand how a population of these heterogeneous and heterogeneously coupled cells activate to work in synchrony is a hot topic in islet physiology research (Pedersen et al., 2013; Benninger and Piston, 2014; Benninger et al., 2014; Markovič et al., 2015; Cappon and Pedersen, 2016). Motivated by complexity science approaches, several studies have investigated beta cell responses in terms of phase transition behavior (Hraha et al., 2014; Loppini et al., 2014; Stamper et al., 2014). In this vein, in our recent study we demonstrated that under physiological circumstances, the initial response to glucose is characterized by a power-law probability distribution of Ca²⁺ wave sizes, which can be maintained in the long run by periodic stimulation, but changes to supercriticality upon constant stimulation, thereby demonstrating empirically the fingerprints and basic preconditions of critical behavior (Gosak et al., 2017).

In the present work, we extend our preceding research to supraphysiological glucose concentrations that are usually used in experiments, but accompany pathophysiological states in vivo. It turns out that higher glucose levels evoke more rapid and qualitatively different beta cell responses when compared to physiological levels of stimulation. To assess the measured non-trivial and rich dynamical patterns, we propose a phenomenological model of coupled excitable cells that accounts for the observed physiological as well as pathophysiological behavior and encompasses both beta cell signaling specifics and heterogeneity. Moreover, in contrast to the exhaustive computational models with many parameters, our phenomenological modeling approach made it easier to empirically explore the necessary ingredients and physiological determinants that ensure a good overlap between experimental and computational results.

Materials and Methods

Multicellular Calcium Imaging in Pancreatic Tissue Slices

Acute pancreatic tissue slices were prepared as described previously (Speier and Rupnik, 2003; Stožer et al., 2013a). Briefly, low-melting point agarose (1.9% V/V) was injected into the proximal common bile duct that was clamped distally at the major duodenal papilla. Retrograde inflow of the agarose served, once cooled, as a scaffold for subsequent tissue cutting into 140 μm thick slices on a vibratome (VT 1000 S, Leica). Staining with the calcium sensitive dye Oregon Green 488 BAPTA-1 AM [6 μM final concentration, 0.03% Pluronic F-127 (w/v), and 0.12% dimethylsulphoxide (v/v) dissolved in HEPES-buffered saline at RT, consisting of (in mM) 150 NaCl, 10 HEPES, 6 glucose, 5 KCl, 2 CaCl₂, 1 MgCl₂; titrated to pH = 7.4 using 1 M NaOH], was performed for 50 min at RT. Confocal imaging of calcium dynamics was performed on the Leica TCS SP5 AOBS Tandem II upright confocal system (20x HCX APO L water immersion objective, NA 1.0) and Leica TCS SP5 DMI6000 CS inverted confocal system (20X HC PL APO water/oil immersion objective, NA 0.7), utilizing a perifusion system filled with extracellular solution, consisting of (in mM) 125 NaCl, 26 NaHCO₃, 6 lactic acid, 3 myo-inositol, 2.5 KCl, 2 Na-pyruvate, 2 CaCl₂, 1.25 NaH₂PO₄, 1 MgCl₂, 0.5 ascorbic acid and added either substimulatory 6 mM or stimulatory 8 or 12 mM glucose. The calcium dye was excited at 488 nm via an argon laser line and the emitted fluoresces was detected in the range of 500–700 nm by a Leica HyD detector. Time series were acquired at 10 Hz (512 × 512 pixels). Further analysis was done off-line by manually selecting ROIs corresponding to beta cells. The exported time series of the F/F₀ ratio were then further processed as explained in the continuation.

Computational Model

Phenomenological Single Cell Model

We utilized a phenomenological model to describe the dynamics of the electrically excitable beta cells. In particular, we made use of a two-dimensional iterated map proposed by Rulkov (2002):

$$u_i(t+1)={\mathrm{\alpha }}_i(t)\mathrm{/}(1+u_i{(t)}^2)+v_i(t)+g_i\sum _j{\mathrm{\epsilon }}_{ij}(u_j-u_i)+\mathrm{\beta }{\xi }_i(t),(1)$$

$$v_i(t+1)=v_i(t)-\mathrm{\sigma }{u}_i(t)-\mathrm{\chi },(2)$$

where u_i(t) and v_i(t) are the slow and the fast dynamical variables for the i-th cell, respectively, and are considered as dimensionless variables, t is the discrete time index, α_i, χ, and σ are systems parameters, and β = 0.0045 defines the strength of Gaussian noise with zero mean and unit variance that accounts for stochasticity in beta cell dynamics. The fast variable u_i(t) describes the dynamics of the membrane potential of the cell, whereas the slow variable v_i(t) reflects the gating variable. Although this is an abstract and simple mathematical model, it mimics well the basic principles of more complex cellular behaviors that are observed in different cell types, including beta cells.

More specifically, in terms of metabolic changes, electrical activity, [Ca²⁺]_IC dynamics, and insulin secretion, beta cells within islets respond to stimulation by glucose in two phases. During the first phase which is transient, they show elevated levels of intracellular ATP and NAD(P)H, followed electrically by very fast bursting or continuous bursting. Bursts are periods of very fast depolarizations called spikes that last a few seconds and continuous bursting consists of an uninterrupted set of spikes at a frequency around 10 Hz. At the level of [Ca²⁺]_IC dynamics, this first phase consists of a transient increase in [Ca²⁺]_IC, during the ascending part of which a few fast [Ca²⁺]_IC oscillations may be present, reflecting fast burst before continuous bursting. It should be pointed out that at present the very fast spikes cannot be resolved in [Ca²⁺]_IC imaging. This first phase lasts a few minutes and at the level of hormone output overlaps with the first phase of insulin secretion. During the second phase, continuous bursting and the accompanying transient increase in [Ca²⁺]_IC change to regular bursting and corresponding fast [Ca²⁺]_IC oscillations. At the level of metabolism and hormone secretion, NADPH and ATP are elevated during this period and insulin secretion shows a stable second phase (Henquin and Meissner, 1984; Gilon and Henquin, 1992; Li et al., 2013, 2014; Gilon et al., 2014; Skelin Klemen et al., 2017; Rorsman and Ashcroft, 2018). It should be noted that during this second phase, insulin is also secreted in bursts synchronized with electrical bursts and fast [Ca²⁺]_IC oscillations (Gilon et al., 1993; Bergsten et al., 1994). This fast electrical, [Ca²⁺]_IC, and secretory pattern is superimposed on a slower set of oscillations in ATP, membrane potential, [Ca²⁺]_IC, insulin secretion, and some other parameters, which has been reviewed in detailed elsewhere (Satin et al., 2015). A further layer of complexity to this behavior comes from the fact that it is glucose-dependent. In higher glucose, the frequency or duration of bursts and correspondingly the fast [Ca²⁺]_IC oscillations increase, such that the active time and insulin secretion increase. Noteworthy, it seems that the underlying slow pattern is not glucose-dependent (Satin et al., 2015; Gosak et al., 2017; Skelin Klemen et al., 2017).

Most importantly, modeling all of the above aspects of beta cell responses to glucose requires the use of realistic biophysical models. However, in this study we focused only on fast [Ca²⁺]_IC dynamics which can be satisfactorily captured by the phenomenological model employed here. Most importantly, since we wanted to study the effects of various types of heterogeneities in a network of coupled beta cells, in comparison with a genuine biophysical model, a phenomenological description of the complex cellular dynamics is not only beneficial in terms of numerical efficiency, but also enables exploration of the system with very few free parameters.

The Rulkov map displays a variety of dynamics depending on the parameter choice, as extensively investigated in the past (Rulkov, 2002; Ibarz et al., 2011; Markovič et al., 2012). To better understand the dynamical phases that occur in our study, we performed a stability analysis. For χ = σ the fixed point equals u^* = −1 and v^* = − 1 − α/2. If α < 2 the steady state is stable and for α < 1.86 the fixed point is asymptotically stable, since the both eigenvalues, λ₁ and λ₂ have only real parts and their absolute value is less than 1 (see Figure 1A). For the values 1.86 < α < 2, the fixed point is still stable (|Re(λ₁)| < 1 and |Re(λ₂)| < 1), but the eigenvalues are complex making the fixed point a spiral sink. For the values of bifurcation parameter 2 < α ≤ 4, the solution (u^*,v^*) becomes unstable and the system exhibits sustained periodic pulses, chaotic bursts of pulses and sustained chaotic pulsing. In our study we focused on the region 1.86 < α < 2, where the steady state is excitable and oscillations can be induced by noise and/or heterogeneity. This is shown in the bifurcation diagram in Figure 1B. For the chosen noise level, oscillations occur at α > 1.93. Noteworthy, with increasing α the excitability level and the cellular activity increase as well. Therefore, increasing α in our model emulates the decrease in glucose-induced K_ATP-channel conductance, the main trigger of beta cells in realistic models (Stamper and Wang, 2019). The temporal behavior of our single-cell phenomenological model is visualized in Figure 1C, for different values the bifurcation parameter α.

FIGURE 1

Figure 1. Dynamical features of the Rulkov map. (A) Real and imaginary eigenvalues λ₁ and λ₂ for different values of the control parameter α. (B) Bifurcation diagram of the fast variable with (blue) and without (green) added noise. (C) Traces of the fast variable with (blue) and without (green) noise for different values of excitability levels α.

Intercellular Coupling Model

The sum in Eq. (1) signifies the electrical coupling and it runs over all cells, whereby ε_ij = 1 if the unit i is coupled to unit j, whilst otherwise ε_ij = 0. g_i is the coupling constant. The structure of the intercellular coupling between beta cells was modeled by the random geometric graph model (Penrose, 2007). First, all N=200 cells were arranged randomly in a unit square with a prescribed minimal possible distance (0.04) to ensure a more homogeneous and realistic spatial distribution of cells. Then, the i-th and the j-th cell were considered to be connected, i.e., ε_ij = 1, if their physical distance was less than $r_{ij}=\sqrt{⟨k⟩/\left(N\mathrm{\pi }\right)}$, where ⟨k⟩ = 6 signifies the average number of connections per cell. A typical intercellular network structure is shown in Figure 2A.

FIGURE 2

Figure 2. Features of the phenomenological model of beta cell population. (A) A typical simulated beta cell network architecture. Red dots denote individual cells and the arrows depict intercellular electrical coupling. (B) Simulated time course of beta cell excitability rate after switching from substimulatory to stimulatory (Δα = 0.08, A = 0.45, B = 0.0004, $T_{\mathrm{m},i}^{(8)}=35000$, t_s = 10000) and suprastimulatory (Δα = 0.09, A = 0.70, B = 0.0008, $T_{\mathrm{m},i}^{(12)}=20000$, t_s = 10000) conditions.

Modeling the Temporal Responses to Stimulation

To simulate the progressive recruitment of beta cells after switching from substimulatory to stimulatory or suprastimulatory glucose concentrations, we introduced a time-dependent function for the parameter α_i that reflects the cellular excitability level:

$${\mathrm{\alpha }}_i(t)={\mathrm{\alpha }}_0+\mathrm{\Delta }{\mathrm{\alpha }}_i\left[A\left(t-t_s\right)B{e}^{-B\left(t-t_s\right)+1}+\frac{{\left(t-t_s\right)}^2}{{\left(t-t_s\right)}^2+{\left(T_{\mathrm{m},i}-t_s\right)}^2}\right],(3)$$

In this manner, we took into account the delay due to glucose metabolism. In Eq. (3) α₀ = 1.90 is the basal substimulatory level of excitability with no activity, Δα is the amplitude of the increased excitability for the i-th cell provoked by increased glucose concentration, and t_s is the initial time before the cells respond to stimulation. The first term within the brackets on the right side of Eq. (3) stands for the initial and the second term for the successive beta cell activations. We implemented such a biphasic and glucose-dependent response to account for the previously observed biphasic and glucose-dependent behavior of beta cells in terms of their metabolic, electrical, [Ca²⁺]_IC, and secretory response described above. The parameters A ∈ [0,1] and B signify the glucose-dependent amplitude and decay rate of the first activation. The parameter T_m,i specifies the temporal scale of the final activation, i.e., elevation in the level of excitability. On the basis of experimental results we hypothesized that under lower and physiological stimulatory conditions the amplitude of the first response and the decay rate are lower (A = 0.45 and B = 0.0004) than under high and supraphysiological stimulatory levels (A = 0.7 and B = 0.0008). Moreover, to account for the longer activation phase observed under 8 mM glucose in comparison to 12 mM stimulation, we set the half-activation times to $T_{\mathrm{m},i}^{(8)}=35000$ and $T_{\mathrm{m},i}^{(12)}=20000$. Finally, to resemble a higher intrinsic beta cell activity under higher stimulation, we set the parameters when emulating the behavior under 8 mM glucose to Δα = 0.08 and σ = χ = 0.001 and to Δα = 0.09 and σ = χ = 0.0012 when emulating the behavior under 12 mM glucose. It should be noted that these small changes in the parameters σand χ have an insignificant effect on bifurcation behavior reported in Figure 1. Temporal traces of simulated excitability rates when switching to stimulatory and suprastimulatory conditions are shown in Figure 2B. Since the parameter α regulates the cellular activity (see Figure 1), by this means a stimulation-specific temporal recruitment of beta cells is modeled.

Heterogeneity of Beta Cells

Previous studies have suggested an extensive heterogeneity among β cells due to differences in topography, cell sizes, functional maturity, channel densities, intercellular coupling, rates of glucose metabolism, membrane potential changes, [Ca²⁺]_IC oscillations, granule content, and secretory capacity, to name only a few examples (MacDonald and Rorsman, 2006; Benninger and Piston, 2014; Bader et al., 2016; Roscioni et al., 2016; Pipeleers et al., 2017; Skelin Klemen et al., 2017; Benninger and Hodson, 2018; Nasteska and Hodson, 2018). To robustly account for the abovementioned differences in glucose sensitivity and metabolism, electrical excitability and [Ca²⁺]_IC signals, as well as intercellular coupling strength, we introduced in our phenomenological model heterogeneity three crucial aspects of cellular signaling: (i) stimulation-induced temporal change in excitability (parameter T_m,i), (ii)stimulation-dependent increase in excitability (parameter Δα_i), and (iii) intercellular coupling strength (parameter g_i). All three parameters were assumed to follow a truncated normal distribution with a relative standard deviation of 30% and a cut-off of 90%. The three types of cellular heterogeneity are schematically visualized in Figure 3.

FIGURE 3

Figure 3. A schematic representation of all three types of heterogeneities in the model. (A) electrical excitability Δα_i, (B) beta cell metabolism T_m,i, (C) electrical coupling g_i. Widths of connections reflect the coupling strength.

Processing of Time Series and Activity Pattern Classification

Time series of individual cells obtained from experimental recordings were first accordingly processed to achieve a coherent and accurate binarization. The main task of this pre-processing step is to level and smooth the individual time series, remove noise, and firmly extract the fast component of Ca²⁺ oscillations. To this purpose, we utilized a band pass filter, whereby the frequency band of interest was determined by visual assessment. The filtered signal was then additionally smoothed with standard sliding window algorithm, with a window size of four frames (Yaroslavsky et al., 2001). Preprocessing of computationally obtained traces was not required. In continuation the time series from experiments and simulations will be referred to as x and the corresponding value at (discrete) time t as x(t). The following binarization procedure was based on the: (i) standard deviation std [x], (ii) first derivative of time series x′⁢(t), and (iii) standard deviations of its the first derivative std [x′]. First, we have defined the potential onset and ending times of individual Ca²⁺ spikes by searching for local extremes in the first derivative. More precisely, our algorithm searches for local maxima’s, which satisfy the condition x′ (t) > 1.5std [x′]. The time, at which the local maxima is found is then the potential onset of Ca²⁺ spikes, t_START. In a time forward direction, we than seek the first local minima that satisfies the condition x′ (t) < −1.5std [x′] and store its time of occurrence, t_END. Lastly we test, if the local maximum of the Ca²⁺ signal within the time interval t ∈ [t_START,t_END] satisfies the condition x(t) > 1.5std [x]. The corresponding binary time series takes on a value of 1 in all time intervals t ∈ [t_START,t_END], where the three conditions are satisfied, whilst otherwise the value is 0.

Afterward, we used binarized time series and the physical positions of individual cells to merge spatially and temporally synchronized events into clusters by performing the space-time cluster analysis, as proposed by Jung (1997) and Jung et al. (1998). In brief, we combined the information about the positions of cells and their binary traces into a space-time cube (STC). In this STC we defined a cubic region of interest (STC-ROI) in which we search for active cells. If two cells in neighboring STC were simultaneously active, they were considered to belong to the same cluster. In other words, we traced the course of the wave from cell-to-cell and if the nearby cells became activated within a short time period and if they were close enough, the given activation was considered as one individual cluster with size p. By this means, an individual STC contains the information about the number of cells that were activated in a given excitation wave, as well as about the temporal extent of the given event, as described previously (Gosak et al., 2017). The spatial side-length of the STC was determined as the average distance to 6 closest neighbors (∼25 μm in experiments and ∼0.12 in simulations). The temporal side-length was determined empirically, so that a firm of discrimination of individual waves was attained. To quantify the spatiotemporal activity patterns in experiments and simulations, we calculated the distribution of cluster sizes N(p) for different stimulation protocols and activity phases. Finally, the results were fitted with a power-law function to qualitatively evaluate the nature of the distribution, i.e., critical vs. supercritical behavior. In particular, by visually assessing deviations from the power-law distribution in the form of an excess of global events we determine the supercritical nature of the spatio-temporal activity, whereas a close-to-power-law behavior implies critical-like behavior, as suggested previously (Levina et al., 2007; Friedman et al., 2012).

Results

First, we present experimentally measured beta cells activity after switching from substimulatory to stimulatory and suprastimulatory levels of glucose. Then, we show the results of our computational model of interconnected excitable cells, which was designed to mimic the activity patterns observed in experiments under physiological as well as under supraphysiological levels of stimulation. Data either from experiments or from simulations were handled in the same manner, in order to provide foundation for further characterization of the spatiotemporal activity.

Experimental Results

To record beta cell responses to glucose stimulation, we used multicellular confocal imaging on acute tissue slices as described in Materials and methods. We stimulated islets with two glucose concentrations: one commonly observed in vivo, i.e., 8 mmol/l, and one measured in conditions of stress or glucose intolerance, i.e., 12 mmol/l. We termed the former physiological and the later supraphysiological concentration. Following either stimulus, beta cells exhibited a two-phase response: (i) an activation phase, characterized by a transient increase in [Ca²⁺]_IC and presence of fast oscillations, during which beta cells were gradually recruited, and (ii) a subsequent plateau phase, characterized by repeated and more regular oscillations of now fully recruited beta cells (Figures 4A, 5A). Heterogeneity of beta cells responses was reflected in the time window during which cells activate within an islet. These intervals differ for the two stimulatory concentrations: about 600 s (100 s < t < 700 s, Figure 4A) for the physiological and about 300 s (150 s < t < 450 s, Figure 5A) for the supraphysiological concentration. To surpass the qualitative description of the two phases, we looked for collective spatiotemporal behavior of beta cells. To this aim, we meticulously detected Ca²⁺ waves during both phases, and plotted them as individual events in space-time for better visualization. While being stimulated with the physiological concentration, the activation phase exhibited very heterogeneous spatiotemporal behavior, one that resulted in calcium waves of very different sizes (Figure 4B). The following plateau phase evoked a more regular pattern of [Ca²⁺]_IC oscillations, with prevailing global intercellular calcium waves that encompassed often the majority of the cells within an islet (Figure 4D). However, the activation/plateau pattern changed during supraphysiological stimulation. Majority of the cells responded with a rapid burst of oscillatory activity followed by brief refractory period during the activation phase (Figure 5B). The subsequent plateau phase was dominated by global intercellular calcium waves (Figure 5D). To be able to quantify the former description, we determined the distribution P(s) of relative wave sizes s and plotted it on log-log scale (Figures 4C,E, 5C,E). While comparing P(s) for the two concentrations in question, we observed that in the lower concentration the P(s) in the activation phase followed the power law, whereas the plateau phase was again dominated by global waves. Such switching in behavior from the critical to the supercritical was not observed in the higher stimulatory concentration, during which the behavior was locked to the supercritical during both phases. Namely, the activation phase under supraphysiological concentrations was too rapid, exhibited a huge activation burst, and lacked the progressive recruitment of cells that featured an emergent behavior with very heterogeneous Ca²⁺ waves.

FIGURE 4

Figure 4. Experimentally measured beta cell responses after stimulation with 8 mM glucose. (A) Three characteristic Ca²⁺ traces and the raster plot of binarized Ca²⁺ activity of all cells in the islet. The orange dotted line indicates the fraction of active cells within the given time-window that was slid throughout the recording. (B,D) 3D raster plots showing the Ca²⁺ activity waveforms for selected intervals for the activation (B) and plateau (D) phase. Colors denote specific Ca²⁺ events. Gray dots on the y, z plane stand for coordinates of cells. (C,E) The distributions of Ca²⁺ wave sizes for the activation (C) and plateau (E) phase. The gray dashed line indicates the power-law fit. The slope in the critical-like activation phase is –1.69.

FIGURE 5

Figure 5. Experimentally measured beta cell responses after stimulation with 12 mM glucose. (A) Three characteristic Ca²⁺ traces and the raster plot of binarized signals of Ca²⁺ oscillations in all cells in the islet. The orange dotted line indicates the fraction of active cells within the given time-window that was slid throughout the recording. (B,D) 3D raster plots showing the Ca²⁺ waves for selected intervals for the activation (B) and plateau (D) phase. Colors denote specific Ca²⁺ events. Gray dots on the y, z plane stand for coordinates of cells. (C,E) The distributions of Ca²⁺ wave sizes for the activation (C) and plateau (E) phase. The gray dashed line indicates the power-law fit.

Computational Results

We developed a phenomenological model of coupled excitable cells with the aim to explore the prerequisites and mechanisms that lead to complex dynamical behavior observed in experiments. The minimalistic map-based description of excitable dynamics mimics the activity of beta cells. The stimulation was modeled as a heterogeneous and delayed increase in the excitability level, whereby a higher increase was used when supraphysiological stimulation was simulated (see section “Materials and Methods”). To further account for beta cell heterogeneity, we additionally included cell-to-cell variability in the absolute levels of excitability and in the intercellular coupling strength. In this case, we obtained a good qualitative agreement with experimental findings. The results are presented in Figures 6, 7 for the simulation of physiological and supraphysiological stimulations, respectively.

FIGURE 6

Figure 6. Simulated beta cell responses after switching from a substimulatory to stimulatory levels of stimulation, i.e., from 6 to 8 mM glucose. (A) Three characteristic traces of simulated cellular dynamics and the raster plot of binarized cellular activity. The orange dotted line indicates the fraction of active cells within the given time-window that was slid throughout the simulation. (B,D) 3D raster plots showing the excitation waves for selected intervals for the activation (B) and plateau (D) phase. Colors denote individual waves. (C,E) The distributions of excitation wave sizes for the activation (C) and plateau (E) phase. The gray dashed line indicates the power-law fit. The slope in the critical-like activation phase is –1.64.

FIGURE 7

Figure 7. Simulated beta cell responses after switching from a substimulatory to a suprastimulatory level of stimulation, i.e., from 6 to 12 mM glucose. (A) Three characteristic traces of simulated cellular dynamics and the raster plot of binarized cellular activity. The orange dotted line indicates the fraction of active cells within the given time-window that was slid throughout the simulation. (B,D) 3D raster plots showing the excitation waves for selected intervals for the activation (B) and plateau (D) phase. Colors denote individual waves. (C,E) The distributions of excitation wave sizes for the activation (C) and plateau (E) phase. The gray dashed line indicates the power-law fit.

Regardless of the stimulation level, we observed a biphasic response. Most importantly, the activation phase under lower stimulation levels was rather long and exhibited waves of various sizes and many of them were confined to small sub-regions of the islet. This is a result of very heterogeneous responses to stimulation, which in turn led to non-trivial and self-organized dynamical patterns. As in experimental measurement, the distribution of spatiotemporal cluster sizes was found to roughly follow a power-law (Figure 6C), which pinpoints toward a transient phase of critical dynamics. In contrast, the activation phase under supraphysiological conditions was shorter and the wave sizes were larger and more homogeneous, similarly as in the experiment. Consequently, the distribution deviates from the pure power-law behavior, mostly on account of an excess of larger excitation events (Figure 7C). The second plateau phase was qualitatively very similar in both scenarios. In both cases the spatiotemporal activity was dominated by global waves, thereby indicating supercritical behavior (Figures 6E, 7E). However, the waves were found to be more frequent and coherent under higher stimulation levels. This resulted to a large extent due to higher excitability levels, which made the cells operate in an even more ordered regime in which stochasticity is less pronounced.

A good agreement between experimental and computational results was obtained only if all three types of heterogeneities, i.e., in excitability level, in the delayed responses to stimulation, and in intercellular coupling strengths were implemented simultaneously. To test the necessity of such a multilayered heterogeneity, we systematically performed simulations with physiological and supraphysiological stimulations without considering one of the particular heterogeneities. Results are presented in Figure 8. It can be seen qualitatively that without any of the heterogeneities the simulations do not match well with experimental results. The most obvious difference occurs in the activation phase after the initial activation, where especially in the case of physiological stimulation diverse wave sizes were observed if all three types of heterogeneities were considered. Here, on the other hand, the dynamics after the initial activation and before the system shifts to the plateau phase, is very inactive and lacks on an emergent transitory phase with progressive recruitment of cells. On the contrary, the plateau phase seems to be weakly affected by the lack of any type of variability and even if one of the heterogeneities is missing, the system behavior very similar as in control simulations. This behavior is somehow expected, since after (probably unphysiological) prolonged stimulation all cells get very excitable and placed in the supercritical regime. Also, the heterogeneity imposed by variability in metabolism diminishes. For a more quantitative evaluation we present in Table 1 the relative activity time in both phases for the experimental data, for the control simulations with 30% variability in all three types cellular heterogeneity, and for simulations without one particular aspect of heterogeneity. The results indicate that indeed the interplay between all three types of heterogeneities is necessary to firmly reproduce the experimentally observed behavior, although it seems that variability in metabolism is the most important determinant, whereas the heterogeneity in the coupling appears to be the least important.

FIGURE 8

Figure 8. Simulated beta cell responses without particular types of cellular heterogeneities. Behavior without heterogeneity in the intrinsic excitability level (A,B), without heterogeneity in the delayed responses to stimulation (C,D), and without heterogeneity in intercellular coupling (E,F), for both physiological (A,C,E) and supraphysiological (B,D,F) stimulation levels.

TABLE 1

Table 1. Activity time during physiological and supraphysiological stimulations.

Finally, after determining that all three types of heterogeneities are required, we explored the impact of their level on the spatio-temporal activity. Figure 9 features the results. It can be observed that no or low degrees of heterogeneities fail to firmly reproduce experimental findings. In case of physiological stimulation, the initial activation of cells is missing and the cells respond much later without a progressive recruitment characterized by excitation waves of different sizes. Moreover, also the emulated supraphysiological stimulation differs if the level of heterogeneity is too low. Especially the dynamical phase after the initial activation is in this case very quiet, in contrast to the experiment and simulations with a higher degree of cellular variability, where a certain fraction of cells oscillates in this intermediate regime before the switch to the plateau phase. A quantitative assessment of this observation is presented in Table 1. It can be seen that with increasing levels of cell-to-cell variability the behavior in simulations becomes more similar to experimental results, although especially in the case with physiological stimulation 20% heterogeneity is not sufficient to achieve good consistency.

FIGURE 9

Figure 9. Simulated beta cell responses for different degrees of heterogeneities. Simulating the behavior without (A,B), with 10% (C,D), and 20% (E,F) of the used degree of heterogeneities for physiological (A,C,E) and supraphysiological (B,D,F) stimulation levels.

Discussion

Information processing in living organisms is orchestrated by large networks of interacting cells. In many cases, the dynamics of these networks is guided by the activation of one or a few elements, which in turn provokes the triggering of other elements, thereby leading to avalanches of activity that propagate through the system. Such emergent behavior is associated with self-organization and very often with critical dynamics resulting in a power-law distribution of the spatial and/or temporal extent of activity profiles. This scenario is particularly appealing for excitable systems, such as neuronal networks (Plenz and Thiagarajan, 2007; Hesse and Gross, 2014; Muñoz, 2018) or excitable cells and tissues (Lopez et al., 2012; Nivala et al., 2012; Gosak et al., 2017). Typically, critical dynamics emerges at the transition between randomness (subcritical dynamics) and order (supercritical dynamics). Variations of dynamical regimes can be induced by changes of global parameters, such as the excitability level, which reflects stimulus intensity. However, in this case a power-law behavior would only be expected in a narrow parameter space in the proximity of a phase transition point. Previous research has underlined the activity-dependent synaptic plasticity, heterogeneity, and hierarchical network organization as plausible mechanisms to overwhelm this drawback. Namely, these biological determinants were found to facilitate scale-free behavior and drive neuronal networks toward the critical state (Levina et al., 2007; Rubinov et al., 2011; Moretti and Muñoz, 2013). These self-organization mechanisms make the oscillators hover around the critical point, which are therefore able to generate effective scale invariance across quite a few scales. The phenomenon is often termed as self-organized quasi-criticality (Muñoz, 2018).

In the present study, we suggest a new mechanism that realistic excitable systems might exploit for expanding the operation in a critical-like regime. As the level of cellular excitability (parameter α) increases with time, the spatio-temporal activity switches from an inactive to an active phase (Osipov et al., 2007). If excitable oscillators are homogeneous, critical behavior is expected only at the phase transition point. In our setting, where the control parameter increases with time, criticality would therefore be observed for a very short transient period of time. However, the combination of cell-to-cell variability and heterogeneously delayed increases in excitability levels substantially broaden this regime of critical-like behavior. As a result, a rather long transient dynamical phase emerges with heterogeneous wave sizes, the distribution of which closely follows a power law. Even though this critical phase is only temporary and later followed by a supercritical phase dominated by global excitatory events, critical-like behavior persists for substantial periods of time. It should be emphasized that due to the transient nature of the scale-invariant activation phase, the observed power law is only an estimation, since much higher number of events would be required to confirm a pure power-law behavior. Moreover, the present study only provides an empirical observation of transient critical-like dynamical state, and further theoretical efforts are needed to uncover the exact mechanisms for the emergence of critical dynamics in such heterogeneous systems with delayed feedbacks.

In real-life settings, it is quite common to have transient or oscillatory stimulation patterns instead of a long-lasting permanent stimulation. In neurons (Buzsaki, 2004; Schroeder and Lakatos, 2009), in pancreatic islets (Pedersen et al., 2013; Satin et al., 2015), and in cardiac myocytes (O’Rourke et al., 1994), the excitation dynamics is governed by basal variations in intrinsic excitability, for instance due to oscillations in hormone or nutrient concentrations. Notably, in our previous study we have shown that such an oscillatory entraining might be a key toward persisting criticality in pancreatic beta cells (Gosak et al., 2017). Apparently, providing proper transitory conditions for heterogeneous excitable elements that exhibit delayed and variable responses to stimulation is a viable route to scale-free behavior. Because of these heterogeneities, changeable and confined regions with elevated excitability emerge from which the excitation waves are triggered. In the activation phase, the waves are typically triggered from cells whose excitability level increased faster, whereas in the plateau phase no specific patterns can be inferred. Moreover, the range of waves in the activation phase depends on the coupling and on the variable excitable state of surrounding elements. This in turn leads to emergent behavior with very heterogeneous spatiotemporal patterns. Most importantly, the critical-like activation phase is only possible to achieve if the stimulation level is not too high. Namely, in case of supraphysiological stimulation levels, the transition to the supercritical state is too abrupt and possibly accompanied by processes that do not occur under physiological conditions, at least not to a notable extent. Consequently, excitable cells are not able to self-organize into a scale-invariant dynamical state, as is the case in physiological stimulation conditions. Cell-to-cell differences are always present to some degree in any cell population, impacting the signaling processes in various tissues and settings (Muotri and Gage, 2006; Marhl et al., 2010; Paszek et al., 2010). Notably, cellular heterogeneity is more than a nuisance and often serves a biological function or contains meaningful information (Altschuler and Wu, 2010). In islet research, beta cell heterogeneity has been one of the key issues for decades and is becoming increasingly popular, particularly in the context of subpopulations (Bader et al., 2016; Avrahami et al., 2017). Thus, the central concept of our model, i.e., multiform beta cell heterogeneity, has a long tradition. In 1987, Pipeleers defined it on grounds of structural, functional, and replicative differences between beta cells. More specifically, he pointed out differences between beta cells in terms of contact with other types of endocrine cells, in gap junctional coupling, and cellular hormone content, in their ability and sensitivity to mount a response to glucose, and in their replicative potential. Ahead of time, he argued that altered beta cell heterogeneity may turn out to be important in development of diabetes and in islet transplantation (Pipeleers, 1987). In two updates shortly thereafter, he provided evidence for functional differences between beta cells in their rates of glucose-induced insulin synthesis and secretion that were attributed to differences in thresholds for glucose utilization and oxidation. Notably, already at that time the idea was put forward that cellular heterogeneity crucially determines dose-dependence of the beta cell response to glucose due to recruitment of beta cells into an active state by increasing glucose concentrations and that elevated levels of glucose can decrease the extent of heterogeneity. Additionally, he proposed that heterogeneity is not just an experimental artifact observed in dispersed beta cells, but at work also in intact tissue (Pipeleers, 1992; Pipeleers et al., 1994). Following the advent of new molecular markers and the omics approaches, in the last decades the concept of beta cell heterogeneity has been further supported at the transcriptomic and proteomic level, together with novel findings that heterogeneity affects beta cell proliferation and survival, as well as their stimulus-secretion coupling, from glucose metabolism and Ca²⁺ signaling to insulin secretion (Benninger and Piston, 2014; Roscioni et al., 2016; Avrahami et al., 2017; Gutierrez et al., 2017; Pipeleers et al., 2017; Benninger and Hodson, 2018). Finally, it has been proposed that the lack of beta cell heterogeneity may importantly contribute to islet failure in diabetes (Johnston et al., 2016; Pipeleers et al., 2017; Skelin Klemen et al., 2017; Benninger and Hodson, 2018; Nasteska and Hodson, 2018). One major drawback, common to most recent work, is the use of dispersed beta cells. It therefore remains to be investigated to what extent the heterogeneity described thus far is translationally relevant in the tissue context or even in vivo (Carrano et al., 2017; Benninger and Hodson, 2018; Gosak et al., 2018; Nasteska and Hodson, 2018).

By employing the tissue slice approach, we studied a large number of coupled beta cells in their normal tissue environment. Our experimental and modeling results intersect with both the original and more recent findings on beta cell heterogeneity at several points and provide some new ideas. First, during the activation phase, differences in glucose sensitivity were observed between different cells within the same islet and these differences were larger in lower glucose (8 mM). In other words, in lower glucose, gradual recruitment of beta cells into an active state, brought about by local [Ca²⁺]_IC waves displaying critical behavior seems to be a major feature of the islet response to constant stimulation. In contrast, in high glucose (12 mM), recruitment is less well pronounced due to early global [Ca²⁺]_IC waves showing supercritical behavior. We wish to speculate that recruitment (Stožer et al., 2013a; Gosak et al., 2017) and local [Ca²⁺]_IC waves (Benninger et al., 2014; Westacott et al., 2017) in islets, as opposed to dispersed cells or clusters of cells (Jonkers and Henquin, 2001), have received less attention in the scientific community due to use of high stimulatory glucose concentrations. Second, further aspects of adaptation in the response to higher glucose are the shorter average delay to activation, shorter activation phase, and higher activity during the plateau phase. Since it has been shown recently that in addition to different pools of granules, the triggering Ca²⁺ signal importantly shapes the biphasic insulin secretion in response to constant stimulation by glucose (Pedersen et al., 2019), our findings shall importantly inform future models of beta cell insulin secretion. Additionally, the behavior of beta cells during the activation phase, when they are not functioning in synchrony with other cells, could be compared with their properties during the plateau phase to more precisely establish the relationships between their different roles. More specifically, such comparison could help answer the question whether the cells that activate first are also the ones that initiate global [Ca²⁺]_IC waves and possess the most functional connections, i.e., function as hubs, during the plateau phase (Stožer et al., 2013b; Johnston et al., 2016; Westacott et al., 2017). Third, using our phenomenological model, we found that all three types of heterogeneities, the choice of which is further substantiated in the following section, are necessary and sufficient to reproduce the experimentally observed behavior. This of course does not exclude the possibility that additional aspects of heterogeneity exist in reality and further modulate beta cell responses.

In comparison with our previous study, we used here a very simple phenomenological model to reproduce the experimentally observed non-trivial activity patterns in islets (Gosak et al., 2017). Such minimalistic modeling approaches have of course limitations, since they do not allow for any mechanistic insights into physiological processes and signaling pathways. At the same time, they offer several advantages. They are numerically very efficient, and most importantly, they contain a small number of parameters whose roles are rather clear, which makes it easier to explicitly study particular aspects of cellular heterogeneity. In contrast, realistic and multi-component cellular models exhibit many parameters that in general affect several aspects of signalization, which hinders a systematic and definitive examination of their particular influences on cellular behavior. Finally, it should be noted that the majority of existing comprehensive beta cell models were mainly focused on the activity on the plateau phase, whereas modeling of collective cellular activations after switching from substimulatory to (supra)stimulatory glucose received very little attention. The mechanisms that govern such stimulus-dependent activation are also understood incompletely. Phenomenological modeling is therefore beneficial in this respect, as long as the empirical description of the processes ensures good agreement between modeling and experimental results. In particular, the time lags and temporal evolution of the excitability level [see eq. (3)] are plausible processes that can be qualitatively linked with previous experimental observations, such as differences in metabolic sensitivity to glucose and the following electrical and [Ca²⁺]_IC responses (Stožer et al., 2013a,b; Benninger et al., 2014; Farnsworth and Benninger, 2014; Johnston et al., 2016). Moreover, we decided to include heterogeneity in intercellular coupling due to the extensive experimental support demonstrating its importance in both normal and pathological islet functioning (Hodson et al., 2013; Farnsworth et al., 2014, 2016; Johnston et al., 2016; Skelin Klemen et al., 2017). However, here we focused only on the fast [Ca²⁺]_IC oscillations and future studies will conceivably need to include additional aspects of heterogeneity to provide a comprehensive and realistic beta cell model capable of describing other components of the [Ca²⁺]_IC pattern, other parameters in the stimulus-secretion cascade, as well as responses to different levels of stimulation and different secretagogues.

Multiscale and multidimensional heterogeneity represent a viable route to critical-like behavior for a substantial period of time. Specifically, the dynamical transition between the inactive and active state occurs in a rather broad temporal interval, especially in the case of physiological levels of stimulation. In other words, as the glucose level increases the activation of cells is not abrupt. Before switching to the dynamical state with global Ca²⁺ events, a transient period of very heterogeneous wave sizes is observed, which implies a critical-like behavior, since the system as a whole bypasses the critical point rather slowly. On the other hand, supraphysiological high stimulation levels lack on such progressive recruitment of cells and lead to a rapid transition to a fully active state. This might be crucial for healthy physiological functioning of pancreatic islets and potentially of other biological tissues as well. Namely, for other biological tissues, there is a large body of evidence indicating severe pathophysiological consequences of an abrupt collective transitions to hyper-regulated synchronous tissue responses. An overview was given by Trefois et al. (2015), showing that critical transitions are identified as early warning signals for the onset of different pathologies ranging from microbiome dysregulations to irritable bowel syndrome, asthma, pulmonary disease, depression, type 1 and type 2 diabetes, inflammation, start and termination of epileptic seizures, cancer, and cardiovascular events.

Further investigations are needed to understand the onset of pathological supercritical behavior in more detail. The molecular and cellular mechanisms are still obscure; however, the results of our study, although only empirical, give at least a hint to an improved methodology, opening a new dimension in studying the (premature) onset of supercriticality by looking at the extent of cell heterogeneity. Some preliminary studies in our lab show that beta cell responses in terms of [Ca²⁺]_IC signals shall also be correlated with other aspects of heterogeneity to get a more complete picture about the mechanisms that make some cells more responsive to glucose and to find out whether this is a stable property or something that changes with time and on which temporal scale. In addition, different concentrations of glucose and additional stimulation protocols, as well other secretagogues shall be used in future studies. Moreover, the general extent of heterogeneity in the islets, as well as the properties of individual cells shall be investigated in mouse models of diabetes and in human islets from normal and diabetic donors to more clearly define the changes under pathological conditions and suggest targets for treatment (Benninger and Hodson, 2018; Stožer et al., 2019). In general, from the viewpoint of clinical approaches, the understanding of critical transitions might help us develop therapies that are more effective. From the viewpoint of preventive health care, an improved understanding of the pathological premature transitions to supercriticality could help us identify and characterize some early warning signals predicting the upcoming pathological transitions. Finally, beyond the preventive and therapeutic role Bargaje et al. (2017) showed that an increase in cell heterogeneity of stem cells just before the critical transition correlates with a branching point on the trajectory of cell fate. This represents a useful tool for forecasting the cell fate outcomes and can be used for optimizing the differentiation protocols in order to obtain desired cell populations, which opens a completely new dimension of bioengineering in the future.

Data Availability

The datasets generated for this study are available on request to the corresponding author.

Ethics Statement

The study was conducted in strict accordance with all national and the European recommendations pertaining to care and work with experimental animals, and all efforts were made to minimize suffering of animals. The protocol was approved by the Veterinary Administration of the Republic of Slovenia (Permit Number: U34401-12/2015/3).

Author Contributions

AS, RM, JD, MP, MSR, MM, and MG designed the study, carried out the research, and wrote the manuscript.

Funding

The authors acknowledge the financial support from the Slovenian Research Agency (Research Core Funding Nos. P1-0055, P1-0403, and P3-0396 and Research Projects Nos. J3-9289, N3-0048, I0-0029, J1-7009, J7-7226, J4-9302, J1-9112, and P1-0403) and the Austrian Science Fund/Fonds zur Förderung der Wissenschaftlichen Forschung (Bilateral Grant I3562-B27 to MSR).

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.

References

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

Altschuler, S. J., and Wu, L. F. (2010). Cellular heterogeneity: do differences make a difference? Cell 141, 559–563. doi: 10.1016/j.cell.2010.04.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Aon, M. A., Cortassa, S., Akar, F. G., and O’Rourke, B. (2006). Mitochondrial criticality: a new concept at the turning point of life or death. Biochim. Biophys. Acta 1762, 232–240. doi: 10.1016/j.bbadis.2005.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Aon, M. A., Cortassa, S., and O’Rourke, B. (2004). Percolation and criticality in a mitochondrial network. Proc. Natl. Acad. Sci. U.S.A. 101, 4447–4452. doi: 10.1073/pnas.0307156101

PubMed Abstract | CrossRef Full Text | Google Scholar

Avrahami, D., Klochendler, A., Dor, Y., and Glaser, B. (2017). Beta cell heterogeneity: an evolving concept. Diabetologia 60, 1363–1369. doi: 10.1007/s00125-017-4326-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Bader, E., Migliorini, A., Gegg, M., Moruzzi, N., Gerdes, J., Roscioni, S. S., et al. (2016). Identification of proliferative and mature β-cells in the islets of Langerhans. Nature 535, 430–434. doi: 10.1038/nature18624

PubMed Abstract | CrossRef Full Text | Google Scholar

Bak, P. (1996). How Nature Works. New York, NY: Springer New York, doi: 10.1007/978-1-4757-5426-5421

CrossRef Full Text | Google Scholar

Bargaje, R., Trachana, K., Shelton, M. N., McGinnis, C. S., Zhou, J. X., Chadick, C., et al. (2017). Cell population structure prior to bifurcation predicts efficiency of directed differentiation in human induced pluripotent cells. Proc. Natl. Acad. Sci. U.S.A. 114, 2271–2276. doi: 10.1073/pnas.1621412114

PubMed Abstract | CrossRef Full Text | Google Scholar

Barral, J., and D Reyes, A. (2016). Synaptic scaling rule preserves excitatory–inhibitory balance and salient neuronal network dynamics. Nat. Neurosci 19, 1690–1696. doi: 10.1038/nn.4415

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. doi: 10.1523/JNEUROSCI.23-35-11167.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Benninger, R. K. P., Dorrell, C., Hodson, D. J., and Rutter, G. A. (2018). The Impact of Pancreatic Beta Cell Heterogeneity on Type 1 Diabetes Pathogenesis. Curr. Diab. Rep. 18:112. doi: 10.1007/s11892-018-1085-1082

PubMed Abstract | CrossRef Full Text | 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

PubMed Abstract | CrossRef Full Text | Google Scholar

Benninger, R. K. P., and Hodson, D. J. (2018). New understanding of β-Cell heterogeneity and in situ islet function. Diabetes 67, 537–547. doi: 10.2337/DBI17-0040

PubMed Abstract | 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

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 Ca2+ wave dynamics. Biophys. J. 107, 2723–2733. doi: 10.1016/j.bpj.2014.10.048

PubMed Abstract | CrossRef Full Text | Google Scholar

Bergsten, P., Grapengiesser, E., Gylfe, E., Tengholm, A., and Hellman, B. (1994). Synchronous oscillations of cytoplasmic Ca2+ and insulin release in glucose-stimulated pancreatic islets. J. Biol. Chem. 269, 8749–8753.

PubMed Abstract | Google Scholar

Berridge, M. J., Lipp, P., and Bootman, M. D. (2000). The versatility and universality of calcium signalling. Nat. Rev. Mol. Cell Biol. 1, 11–21. doi: 10.1038/35036035

PubMed Abstract | CrossRef Full Text | Google Scholar

Bertram, R., Sherman, A., and Satin, L. S. (2007). Metabolic and electrical oscillations: partners in controlling pulsatile insulin secretion. Am. J. Physiol. Metab. 293, E890–E900. doi: 10.1152/ajpendo.00359.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Bialek, W., Cavagna, A., Giardina, I., Mora, T., Silvestri, E., Viale, M., et al. (2012). Statistical mechanics for natural flocks of birds. Proc. Natl. Acad. Sci. U.S.A. 109, 4786–4791. doi: 10.1073/pnas.1118633109

PubMed Abstract | CrossRef Full Text | Google Scholar

Bosco, D., Haefliger, J.-A., and Meda, P. (2011). Connexins: key mediators of endocrine function. Physiol. Rev. 91, 1393–1445. doi: 10.1152/physrev.00027.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Brochini, L., de Andrade Costa, A., Abadi, M., Roque, A. C., Stolfi, J., and Kinouchi, O. (2016). Phase transitions and self-organized criticality in networks of stochastic spiking neurons. Sci. Rep. 6:35831. doi: 10.1038/srep35831

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsaki, G. (2004). Neuronal oscillations in cortical networks. Science 304, 1926–1929. doi: 10.1126/science.1099745

PubMed Abstract | CrossRef Full Text | Google Scholar

Calabrese, A., Zhang, M., Serre-Beinier, V., Caton, D., Mas, C., Satin, L. S., et al. (2003). Connexin 36 controls synchronization of Ca2+ oscillations and insulin secretion in MIN6 cells. Diabetes Metab. Res. Rev. 52, 417–424. doi: 10.2337/diabetes.52.2.417

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:053103. doi: 10.1063/1.4949020

PubMed Abstract | CrossRef Full Text | Google Scholar

Carrano, A. C., Mulas, F., Zeng, C., and Sander, M. (2017). Interrogating islets in health and disease with single-cell technologies. Mol. Metab. 6, 991–1001. doi: 10.1016/j.molmet.2017.04.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Carvalho, C. P. F., Oliveira, R. B., Britan, A., Santos-Silva, J. C., Boschero, A. C., Meda, P., et al. (2012). Impaired β-cell-β-cell coupling mediated by Cx36 gap junctions in prediabetic mice. Am. J. Physiol. Metab. 303, E144–E151. doi: 10.1152/ajpendo.00489.2011

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

di Santo, S., Villegas, P., Burioni, R., and Muñoz, M. A. (2018). Landau–Ginzburg theory of cortex dynamics: scale-free avalanches emerge at the edge of synchronization. Proc. Natl. Acad. Sci. U.S.A. 115, E1356–E1365. doi: 10.1073/pnas.1712989115

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–16. doi: 10.1371/journal.pone.0082374

PubMed Abstract | CrossRef Full Text | Google Scholar

Ellis, G. F. R., and Kopel, J. (2019). The dynamical emergence of biology from physics: branching causation via biomolecules. Front. Physiol. 9:1966. doi: 10.3389/fphys.2018.01966

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., Hemmati, A., Pozzoli, M., and Benninger, R. K. P. (2014). Fluorescence recovery after photobleaching reveals regulation and distribution of connexin36 gap junction coupling within mouse islets of Langerhans. J. Physiol. 592, 4431–4446. doi: 10.1113/jphysiol.2014.276733

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

Friedman, N., Ito, S., Brinkman, B. A. W., Shimono, M., DeVille, R. E. L., Dahmen, K. A., et al. (2012). Universal critical dynamics in high resolution neuronal avalanche data. Phys. Rev. Lett. 108:208102. doi: 10.1103/PhysRevLett.108.208102

PubMed Abstract | CrossRef Full Text | Google Scholar

Furusawa, C., and Kaneko, K. (2012). Adaptation to optimal cell growth through self-organized criticality. Phys. Rev. Lett. 108:208103. doi: 10.1103/PhysRevLett.108.208103

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerdes, J. M., Christou-Savina, S., Xiong, Y., Moede, T., Moruzzi, N., Karlsson-Edlund, P., et al. (2014). Ciliary dysfunction impairs beta-cell insulin secretion and promotes development of type 2 diabetes in rodents. Nat. Commun. 5:5308. doi: 10.1038/ncomms6308

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilon, P., Chae, H.-Y., Rutter, G. A., and Ravier, M. A. (2014). Calcium signaling in pancreatic β-cells in health and in Type 2 diabetes. Cell Calcium 56, 340–361. doi: 10.1016/j.ceca.2014.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | Google Scholar

Gilon, P., Shepherd, R. M., and Henquin, J. C. (1993). Oscillations of secretion driven by oscillations of cytoplasmic Ca2+ as evidences in single pancreatic islets. J. Biol. Chem 268, 22265–22268.

PubMed Abstract | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Gosak, M., Stožer, A., Markovič, R., Dolenšek, J., Perc, M., Rupnik, M. S., et al. (2017). Critical and supercritical spatiotemporal calcium dynamics in beta cells. Front. Physiol. 8:1106. doi: 10.3389/fphys.2017.01106

PubMed Abstract | CrossRef Full Text | Google Scholar

Gutierrez, G. D., Gromada, J., and Sussel, L. (2017). Heterogeneity of the pancreatic beta cell. Front. Genet. 8:22. doi: 10.3389/fgene.2017.00022

PubMed Abstract | CrossRef Full Text | Google Scholar

Hahn, G., Ponce-Alvarez, A., Monier, C., Benvenuti, G., Kumar, A., Chavane, F., et al. (2017). Spontaneous cortical activity is transiently poised close to criticality. PLoS Comput. Biol. 13:e1005543. doi: 10.1371/journal.pcbi.1005543

PubMed Abstract | CrossRef Full Text | Google Scholar

Haimovici, A., Tagliazucchi, E., Balenzuela, P., and Chialvo, D. R. (2013). Brain organization into resting state networks emerges at criticality on a model of the human connectome. Phys. Rev. Lett. 110:178101. doi: 10.1103/PhysRevLett.110.178101

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamelin, R., Allagnat, F., Haefliger, J.-A., and Meda, P. (2009). Connexins, diabetes and the metabolic syndrome. Curr. Protein Pept. Sci. 10, 18–29. doi: 10.2174/138920309787315167

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., and Meissner, H. P. (1984). Significance of ionic fluxes and changes in membrane potential for stimulus-secretion coupling in pancreatic B-cells. Experientia 40, 1043–1052. doi: 10.1007/BF01971450

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

Hobbs, J. P., Smith, J. L., and Beggs, J. M. (2010). Aberrant neuronal avalanches in cortical tissue removed from juvenile epilepsy patients. J. Clin. Neurophysiol. 27, 380–386. doi: 10.1097/WNP.0b013e3181fdf8d3

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

Hraha, T. H., Westacott, M. J., Pozzoli, M., Notary, A. M., McClatchey, P. M., and Benninger, R. K. P. (2014). 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 | Google Scholar

Ibarz, B., Casado, J. M., and Sanjuán, M. A. F. (2011). Map-based models in neuronal dynamics. Phys. Rep. 501, 1–74. doi: 10.1016/j.physrep.2010.12.003

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 Ca2+ in Islet cell clusters show that glucose rapidly recruits β-Cells and gradually increases the individual cell response. Diabetes Metab. Res. Rev. 50, 540–550. doi: 10.2337/diabetes.50.3.540

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. doi: 10.1152/jn.1998.79.2.1098

PubMed Abstract | CrossRef Full Text | 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-62156

PubMed Abstract | CrossRef Full Text | Google Scholar

Karin, O., and Alon, U. (2017). Biphasic response as a mechanism against mutant takeover in tissue homeostasis circuits. Mol. Syst. Biol. 13:933. doi: 10.15252/msb.20177599

PubMed Abstract | CrossRef Full Text | Google Scholar

Karin, O., Swisa, A., Glaser, B., Dor, Y., and Alon, U. (2016). Dynamical compensation in physiological circuits. Mol. Syst. Biol. 12:886. doi: 10.15252/msb.20167216

PubMed Abstract | CrossRef Full Text | Google Scholar

Khaluf, Y., Ferrante, E., Simoens, P., and Huepe, C. (2017). Scale invariance in natural and artificial collective systems: a review. J. R. Soc. Interface 14:20170662. doi: 10.1098/rsif.2017.0662

PubMed Abstract | CrossRef Full Text | 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

Konstantinova, I., Nikolova, G., Ohara-Imaizumi, M., Meda, P., Kuc̆era, 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

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Shuai, H. Y., Gylfe, E., and Tengholm, A. (2013). Oscillations of sub-membrane ATP in glucose-stimulated beta cells depend on negative feedback from Ca2+. Diabetologia 56, 1577–1586. doi: 10.1007/s00125-013-2894-2890

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L., Trifunovic, A., Kohler, M., Wang, Y., Petrovic Berglund, J., Illies, C., et al. (2014). Defects in β-Cell Ca2+ dynamics in age-induced diabetes. Diabetes Metab. Res. Rev 63, 4100–4114. doi: 10.2337/db13-1855

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

Loppini, A., Capolupo, A., Cherubini, C., Gizzi, A., Bertolaso, M., Filippi, S., et al. (2014). On the coherent behavior of pancreatic beta cell clusters. Phys. Lett. A 378, 3210–3217. doi: 10.1016/J.PHYSLETA.2014.09.041

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

MacDonald, P. E., and Rorsman, P. (2006). Oscillations, intercellular coupling, and insulin secretion in pancreatic β cells. PLoS Biol. 4:e49. doi: 10.1371/journal.pbio.0040049

PubMed Abstract | CrossRef Full Text | Google Scholar

Marhl, M., Gosak, M., Perc, M., and Roux, E. (2010). Importance of cell variability for calcium signaling in rat airway myocytes. Biophys. Chem. 148, 42–50. doi: 10.1016/j.bpc.2010.02.006

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. doi: 10.1038/srep07845

PubMed Abstract | CrossRef Full Text | Google Scholar

Markovič, R., Gosak, M., and Marhl, M. (2012). “The role of topological features of intercellular communication networks by the synchronization of cellular oscillators,” in Let’s Face Chaos through Nonlinear Dynamics, eds M. Robnik and V. G. Romanovski (Maribor: AIP Conference Proceedings), 256–267. doi: 10.1063/1.4745588

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

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. doi: 10.1038/srep10578

PubMed Abstract | CrossRef Full Text | 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

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Muñoz, M. A. (2018). Colloquium: criticality and dynamical scaling in living systems. Rev. Mod. Phys. 90:031001. doi: 10.1103/RevModPhys.90.031001

CrossRef Full Text | Google Scholar

Muotri, A. R., and Gage, F. H. (2006). Generation of neuronal variability and complexity. Nature 441, 1087–1093. doi: 10.1038/nature04959

PubMed Abstract | CrossRef Full Text | Google Scholar

Nasteska, D., and Hodson, D. J. (2018). The role of beta cell heterogeneity in islet function and insulin release. J. Mol. Endocrinol. 61, R43–R60. doi: 10.1530/JME-18-11

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

Nonnenmacher, M., Behrens, C., Berens, P., Bethge, M., and Macke, J. H. (2017). Signatures of criticality arise from random subsampling in simple population models. PLoS Comput. Biol. 13:e1005718. doi: 10.1371/journal.pcbi.1005718

PubMed Abstract | CrossRef Full Text | Google Scholar

Nykter, M., Price, N. D., Aldana, M., Ramsey, S. A., Kauffman, S. A., Hood, L. E., et al. (2008). Gene expression dynamics in the macrophage exhibit criticality. Proc. Natl. Acad. Sci. U.S.A. 105, 1897–1900. doi: 10.1073/pnas.0711525105

PubMed Abstract | CrossRef Full Text | Google Scholar

O’Rourke, B., Ramza, B., and Marban, E. (1994). Oscillations of membrane current and excitability driven by metabolic oscillations in heart cells. Science 265, 962–966. doi: 10.1126/science.8052856

PubMed Abstract | CrossRef Full Text | Google Scholar

Osipov, G. V., Kurths, J., and Zhou, C. (2007). Synchronization in Oscillatory Networks. Berlin: Springer Berlin Heidelberg, doi: 10.1007/978-3-540-71269-71265

CrossRef Full Text | Google Scholar

Pasquale, V., Massobrio, P., Bologna, L. L., Chiappalone, M., and Martinoia, S. (2008). Self-organization and neuronal avalanches in networks of dissociated cortical neurons. Neuroscience 153, 1354–1369. doi: 10.1016/j.neuroscience.2008.03.050

PubMed Abstract | CrossRef Full Text | Google Scholar

Paszek, P., Ryan, S., Ashall, L., Sillitoe, K., Harper, C. V., Spiller, D. G., et al. (2010). Population robustness arising from cellular heterogeneity. Proc. Natl. Acad. Sci. U.S.A. 107, 11644–11649. doi: 10.1073/pnas.0913798107

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 Ca2+ Entrainment in Pancreatic Islets by Oscillatory Glucose. Biophys. J. 105, 29–39. doi: 10.1016/j.bpj.2013.05.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Pedersen, M. G., Tagliavini, A., and Henquin, J.-C. (2019). Calcium signaling and secretory granule pool dynamics underlie biphasic insulin secretion and its amplification by glucose: experiments and modeling. Am. J. Physiol. Metab. 316, E475–E486. doi: 10.1152/ajpendo.00380.2018

PubMed Abstract | CrossRef Full Text | Google Scholar

Penrose, M. (2007). Random Geometric Graphs, Vol. 5. Oxford: Oxford University Press, doi: 10.1093/acprof:oso/9780198506263.001.0001

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-82207

PubMed Abstract | CrossRef Full Text | Google Scholar

Pipeleers, D. (1987). The biosociology of pancreatic B cells. Diabetologia 30, 277–291. doi: 10.1007/bf00299019

CrossRef Full Text | Google Scholar

Pipeleers, D., De Mesmaeker, I., Robert, T., and Van Hulle, F. (2017). Heterogeneity in the beta-cell population: a guided search into its significance in pancreas and in implants. Curr. Diab. Rep. 17:86. doi: 10.1007/s11892-017-0925-929

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.

Google Scholar

Pipeleers, D. G. (1992). Heterogeneity in pancreatic β-cell population. Diabetes Metab. Res. Rev 41, 777–781. doi: 10.2337/diab.41.7.777

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

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

Ravier, M. A., Guldenagel, M., Charollais, A., Gjinovci, A., Caille, D., Sohl, G., et al. (2005). Loss of Connexin36 channels alters β-Cell coupling, islet synchronization of glucose-induced Ca2+ and insulin oscillations, and basal insulin release. Diabetes Metab. Res. Rev. 54, 1798–1807. doi: 10.2337/diabetes.54.6.1798

PubMed Abstract | CrossRef Full Text | Google Scholar

Rorsman, P., and Ashcroft, F. M. (2018). Pancreatic β-Cell electrical activity and insulin secretion: of mice and men. Physiol. Rev. 98, 117–214. doi: 10.1152/physrev.00008.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Rorsman, P., and Braun, M. (2013). Regulation of insulin secretion in human pancreatic islets. Annu. Rev. Physiol. 75, 155–179. doi: 10.1146/annurev-physiol-030212-183754

PubMed Abstract | CrossRef Full Text | Google Scholar

Roscioni, S. S., Migliorini, A., Gegg, M., and Lickert, H. (2016). Impact of islet architecture on β-cell heterogeneity, plasticity and function. Nat. Rev. Endocrinol. 12, 695–709. doi: 10.1038/nrendo.2016.147

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

Rulkov, N. F. (2002). Modeling of spiking-bursting neural behavior using two-dimensional map. Phys. Rev. E 65:041922. doi: 10.1103/PhysRevE.65.041922

PubMed Abstract | CrossRef Full Text | Google Scholar

Sasai, Y. (2013). Cytosystems dynamics in self-organization of tissue architecture. Nature 493, 318–326. doi: 10.1038/nature11859

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

Schroeder, C. E., and Lakatos, P. (2009). Low-frequency neuronal oscillations as instruments of sensory selection. Trends Neurosci. 32, 9–18. doi: 10.1016/j.tins.2008.09.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Shew, W. L., Clawson, W. P., Pobst, J., Karimipanah, Y., Wright, N. C., and Wessel, R. (2015). Adaptation to sensory input tunes visual cortex to criticality. Nat. Phys. 11, 659–663. doi: 10.1038/nphys3370

PubMed Abstract | CrossRef Full Text | Google Scholar

Shew, W. L., and Plenz, D. (2013). The functional benefits of criticality in the cortex. Neurosci. 19, 88–100. doi: 10.1177/1073858412445487

PubMed Abstract | CrossRef Full Text | 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 0:e1342022. doi: 10.1080/19382014.2017.1342022

PubMed Abstract | CrossRef Full Text | Google Scholar

Skupin, A., Kettenmann, H., Winkler, U., Wartenberg, M., Sauer, H., Tovey, S. C., et al. (2008). How does intracellular Ca2+ oscillate: by chance or by the clock? Biophys. J. 94, 2404–2411. doi: 10.1529/biophysj.107.119495

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-1099

PubMed Abstract | CrossRef Full Text | Google Scholar

Squires, P. E., Persaud, S. J., Hauge-Evans, A. C., Gray, E., Ratcliff, H., and Jones, P. M. (2002). Co-ordinated Ca2+-signalling within pancreatic islets: does β-cell entrainment require a secreted messenger. Cell Calcium 31, 209–219. doi: 10.1016/S0143-4160(02)00034-39

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:012719. doi: 10.1103/PhysRevE.89.012719

PubMed Abstract | CrossRef Full Text | Google Scholar

Stamper, I. J., and Wang, X. (2019). Integrated multiscale mathematical modeling of insulin secretion reveals the role of islet network integrity for proper oscillatory glucose-dose response. J. Theor. Biol. 475, 1–24. doi: 10.1016/j.jtbi.2019.05.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Stoop, R., and Gomez, F. (2016). Auditory power-law activation avalanches exhibit a fundamental computational ground state. Phys. Rev. Lett. 117:038102. doi: 10.1103/PhysRevLett.117.038102

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

Stožer, A., Hojs, R., and Dolenšek, J. (2019). Beta cell functional adaptation and dysfunction in insulin resistance and the role of chronic kidney disease. Nephron doi: 10.1159/000495665 [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

Tagliazucchi, E., Chialvo, D. R., Siniatchkin, M., Amico, E., Brichant, J.-F., Bonhomme, V., et al. (2016). Large-scale signatures of unconsciousness are consistent with a departure from critical dynamics. J. R. Soc. Interface 13:20151027. doi: 10.1098/rsif.2015.1027

PubMed Abstract | CrossRef Full Text | Google Scholar

Tamayo, P., Slonim, D., Mesirov, J., Zhu, Q., Kitareewan, S., Dmitrovsky, E., et al. (1999). Interpreting patterns of gene expression with self-organizing maps: methods and application to hematopoietic differentiation. Proc. Natl. Acad. Sci. U.S.A. 96, 2907–2912. doi: 10.1073/pnas.96.6.2907

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

Timme, N. M., Marshall, N. J., Bennett, N., Ripp, M., Lautzenhiser, E., and Beggs, J. M. (2016). Criticality maximizes complexity in neural tissue. Front. Physiol. 7:425. doi: 10.3389/fphys.2016.00425

PubMed Abstract | CrossRef Full Text | Google Scholar

Tkačik, G., Mora, T., Marre, O., Amodei, D., Palmer, S. E., Berry, M. J., et al. (2015). Thermodynamics and signatures of criticality in a network of neurons. Proc. Natl. Acad. Sci. U.S.A. 112, 11508–11513. doi: 10.1073/pnas.1514188112

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

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

Yang, B. T., Dayeh, T. A., Kirkpatrick, C. L., Taneera, J., Kumar, R., Groop, L., et al. (2011). Insulin promoter DNA methylation correlates negatively with insulin gene expression and positively with HbA1c levels in human pancreatic islets. Diabetologia 54, 360–367. doi: 10.1007/s00125-010-1967-1966

PubMed Abstract | 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), 155–169. doi: 10.1117/12.424970

CrossRef Full Text | Google Scholar

Zare, M., and Grigolini, P. (2013). Criticality and avalanches in neural networks. Chaos Soliton. Fract. 55, 80–94. doi: 10.1016/j.chaos.2013.05.009

CrossRef Full Text | Google Scholar