Random Matrix Analysis of Ca2+ Signals in β-Cell Collectives

Even within small organs like pancreatic islets, different endocrine cell types and subtypes form a heterogeneous collective to sense the chemical composition of the extracellular solution and compute an adequate hormonal output. Erroneous cellular processing and hormonal output due to challenged heterogeneity result in various disorders with diabetes mellitus as a flagship metabolic disease. Here we attempt to address the aforementioned functional heterogeneity with comparing pairwise cell-cell cross-correlations obtained from simultaneous measurements of cytosolic calcium responses in hundreds of islet cells in an optical plane to statistical properties of correlations predicted by the random matrix theory (RMT). We find that the bulk of the empirical eigenvalue spectrum is almost completely described by RMT prediction, however, the deviating eigenvalues that exist below and above RMT spectral edges suggest that there are local and extended modes driving the correlations. We also show that empirical nearest neighbor spacing of eigenvalues follows universal RMT properties regardless of glucose stimulation, but that number variance displays clear separation from RMT prediction and can differentiate between empirical spectra obtained under non-stimulated and stimulated conditions. We suggest that RMT approach provides a sensitive tool to assess the functional cell heterogeneity and its effects on the spatio-temporal dynamics of a collective of beta cells in pancreatic islets in physiological resting and stimulatory conditions, beyond the current limitations of molecular and cellular biology.

Introduction

Pancreatic islets are collectives of endocrine cells. Based on the end hormone that these cells exocytose in a Ca²⁺-dependent manner after being stimulated, several types of cells have been described to compose an islet, with 3 dominant cell types: alpha, beta and delta (Briant et al., 2017). Islets from different parts of pancreas can contain different fractions of each of these cell types, but the bulk cellular mass in a typical islet is in an non-diabetic organism composed of a collectives of insulin secreting beta cells (Dolenšek et al., 2015; Rorsman and Ashcroft, 2017). Early studies assumed these beta cell collectives to be a rather homogeneous population of cells, however, the subsequent functional analyses have revealed a remarkable degree of heterogeneity even in dissociated beta cells in culture. The beta cells were found to differ in a number of physiological parameters, among others in glucose sensitivity and Ca²⁺ oscillation pattern (Zhang et al., 2003), electrical properties (Misler et al., 1986), redox states (Kiekens et al., 1992), or pattern on cAMP oscillations (Dyachok et al., 2006). These early quests (Pipeleers, 1992) have been mostly searching for morphological, physiological and molecular features that would presumably satisfy at least 3 criteria: (a) entitle special roles for individual cells within the collectives, (b) remain valid even after cell dissociation, and (c) enable to trace embryonic and postnatal development as well as changes during pathogeneses of different forms of diabetes. Recent onset of efficient high-throughput analyses has catapulted these approaches on mostly dissociated cells to a completely new level and enabled identification of a multitude of functional and non-functional subpopulations with their functional characteristics, gene and protein expression, incidences and diabetes-related changes (for a recent review refer to Benninger and Hodson, 2018). A major limitation of these analyses is that the subpopulations described in different studies have relatively little in common and currently their translational relevance is weak. This is, however, not surprising since these approaches primarily deal with sample averages and present merely a small number of discrete snapshots of a very dynamic complex activity. Nevertheless, they turned out to be extremely useful in initial attempts to construct a pseudotime map of the sequence of events in pancreatic endocrine system (Damond et al., 2019) and its interaction with the immune system (Wang et al., 2019) during the progression of type 1 diabetes mellitus (T1D) in humans. Still, better tools are needed to assess the general rules underlying functional heterogeneity in a real-time spatio-temporal dynamics within collectives of beta cells or collectives of any other cell types in an organism. Instead of an ultra-reductionist approach to precisely associate individual molecular markers yielding more or less random correlations in respect to functional heterogeneity, a network of cells is first segmented to their level of functional influence expressed in high deviating eigenvalues to lead a more efficient further discovery of few collective parameters which may also have an identifiable molecular signature. Pancreatic islets have been described as broad scale network (Stožer et al., 2013b) and the search for important principal components is well justified.

The fact is that most of what we know about pancreatic beta cells has been gained by studying dissociated beta cells in cell culture. Therefore, even phenomena that can only be observed in isolated groups of electrically coupled beta cells, like electrical activity (Rorsman and Trube, 1986) or cytosolic Ca²⁺ oscillation, are currently still mostly modeled within the framework of a single cell excitability (Sherman et al., 1988; Bertram et al., 2018). However, each beta cell interacts with several immediate and distal neighboring cells in a pancreatic islet, implicating high-ordered interactions between a large number of elements. Therefore, there is a rich exchange of signals within such a beta cell collective, both through direct cell-cell coupling (Bavamian et al., 2007) as well as through paracrine signaling (Caicedo, 2013; Capozzi et al., 2019). Such an organization necessarily yields a complex response patterns of cell activity after stimulation with physiological or pharmacological stimuli. Until recently, the richness of the aforementioned cell-cell interactions also could not be experimentally documented. However, recent technological advancements made it possible to use various optical tools to address these issues (Frank et al., 2018). For example, the functional multicellular imaging (fMCI) enabled completely new insights into our understanding of a beta cell in an islet as a biological network (Dolenšek et al., 2013; Stožer et al., 2013a,b). The dynamics of a measurable physiological parameter can namely be recorded in hundreds of beta cells within their intact environment (Speier and Rupnik, 2003; Marciniak et al., 2014) simultaneously. The measured oscillatory cytosolic Ca²⁺ concentration changes, which are required to drive insulin release turned out to be a practical tool to trace cellular activity and fundamental to study their interactions in such big collectives (Dolenšek et al., 2013; Stožer et al., 2013a,b). With the use of the tools of statistical physics we and others reconstructed, for example the complex network topologies in beta cell activation, activity and deactivation during transient glucose challenges (Stožer et al., 2013b; Gosak et al., 2015; Markovič et al., 2015; Johnston et al., 2016). As in some other, previously analyzed biological systems, also for the pancreatic islets, the minimal model incorporating pairwise interactions provides accurate predictions about the collective effects (Schneidman et al., 2006; Korošak and Slak Rupnik, 2018).

