Collective Sensing of β-Cells Generates the Metabolic Code

Major part of a pancreatic islet is composed of β-cells that secrete insulin, a key hormone regulating influx of nutrients into all cells in a vertebrate organism to support nutrition, housekeeping or energy storage. β-cells constantly communicate with each other using both direct, short-range interactions through gap junctions, and paracrine long-range signaling. However, how these cell interactions shape collective sensing and cell behavior in islets that leads to insulin release is unknown. When stimulated by specific ligands, primarily glucose, β-cells collectively respond with expression of a series of transient Ca²⁺ changes on several temporal scales. Here we reanalyze a set of Ca²⁺ spike trains recorded in acute rodent pancreatic tissue slice under physiological conditions. We found strongly correlated states of co-spiking cells coexisting with mostly weak pairwise correlations widespread across the islet. Furthermore, the collective Ca²⁺ spiking activity in islet shows on-off intermittency with scaling of spiking amplitudes, and stimulus dependent autoassociative memory features. We use a simple spin glass-like model for the functional network of a β-cell collective to describe these findings and argue that Ca²⁺ spike trains produced by collective sensing of β-cells constitute part of the islet metabolic code that regulates insulin release and limits the islet size.

1. Introduction

Endocrine cells in vertebrates act both as coders and decoders of metabolic code (Tomkins, 1975) that carries information from primary endocrine sensors to target tissues. In endocrine pancreas, energy-rich ligands provide a continuous input to a variety of specific receptor proteins on and in individual β-cells and initiate signaling events in and between these cells (Henquin, 2009). In an oversimplified medical physiology textbook interpretation, glucose is transported into a β-cell through facilitated diffusion, is phosphorylated and converted within a metabolic black box to ATP, leading to closure of K_ATP channels, cell membrane depolarization and activation of voltage-activated calcium channels (VACCs), followed by a rise in cytosolic Ca²⁺ to a micromolar range and triggering of SNARE-dependent insulin release (Ashcroft and Rorsman, 1989). However, glucose may influence β-cells signaling through several additional routes. There may be alternative glucose entry routes, like for example active Na-glucose cotransport (Tomita, 1976; Trautmann and Wollheim, 1987), alternative calcium release sites, like ryanodine (Islam, 2002) and IP₃ receptors (Lang, 1999) or glucose may directly activate the sweet taste receptor and initiate signaling (Henquin, 2012), to name a few. Activation of a β-cell on a single cell level therefore likely involves triggering of a variety of elementary Ca²⁺ events (Berridge et al., 2000), which interfere in space and time into a unitary β-cell Ca²⁺ response to support Ca²⁺-dependent insulin release. This Ca²⁺-dependent insulin release can be further modulated by activation of different protein phosphorylation/dephosphorylation patterns (PKA, PKC, Cdk5, etc.) (Mandic et al., 2011; Skelin and Rupnik, 2011) or other protein modifications (Paulmann et al., 2009) to either reduce or increase the insulin output.

One of the important features of the sensory collectives is the optimization of the spatial relations between its elements to maximize the precision of sensing (Fancher and Mugler, 2017; Saakian, 2017). In islets of Langerhans, β-cells dwell as morphologically well defined cellulo-social collectives. These ovoid microorgans are typically not longer than 150 μm. The relatively small and constant pancreatic islet size is an intriguing feature in vertebrate biology. The size distribution of islets is comparable in humans, rodents and wider within different vertebrate species, irrespective of evident differences in overall body and pancreas size as well as total β-cell mass (Kim et al., 2009; Dolenšek et al., 2015). In mice, islet sizes range between 50 and 600 μm, with a median values below 150 μm (Lamprianou et al., 2011). To accommodate differences in the body size, there is nearly a linear relationship between the total number of similarly sized islets and body mass across different vertebrate species (Montanya et al., 2000; Bouwens and Rooman, 2005). However, why are islets so conserved in size is unknown.

All β-cells within an islet collective represent a single functional unit, electrically and chemically coupled network, with gap junction proteins, Connexins 36 (Cx36) (Bavamian et al., 2007), for short-range interactions and with paracrine signaling (Caicedo, 2013) for long-range interactions between cells. The unitary cell response in one β-cell influences the formation of similar responses in neighboring β-cells and contributes to coordination of a large number of β-cells (Cigliola et al., 2013; Stožer et al., 2013a). Explorations of these functional β-cell networks, constructed from thresholded pairwise correlations of Ca²⁺ imaging signals (Stožer et al., 2013b; Markovič et al., 2015; Johnston et al., 2016; Gosak et al., 2017a), showed that strongly correlated subsets of β-cell collective organize into modular, broad-scale networks with preferentially local correlations reaching up to 40 μm (Markovič et al., 2015), but understanding of mechanisms that lead to these strongly correlated networks states in β-cell populations is still lacking. We argue that β-cells sense, compute and respond to information as a collective, organized in a network similar to sensory neuron populations (Schneidman et al., 2006; Tkačik and Bialek, 2016), and not as a set of independent cells strongly coupled only when stimulation is high enough.

