ORIGINAL RESEARCH article
Heterogeneity and Delayed Activation as Hallmarks of Self-Organization and Criticality in Excitable Tissue
- 1Faculty of Medicine, University of Maribor, Maribor, Slovenia
- 2Faculty of Natural Sciences and Mathematics, University of Maribor, Maribor, Slovenia
- 3Faculty of Education, University of Maribor, Maribor, Slovenia
- 4Faculty of Energy Technology, University of Maribor, Krško, Slovenia
- 5Center for Applied Mathematics and Theoretical Physics, University of Maribor, Maribor, Slovenia
- 6Complexity Science Hub Vienna, Vienna, Austria
- 7Institute of Physiology and Pharmacology, Medical University of Vienna, Vienna, Austria
- 8Alma Mater Europaea – ECM, Maribor, Slovenia
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.
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 Ca2+ waves. Jung et al. (1998) reported a power law distribution of noise-induced spiral Ca2+ wave sizes in cultured networks of astrocytes. Heavy-tailed distributions and an avalanche-like behavior have also been observed in intracellular Ca2+ signalization in cardiac myocytes (Nivala et al., 2012) and in immature oocytes (Lopez et al., 2012). In both studies Ca2+ waves in individual cells resulted from random local Ca2+ events, reflecting small Ca2+ release events from individual channels or a cluster of channels, which can occasionally integrate to global events reflecting a whole-cell Ca2+ signal (Berridge et al., 2000). As the localized subcellular Ca2+ 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 Ca2+ 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 Ca2+ concentration ([Ca2+]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 [Ca2+]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 [Ca2+]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 Ca2+ 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 CaCl2, 1 MgCl2; 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 NaHCO3, 6 lactic acid, 3 myo-inositol, 2.5 KCl, 2 Na-pyruvate, 2 CaCl2, 1.25 NaH2PO4, 1 MgCl2, 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/F0 ratio were then further processed as explained in the continuation.
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):
where ui(t) and vi(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 ui(t) describes the dynamics of the membrane potential of the cell, whereas the slow variable vi(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, [Ca2+]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 [Ca2+]IC dynamics, this first phase consists of a transient increase in [Ca2+]IC, during the ascending part of which a few fast [Ca2+]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 [Ca2+]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 [Ca2+]IC change to regular bursting and corresponding fast [Ca2+]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 [Ca2+]IC oscillations (Gilon et al., 1993; Bergsten et al., 1994). This fast electrical, [Ca2+]IC, and secretory pattern is superimposed on a slower set of oscillations in ATP, membrane potential, [Ca2+]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 [Ca2+]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 [Ca2+]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, λ1 and λ2 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)| < 1 and |Re(λ2)| < 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 KATP-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. Dynamical features of the Rulkov map. (A) Real and imaginary eigenvalues λ1 and λ2 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. gi 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 , where ⟨k⟩ = 6 signifies the average number of connections per cell. A typical intercellular network structure is shown in Figure 2A.
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, , ts = 10000) and suprastimulatory (Δα = 0.09, A = 0.70, B = 0.0008, , ts = 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:
In this manner, we took into account the delay due to glucose metabolism. In Eq. (3) α0 = 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 ts 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, [Ca2+]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 Tm,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 and . 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, [Ca2+]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 [Ca2+]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 Tm,i), (ii)stimulation-dependent increase in excitability (parameter Δαi), and (iii) intercellular coupling strength (parameter gi). 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. A schematic representation of all three types of heterogeneities in the model. (A) electrical excitability Δαi, (B) beta cell metabolism Tm,i, (C) electrical coupling gi. 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 Ca2+ 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 Ca2+ 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 Ca2+ spikes, tSTART. 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, tEND. Lastly we test, if the local maximum of the Ca2+ signal within the time interval t ∈ [tSTART,tEND] satisfies the condition x(t) > 1.5std [x]. The corresponding binary time series takes on a value of 1 in all time intervals t ∈ [tSTART,tEND], 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).
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.
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 [Ca2+]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 Ca2+ 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 [Ca2+]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 Ca2+ waves.
Figure 4. Experimentally measured beta cell responses after stimulation with 8 mM glucose. (A) Three characteristic Ca2+ traces and the raster plot of binarized Ca2+ 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 Ca2+ activity waveforms for selected intervals for the activation (B) and plateau (D) phase. Colors denote specific Ca2+ events. Gray dots on the y, z plane stand for coordinates of cells. (C,E) The distributions of Ca2+ 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. Experimentally measured beta cell responses after stimulation with 12 mM glucose. (A) Three characteristic Ca2+ traces and the raster plot of binarized signals of Ca2+ 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 Ca2+ waves for selected intervals for the activation (B) and plateau (D) phase. Colors denote specific Ca2+ events. Gray dots on the y, z plane stand for coordinates of cells. (C,E) The distributions of Ca2+ wave sizes for the activation (C) and plateau (E) phase. The gray dashed line indicates the power-law fit.
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. 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. 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. 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.
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. 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.
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 Ca2+ 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 [Ca2+]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 [Ca2+]IC waves showing supercritical behavior. We wish to speculate that recruitment (Stožer et al., 2013a; Gosak et al., 2017) and local [Ca2+]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 Ca2+ 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 [Ca2+]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 [Ca2+]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 [Ca2+]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 [Ca2+]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 Ca2+ 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 [Ca2+]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.
The datasets generated for this study are available on request to the corresponding author.
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).
AS, RM, JD, MP, MSR, MM, and MG designed the study, carried out the research, and wrote the manuscript.
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.
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
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
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
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
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
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
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
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
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
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
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.
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
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
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
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
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
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
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
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
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
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
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
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.
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
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
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
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
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
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
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
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
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
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
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
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
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.1922.214.171.1248
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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].
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
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
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
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
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
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
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
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
Keywords: excitable cells, self-organized criticality, beta cells, calcium imaging, computational model, cellular heterogeneity, activation delay
Citation: Stožer A, Markovič R, Dolenšek J, Perc M, Marhl M, Slak Rupnik M and Gosak M (2019) Heterogeneity and Delayed Activation as Hallmarks of Self-Organization and Criticality in Excitable Tissue. Front. Physiol. 10:869. doi: 10.3389/fphys.2019.00869
Received: 20 February 2019; Accepted: 21 June 2019;
Published: 05 July 2019.
Edited by:Francisco Monroy, Complutense University of Madrid, Spain
Reviewed by:Nir Shachna Gov, Weizmann Institute of Science, Israel
Krasimira Tsaneva-Atanasova, University of Exeter, United Kingdom
Copyright © 2019 Stožer, Markovič, Dolenšek, Perc, Marhl, Slak Rupnik and Gosak. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.