Along these lines we have recently shown that beta cell collectives work as a broad-scale complex networks (Stožer et al., 2013b; Markovič et al., 2015; Gosak et al., 2017a), sharing similarities in global statistical features and structural design principles with internet and social networks (Milo et al., 2002; Barabási and Márton Pósfai, 2016; Daniels et al., 2016; Perc et al., 2017; Duh et al., 2018). In addition to complex network description when strong cell-cell interaction are primarily taken into account, the analyses of weak pairwise interaction enabled us to use a spin glass model (Korošak and Slak Rupnik, 2018), as well as the assessment of self organized criticality (Bak, 2013; Marković and Gros, 2014; Gosak et al., 2017b; Stožer et al., 2019), also often found in biological samples (Schneidman et al., 2006). The important result from these functional studies is that a faulty pattern of hormone release due to deviating numbers of individual cell types or changes in their function lead to one of the forms of a large family of metabolic diseases called diabetes mellitus (American Diabetes Association, 2014; Pipeleers et al., 2017; Skelin Klemen et al., 2017; Nasteska and Hodson, 2018; Capozzi et al., 2019).

The basic object that we study here is the correlation matrix C with elements computed from measured Ca²⁺ signals:

$\begin{array}{ll}Corr(y_i,y_j)=C_{ij}=\frac{\langle y_i{y}_j\rangle -\langle y_i\rangle \langle y_j\rangle }{{\sigma }_i{\sigma }_j},& (1)\end{array}$

where y_i(t) is the i-th time series of Ca²⁺ signal out of N signals measured simultaneously in a collective of pancreatic beta cells.

Random matrix theory (Guhr et al., 1998; Mehta, 2004) (RMT) is concerned with statistical properties of matrices with random elements. Applying RMT to correlation matrices, we study the spectrum of the correlation matrix C given by the set of its eigenvalues λ_n:

$\begin{array}{ll}\text{C}{\text{u}}_n={\lambda }_n{\text{u}}_n,& (2)\end{array}$

where u_n are the corresponding eigenvectors.

Statistical properties of the spectra of random correlation matrices for N uncorrelated time series with M random elements where q = N/M is finite in the limit N, M → ∞ are known analytically (Marchenko and Pastur, 1967; Bun et al., 2017). The eigenvalue probability density is:

$\begin{array}{ll}\rho (\lambda )=\frac{1}{2\pi q\lambda }\sqrt{({\lambda }_{+}-\lambda )(\lambda -{\lambda }_{-})},& (3)\end{array}$

where the spectral boundaries are:

$\begin{array}{ll}{\lambda }_{\pm }={\left(\sqrt{q}\pm 1\right)}^2& (4)\end{array}$

When the spectrum of the correlation matrix is unfolded (Guhr et al., 1998) by mapping eigenvalues λ_k → ξ_k so that the probability density of the unfolded eigenvalues is constant ρ(ξ) = 1, the RMT predicts that the distribution P(s) of nearest neighbor spacings s_k = ξ_k+1 − ξ_k is approximately given by the Wigner surmise (Mehta, 2004):

$\begin{array}{ll}P(s)=\frac{\pi }{2}sexp\left(-\frac{\pi }{4}{s}^2\right).& (5)\end{array}$

Possible pair correlations in the eigenvalue spectrum on scales larger than nearest neighbors can be revealed with the use of variance of n_ξ(L), the number of eigenvalues in the interval of length L around eigenvalue ξ. This number variance (Mehta, 2004) is given by:

$\begin{array}{ll}{\Sigma }^2(L)=\langle {(n_{\xi }(L)-L)}^2\rangle .& (6)\end{array}$

If the eigenvalue spectrum is poissonian the number variance is Σ²(L) ~ L, while real, symmetric random matrices exhibit correlated spectra for which RMT predicts Σ²(L) ~ logL (Mehta, 2004).