Here we reanalyze pairwise correlations of Ca²⁺ spike trains (unitary β-cell responses on the shortest temporal scale) in β-cell collective recorded in fresh pancreatic tissues slice under changing glucose stimulation conditions (6 mM subthreshold–8 mM stimulatory) using methodological approaches previously described (Stožer et al., 2013b; Markovič et al., 2015; Gosak et al., 2017a,b). We specifically look at weak correlations between β-cells which we found to be widely spread across the islet (Azhar and Bialek, 2010). Guided by the use of statistical physics models in describing populations of neurons (Schneidman et al., 2006; Tkacik et al., 2009), we use a simple spin glass model for Ca²⁺ β-cells activity and show that it well captures the features observed in the measured data. In a way, we recognize this efficiency of simple models in both neuronal and endocrine cell collectives as one manifestation of the “beauty in function” (Rasmussen, 1970).

2. Spin Model of a β-Cell Collective

Spin models have been borrowed from statistical physics to describe the functional behavior of large, highly interconnected systems like sensory neurons (Schneidman et al., 2006; Tkacik et al., 2009; Tkačik et al., 2014), immune system (Parisi, 1990), protein interactions (Bryngelson and Wolynes, 1987), financial markets (Bornholdt, 2001; Krawiecki et al., 2002), and social interactions between mammals (Daniels et al., 2016, 2017).

The model of the islet consist of N cells; at time t each of the cells can be in one of two states, spiking or silent, represented by a spin variable S_i(t) = ±1, (i = 1, …, N). The discrete time steps in model computations correspond to 2 s binning size of the Ca²⁺ data. The effective field E_i of the i-th cell has two contributions: one from the cell interacting with all other cells with interaction strength J_ij, and one from external field h_i. We assume that interactions extend over the whole system.

$\begin{array}{rc}{E}_i(t)=h_i(t)+{\displaystyle \sum _{j=1}^N}{J}_{ij}{S}_j(t)& (1)\end{array}$

At the next moment (t + 1) each cell updates its state S_i(t) with the probability p to S_i(t + 1) = +1 and with the probability 1 − p to S_i(t + 1) = −1. The probability p depends on the effective field E_i that the i-th cell senses:

$\begin{array}{rc}p=\frac{1}{1+exp(-2E_i)}.& (2)\end{array}$

The interaction strength J_ij is a fluctuating quantity with contributions from amplitude J common to all links and from the pairwise contributions with amplitude I (Krawiecki et al., 2002): J_ij = Jλ(t)+Iν_ij(t). Here are the fluctuations λ(t) and ν_ij(t) random variables uniformly distributed in the interval [−1, 1]. The external field h_i(t) = η(t) is also a random variable, uniformly distributed in the interval η(t) = [h_min, h_max]. In the mean-field approximation the average state of the system $m(t)=\frac{1}{N}\sum _j{S}_j$, evolves with time according to Krawiecki et al. (2002):

$\begin{array}{rc}m(t+1)=tanh(J\lambda (t)m(t)+h_{mf}(t)),& (3)\end{array}$

where we h_mf = η(t)/N. In the Results section below we demonstrate that the model describes the important features observed in the data well. In all computations we used a model with N = 200 spins, and we set, following the original model (Krawiecki et al., 2002), the pairwise interaction amplitude to I = 2J. The values of the remaining three free parameters of the model, J, h_max and h_min, were chosen to fit the model computations to the qualitative features of the Ca²⁺ data as described in the next section.

3. Results

The functional multicellular imaging (fMCI) records a full temporal activity trace for every cell in an optical plane of an islet from which meaningful quantitative statements about the dynamics of unitary Ca²⁺ responses and information flow in the β-cell collective are possible (Dolenšek et al., 2013; Stožer et al., 2013a). Briefly, after the stimulation with increased glucose level, first asynchronous Ca²⁺ transients appear, followed by a sustained plateau phase with oscillations on different temporal scales, from slow oscillations (100–200 s) to trains of fastest Ca²⁺ spikes (1–2 s). As the relation between the rate of insulin release and cytosolic Ca²⁺ activity shows saturation kinetics with high cooperativity (Skelin and Rupnik, 2011), the insulin release probability is significantly increased during these Ca²⁺ spikes.

Initially, fMCI has been done at the glucose concentrations much higher than those at which β-cells usually operate. The main reason for this was to ensure comparability of the results with the mainstream research in the field using mostly biochemical approaches. At 16 mM glucose, a collective of β-cells responds in a fast, synchronized, and step-like manner. Therefore the first interpretation has been that gap junction coupling between neighboring β-cells presents a major driving force for the β-cell activation and inhibitory dynamics (Hraha et al., 2014; Markovič et al., 2015). Accordingly, the removal of Cx36 proteins does cause hyperinsulinemia at resting glucose levels and blunted responses to stimulatory glucose concentration (Speier et al., 2007). Since β-cells in fresh pancreatic tissue slices are sensitive to physiological concentration of glucose (6–9 mM) (Speier and Rupnik, 2003), we here focused on this less explored concentration range. We looked at the spiking part of the Ca²⁺ imaging signals for which it has been previously shown to contain enough information to allow reconstruction of functional cell networks (Stetter et al., 2012).