Previous work using RMT in different systems, e.g., on statistics of energy levels of complex quantum systems (Guhr et al., 1998; Mehta, 2004) or correlations in financial markets (Plerou et al., 2002) identified that a bulk of the eigenvalue spectrum agreed with RMT predictions, which suggested a large degree of randomness in the measured cross-correlations in these systems. Only a small fraction of typically a few percent of eigenvalues were found to deviate from universal RMT predictions and were instrumental to identify system specific, non-random properties of the observed systems and yielding key information about the underlying interactions. Biological systems are often complex with large number of interacting parts, high dimensional systems that are basically black boxes “in which a large number of particles are interacting according to unknown laws” (Dyson, 1962). One way to approach such high dimensional systems is to look at the spectrum of the covariance matrix and try to find few principal components that describer most of system variance. This method, principal component analysis (PCA), has been suggested as “‘hypothesis generating’ tool creating a statistical mechanics frame for biological systems modeling” (Giuliani, 2017). PCA works best for systems where we can find few eigenvalues in the covariance spectrum well separated from the bulk and the system can be described in low dimensional space. Usually, there is no clear separation in the eigenvalue spectrum and other methods such as RMT or methods using renormalization group approach (Bradde and Bialek, 2017) are more suitable. In biological systems, RMT has been used to filter out critical correlations from data-rich samples in genomic, proteomic and metabolomic analyses (Luo et al., 2006; Agrawal et al., 2014), as well as in brain activity measured by EEG (Šeba, 2003) and dynamic brain networks constructed from fMRI data (Wang et al., 2015, 2016). While eigenvalue spacing distributions showed agreement with RMT predictions, the number variance distributions often display deviations pointing to the physiologically relevant reduction in correlated eigenvalues fluctuations with partially decoupled components transiting toward Poisson distribution (Šeba, 2003). Such transitions have also been used as an objective approach for the identification of functional modules within large complex biological networks (Luo et al., 2006). Additionally, as for protein-protein interactions in different species, these latter deviation from RMT predictions has been interpreted as an evidence to support the prevalence of non-random mutations in biological systems (Agrawal et al., 2014).

In this paper we used the RMT approach to test the cross-correlations in the cytosolic Ca²⁺ oscillations under non-stimulated and glucose stimulated conditions. We demonstrate that statistical properties of cross-correlations based on functional multicellular imaging data follows those predicted by RMT, with both high- and low-end deviating eigenvalues, suggesting local as well as global modes driving this correlation in functional islet. In addition, our results show that the long range correlations in eigenvalue spectrum deviate in a stimulus dependent manner.

Dataset Description

We define beta cell collective activity to sense nutrients and produce metabolic code as the relevant constraining context for the physical outcomes of analysis (Ellis and Kopel, 2018; Korošak and Slak Rupnik, 2018). Our data consist of a typical Ca²⁺ activity recorded by multicellular imaging on an islet of Langerhans of fresh mouse pancreatic slice. All data analysis has been performed using custom scripts in Python 3.5 software and customized scripts (RMThreshold) in R software. We used raw data for each calcium signal, but we detrended the signals to remove possible sources of spurious correlations due to systematic slow variations caused by the imaging technique. A common problem in the analysis of fMCI Ca²⁺ signals in living tissue is selection of regions of interest corresponding to a true signal originating from a cytosolic area of an individual cell and not two or more neighboring cells. In practice the reproducibility of the results depends on the level of experience of the operator to subjectively recognize structure from the patterns of activity. While we are primarily interested in the activity of a large population of cells, their interactions/correlations and their collective response, it is crucial that this signals originate from regions of interest that correspond to individual cells. Collectives of cells, like beta cells in the islet of Langerhans are densely packed structures, where extracellular space and the cell membrane represent a relatively small to negligible cross-section area on the image of two-dimensional optical section obtained by confocal microscopy. Therefore we decided to avoid the aforementioned subjectivity issue by the random sampling of pixel level signals in the recorded time series. For this analysis we randomly selected N = 4000 signals out of the complete dataset of 256 × 256 signals each M = 23,820 timesteps long (about 40 min recordings at 10 Hz resolution). Glucose concentration was changed during the recording from 6 mM (lasting about first 5,000 timesteps) to 8 mM and back to 6 mM (approx. last 5,000 timesteps) near the end of experiment.

The source of correlations in a cell population where the terminal action is a calcium-dependent process (e.g., exocytosis of insulin in beta cells) are the individual events in a form of plasma membrane ion channel or transporter activity, internal membrane ion channel or transporter activity, as well as calcium leak from activated immediate neighboring beta cell (Berridge et al., 2000). The correlations between the activities of beta cells depend strongly upon the glucose concentration (Dolenšek et al., 2013; Markovič et al., 2015), however in the physiological plasma glucose range (6–9 mM), most correlations are weak (Korošak and Slak Rupnik, 2018), so that the probability of detecting co-activation basically equals the product of the probabilities of activities of individual beta cells. The correlations are statistically significant for almost all pairs of immediate neighbors.

Results

The distribution of correlation coefficients reveals that most of the correlations between the pairs of Ca²⁺ signals are weak, but there is also non-negligible contribution of highly correlated pairs of signals (Figure 1, left, black outer line). We also checked the sampling procedure by comparing the computed distribution of distances between pairs of randomly chosen points from 256 × 256 image square to the analytical probability distribution of distances between two random points in a square (Philip, 2007) (Figure 1, right). We found a perfect match between the distance distribution computed from data and the theoretical distribution, confirming that our random sampling of data points was non-biased.

FIGURE 1

Figure 1. (Left) Correlation coefficient distribution of N = 4000 Ca²⁺ signals randomly chosen from experimental dataset (outer, black line). Dashed blue line is the distribution of correlations between residuals after removing the influence of the largest eigenvalue in signals. (Right) Distribution of distances between signals randomly chosen from Ca²⁺ imaging data. Black line is the experimental data, thinner red line is the theoretical distribution of distances between random points in a square (Philip, 2007).

Guided by the observed non-gaussian nature of correlation distribution (Figure 1, left) we explored a detailed structure of the correlation matrix, since distribution of correlation coefficients only partially hints to the nature of cell to cell coordination. To this end we computed the eigenvalues and eigenvectors of the correlation matrix (Equation 2) and compared the obtained eigenspectrum with the RMT prediction. In Figure 2 (top left) we show the distribution of eigenvalues that belong to the empirical correlation matrix (black trace) and the RMT prediction (red line) given by Equation (3). While most of the eigenvalues falls within the limits λ_± of the RMT spectrum, there are also significant deviations from RMT prediction. We found the largest empirical eigenvalue λ_max two orders of magnitude away from the upper limit of the RMT spectrum, and also a part of the empirical spectrum that extends below the lower RMT limit. To see if the deviations from the RMT are inherent to the measured Ca²⁺ signals, we prepared a surrogate dataset by randomly shuffling each signal's time series. We then computed the correlation matrix and its eigenvalue spectrum from randomized surrogate dataset. As shown in Figure 2 (top right), the match between the eigenvalue distribution of randomized dataset and RMT is perfect.

FIGURE 2

Figure 2. (Top left) Probability distribution of eigenvalues of the empirical correlation matrix for N = 4,000 randomly picked signals (black solid line) compared to the distribution of eigenvalues of random correlation matrix of the same size (red solid line). (Top right) Probability distribution of eigenvalues of the surrogate correlation matrix constructed from shuffled empirical values (black solid line) compared to the random correlation matrix of the same size (red solid line). (Bottom left) Inverse participation ratio plot of eigenvalues showing a random band matrix structure of C with large IPR values at both edges of the eigenvalue spectrum. The dashed vertical lines show RMT bounds. The IPR spectrum for randomized correlations is shown in red. (Bottom right) The fractional rank plot of the entire spectrum of eigenvalues (open dots). For comparison we added the same plot of eigenvalues of correlation matrix computed from randomized data (open squares). The full line shows the fractional rank plot of eigenvalue spectra obtained from distribution given by Equation (3). The shape of the distribution of large eigenvalues points to a scaling relationship.

Previous RMT analysis of stock correlations in financial markets consistently showed (Laloux et al., 1999; Plerou et al., 1999, 2002) that the distribution of components of the eigenvector u_max corresponding to largest eigenvalue λ_max strongly deviates from Gaussian form, suggesting that this mode reflects the collective response of the system to the stimuli. In our case this corresponds to collective response of beta-cells to glucose stimulus. In a linear statistical model for Ca²⁺ signals, we model the response common to all beta-cell with Y(t) and the signals are expressed as:

$\begin{array}{ll}{y}_i(t)=a_i+b_{i}Y(t)+\delta y_i(t),& (7)\end{array}$

where δy_i(t) is the residual part of each signal. Coefficients a_i, b_i are obtained by regression. Following Plerou et al. (2002) we approximated the common response Y(t) with the projection of all signals on the largest eigenvector:

$\begin{array}{ll}{y}_{max}(t)={\displaystyle \sum _{i=1}^N}{u}_i({\lambda }_{max})y_i(t),& (8)\end{array}$

where u_{max, i} is the i-th component of the eigenvector corresponding to largest eigenvalue λ_max. To see the influence of the collective response to the distribution of correlation coefficents, we computed using Y = y_max the residuals δy_i(t) for all N signals and their correlations C_{res(i, j)} = Corr(δy_i, δy_j). The dashed blue line (inner trace) on Figure 1 (left) shown the distribution of C_res and reveals that the collective response predominantly contributes to large correlations.

To test further if the largest eigenvalue and the corresponding eigenvector capture the collective calcium response we compared the average signal $\overline{y}(t)=1/N\sum _j{y}_j(t)$ with y_max. The correlation between signals projected on the largest eigenvalue mode and mean signal was high: $Corr(y_{max},\overline{y})\approx 0.8$, confirming the expectation that the largest eigenvalue represents collective effect. Similarly, we checked how similar are the signals corresponding to the bulk RMT eigenvalues:

$\begin{array}{ll}{y}_{bulk,i}(t)={\displaystyle \sum _{j=1}^N}{u}_j({\lambda }_i)y_j(t),& (9)\end{array}$

where λ_i is the eigenvalues from the RMT interval [λ₊, λ₋]. The computed correlation between signals projected on bulk eigenvectors and mean signal, averaged over all signals was $<Corr(y_{bulk,i},\overline{y})>=-0.0044\pm 0.0047$ suggesting no correlation between the mean signal and signals coming from the bulk RMT regime. To further characterize the eigenvector structure of the empirical Ca²⁺ correlation matrix, we looked at the inverse participation ratio (IPR) of eigenvector u(λ) corresponding to eigenvalue λ defined as (Plerou et al., 1999, 2002):