For the present analysis we reused a dataset of individual Ca²⁺-dependent events from N = 188 ROIs with known positions from the central part of the fresh rodent pancreatic oval shaped islet (370 um in length and 200 um wide), representing β-cells, recorded with fMCI technique at 10 Hz over period of 40 min (for methodological details see Stožer et al., 2013b; Markovič et al., 2015; Gosak et al., 2017a,b). During the recording the glucose concentration in the solution filling the recording chamber has been increased from 6–8 mM, reaching equilibration at around 200 s after the start of the experiment, and then decreased to initial concentration near the end of experiment at around 2,000 s (dashed red lines in Figures 4, 5 represent points where glucose levels were completely equilibrated in the recording chamber). We applied ensemble empirical mode decomposition (Luukko et al., 2016) on recorded traces to isolate the Ca²⁺ spiking component of the signal. Finally, based on previous experiments in our laboratory, we binarized the signals using 2 s wide bins (Figure 1, left panel) and obtained binary spike trains S_j(t)±1, (j = 1…N), of β-cells' Ca²⁺ activity, each cell represented as a spin. As can be seen from the Figure 1 the chosen bin width adequately describe the unitary events seen in the calcium traces. An example of spiking dynamics of 30 randomly chosen spins is shown as a raster plot in the right panel of Figure 1.

FIGURE 1

Figure 1. (Left) A Ca²⁺ trace showing a short train of spikes after ensemble empirical mode decomposition with overlaid binary form with 2 s wide bins. (Right) Spin raster plot of 30 randomly picked β-cells.

Statistical methods based on mostly pairwise correlations between neurons populations have been successfully used in predicting spiking patterns in cell populations (Schneidman et al., 2006; Tkacik et al., 2009; Tkačik et al., 2014; Ferrari et al., 2017). It may seem surprising that models with first and second-order correlation structure work not only when the cell activity is very sparse so the correlations could be described by perturbation theory (Roudi et al., 2009), but can reproduce the statistics of multiple co-spiking activity (Barton and Cocco, 2013; Merchan and Nemenman, 2016; Ferrari et al., 2017). We computed truncated correlations

$\begin{array}{rc}c(i,j)=\langle S_i{S}_j\rangle -\langle S_i\rangle \langle S_j\rangle & (4)\end{array}$

for all pairs of cells. The pairwise correlations found are mostly weak with the distribution shown in Figure 2 (left panel), but they extend widely over the distances up to 170 μm across the islet, which is larger than an average vertebrate islet size (Figure 2, right panel). At distances larger than 170 μm the correlations decrease sharply toward zero. Such weak and long-ranging pairwise correlations could be the root of criticality and of strongly correlated network states in biological systems (Schneidman et al., 2006; Azhar and Bialek, 2010; Mora and Bialek, 2011; Tkačik et al., 2015).

FIGURE 2

Figure 2. (Left) Distribution of pairwise correlations of β-cell collective computed from Ca²⁺ imaging spiking signals. (Right) Pair correlations distribution over distance. Weak correlations extend over the whole system up to 170 μm. Black line shows the average values of correlations at particular cell-cell distances.

To check for the existence of strongly correlated states in weakly correlated β-cell collective we computed probability distributions P_N(K) of K simultaneously spiking cells in groups of N = 10, 20, 30 cells. Here, we used the entire dataset, the low and the high glucose concentration parts, from which we sampled cells signals. While the P_N(K) of randomly reshuffled spike trains expectedly follows Poisson distribution (left panel in Figure 3, black crosses and dashed line for N = 10 spins), the observed co-spiking probabilities are orders of magnitude higher (diamonds in left panel of Figure 3 for N = 10 spins) than corresponding probabilities in groups of independent spins. The statistics of these co-spiking events were described by an exponential distribution (Schneidman et al., 2006), by finding the effective potential (Tkačik et al., 2013, 2014) matching the observed P_N(K) and adding it to the hamiltonian of the model, or by using beta-binomial distribution (Nonnenmacher et al., 2017) P_N(K) = C(N, K)B(α + K, β + N − K)/B(α, β) where C(N, K) is binomial coefficient and B(α, β) is the beta function.

FIGURE 3

Figure 3. (Left) Probability distributions of K cells among N spiking simultaneously. Randomly shuffled spike trains (black crosses, N = 10) with dashed line - Poisson distribution; N = 10 (diamonds), N = 20 (squares), N = 30 (open dots); model (red pluses + red dashed line with J = 2.0 used for the entire dataset, N_spins = 200 spins, h_min = −2.65, h_max = −1.65), beta-binomial model (Nonnenmacher et al., 2017) (black dashed line, α = 0.38, β = 11.0); (Right) Scaling of mean field return: open dots - data, red pluses - mean field approximation from the spin model of β-cells computed with J = 2.0, h_mf = η(t)/N. Dashed line P(G) ~ G^−2.0

We next run the spin model of 200 β-cells and then sampled the computed spike trains to obtain P_N(K) from the model for N = 10, 20, 30. Despite its simple structure, the model matches order of magnitude of the observed P_N(K) well when we set the interaction strength at J = 2.0, as shown in the left panel of Figure 3 (red pluses and red dashed line), particularly for larger K values. In the model here we did not treat the low and the high glucose concentration part separately, we used J = 2.0 for the entire dataset. For comparison, we also show how the beta-binomial model fits to the observed data using the parameters α = 0.38, β = 11.0 in all N = 10, 20, 30 cases. These values are also close to the best-fitting parameters (α = 0.38, β = 12.35) to the simulated and observer correlated neural population activity data as reported in Nonnenmacher et al. (2017).