$\begin{array}{ll}I(\lambda )={\displaystyle \sum _j^N}{(u_j(\lambda ))}^4.& (10)\end{array}$

The value of 1/I(λ) reflects the number of nonzero eigenvector components: if an eigenvector consist of equal components $u{(\lambda )}_i=1/\sqrt{(N)}$ then 1/I(λ) = N, in other extreme case 1/I(λ) = 1 when an eigenvector has one component equal to 1 and all others are zero.

Figure 2 (lower left) shows the computed values of IPR for all eigenvectors as function of corresponding eigenvalue. The red datapoints are the IPR data computed for the surrogate, randomized timeseries data for which we found 1/I ~ N as expected. We found similarly values 1/I ~ N for the largest eigenvalues of the empirical spectrum (black datapoints) suggesting that to this eigenvectors almost all signals contribute. Deviations from flat RMT prediction at the edges of the RMT spectrum ([λ₊, λ₋] interval, vertical dashed lines) with large I(λ) values suggests that these states are localized with only a few signals contributing. This points to a complex structure of the empirical correlation matrix C with coexisting extended and localized eigenvectors similar to one found in correlations in financial markets (Plerou et al., 1999, 2002). In addition, as shown in Figure 2 (lower right, open dots), we observe a scaling behavior in rank-ordered plot of eigenvalues of empirical correlation matrix that has been connected with a fixed point in renormalization group sense (Bradde and Bialek, 2017; Meshulam et al., 2018). For comparison, we plot also the rank-ordered eigenvalues of randomized data (open squares) and RMT prediction based on eigenvalue density given by Equation (3) (full line) which perfectly describes the randomized dataset. The observed scaling of eigenvalues hints toward the critical behavior that was conjectured for beta-cell collective at the transition from glucose non-stimulated to stimulated phase (from 6 mM to 8 mM) (Gosak et al., 2017b; Stožer et al., 2019).

To explore the statistical differences of signals in non-stimulated and stimulated phase, we separated the original data into two groups of N signals each with M = 10⁴ timesteps corresponding to response to 6 and 8 mM glucose stimuli. For each group we computed the unfolded eigenvalue spectra and also for randomized data. The results for the nearest-neighbor spacing and number variances are shown in Figure 3. For nearest-neighbor spacing distribution we find a good agreement with the RMT prediction both, for non-stimulatory and stimulatory conditions, as well as shuffled stimulated data. All three datasets are well described with the Wigner surmise (Equation 5), so nearest-neighbor spacing does not seem to be sensitive to stimuli changes. On the other hand, however, the number variance is sensitive to stimuli change already during physiological stimulation of the beta-cell collective. The random matrix prediction is in this case valid for shuffled stimulated data only (Figure 3, right).

FIGURE 3

Figure 3. (Left) Nearest-neighbor spacing distribution of empirical correlation matrix eigenvalues of calcium signals in non-stimulatory and stimulatory regime. Open squares: shuffled, randomized data; open dots: 8 mM glucose; crosses 6 mM glucose; full line Wigner surmise (Equation 5). (Right) Number variance of eigenvalue spectra of calcium signals. Open squares: shuffled, randomized data; open dots: 8 mM glucose; crosses 6 mM glucose; full line RMT prediction Σ²(L) = 1/π²(log(2πL) + 1 + γ – π²/8) (Mehta, 2004); dashed line Poissonian limit Σ²(L) = L.

Discussion

The unique spatio-temporal resolution of functional multicellular imaging and sensitivity of advanced statistical approaches for a plethora of different modes of complex network scenarios and levels of criticality, makes these approaches a method of choice to assess the nature of cell-cell interactions under different stimulation conditions. At the same time it enables us to test the validity of experimental designs for study of beta cell function, primarily in the domain of stimulation strength and dynamics. We suggest that without such validation the most critical events in the activation chain within the beta cell collective have been and shall be overlooked or misinterpreted (Stožer et al., 2019). The predominant use of supraphysiological glucose concentration can namely severely deform the relatively slow beta cell recruitment in a collective at physiological glucose concentrations (Stožer et al., 2013a; Gosak et al., 2017a) and miss the typically segregated network clusters of Ca²⁺ events (Benninger and Piston, 2014; Markovič et al., 2015; Westacott et al., 2017), turning critical behavior into disruptive supracritical activity (Gosak et al., 2017b; Stožer et al., 2019). Only under rather narrow physiological conditions it shall be possible to extract the fine structure of cell-cell interactions causing long-term and efficient cell collaboration with the collective. Breaking apart this delicate structure of cell-cell interaction does result in a massive activity, which can be readily described by tools of statistical physics, but this activity does not necessarily serve its physiological or biological purpose (Ellis and Kopel, 2018).

A common denominator of the previous attempts to categorize different beta cell types points to some metabolic and secretory features that can be either reproduced between different classifications or not. Usually there exist a bulk of one subtype and one or more less frequent subtypes (Benninger and Hodson, 2018). These less frequent subtypes can nevertheless have important regulatory roles that may not be immediately apparent. This issue is particularly critical if the frequency of a beta cell subtype represents only a couple of percent of the entire beta cell population in an islet. Along these lines there have been some indications regarding the beta cell subtypes that can serve as pacemakers or hubs within a dynamic islet cell network (Johnston et al., 2016; Lei et al., 2018), however due to the nature of complexity of network features, we may still be short of evidence for definitive conclusions. The full description of heterogeneity of endocrine cells within an islet, ultimately producing a adequate release of hormones is therefore still lacking. In trying to grasp this complexity, it is important to take into account interaction of beta cell collectives with other cell types in and around an islet, like glucagon-secreting alpha cells (Svendsen et al., 2018; Capozzi et al., 2019) or somatostatin secreting delta cells (Rorsman and Huising, 2018) as well as neurons and glial cells (Meneghel-Rozzo et al., 2004), but also endothelial, immune cells (Damond et al., 2019), as well as acinar and ductal cells (Bertelli and Bendayan, 2005).

Random matrix theory is a fitting mathematical framework which provides powerful analytical tools to separate cell-cell interactions happening by chance from those produced by specific coordinated interactions after a changed chemical composition of cell's surrounding. In the financial sector, adequate asset allocation and portfolio risk-estimation can lead to a higher profit and is therefore clear why it makes sense to invest time into cross-correlation analyses (Plerou et al., 2002). But what would be the gain of knowing that randomness of cell-cell correlation matrices is physiologically regulated? Firstly, we suggest that the analysis of the universal properties of empirical cross-correlations is a valuable tool to identify distinct types and further subtypes of endocrine cells within an islet through their non-local and local effects. The largest eigenvalue of C namely represents the influence of non-local modes common to all measured Ca²⁺ fluctuations. Other large eigenvalues can be used to address cross-correlations between cells of the same type, cells with specific functions in the collective or that these cells reside in topologically similar area of the islet. Quantifying correlations between different beta cells in an islet is therefore an exciting scientific effort that can help us understand cell communities as a complex dynamical system, estimate the amount of factors ruling the system or potential presence of a stress situation (Gorban et al., 2010).

The large values of inverse participation ratio (IPR) (Figure 2, bottom left) compared to the IPR values in the bulk, indicate that only a few cells contribute to these eigenstates with eigenvalues at the edges of the RMT bulk spectrum. In contrast, all cells contribute to the eigenstates corresponding to the largest eigenvalues. This means that we find delocalized states for the largest eigenvalues and localized states as we move toward the RMT edge of the spectrum. Similar findings were recently reported in RMT analysis of single-cell sequencing data (Aparicio et al., 2018), where the spectrum of covariance matrix of single-cell genomic data followed RMT predictions with deviations at the bulk edge. The localized states at the edge of the bulk spectrum were connected with the true biological signal.

Our results show that the number variance reflecting the correlation between subsequent eigenvalues (a measure for long range correlations in eigenvalue spectrum) follows the RMT predictions up to a certain distance L, however at larger distances it starts to deviate in a stimulus dependent manner, suggesting structural features in the beta cell network. Transitions between Poissonian and GOE statistics in biological systems have been previously described during the process of either integration or segregation of complex biological networks, showing various degrees of long range correlations at various physiological conditions (Luo et al., 2006). This understanding has a vital practical value since it can help decipher different roles that beta cells can play in a collective and to further validate the importance, if any, of previously defined and continuously appearing novel molecular markers of beta cell heterogeneity (Benninger and Hodson, 2018; Damond et al., 2019; Wang et al., 2019). An advanced knowledge about the dynamic properties of the functional cell types will shed a new light into understanding of physiological regulation of insulin release and the assessment of perils of stimulation outside of the physiological range. Furthermore, it can help us elucidate the mechanisms on how this function changes during the pathogenesis of different forms of diabetes mellitus and lead us to novel approaches of therapy planning and prevention. And finally, it can help us understand the general principles ruling the interactions in collectives of other cell types.

Data Availability Statement

The raw data used in this research is available upon request to the corresponding author.

Author Contributions

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

Conflict of Interest

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 DK and MS, acknowledge the financial support from the Austrian Science Fund/Fonds zur Förderung der Wissenschaftlichen Forschung [Bilateral grants I3562-B27 and I4319-B30 (to MS)] and the Slovenian Research Agency (research core funding, No. P3-0396), as well as research project, No. N3-0048).

References

Agrawal, A., Sarkar, C., Dwivedi, S. K., Dhasmana, N., and Jalan, S. (2014). Quantifying randomness in protein–protein interaction networks of different species: a random matrix approach. Physica A 404, 359–367. doi: 10.1016/j.physa.2013.12.005

CrossRef Full Text | Google Scholar

American Diabetes Association (2014). Diagnosis and classification of diabetes mellitus. Diabetes Care 37(Suppl. 1):S81–S90. doi: 10.2337/dc14-S081

CrossRef Full Text

Aparicio, L., Bordyuh, M., Blumberg, A. J., and Rabadan, R. (2018). Quasi-universality in single-cell sequencing data. arXiv:1810.03602. doi: 10.1101/426239

CrossRef Full Text | Google Scholar

Bak, P. (2013). How Nature Works: The Science of Self-Organized Criticality. New York, NY: Springer Science & Business Media.

Google Scholar