The microscopic model of interacting spins with interactions randomly varying in time (Krawiecki et al., 2002), adopted here to describe interacting β-cell collective, exhibits scaling of price fluctuations (Bornholdt, 2001) observed in financial markets (Gopikrishnan et al., 1999) and on-off intermittency with attractor bubbling dynamics of average price (Krawiecki et al., 2002). Following this idea, we looked at the logarithmic return of average state of β-cell collective at time t (Bornholdt, 2001): G(t) = log(m(t)) − log(m(t − 1)). As presented in the right panel of Figure 3, the distribution P(G) (of positive G values) can indeed be approximated with a scaling law: P(G) ~ G^−γ with γ = 2.0. There is an analytical relationship (Krawiecki et al., 2002) between J and exponent γ of the distribution of amplitudes of the return of the mean field, J = γ^1/(γ−1), which gives J = 2.0 for γ = 2.0. We used this as a consistency check between the model computations and mean field approximation. Computing the average state with the Equation (3) of the model, we can reproduce the observed distribution by setting on the interaction strength to J = 2.0 at t_on = 400 s and off to J = 0 at t_off = 2,200 s. The amplitude of the interaction J is consistent with the computation of the co-spiking probability.

In Figure 4 we show the plots of both, observed and computed, returns of average state of interacting β-cells for comparison. The glucose concentration was changed during the experiment in a stepwise manner: from 6 to 8 mM at the beginning and back to 6 mM near the end of recording period. The effect of both changes is nicely visible in the G(t) plot (upper panel, Figure 4) where the on-off intermittent dynamics of the average state starts around t_on = 400 s and lasts until around t_off = 2, 200 s in the experiment. Both observed events are delayed with respect to the times of glucose concentration change due to the asynchronous Ca²⁺ transients (Stožer et al., 2013a). We expect that the response of β-cell collective to the stimulus increase must be visible in the variance of average state Var(m) which is in Ising-like model we are using here equal to susceptibility of the system χ = Var(m) = 〈m²〉−〈m〉². We used the low glucose concentration part (6 mM) of the data to estimate the boundaries of the external field interval [h_max, h_min] to describe the first part of susceptibility. Using the maximal and minimal spiking rates of cells (m_max, m_min) in 6 mM glucose from the data and the mean-field approximation with J = 0 corresponding to the non-stimulatory glucose regime we have $\left[h_{max},h_{min}\right]=\left[{tanh}^{-1}(m_{max}),{tanh}^{-1}(m_{min})\right]=\left[-1.65,-2.65\right]$. In upper panel of Figure 5 (open black dots) we show the plot of susceptibility as a function of recording time, focusing around the transition to increased glucose concentration during the experiment. There is a sharp increase of susceptibility at around t_on, the same time the on-off intermittency starts to appear in G(t). Using mean field approximation of the spin model Equation (3) for computation of susceptibility (averaged over many runs) and setting J = 0 for t<t_on and J = 2.0 for t>t_on we can well describe the observed evolution of susceptibility and capture the rapid onset of increased sensibility of the islet (red line in upper part of Figure 5).

FIGURE 4

Figure 4. (Upper) Observed logarithmic return of the average state of of β-cell collective G(t), (Lower) logarithmic return of the average state computed from the model with J = 2.0 for t_on < t < t_off, denoted with vertical dashed lines in figures. Dashed red lines represent points where glucose levels were completely equilibrated in the recording chamber.

FIGURE 5

Figure 5. (Upper) Susceptibility of β-cell collective around transition to stimulatory glucose level. Open dots are the experimental data, red line shows the result of the mean field computations with J = 2.0 onset at t = t_on(blackdashedline). Dashed red line represent the point where glucose level completely equilibrated during the 6–8 mM transition in the recording chamber. (Lower) Normalized conditional entropy. Open dots are experimental data at 8 mM glc, open squares at 6 mM glc. Red pluses show the results of the spin model computations with N_spins = 200 spins, and the parameters: h_min = −2.65, h_max = −1.65, J = 2.0 for the upper, and J = 0 for lower the lower part.

Pairwise correlation structure enables error-correction features of population coding in neural systems (Schneidman et al., 2006). To check for memory-like or error-correcting properties in islets, we use the conditional entropy H(S_i|S), the measure for the information we need to determine the state of N-th cell (i.e., spiking or not) if we know the states of N − 1 cells (S = S_j ≠ i) in a group of N cells. If the state of the N-th cell is completely determined by other N − 1 cells, the conditional entropy is zero H(S_i|S) = 0 and the error correction is perfect. When S_j are independent random states, the conditional entropy equals the entropy of the N-th cell H(S_i).

We computed the quantity 1 − H(S_i|S)/H(S_i) (normalized mutual information) as a function of number of cells (for small groups of cells) and extrapolate the trend toward the limit H(S_i|S) = 0 that determines the critical number of cells, N_c, needed to predict the state of another cell. As seen in the lower panel of Figure 5, the predictability is a glucose-dependent parameter. With non-stimulatory glucose concentration, the complete set of data is required for predictions, whereas at 8 mM glucose we find that order of magnitude smaller number of measured cells are needed to predict the states of other cells.

4. Discussion