Barabási, A.-L., and Márton Pósfai. (2016). Network Science. Cambridge: Cambridge University Press.

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. Diab. Obes. Metab. 9, 118–132. doi: 10.1111/j.1463-1326.2007.00780.x

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

Bertelli, E., and Bendayan, M. (2005). Association between endocrine pancreas and ductal system. more than an epiphenomenon of endocrine differentiation and development? J. Histochem. Cytochem. 53, 1071–1086. doi: 10.1369/jhc.5R6640.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Bertram, R., Satin, L. S., and Sherman, A. S. (2018). Closing in on the mechanisms of pulsatile insulin secretion. Diabetes 67, 351–359. doi: 10.2337/dbi17-0004

PubMed Abstract | CrossRef Full Text | Google Scholar

Bradde, S., and Bialek, W. (2017). PCA meets RG. J. Stat. Phys. 167, 462–475. doi: 10.1007/s10955-017-1770-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Briant, L. J., Zhang, Q., Vergari, E., Kellard, J. A., Rodriguez, B., Ashcroft, F. M., et al. (2017). Functional identification of islet cell types by electrophysiological fingerprinting. J. R. Soc. Interf. 14:20160999. doi: 10.1098/rsif.2016.0999

PubMed Abstract | CrossRef Full Text | Google Scholar

Bun, J., Bouchaud, J.-P., and Potters, M. (2017). Cleaning large correlation matrices: tools from random matrix theory. Phys. Rep. 666, 1–109. doi: 10.1016/j.physrep.2016.10.005

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

Capozzi, M. E., Svendsen, B., Encisco, S. E., Lewandowski, S. L., Martin, M. D., Lin, H., et al. (2019). β-cell tone is defined by proglucagon peptides through cyclic amp signaling. JCI Insight. 7:126742. doi: 10.1172/jci.insight.126742

CrossRef Full Text

Damond, N., Engler, S., Zanotelli, V. R., Schapiro, D., Wasserfall, C. H., Kusmartseva, I., et al. (2019). A map of human type 1 diabetes progression by imaging mass cytometry. Cell Metab. 29, 755–768.e5. doi: 10.1016/j.cmet.2018.11.014

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

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

Duh, A., Slak Rupnik, M., and Korošak, D. (2018). Collective behavior of social bots is encoded in their temporal twitter activity. Big Data 6, 113–123. doi: 10.1089/big.2017.0041

PubMed Abstract | CrossRef Full Text | Google Scholar

Dyachok, O., Isakov, Y., Sågetorp, J., and Tengholm, A. (2006). Oscillations of cyclic amp in hormone-stimulated insulin-secreting β-cells. Nature 439:349. doi: 10.1038/nature04410

PubMed Abstract | CrossRef Full Text | Google Scholar

Dyson, F. J. (1962). Statistical theory of the energy levels of complex systems. I. J. Math. Phys. 3, 140–156. doi: 10.1063/1.1703773

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Frank, J. A., Broichhagen, J., Yushchenko, D. A., Trauner, D., Schultz, C., and Hodson, D. J. (2018). Optical tools for understanding the complexity of β-cell signalling and insulin release. Nat. Rev. Endocrinol. 14, 721–737. doi: 10.1038/s41574-018-0105-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Giuliani, A. (2017). The application of principal component analysis to drug discovery and biomedical data. Drug Discov. Today 22, 1069–1076. doi: 10.1016/j.drudis.2017.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Gorban, A. N., Smirnova, E. V., and Tyukina, T. A. (2010). Correlations, risk and crisis: from physiology to finance. Physica A 389, 3193–3217. doi: 10.1016/j.physa.2010.03.035

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. 24, 118–135. doi: 10.1016/j.plrev.2017.11.003

CrossRef Full Text | Google Scholar

Gosak, M., Stozer, A., Markovic, R., Dolensek, 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

Guhr, T., Müller-Groeling, A., and Weidenmüller, H. A. (1998). Random-matrix theories in quantum physics: common concepts. Phys. Rep. 299, 189–425. doi: 10.1016/S0370-1573(97)00088-4

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

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Korošak, D., and Slak Rupnik, M. (2018). Collective sensing of β-cells generates the metabolic code. Front. Physiol. 9:31. doi: 10.3389/fphys.2018.00031

PubMed Abstract | CrossRef Full Text | Google Scholar

Laloux, L., Cizeau, P., Bouchaud, J.-P., and Potters, M. (1999). Noise dressing of financial correlation matrices. Phys. Rev. Lett. 83:1467. doi: 10.1103/PhysRevLett.83.1467

CrossRef Full Text | Google Scholar

Lei, C.-L., Kellard, J. A., Hara, M., Johnson, J. D., Rodriguez, B., and Briant, L. J. (2018). Beta-cell hubs maintain ca2+ oscillations in human and mouse islet simulations. Islets 10, 151–167. doi: 10.1080/19382014.2018.1493316

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, F., Zhong, J., Yang, Y., Scheuermann, R. H., and Zhou, J. (2006). Application of random matrix theory to biological networks. Phys. Lett. A 357, 420–423. doi: 10.1016/j.physleta.2006.04.076

CrossRef Full Text | Google Scholar

Marchenko, V. A., and Pastur, L. A. (1967). Distribution of eigenvalues for some sets of random matrices. Matematicheskii Sbornik 114, 507–536.