Pancreatic β-cell continuously intercepts a variety of energy-rich or signaling ligands using the whole spectrum of specific receptors on the cell membrane, as well as in metabolic and signaling pathways within the cell. The cell converts these signals into a binary cellular code, for example a train of Ca²⁺ spikes, which drive insulin release that fits current physiological needs of the body. This allow already a single cell to sense its chemical environment with extraordinary, often diffusion limited precision (Bialek and Setayeshgar, 2005), however, judging by their heterogeneous secretory behavior in cell culture, the precision of sensing among the individual β-cells is quite diverse (Hiriart and Ramirez-Medeles, 1991). Recent experimental evidence and modeling have shown that cell collectives sense better compared to an individual cell. The precise mechanism of this collective sensing improvement depends on cell-cell communication type, which can be short-range with direct cell contacts or long-range with paracrine signaling (Fancher and Mugler, 2017; Saakian, 2017). Furthermore, also long-range interaction have its finite reach which can poise a limit to the cell collective size and therefore determines its optimal as well as maximal size. As mentioned in the Introduction here, it is intriguing how well conserved the pancreatic islet size is in vertebrates of dramatically different body dimensions (Montanya et al., 2000). In a single vertebrate organism the size of the islets can be bigger that 150 um, but functional studies revealed that the islets bigger than 200 um secrete 50% less insulin after glucose stimulation (Fujita et al., 2011). These functional differences between small and large islets have been partially attributed to diffusion barriers for oxygenation and nutrition, limiting the survival of core β-cells in bigger islets after isolation. However, reducing these diffusion barriers had no influence on insulin secretory capacity (Williams et al., 2010) suggesting that other factors, like diffusion of paracrine signaling molecules (Caicedo, 2013), could limit the collective β-cell function in bigger islets. This dominance of a long-range information flow, likely limited to some physical constraints, has indicated the use of the mathematical equivalency with spin glass-like systems (Tkačik and Bialek, 2016).

We strongly believe that advanced complex network analysis based on strong short-range correlations can continue to provide valuable information regarding the β-cell network topologies, network on network interactions and describe the functional heterogeneity of individual β-cells (Gosak et al., 2015, 2017a; Markovič et al., 2015; Johnston et al., 2016). However, the main goal of the present study was to determine the influence of weak long-range correlations between pairs of β-cells on the probability of activation of single β-cells. Recently has been shown that it suffice to use pairwise correlations to quantitatively describe the collective behavior of cell collectives (Merchan and Nemenman, 2016). The typically small values of pairwise correlation coefficients with the median values below 0.02, would intuitively be ignored and β-cells described as if they act independently, however in larger populations of cells this assumption completely fails (Schneidman et al., 2006). In fact, at physiological stimulatory glucose levels between 6 and 9 mM, β-cell collectives are entirely dominated by weak average pairwise correlations (Figure 2). Nevertheless, this is the glucose concentration range, where β-cells are most responsive to the nutrient to, as a collective, compute their activity state and pulsatile insulin release, and to meet the organismal needs between the environmental and behavioral extremes of food shortage and excess (Schmitz et al., 2008)?

Based on the range of the calculated weak pairwise correlations of up to 170 um (Figure 2), we predict that β-cells collective falls into a category of sparse packed tissues with dominant paracrine interactions and that cell-cell distances contribute to optimal sensing and functional response in creating the metabolic code governing the release of insulin. It remains unclear whether and how the position of β-cells within an islet is controllable. As many other cells, β-cells are polarized and possess a primary cilium (Gan et al., 2017), which should have a primary role in sensory function, i.e., insulin sensing in paracrine signaling (Doğaner et al., 2016), and not in cell motility. It is quite interesting though, that the ciliopathies are highly associated with reduced β-cell function and increased susceptibility to diabetes mellitus (Gerdes et al., 2014). Future experiments are required to test for the possible motility of β-cells within the islet to adopt an optimal separation of key sensitive β-cells. To further extrapolate the collective sensing idea, it is also possible that the diffuse arrangement of a collective of islets within different parts of pancreas, which are exposed to different vascular inputs (Dolenšek et al., 2015), serves to optimize nutrient sensing experience, yet on a higher organizational level, providing a topological information regarding the nutrient levels in different parts of the gastrointestinal tract. The nature and level of interactions between individual islets in the pancreas are currently also unknown.

As in retinal neuron networks, β-cells encode information about the presence of energy-rich nutrients into sequences of intermittent Ca²⁺ spikes. In a natural setting of sensory neural networks with stimuli derived from a space with very high dimensionality the coding seems challenging and interpretations require some strong assumptions (Tkačik et al., 2014). We currently do not understand the input dimensionality of a typical ligand mixture around the β-cells, we simply assume it is not high. As in retinal networks (Schneidman et al., 2006; Tkačik et al., 2014), the predictability regarding the functional state of individual β-cells is defined by the network and not the chemical environment. This suggests that the sensory information at physiological glucose levels is substantially redundant. It is likely that the nutrient mixture presents a noisy challenge for the information transfer which is typical for biological system. But why study the insulin release pattern or the metabolic code? The β-cell network possess associative or error-correcting properties (Figure 5), so this idea from the sensory neuron networks can be generalized also to populations of endocrine cells (Schneidman et al., 2006), which may again influence the optimal islet size and suggest the presence of functional subunits within the islet that could adapt, for example, to changing environment in a dynamic fashion. Furthermore, error-correction properties are glucose dependent and can be physiologically modulated (Figure 5). The trains of Ca²⁺ spikes at constant glucose stimulation (8 mM) are inhomogeneous, display on-off intermittency (Figure 4) and scaling of log returns of average state (Figure 3) analog to models of financial time series (Krawiecki et al., 2002). For the spin glass approach we also postulate that the sources of stochasticity in an islet collective are various. On one hand, the β-cells make decisions on activation under the influence of the external environment and other β-cells. Second, also the time-dependent interaction strength among β-cells is random, which could reflect their socio-cellular communication network and indicate that the external environment can be sensed differently between different β-cells in an islet (Gosak et al., 2017b).