Google Scholar

Marciniak, A., Cohrs, C. M., Tsata, V., Chouinard, J. A., Selck, C., Stertmann, J., et al. (2014). Using pancreas tissue slices for in situ studies of islet of langerhans and acinar cell biology. Nat. Protoc. 9:2809. doi: 10.1038/nprot.2014.195

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Mehta, M. L. (2004). Random Matrices, Vol. 142. Amsterdam: Elsevier.

Google Scholar

Meneghel-Rozzo, T., Rozzo, A., Poppi, L., and Rupnik, M. (2004). In vivo and in vitro development of mouse pancreatic β-cells in organotypic slices. Cell Tissue Res. 316, 295–303. doi: 10.1007/s00441-004-0886-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Meshulam, L., Gauthier, J. L., Brody, C. D., Tank, D. W., and Bialek, W. (2018). Coarse–graining and hints of scaling in a population of 1000+ neurons. arXiv:1812.11904v1.

Google Scholar

Milo, R., Shen-Orr, S., Itzkovitz, S., Kashtan, N., Chklovskii, D., and Alon, U. (2002). Network motifs: simple building blocks of complex networks. Science 298, 824–827. doi: 10.1126/science.298.5594.824

PubMed Abstract | CrossRef Full Text | Google Scholar

Misler, S., Falke, L. C., Gillis, K., and McDaniel, M. L. (1986). A metabolite-regulated potassium channel in rat pancreatic b cells. Proc. Natl. Acad. Sci. U.S.A. 83, 7119–7123. doi: 10.1073/pnas.83.18.7119

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Perc, M., Jordan, J. J., Rand, D. G., Wang, Z., Boccaletti, S., and Szolnoki, A. (2017). Statistical physics of human cooperation. Phys. Rep. 687, 1–51. doi: 10.1016/j.physrep.2017.05.004

CrossRef Full Text | Google Scholar

Philip, J. (2007). The Probability Distribution of the Distance Between Two Random Points in a Box. Stockholm: KTH Mathematics, Royal Institute of Technology.

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Plerou, V., Gopikrishnan, P., Rosenow, B., Amaral, L. A. N., Guhr, T., and Stanley, H. E. (2002). Random matrix approach to cross correlations in financial data. Phys. Rev. E 65:066126. doi: 10.1103/PhysRevE.65.066126

PubMed Abstract | CrossRef Full Text | Google Scholar

Plerou, V., Gopikrishnan, P., Rosenow, B., Amaral, L. A. N., and Stanley, H. E. (1999). Universal and nonuniversal properties of cross correlations in financial time series. Phys. Rev. Lett. 83:1471. doi: 10.1103/PhysRevLett.83.1471

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Rorsman, P., and Huising, M. O. (2018). The somatostatin-secreting pancreatic δ-cell in health and disease. Nat. Rev. Endocrinol. 14, 404–414. doi: 10.1038/s41574-018-0020-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Rorsman, P., and Trube, G. (1986). Calcium and delayed potassium currents in mouse pancreatic beta-cells under voltage-clamp conditions. J. Physiol. 374, 531–550. doi: 10.1113/jphysiol.1986.sp016096

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

Šeba, P. (2003). Random matrix analysis of human EEG data. Phys. Rev. Lett. 91:198104. doi: 10.1103/PhysRevLett.91.198104

PubMed Abstract | CrossRef Full Text | Google Scholar

Sherman, A., RiNZEL, J., and Keizer, J. (1988). Emergence of organized bursting in clusters of pancreatic beta-cells by channel sharing. Biophys. J. 54, 411–425. doi: 10.1016/S0006-3495(88)82975-8

PubMed Abstract | CrossRef Full Text | Google Scholar

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

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

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

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Svendsen, B., Larsen, O., Gabe, M. B. N., Christiansen, C. B., Rosenkilde, M. M., Drucker, D. J., et al. (2018). Insulin secretion depends on intra-islet glucagon signaling. Cell Rep. 25, 1127–1134. doi: 10.1016/j.celrep.2018.10.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, R., Wang, L., Yang, Y., Li, J., Wu, Y., and Lin, P. (2016). Random matrix theory for analyzing the brain functional network in attention deficit hyperactivity disorder. Phys. Rev. E 94:052411. doi: 10.1103/PhysRevE.94.052411

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, R., Zhang, Z.-Z., Ma, J., Yang, Y., Lin, P., and Wu, Y. (2015). Spectral properties of the temporal evolution of brain network structure. Chaos 25:123112. doi: 10.1063/1.4937451

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y. J., Traum, D., Schug, J., Gao, L., Liu, C., Atkinson, M. A., et al. (2019). Multiplexed in situ imaging mass cytometry analysis of the human endocrine pancreas and immune system in type 1 diabetes. Cell Metab. 29, 769–783.e4. doi: 10.1016/j.cmet.2019.01.003

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, M., Goforth, P., Bertram, R., Sherman, A., and Satin, L. (2003). The ca2+ dynamics of isolated mouse β-cells and islets: implications for mathematical models. Biophys. J. 84, 2852–2870. doi: 10.1016/S0006-3495(03)70014-9

PubMed Abstract | CrossRef Full Text | Google Scholar