Biological systems seem to poise themselves at criticality, with a major advantage of enhanced reactivity to external perturbations (Mora and Bialek, 2011). Often a limited number of individual functional entities, cell or groups of cells as found in pancreatic islets, appeared to be limiting to address criticality. However, it has been recently demonstrated that even in biological systems with small number of interacting entities one can operationally define criticality and observe changes in robustness and sensitivity of adaptive collective behavior (Daniels et al., 2017). Our results suggest that β-cells collective within the islet sits near its critical point and we could determine the susceptibility in the islet. Stimulatory glucose concentration (8 mM) has been decreasing distance to criticality by increasing sensitivity (Figure 5). Smaller distance to criticality at unphysiologically high glucose levels has its possible adverse consequences in a phenomenon called critical slowing down as the system takes more and more time to relax as it comes nearer to the critical point (Mora and Bialek, 2011). Our preliminary results show that at very strong stimulation (i.e., glucose levels above 12 mM) the whole system freezes into a certain state where short-term interaction take over enabling global phenomena within the islets, e.g., Ca²⁺ waves (Stožer et al., 2013b) requiring progressively longer periods to relax to baseline with increasing glucose concentrations.

Further work will be needed to exploit at what circumstances deviations in islet size can contribute to islet malfunction and pathogeneses of different forms of diabetes mellitus. Until recently it has been thought that insulin release is no longer functional in type 1 diabetes mellitus. We now know that even in type 1 diabetic patients small and functional collectives of β-cells persist in the pancreata of these patients even decades after the diagnosis (Faustman, 2014). On the other hand, the β-cells mass in an islet can be increased in type 2 diabetic patients in the initial phases after the diagnosis (Rahier et al., 2008) or in animal models (Daraio et al., 2017) and can only be reduced in the later phases (Rahier et al., 2008). The detailed relations between the reduced or increased insulin release, changed islet size and therefore changed circumstances for paracrine signaling in disturbed collective nutrient sensing and during the aforementioned pathogeneses of diabetes mellitus remain to be established.

Author Contributions

All authors listed, have made substantial, direct and intellectual contribution to the work, and approved it for publication.

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.

Acknowledgments

The authors acknowledge the financial support from the Slovenian Research Agency (research core funding, No. P3-0396), as well as research project, No. N3-0048).

References

Ashcroft, F. M., and Rorsman, P. (1989). Electrophysiology of the pancreatic β-cell. Prog. Biophys. Mol. Biol. 54, 87–143. doi: 10.1016/0079-6107(89)90013-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Azhar, F., and Bialek, W. (2010). When are correlations strong? arXiv preprint arXiv:1012.5987.

Google Scholar

Barton, J., and Cocco, S. (2013). Ising models for neural activity inferred via selective cluster expansion: structural and coding properties. J. Stat. Mech. Theory Exp. 2013:P03002. doi: 10.1088/1742-5468/2013/03/P03002

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

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

Bialek, W., and Setayeshgar, S. (2005). Physical limits to biochemical signaling. Proc. Natl. Acad. Sci. U.S.A. 102, 10040–10045. doi: 10.1073/pnas.0504321102

PubMed Abstract | CrossRef Full Text | Google Scholar

Bornholdt, S. (2001). Expectation bubbles in a spin model of markets: intermittency from frustration across scales. Int. J. Mod. Phys.C 12, 667–674. doi: 10.1142/S0129183101001845

CrossRef Full Text | Google Scholar

Bouwens, L., and Rooman, I. (2005). Regulation of pancreatic beta-cell mass. Physiol. Rev. 85, 1255–1270. doi: 10.1152/physrev.00025.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Bryngelson, J. D., and Wolynes, P. G. (1987). Spin glasses and the statistical mechanics of protein folding. Proc. Natl. Acad. Sci. U.S.A. 84, 7524–7528. doi: 10.1073/pnas.84.21.7524

PubMed Abstract | CrossRef Full Text | Google Scholar

Caicedo, A. (2013). Paracrine and autocrine interactions in the human islet: more than meets the eye. Semin. Cell Dev. Biol. 24, 11–21. doi: 10.1016/j.semcdb.2012.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Cigliola, V., Chellakudam, V., Arabieter, W., and Meda, P. (2013). Connexins and β-cell functions. Diabetes Res. Clin. Pract. 99, 250–259. doi: 10.1016/j.diabres.2012.10.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Daniels, B. C., Ellison, C. J., Krakauer, D. C., and Flack, J. C. (2016). Quantifying collectivity. Curr. Opin. Neurobiol. 37, 106–113. doi: 10.1016/j.conb.2016.01.012

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Daraio, T., Bombek, L. K., Gosak, M., Valladolid-Acebes, I., Klemen, M. S., Refai, E., et al. (2017). Snap-25b-deficiency increases insulin secretion and changes spatiotemporal profile of ca2+ oscillations in β cell networks. Sci. Rep. 7:7744. doi: 10.1038/s41598-017-08082-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Doğaner, B. A., Yan, L. K., and Youk, H. (2016). Autocrine signaling and quorum sensing: extreme ends of a common spectrum. Trends Cell Biol. 26, 262–271. doi: 10.1016/j.tcb.2015.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Fancher, S., and Mugler, A. (2017). Fundamental limits to collective concentration sensing in cell populations. Phys. Rev. Lett. 118:078101. doi: 10.1103/PhysRevLett.118.078101

PubMed Abstract | CrossRef Full Text | Google Scholar

Faustman, D. L. (2014). Why were we wrong for so long? the pancreas of type 1 diabetic patients commonly functions for decades. Diabetologia 57, 1–3. doi: 10.1007/s00125-013-3104-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferrari, U., Obuchi, T., and Mora, T. (2017). Random versus maximum entropy models of neural population activity. Phys. Rev. E 95:042321. doi: 10.1103/PhysRevE.95.042321

PubMed Abstract | CrossRef Full Text | Google Scholar

Fujita, Y., Takita, M., Shimoda, M., Itoh, T., Sugimoto, K., Noguchi, H., et al. (2011). Large human islets secrete less insulin per islet equivalent than smaller islets in vitro. Islets 3, 1–5. doi: 10.4161/isl.3.1.14131

PubMed Abstract | CrossRef Full Text | Google Scholar

Gan, W. J., Zavortink, M., Ludick, C., Templin, R., Webb, R., Webb, R., et al. (2017). Cell polarity defines three distinct domains in pancreatic β-cells. J. Cell Sci. 130, 143–151. doi: 10.1242/jcs.185116

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

Gopikrishnan, P., Plerou, V., Amaral, L. A. N., Meyer, M., and Stanley, H. E. (1999). Scaling of the distribution of fluctuations of financial market indices. Phys. Rev. E 60:5305. doi: 10.1103/PhysRevE.60.5305

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Gosak, M., Markovič, R., Dolenšek, J., Rupnik, M. S., Marhl, M., Stožer, A., et al. (2017a). Network science of biological systems at different scales: a review. Phys. Life Rev. doi: 10.1016/j.plrev.2017.11.003. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

Gosak, M., Stožer, A., Markovič, R., Dolenšek, J., Perc, M., Slak Rupnik, M., et al. (2017b). 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

Henquin, J.-C. (2009). Regulation of insulin secretion: a matter of phase control and amplitude modulation. Diabetologia 52:739. doi: 10.1007/s00125-009-1314-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Henquin, J.-C. (2012). Do pancreatic β cells taste nutrients to secrete insulin? Sci. Signal. 5:pe36. doi: 10.1126/scisignal.2003325

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

Islam, M. S. (2002). The ryanodine receptor calcium channel of β-cells. Diabetes 51, 1299–1309. doi: 10.2337/diabetes.51.5.1299

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, A., Miller, K., Jo, J., Kilimnik, G., Wojcik, P., and Hara, M. (2009). Islet architecture: a comparative study. Islets 1, 129–136. doi: 10.4161/isl.1.2.9480

PubMed Abstract | CrossRef Full Text | Google Scholar

Krawiecki, A., Hołyst, J. A., and Helbing, D. (2002). Volatility clustering and scaling for financial time series due to attractor bubbling. Phys. Rev. Lett. 89:158701. doi: 10.1103/PhysRevLett.89.158701

PubMed Abstract | CrossRef Full Text | Google Scholar

Lamprianou, S., Immonen, R., Nabuurs, C., Gjinovci, A., Vinet, L., Montet, X. C., et al. (2011). High-resolution magnetic resonance imaging quantitatively detects individual pancreatic islets. Diabetes 60, 2853–2860. doi: 10.2337/db11-0726

PubMed Abstract | CrossRef Full Text | Google Scholar

Lang, J. (1999). Molecular mechanisms and regulation of insulin exocytosis as a paradigm of endocrine secretion. FEBS J. 259, 3–17. doi: 10.1046/j.1432-1327.1999.00043.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Luukko, P., Helske, J., and Räsänen, E. (2016). Introducing libeemd: a program package for performing the ensemble empirical mode decomposition. Comput. Stat. 31, 545–557. doi: 10.1007/s00180-015-0603-9

CrossRef Full Text | Google Scholar

Mandic, S. A., Skelin, M., Johansson, J. U., Rupnik, M. S., Berggren, P.-O., and Bark, C. (2011). Munc18-1 and munc18-2 proteins modulate β-cell ca2+ sensitivity and kinetics of insulin exocytosis differently. J. Biol. Chem. 286, 28026–28040. doi: 10.1074/jbc.M111.235366

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

Merchan, L., and Nemenman, I. (2016). On the sufficiency of pairwise interactions in maximum entropy models of networks. J. Stat. Phys. 162, 1294–1308. doi: 10.1007/s10955-016-1456-5

CrossRef Full Text | Google Scholar

Montanya, E., Nacher, V., Biarnés, M., and Soler, J. (2000). Linear correlation between beta-cell mass and body weight throughout the lifespan in lewis rats: role of beta-cell hyperplasia and hypertrophy. Diabetes 49, 1341–1346. doi: 10.2337/diabetes.49.8.1341

PubMed Abstract | CrossRef Full Text | Google Scholar

Mora, T., and Bialek, W. (2011). Are biological systems poised at criticality? J. Stat. Phys. 144, 268–302. doi: 10.1007/s10955-011-0229-4

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

Parisi, G. (1990). A simple model for the immune network. Proc. Natl. Acad. Sci. U.S.A. 87, 429–433. doi: 10.1073/pnas.87.1.429

PubMed Abstract | CrossRef Full Text | Google Scholar

Paulmann, N., Grohmann, M., Voigt, J.-P., Bert, B., Vowinckel, J., Bader, M., et al. (2009). Intracellular serotonin modulates insulin secretion from pancreatic β-cells by protein serotonylation. PLoS Biol. 7:e1000229. doi: 10.1371/journal.pbio.1000229

PubMed Abstract | CrossRef Full Text | Google Scholar

Rahier, J., Guiot, Y., Goebbels, R. M., Sempoux, C., and Henquin, J.-C. (2008). Pancreatic β-cell mass in european subjects with type 2 diabetes. Diabetes Obes. Metab. 10, 32–42. doi: 10.1111/j.1463-1326.2008.00969.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Rasmussen, H. (1970). Cell communication, calcium ion, and cyclic adenosine monophosphate. Science 170, 404–412. doi: 10.1126/science.170.3956.404

PubMed Abstract | CrossRef Full Text | Google Scholar

Roudi, Y., Nirenberg, S., and Latham, P. E. (2009). Pairwise maximum entropy models for studying large biological systems: when they can work and when they can't. PLoS Comput. Biol. 5:e1000380. doi: 10.1371/journal.pcbi.1000380

PubMed Abstract | CrossRef Full Text | Google Scholar

Saakian, D. B. (2017). Kinetics of biochemical sensing by single cells and populations of cells. Phys. Rev. E 96:042413. doi: 10.1103/PhysRevE.96.042413

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Schneidman, E., Berry, M. J. II., Segev, R., and Bialek, W. (2006). Weak pairwise correlations imply strongly correlated network states in a neural population. Nature 440:1007. doi: 10.1038/nature04701

PubMed Abstract | CrossRef Full Text | Google Scholar

Skelin, M., and Rupnik, M. (2011). Camp increases the sensitivity of exocytosis to ca 2+ primarily through protein kinase a in mouse pancreatic beta cells. Cell Calcium 49, 89–99. doi: 10.1016/j.ceca.2010.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Speier, S., Gjinovci, A., Charollais, A., Meda, P., and Rupnik, M. (2007). Cx36-mediated coupling reduces β-cell heterogeneity, confines the stimulating glucose concentration range, and affects insulin release kinetics. Diabetes 56, 1078–1086. doi: 10.2337/db06-0232

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Stetter, O., Battaglia, D., Soriano, J., and Geisel, T. (2012). Model-free reconstruction of excitatory neuronal connectivity from calcium imaging signals. PLoS Comput. Biol. 8:e1002653. doi: 10.1371/journal.pcbi.1002653

PubMed Abstract | CrossRef Full Text | Google Scholar

Stožer, A., Dolenšek, J., and Rupnik, M. S. (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

Tkačik, G., and Bialek, W. (2016). Information processing in living systems. Annu. Rev. Condens. Matt. Phys. 7, 89–117. doi: 10.1146/annurev-conmatphys-031214-014803

CrossRef Full Text | Google Scholar

Tkačik, G., Marre, O., Amodei, D., Schneidman, E., Bialek, W., and Berry, M. J. II. (2014). Searching for collective behavior in a large network of sensory neurons. PLoS Comput. Biol. 10:e1003408. doi: 10.1371/journal.pcbi.1003408

PubMed Abstract | CrossRef Full Text | Google Scholar

Tkačik, G., Marre, O., Mora, T., Amodei, D., Berry, M. J. II., and Bialek, W. (2013). The simplest maximum entropy model for collective behavior in a neural network. J. Stat. Mech. Theory Exp. 2013:P03011.

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

Tkacik, G., Schneidman, E., Berry, I., Michael, J., and Bialek, W. (2009). Spin glass models for a network of real neurons. arXiv preprint arXiv:0912.5409.

Google Scholar

Tomita, T. (1976). Phlorizin: its effect on glucose-induced insulin secretion and protection against the alloxan effect in isolated islets. FEBS Lett. 65, 140–143. doi: 10.1016/0014-5793(76)80465-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Tomkins, G. M. (1975). The metabolic code. Science 189, 760–763. doi: 10.1126/science.169570

PubMed Abstract | CrossRef Full Text | Google Scholar

Trautmann, M. E., and Wollheim, C. B. (1987). Characterization of glucose transport in an insulin-secreting cell line. Biochem. J. 242, 625–630. doi: 10.1042/bj2420625

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, S. J., Huang, H.-H., Kover, K., Moore, W. V., Berkland, C., Singh, M., et al. (2010). Reduction of diffusion barriers in isolated rat islets improves survival, but not insulin secretion or transplantation outcome. Organogenesis 6, 115–124. doi: 10.4161/org.6.2.10373

CrossRef Full Text | Google Scholar