Skip to main content
Advertisement
  • Loading metrics

Network representation of multicellular activity in pancreatic islets: Technical considerations for functional connectivity analysis

  • Marko Šterk,

    Roles Conceptualization, Formal analysis, Investigation, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing

    Affiliations Faculty of Natural Sciences and Mathematics, University of Maribor, Maribor, Slovenia, Faculty of Medicine, University of Maribor, Maribor, Slovenia

  • Yaowen Zhang,

    Roles Data curation, Investigation

    Affiliation Department of Pediatrics, Department of Bioengineering, University of California, San Diego, La Jolla, California, United States of America

  • Viljem Pohorec,

    Roles Data curation, Investigation

    Affiliation Faculty of Medicine, University of Maribor, Maribor, Slovenia

  • Eva Paradiž Leitgeb,

    Roles Data curation, Investigation, Writing – review & editing

    Affiliation Faculty of Medicine, University of Maribor, Maribor, Slovenia

  • Jurij Dolenšek,

    Roles Investigation, Writing – review & editing

    Affiliations Faculty of Natural Sciences and Mathematics, University of Maribor, Maribor, Slovenia, Faculty of Medicine, University of Maribor, Maribor, Slovenia

  • Richard K. P. Benninger,

    Roles Investigation, Resources, Writing – review & editing

    Affiliations Department of Bioengineering, Barbara Davis Center for Diabetes, Aurora, Colorado, United States of America, Barbara Davis Center for Childhood Diabetes, University of Colorado Anschutz Medical Campus, Aurora, Colorado, United States of America

  • Andraž Stožer ,

    Roles Conceptualization, Funding acquisition, Investigation, Validation, Writing – review & editing

    andraz.stozer@um.si (AS); vkravets@ucsd.edu (VK); marko.gosak@um.si (MG)

    Affiliation Faculty of Medicine, University of Maribor, Maribor, Slovenia

  • Vira Kravets ,

    Roles Conceptualization, Funding acquisition, Investigation, Resources, Validation, Writing – review & editing

    andraz.stozer@um.si (AS); vkravets@ucsd.edu (VK); marko.gosak@um.si (MG)

    Affiliations Department of Pediatrics, Department of Bioengineering, University of California, San Diego, La Jolla, California, United States of America, Department of Bioengineering, Jacobs School of Engineering, University of California, San Diego, La Jolla, California, United States of America

  • Marko Gosak

    Roles Conceptualization, Funding acquisition, Investigation, Methodology, Software, Supervision, Validation, Writing – original draft

    andraz.stozer@um.si (AS); vkravets@ucsd.edu (VK); marko.gosak@um.si (MG)

    Affiliations Faculty of Natural Sciences and Mathematics, University of Maribor, Maribor, Slovenia, Faculty of Medicine, University of Maribor, Maribor, Slovenia, Alma Mater Europaea, Maribor

Abstract

Within the islets of Langerhans, beta cells orchestrate synchronized insulin secretion, a pivotal aspect of metabolic homeostasis. Despite the inherent heterogeneity and multimodal activity of individual cells, intercellular coupling acts as a homogenizing force, enabling coordinated responses through the propagation of intercellular waves. Disruptions in this coordination are implicated in irregular insulin secretion, a hallmark of diabetes. Recently, innovative approaches, such as integrating multicellular calcium imaging with network analysis, have emerged for a quantitative assessment of the cellular activity in islets. However, different groups use distinct experimental preparations, microscopic techniques, apply different methods to process the measured signals and use various methods to derive functional connectivity patterns. This makes comparisons between findings and their integration into a bigger picture difficult and has led to disputes in functional connectivity interpretations. To address these issues, we present here a systematic analysis of how different approaches influence the network representation of islet activity. Our findings show that the choice of methods used to construct networks is not crucial, although care is needed when combining data from different islets. Conversely, the conclusions drawn from network analysis can be heavily affected by the pre-processing of the time series, the type of the oscillatory component in the signals, and by the experimental preparation. Our tutorial-like investigation aims to resolve interpretational issues, reconcile conflicting views, advance functional implications, and encourage researchers to adopt connectivity analysis. As we conclude, we outline challenges for future research, emphasizing the broader applicability of our conclusions to other tissues exhibiting complex multicellular dynamics.

Author summary

Islets of Langerhans, multicellular microorgans in the pancreas, are pivotal for whole-body energy homeostasis. Hundreds of beta cells within these networks synchronize to produce insulin, a crucial hormone for metabolic control. Coordinated activity disruptions in these multicellular networks contribute to irregular insulin secretion, a hallmark of diabetes. Recognizing the significance of collective activity, network science approaches have been increasingly applied in islet research. However, variations in experimental setups, imaging techniques, signal processing, and connectivity analysis methods across different research groups pose challenges for integrating findings into a comprehensive picture. Therefore, we present here a systematic analysis of various approaches impacting results in islet activity network representation. We find that methods for constructing functional connectivity maps aren’t critical, but caution is necessary when aggregating data from different islets. Network analysis conclusions are notably influenced by factors such as time series pre-processing, the oscillatory component of signals, and experimental preparation. Despite these challenges, this paper advocates for the adoption of connectivity analysis in future islet research, emphasizing that the insights gained extend beyond pancreatic islets to provide valuable contributions for understanding connectivity in other multicellular systems.

Introduction

Proper insulin secretion and insulin sensitivity of peripheral tissues are crucial in regulating uptake and disposal of energy rich molecules, thereby sustaining metabolic homeostasis [1]. Pancreatic beta cells constitute part of a crucial negative feedback loop, sensing changes in plasma levels of energy-rich nutrients and accordingly adjusting release of insulin into the bloodstream [2]. The cascade of cellular events connecting the changes in plasma nutrient levels with the proper insulin secretion have been studied in detail [310]. The crucial steps in the stimulus-secretion coupling cascade involve an increase in intracellular ATP concentration, closure of ATP-sensitive potassium channels, membrane depolarization, opening of voltage-activated Ca2+ channels and an increase in intracellular calcium concentration ([Ca2+]i), leading ultimately to exocytosis of insulin-containing vesicles. Beta cell response is further modulated by homo- and heterologous cell-to-cell interactions within islets [1114], by autonomous nerve control [1517], and by hormones released by the gut [1820]. Of importance, beta cells display complex oscillatory activity and are intrinsically heterogeneous [8,21,22], with differences observed on molecular [23], morphological [24,25] and functional levels [26,27], and it is only due to strong coupling within the islets that beta cells properly respond to glucose excursions.

In a coupled system of beta cells, glucose stimulation triggers two distinct and qualitatively different phases [8,2732]. The initial response consists of a phasic increase in activity, characterized by membrane depolarization and increase in [Ca2+]i which occurs sooner in higher glucose concentrations [27,33]. In case that the stimulus is still present, a complex tonic activity follows. This second phase is characterized by repetitive membrane potential and [Ca2+]i oscillations, as well as pulses of insulin secretion. These oscillations are not generated randomly among cells. Rather, they are phase-lagged between cells, such that waves of membrane depolarization and [Ca2+]i are formed, spreading from cell to cell from different wave-initiating cells near the islet periphery [13,16,3436]. An increase in glucose concentration is coded as a fractional increase in activity within a time period, termed also relative active time [16,27,33] or duty cycle. The mechanistic substrate for such cohesive functioning of beta cells is intercellular communication via gap junction channels, consisting of connexin 36 (Cx36) [34,3740]. Cx36 provides both metabolic and electrical coupling between spatially organized heterogenous beta cells [34,39,41]. While other mechanism, such as autonomic innervation [42], autocrine and paracrine [43,44] signaling also contribute to cell-cell communication, Cx36 have been shown to play the main role in synchronizing beta cell collectives and maintaining proper insulin secretion [45,46]. Indeed, expression of Cx36 is decreased in diabetic conditions [47,48] leading to desynchronisation of [Ca2+]i oscillations and perturbations in pulsatile insulin secretion [37,40,4952]. Hence, the gap-junctional connections among beta cells are imperative for optimal beta cell function and comprehending their collective dynamics holds significance in elucidating the mechanisms underlying diabetes pathogenesis and its treatment.

Due to their highly heterogeneous nature, the presence of distinct subpopulations, and an ever-changing environment, beta cells display intricate yet coherent intercellular activity patterns [8,53]. Because coordinated intercellular activity is not only crucial for tightly regulated insulin secretion but is also known to be altered in diabetes, researchers are investing considerable effort in describing and studying how collective rhythmicity is established in beta cell populations and how the underlying mechanisms change in disease. In recent years, the emergence of network analyses has provided a promising tool for evaluating data obtained through advanced multicellular imaging, with the goal to objectively characterize collective activity in islets [16,42,5458]. In this approach, individual cells serve as nodes, and their positions correspond to their physical locations within the tissue. The connections between cells reflect functional associations and are determined based on the temporal similarity of the measured cellular dynamics, most often [Ca2+]i activity [56]. The application of network approaches has uncovered a modular organization in the functional beta cell networks, that exhibit greater heterogeneity than anticipated in a gap junction coupled syncytium. The identified indicators of small-worldness and a heavy-tailed degree distribution imply the existence of highly connected cells, called hubs [54,56]. Although their precise function remains somewhat enigmatic, these hubs are believed to represent a subpopulation with distinct attributes that confer upon them an above-average impact on the synchronized behavior [8,13,55,5961]. Furthermore, the collective responses to stimulation and the mediation of intercellular signals were also found to be influenced by other beta cell subpopulations. Specifically, the first responder cells were found crucial in mediating the responses to increasing stimulation during first phase of the islet’s response [57], whilst the wave initiator cells act as triggers of intercellular signals that synchronize the cells [13,34,62], being thereby presumably implicated in the regulation of pulsatile insulin release during the second phase [14,41]. In recent years, advanced methodological approaches, including optogenetics, photopharmacological methods, and RNA sequencing, along with network analyses, have unveiled specific characteristics within these subpopulations [8,13,41,54,63]. Acknowledging their unique attributes and significant contribution to shaping overall islet activity, there is a growing interest in their role in diabetes development [45,54,64].For this reason, it becomes even more important to precisely define these subpopulations, and objectively determine them through network analyses.

Nevertheless, due to variations in experimental preparations, microscopic imaging techniques, the nature of recorded signals, the following signal processing techniques, and the methods for deriving functional connectivity patterns that are employed by different research groups, comparing findings and integrating them into a comprehensive bigger picture becomes challenging even for experts in islet research. Additionally, the introduction of new terminology has further contributed to disputes in data interpretation, as well as to apparent contradictions regarding functional connectivity and the role of different beta cell subpopulations, which can be in part attributed to aforementioned methodological discrepancies [8,13,27,53,6567]. To at least partly address these issues, we present here a systematic analysis of how different experimental designs and computational approaches impact the results obtained from network representations of multicellular islet activity. Specifically, we analyze how the results are affected by different methods used to evaluate coordinated cellular behavior and network construction, different timescales of observed oscillatory calcium activity, different mouse strains used for tissue slice preparation, and the type of experimental preparation (i.e., tissue slices vs. isolated islets). All of the above represents some of the most prevalent genuine variations due to the diverse nature of work, experimental techniques, and the availability of equipment in laboratories worldwide.

Results

The role of different methods for the evaluation of time series similarities

We start by examining the effect of the type of time series similarity measure used to extract functional beta cell networks. We analyzed the beta cell [Ca2+]i responses to glucose stimulation obtained by means of multicellular confocal imaging in acute tissue slices from NMRI mice. The stimulatory glucose concentration was 12 mM and a 15-minute interval of sustained oscillatory activity (i.e., plateau phase) was used for the analysis, as indicated in Fig 1A. Fig 1B shows the extracted functional networks obtained by three different techniques: Pearson correlation coefficient (left panel, red), coactivity coefficient (blue, middle-left panel), and mutual information (purple, middle-right panel). A variable threshold was used so that roughly the same average node degrees (kavg between 8 and 9) were obtained in all three networks, facilitating a robust inter-network comparison. The comparison of methods for constructing networks from similarity matrices is analyzed separately in continuation. The right-most panel of Fig 1B shows calculated network parameters. Upon visual inspection of the networks and their corresponding parameters we can see a high degree of similarity between all displayed networks. Consistent with previous findings, all the networks exhibit high levels of clustering, modularity, and small-worldness [16,27,56,6870]. In Fig 1C and 1D the degree and edge length distributions of the same networks as in Fig 1B are presented, and Fig 1E shows the calculated internetwork similarities (see Methods section for details). All three panels further underline the observed resemblance between the networks extracted from different methods. Furthermore, the similarity in the degree distribution of all three networks suggests comparable variations in the number of functional connections, and, in addition, the level of heterogeneity indicates the presence of hub cells in all three cases.

thumbnail
Fig 1. The role of time series similarity measures on the functional beta cell network characteristics.

(A) Average signals of unprocessed (black line) and fast oscillatory activity (grey line) (upper panel) and the raster plot (lower panel) showing the binarized fast beta cell dynamics of all cells in the slice. (B) Functional networks designed with fixed average network node degrees (kavg ≈ 8) based on three distinct time series similarity measures, and the corresponding network parameters. The correlation method (red) is represented on the left, the coactivity method (blue) in the middle, and the mutual information method (purple) on the right. Colored dots indicate physical locations of cells within islets, while grey lines represent functional connections between them. The table on the right shows the average network node degree (kavg), average clustering coefficient (Cavg), modularity (Q), global efficiency (E), relative largest component (Smax), and the small-world coefficient (SW) for each network. Degree distributions (C) and distribution of functional connection lengths (D) for networks presented in panel A. Boxes on panels (D) determine the 25th and 75th percentile, whiskers denote the 10th and 90th percentile, and the lines within boxes indicate the median values. E) Jaccard internetwork similarity of the three networks extracted from different methods. (F) Pairwise comparison of node degrees in all networks from panel (B). Gray dots represent the node degree of the same cells in the graphs. (G) Relative active time as a function of node degree in networks built with the correlation method (red), coactivity method (blue) and mutual information method (purple). Dots denote average values and bars denote the standard error.

https://doi.org/10.1371/journal.pcbi.1012130.g001

To investigate the level of similarity between networks further, we present the pair-wise relationships between the node degrees of the same cells in all three constructed networks in Fig 1F. A clear relationship can be observed, with cells having a high degree in one network also having a high degree in the other. There is a consistent trend of matching connection numbers between different networks. The highest overlap is noticed between the coactivity and mutual information, as both these methods rely on binarized signals. Furthermore, previous studies have indicated that there is a tendency that the cells with many functional connections exhibit a higher-than-average activity [13,27]. This characteristic was alluded to when describing hub cell [Ca2+]i dynamics as “preceding and outlasting that of the follower cells” in Johnston et al, 2016 [54]. As such, hub cells typically manifest durations of oscillations exceeding the average, which contributes to higher cellular activity. Nevertheless, it is crucial to emphasize that this does not inherently suggest the role of hubs as wave initiators. In other words, the cells with longest oscillations do not necessarily initiate the intercellular waves. Subsequent analyses conducted on more extensive datasets in human [61] and mouse [13] islets revealed a distinct lack of overlap between wave initiators and hubs, although they confirmed that hubs tend to have the longest oscillations. Here, we aimed to confirm whether a comparable relation between relative active time and node degrees can be obtained when employing different techniques to quantify the similarity between [Ca2+]i signals. In Fig 1G the relationships between the relative active time and node degree of the same cells are depicted. For all three methods, very similar trends are observed. Specifically, there is a positive relationship between the relative active times of cells and their corresponding node degrees, which indicates that the extracted functional relationships are roughly independent of the methods used to evaluate synchronous cellular activity.

The role of different techniques used for network construction

Since no significant differences were found among various methods for evaluating intercellular synchronicity, we proceeded with the correlation method in subsequent analyses, which also has the advantage of not requiring signal binarization. In Fig 2 we explore the influence of network construction methods on the functional beta cell network structure and their relationship with the cellular activity parameter. Fig 2A displays functional networks of three different islets (rows) constructed using a fixed similarity threshold (left column), fixed average node degree (middle column), and the multilayer minimum spanning tree (MST) technique (right column). Due to the differences in [Ca2+]i signals, such as differences in overall activity, nature of [Ca2+]i signals (e.g., fast electrical vs. slow metabolic oscillations, see below for more details), presence of noise, etc., the average node degrees of networks constructed with the fixed similarity thresholds can vary greatly (between 4.4 and 17.5 in the three islets analyzed, see red networks in Fig 2A), whereas the average degrees for the other two techniques are fixed around 8 (see blue and purple networks in Fig 2A). In Fig 2B, 2C and 2D, we present the relationships between the relative active times for the data pooled from all three islets. Notably, a positive correlation between node degree and relative active time is inferred only for the variable threshold, i.e., a fixed average degree, and multilayer MST techniques. This tendency of hub cells exhibiting higher-than-average activity is in accordance with previous reports. In contrast, this relationship is apparently blurred for the fixed threshold technique due to variations in the number of connections along with differences in intrinsic activities among different islets, and hence even an opposite trend can be obtained, despite the fact that the relation is positive in all individual islets.

thumbnail
Fig 2. The role of network construction method in determining the beta cell network structure and the relationship to cellular activity parameters.

(A) Networks of three different islets (rows) constructed with three distinct construction methods (columns) with indicated average network node degrees (kavg). Construction methods are fixed correlation threshold (left), fixed average network node degree (middle), and multilayer minimum-spanning tree (right). (B-D) Relative active time of cells as a function of the relative node degree for networks constructed with the fixed correlation threshold method (B; red), fixed average network node degree method (C; blue) and multilayer minimum-spanning tree method (D; purple). Dots indicate the average values of cells within the same degree intervals and the error bars the corresponding SE. Note that the node degrees were normalized to facilitate the comparison of different islets.

https://doi.org/10.1371/journal.pcbi.1012130.g002

To further illustrate the challenge of pooling data from different islets, we present experimental data in Fig 3, where islets were subjected to stimulatory glucose in the first interval and subsequently treated with the GLP-1 receptor agonist exendin-4 (Ex-4) in the second interval. The precise protocol and the average unprocessed Ca2+ signal of all cells is shown in Fig 3A. Previous research has demonstrated that Ex-4 increases the density of functional beta cell networks [71], and our results validate this observation under both fixed threshold Rth and fixed average degree kavg methods (Fig 3B). However, the fixed Rth approach encapsulates substantial heterogeneity among islets, obscuring the effects of Ex-4 stimulation when aggregating data across multiple islets, despite discernible trends at the individual islet level. In such scenarios, employing a fixed kavg technique proves to be a more appropriate option for analyzing network differences induced by the pharmacological agent. Specifically, utilizing a variable threshold to maintain a consistent average degree in the initial interval mitigates inter-islet heterogeneity. Subsequently applying the same threshold in the second interval enables an unbiased assessment of the pharmacological intervention, normalized by the network characteristics observed in the first interval of each respective islet. This not only facilitates a robust comparison of network parameters across different islets but also ensures a more accurate statistical evaluation, as demonstrated in Fig 3C.

thumbnail
Fig 3. The role if different methods for designing and quantitative analysis of functional beta cell networks.

A) Average Ca2+ signals from all cells within a representative islet are depicted, stimulated with 10 mM glucose and 10 mM glucose + the GLP-1 receptor agonist exendin-4 (Ex-4), with specified intervals for network analysis. B) Functional networks were constructed for two intervals (Interval 1: 10 mM glucose only; Interval 2: 10 mM glucose + 20 nM Ex-4) across two different islets, utilizing two distinct network design techniques (fixed Rth = 0.75 and fixed average degree with kavg (Int. 1) = 8). The use of the fixed Rth method resulted in significant variations in network density between the two islets, complicating the comparison of network metrics. Conversely, fixing the average degree in Interval 1 normalized inherent differences in overall coherence of intercellular Ca2+ activity, facilitating the assessment of the pharmacological manipulation effect across different islets. C) To illustrate the issue posed by the fixed Rth method and the resulting high disparities in network densities, we compared the pooling of data from 10 islets subjected to closely matched protocols using both thresholding techniques (i.e., fixed Rth and fixed kavg). While both methods revealed a denser network in Interval 2 in response to Ex-4, these differences were almost completely masked by inter-islet variability inherent in the fixed Rth method. In such scenarios, normalizing the average degree proves to be the superior approach, as it facilitates a robust evaluation of data from different islets. Data used for this analysis is from Ref [71].

https://doi.org/10.1371/journal.pcbi.1012130.g003

It is worth noting, however, that in some cases, the fixed Rth method is the more suitable choice. It is known that with gradually increasing glucose concentration, both activity and intercellular communication levels increase, leading to denser functional networks in these cases [72]. When using fixed kavg or multilayer MST methods, which impose a specific number of connections regardless of the nature of activity, such networks do not differ in the number of connections, which is incorrect in these scenarios. Moreover, under conditions of low stimulation, an unjustifiably large number of connections is obtained due to, for example, a low threshold, which does not reflect correlations in dynamics but rather random associations. Such an example is illustrated in S1 Fig.

The role of the mouse strain used in tissue slice preparation

Laboratory mice are a vital source of islets of Langerhans in beta cell physiology research; however, various laboratories employ various mouse models. Previous research has indicated that there is a considerable phenotypic variation between different mouse strains as well as substrains of the inbred strains [7375] which manifest themselves also in beta cell responses to glucose and [Ca2+]i signalling characteristics [33,76]. For that reason, we investigate here whether the functional beta cell network structure extracted from multicellular [Ca2+]i recordings in tissue slices depends on the mouse strain. To this purpose, we compared the beta cell networks from outbred NMRI mice and inbred C57BL/6J mice. We used the correlation method to evaluate similarity between [Ca2+]i signals and the fixed average degree method (kavg = 8) to construct networks. In all recordings we used a 6-10-6 mM glucose protocol, as presented in Fig 4A. Intervals of 10–20 min sustained activity in the plateau phase were then used for the analysis. In Fig 4B we show typical networks from both strains, which, upon visual inspection, exhibit a rather similar topological organization. To provide a more detailed and quantitative insight, we computed various network metrics from pooled data from multiple islets. Results in Fig 4C, 4D and 4E indicate that the edge length, clustering coefficient, and degree distributions are very similar. Furthermore, the computation of network parameters presented in Fig 4F has revealed that beta cell networks from different mouse strains exhibit a similar degree of functional segregation, efficiency, and small-worldness; none of the results were identified as significant.

thumbnail
Fig 4. The impact of mouse strain used for tissue slice preparation on the parameters of beta cell networks.

(A) Average signals of unprocessed and fast oscillatory activity and the raster plot showing the binarized fast beta cell dynamics of all cells in slices from NMRI mice islets (upper panel, blue) and C57BL/6J mice islets (lower panel, purple). (B) Functional networks derived from representative recordings in islet from NMRI (blue, upper panel) and C57BL/6J (purple, lower panel) mice. (C) Edge length distributions, (D) clustering coefficient distributions, and (E) node degree distributions from a pooled data set from NMRI (blue) and BL6J (purple) mouse recordings. (F) Network parameters for extracted networks from NMRI and BL6J mouse recordings: modularity (left), relative largest component (middle left), global efficiency (middle right), and small-worldness coefficient (right). Dots represent values of individual recordings with horizontal lines indicating median values. Boxes on panels (B) and (C) determine the 25th and 75th percentile, whiskers denote the 10th and 90th percentile and the lines within boxes indicate median values. Data was pooled from islets/cells: 6/779 (NMRI), 6/617 (C57BL/6J). In all recordings, the islets were stimulated with 10 mM glucose and 10–20 minute intervals in the plateau phase were used for the analysis.

https://doi.org/10.1371/journal.pcbi.1012130.g004

The role of different time scales of oscillatory [Ca2+]i activity and time series preparation

Next, we investigate how the type of oscillatory activity and signal preparation impact the functional beta cell network topology. To this purpose, we performed prolonged multicellular imaging in tissue slices from NMRI mice. Fig 5A displays an average [Ca2+]i signal of all cells in a representative islet under stimulation with 8 mM glucose. Three different temporal traces are presented: the unprocessed (i.e., raw recorded) signal (top, red), the filtered slow oscillations (middle, purple), and the filtered fast oscillations (middle, blue).The fast and slow oscillations principally represent the electrical and metabolic activity of cells, respectively [77]. The lower panels in Fig 5A feature raster plots depicting binarized activity of the slow and fast oscillatory component. Notably, both types of oscillatory activity exhibit distinct, regular patterns. In Fig 5B we present correlation-based functional networks constructed with the fixed average degree technique for the three distinct signal types. A visual assessment points out a clear difference between the three extracted networks. The fast oscillation-based network (middle panel) exhibits shorter edge lengths and a more clustered, localized, structure, while the slow oscillation-based network (right panel) shows more long-range edges and a less clustered structure. A quantitative assessment of the networks confirms the observed differences. The slow oscillatory component network is more heterogeneous, less clustered and exhibits longer functional connections (Fig 5C, 5D and 5E). The reason for this is in the type of cellular dynamics the networks encode. The fast oscillations are representative of the electrical activity of cells, which is mediated by gap-junction-driven intercellular waves and thus contributes to the shorter, more clustered network structure which is quite similar to the underlying physical network. On the other hand, the slow component signal is associated with cellular metabolism which is to a greater extent affected by the similarity of intrinsic metabolic characteristics of cells and less by cell-to-cell coupling [56,70,78,79]. Interestingly, the raw-signal-derived functional network appears to be poised in between, which is somehow expected, as it encompasses both types of oscillatory activity. To evaluate the properties of different networks further, we quantified the extracted functional connectivity patterns using conventional network metrics (Fig 5F). The results indicate that the networks derived from different dynamical components have comparable values of the small-world coefficient and the relative largest component, but there are profound differences in modularity and global efficiency. Namely, the fast oscillations-derived network is more segregated and exhibits lower efficiency, primarily due to the less pronounced long-range connections. Moreover, we present in Fig 5G the relationship between the relative active time of cells and their corresponding node degrees in all three types of networks. The tendency of hub cells being the most active is most pronounced in the case of fast oscillations, whereas the relation is less apparent for the raw and slow component. Notably, the latter aligns with recent theoretical predictions [78]. Finally, we assess the similarities between the three networks and present in Fig 5H the pair-wise relationships between the node degrees in different networks. The results indicate that the strongest relation exists between the fast and raw oscillatory signals, while the relationship is the weakest between the fast and the slow component. To investigate this in further detail, we quantified the overlap between different networks, including the hypothesized structural network that was modeled as a geometric network in which nearby cells are connected. From Fig 5I, we can observe a substantial similarity in both inter-network similarity and overlap of hub cells between the unprocessed signal and the signals of both oscillatory components, with a higher level of similarity observed in the fast component. This is expected in signals from slices, as the fast component is very pronounced. However, the key point is that the highest level of similarity between the structural network and the functional networks is obtained from fast oscillations, while the similarity between the structural and slow networks is substantially lower. Similarly, the connection between the fast and slow component-derived networks is relatively low, as previously indicated by the results in Fig 5H. These quantitative results can be further visually assessed with the illustrations in S2 Fig, depicting all four types of networks for all 5 islets included in the analysis. It can be observed that the networks of unprocessed signals and signals of the slow component contain many long-range connections, while those in the fast component network are significantly fewer, making it visually more similar to the structural network. Importantly, fast oscillations may be more strongly determined by slow oscillations, such as in the case of compound oscillations [80,81]. Such an example is depicted in S3 Fig and in this case, the functional network based on slow oscillations rather than fast oscillations is most similar to the functional network based on the raw signal. However, the functional network based on fast oscillations remains the one that is most similar to the structural network (S3B and S3C Fig).

thumbnail
Fig 5. The role of the type of oscillations and signal preparation on the characteristics of functional network topology and the relations to the cellular activity.

(A) Unprocessed (red), fast-component only (blue), and slow-component only (purple) average [Ca2+]i signal of all cells in the islet from acute tissue slice from NMRI mouse. The lower panels display raster plots that represent the binarized activity of the slow and fast oscillatory components. (B) Functional networks designed based on raw cellular signals (left), fast-component only signals (middle), and slow-component only signals (right). Networks were constructed with the fixed average network node degree method (kavg≈8.0) based on time series correlations as the similarity measure. Distribution of node degrees (C), clustering coefficients (D) and functional connection lengths (E) for the three networks presented in panel B. (F) Network parameters extracted from functional connectivity maps derived from different oscillatory components. (G) Relative active time of cells as a function of their corresponding node degrees in networks constructed based on raw signals (red), fast-component only signals (blue) and slow-component only signals (purple). Colored dots represent average values of cells within the same degree intervals and the error bars denote SE. Individual values were normalized by the average value of the relative active time within the given islet so to ease comparison between different islets. (H) The pairwise relationships between node degrees in different networks. The grey dots denote values from individual cells and the black line indicates the linear fit, whereby R2 indicates goodness-of-fit. I) Similarity between different types of networks (left) and the relative overlap of hub cells (right), identified as the top 1/6 of the most connected cells. The structural networks were modeled as equivalent geometric networks, in which nearby cells are deemed connected (see Materials and Methods and S4 Fig). Boxes in panels (D) and (E) determine the 25th and 75th percentile, whiskers denote the 10th and 90th percentile and the horizontal lines within boxes indicate the median values. Dots in panel (F) indicate the values from individual islets and the horizontal line denote the median. Stars denote statistical differences; *p<0.05,**p<0.01. Data presented in panels (F-I) is based on 5 different islets.

https://doi.org/10.1371/journal.pcbi.1012130.g005

Functional connectivity networks in isolated islets

In addition to acute tissue slices, isolated islets play a prevalent role in beta cell research, including in the context of collective activity network analyses. Thus, we proceed with analyzing the nature of multicellular dynamics and the underlying functional networks within islets isolated from C57BL/6J mice. In Fig 6A, we present the responses of a representative isolated islet upon transitioning from 2 mM to 11 mM glucose. The cells exhibit an initial, profound elevation in [Ca2+]i levels, followed by the emergence of coordinated [Ca2+]i oscillations after approximately 8–10 minutes. The raster plots indicate that these oscillations frequently span the entire islet. The functional network extracted from the phase of sustained oscillatory activity, constructed based on time series correlation as the similarity measure along with the fixed average network node degree method, is shown in Fig 6B. The characterization of beta cell networks was based on 5 different isolated islets subjected to the same protocol. In the table shown in Fig 6C the average values of network parameters are provided and Fig 6D shows the pooled degree distributions. We can observe that the topological parameters of networks from isolated islets do not differ much from those in slice-based networks: they are quite modular and exhibit features of small-world networks. However, upon visually evaluating the network illustrated in Fig 6B and considering the characteristics of clustering coefficient (Fig 6E) and functional connection length distributions (Fig 6F), it becomes evident that the networks observed in isolated islets exhibit properties that are more similar to the networks characterized by slow activity in slices. Note that for comparison the data on fast and slow activity-derived networks from slices from C57BL/6J mice are provided separately. For this comparison, the same dataset was used as in Fig 4, where also the stimulatory glucose concertation was similar (i.e., 10 mM). In contrast to fast oscillation-based networks in slices, isolated islet networks manifest a higher efficiency, a reduced modularity, and low clustering coefficient values. Moreover, the distribution of relative connection lengths indicates that there is a larger fraction of long-range connections in isolated islets. All these attributes can be observed in slow oscillation-based networks in slices. Notably, within isolated islets, a discernible trend emerges where cells with an increased number of functional connections consistently demonstrate higher relative active times—reminiscent of the observed behavior in slices—regardless of the temporal aspect (Fig 6G).

thumbnail
Fig 6. Functional beta cell connectivity patterns in isolated islets.

(A) Average [Ca2+]i signal of a representative isolated islet recording with indicated plateau phase for signal analysis (upper panel) and corresponding binarized oscillatory activity of all cells in the recording (lower panel). (B) Extracted functional network based on cellular signals in panel (A) constructed with the fixed average network node degree method with an average network node degree kavg≈8.0. Green dots represent physical locations of cells within the islet and grey lines indicate functional connections between them. (C) Extracted average functional network parameters: average network node degree (kavg), average clustering coefficient (Cavg), modularity (Q), global efficiency (E), average shortest path length (Lavg), small-world coefficient (SW), and relative largest component (Smax). Degree distributions of all extracted functional networks (D), and corresponding distributions of clustering coefficients (E), and relative edge lengths (F). To ease comparison between different islets, the physical lengths of connections were normalized with the average distance to the eight nearest neighbors. Additionally, in panels (D-F) data illustrating network attributes derived from fast and slow activities in slices from C57BL/6J mice are presented for comparison. (G) Relative active time as a function of node degree for all extracted functional networks. Boxes on panels (E-F) determine the 25th and 75th percentile, whiskers denote the 10th and 90th percentile, and the lines within boxes indicate the median values. Dots in panel (G) represent average values and vertical bars denote the standard error. Data for panels (C-G) for isolated islets was pooled from islets/cells: 5/468 and for slices the same dataset was used as in Fig 4 (islets/cells: 6/617). *p<0.05,**p<0.01, ***p<0.001.

https://doi.org/10.1371/journal.pcbi.1012130.g006

Furthermore, to further assess the differences and similarities between beta cell networks from slices and isolated islets and how they relate to different types of oscillatory activity, we present in S3 Fig an analysis of an isolated islet where the fast component of oscillations was relatively well present, which is frequency-wise highly comparable to that in tissue slices. This enabled the separate consideration of individual oscillatory components, and similarly to tissue slices, it was found in this case as well that there is a significant similarity between the structural network and the functional network obtained from the fast component, while the similarity between the slow and structural is considerably lower. It is also worth mentioning that in isolated islets, there is much greater similarity between networks derived from unprocessed signals and the slow component, whereas in tissue slices, there is greater similarity between networks based on unprocessed signals and the fast component. The reason for this is that in tissue slices, fast oscillations are the more dominant type of signal, while in isolated islets, slow oscillatory activity prevails.

Discussion

Functional connectivity analysis is a powerful tool applicable to studying the interactions between different components in a plethora of real-life systems. In recent years, it is becoming increasingly more popular to describe interactions between individual cells, particularly within the islets (for review see [56]). However, due to relatively demanding computational approaches, encompassing both data extraction and subsequent analyses of coordinated functioning, obtaining patterns of functional connectivity is not straightforward and can easily become ambiguous. In neuroscience, where the greatest progress in this field has been made, it has become evident that objectively assessing connectivity patterns is challenged by various objective reasons tied to experimental variations and computational methodologies, such as thresholding techniques [8284], techniques used for data pooling [85,86], number of sensors used to record brain activity [87,88], and the selection of frequency intervals [89,90]. Most importantly, similar issues are witnessed in the network-based analysis of spatiotemporal cellular dynamics in islets. More specifically, different research groups employ diverse experimental techniques and preparations leading to discrepancies in types of oscillatory signals and the multicellular activity is recorded at varying spatial and temporal resolutions. There are also variations in how recordings are preprocessed before network analysis, as well as in the techniques used for the analysis itself. These, along with some terminological discrepancies in the scientific literature, are the primary reasons why we chose to investigate how various factors influence network analyses and their interpretation.

First, we evaluated the role of metrics that are used for the evaluation of synchronized activity between the measured cellular dynamics. We compared three different methods, namely one that is based directly on the recorded [Ca2+]i activity (Pearson’s correlation), and two that are based on binarized time series (coactivity and mutual information). It turned out that irrespective of the method used to quantify synchronous behavior, similar networks are obtained, characterized by small-worldness, modularity, high degree of clustering, a heavy-tailed degree distribution which indicates the presence of hub cells, and a similar relation between the relative active time and the node degree (see Fig 1). Another crucial aspect in the process of extracting functional connectivity maps involves the thresholding of similarity matrices. As highlighted in Figs 2A and 3B, utilizing a fixed threshold can yield significant disparities among different islets, potentially introducing biases into the relations drawn from aggregated data. To mitigate this concern, using a variable threshold and a fixed average degree proves advantageous. With this approach we can firmly evaluate the effect of pharmacological interventions or extract the relations between network and classical physiological parameters when data is pooled from multiple islets, as a variable threshold can mask the inter-islet heterogeneity (see Figs 2C and 3). Specifically in multi-phase experiments, where consecutive intervals have to be analyzed [71,91], application of a variable threshold has proven beneficial, as it overcomes inter-islet variability. For example, by establishing the variable threshold based on the first interval, thereby maintaining a fixed average node degree, one can consistently apply the same threshold to construct networks during the second interval. This normalization procedure facilitates an objective assessment of alterations in islet network structure, despite inherent differences in networks from different islets (Fig 3). It is important to note, however, that this method has a limitation: its fixed average number of connections prevents it from capturing the variations in overall synchronicity that are depicted by the network density. For instance, it is known that an increase in glucose concentration leads to increased and more global spatiotemporal activity, resulting in denser functional networks [27,72]. If a fixed average degree is then employed, these differences become obscured, and in conditions of low stimulation, numerous connections emerge that lack statistical significance. This occurs because, with a low threshold, these connections predominantly signify random associations rather than synchronized activity (S1 Fig). In this study we have also introduced a third option encompassing the construction of functional networks through a multilayered MST. A notable advantage of this method lies in its absence of explicit thresholding, with the singular free parameter being the number of layers, which in turn specify the average degree. Nonetheless, the drawback of the minimum spanning tree method is that it enforces at least one connection to each cell (or more in case of multilayered MST), so that even the cells which are completely desynchronized can have a comparable number of functional connections as an average cell. Therefore, while the method is attractive for its apparent objectivity, its appropriateness diminishes when the signals are rather heterogeneous and if there are subpopulations of cells whose dynamics are weakly or not at all correlated with the rest of the cells (such as those of alpha cells, see S4 Fig). To sum up at this point, the choice of the best method to construct networks is not always straightforward and may depend on the context, i.e., both the experimental protocol and the parameters we want to objectively describe through network analysis. In doing so, we must, of course, be aware of the strengths and weaknesses of different approaches.

In previous studies, variations in glucose-induced [Ca2+]i activity among different mouse strains and substrains have been reported. Compared to outbred NMRI mice, cells from the inbred C57BL/6J and C57BL/6N substrains show a rightward shift in activations and earlier deactivations. In addition, during the plateau phase, the encoding mechanisms to enhance calcium activity in response to glucose differ quantitatively in all three groups [33]. Secretagogues other than glucose also cause [Ca2+]i oscillations to vary greatly [76]. Generally, however, there are similarities between C57BL/6J, C57BL/6N, and NMRI mice in the sense that all three groups showed glucose-dependent activation and deactivation responses, as well as a 3% increase in relative active time per millimole of glucose [33]. Notably, up until now, differences between strains of mice at the level of multicellular activity have not been studied. In this study, we addressed these questions using network analyses and found that the functional networks of islets in different mice are structurally very similar. Apparently, the mechanisms that coordinate fast oscillatory activity across the islets from NMRI or C57BL/6N mice, i.e., gap-junction mediated depolarization and [Ca2+]i waves, are the same and do not differ between mouse strains.

In response to glucose and many other secretagogues, electrical activity, intracellular calcium, and insulin secretion oscillate in synchrony at two different time scales [92,93]. The first are the so-called metabolic or slow oscillations with a frequency of around 0.1–0.2 min-1, and the second the so-called electrical or fast oscillations with a frequency of around 1–5 min-1 [53,94]. Noteworthy, fast oscillations show variations and have the highest frequency rates around the peaks of the slow component and the lowest around the nadirs [77,95,96]. Additionally, the relative active time or duty cycle of the fast component characteristically increases with increasing stimulation, whereas the frequency of slow oscillations remain unaltered [27,33,9698]. According to the recent metronome model of beta cell function, slow oscillations set the pace for insulin pulses, whereas the fast oscillations fine-tune their amplitude [94]. Both slow and fast oscillations are phase-locked between different beta cells within a given islet by means of intercellular waves [14,34,35,55,62,68]. In accordance with this, the average correlation between calcium traces of different cells from the same islet decreases with intercellular distance for both the slow and the fast component, implying that intercellular coupling mediates the synchronicity of both types of oscillations [78,96].

If one constructs and compares functional connectivity maps for the raw signal and both dynamic components separately (Fig 5), the distributions of node degrees do not differ significantly. However, the networks of fast oscillatory activity are more locally clustered and segregated, more modular, and have lower average edge lengths and global efficiency, while the slow oscillations are principally more global, resulting in numerous long-range connections and consequently a more cohesive structure that shows a lower modularity and higher global efficiency. Importantly, for the raw signal, it seems that except for the node degree, most of the network measures are more strongly determined by the slow component [56]. A logical consequence of the abovementioned differences in functional network structure is the finding that there is a relatively weak correlation between the fast and slow network layer [96], implying that different synchronization principles are at work [70,78], and one should not directly compare results of studies relying on fast oscillations with the ones relying on slow oscillations. Importantly, even with the same experimental model, e.g., isolated mouse islets, and set of analytical tools applied to extracting and analyzing [Ca2+]i oscillations, islets with preponderance of fast, mixed or slow oscillations might coexist [99101], and in this case, data should not be simply pooled, since this may obscure relevant biological differences, but analyzed for the two temporal components and for oscillatory phenotypes separately. Extrapolating this reasoning further, the caveats we pointed out in this paragraph should also be kept in mind when comparing experimental traces from different animal models, even when using the same experimental approach and the same set of analytical tools. For instance, the presence and relative importance of fast and slow oscillations may vary between beta cells from zebrafish [55,102], mice [100,103], rats [104,105], sand rats [106,107], pigs [108,109], and humans [61,110], to name only a few. To facilitate interspecies comparison, future studies shall clearly specify the type of oscillations they are addressing. Finally, at present, it is difficult to experimentally compare the relationship between the structural networks of beta cells and their functional counterparts, but modelling studies suggest that the intricate structure of functional beta cell networks based on fast and slow oscillations may be at least partly explained by heterogeneity in beta cell activity and heterogenous intercellular coupling [68,70,78].

Different groups that employ network measures in their analyses typically use different experimental approaches to obtain [Ca2+]i traces. While most groups use cultured isolated islets in combination with CCD camera-based or confocal imaging, some use the acute tissue slices in combination with confocal imaging. To be able to compare findings from different groups or combine them into a coherent bigger picture of islet network properties, these differences also need to be addressed as they are an important possible systematic confounding variable. Essentially, the methodology and experimental setup would not seem to be key parameters if they did not entail differences in the nature of the oscillatory signals. In tissue slices fast or mixed oscillations are more predominant (see Figs 4A or 5A), whereas in isolated islets the slow oscillations are predominant (see Fig 6A). Here, we explicitly demonstrated that the distinct nature of oscillations leads to different functional beta cell networks. While some network properties in fast-derived and slow-derived networks are similar, such as heterogeneity and small-worldness, they fundamentally differ from each other, and the significance of certain subpopulations in one network is therefore not equivalent to that in the other network. Moreover, even if oscillations qualify as fast, in isolated islets, they are typically longer than 10 seconds at concentrations > 10 mM glucose [34,92,99], whereas in slices, they tend to be shorter than 10 seconds [16,27,33,111]. The exact mechanism behind these differences remains to be explained, but in addition to possible differences in ionic composition and the presence of additional secretagogues in the extracellular fluid that can affect the patterns of oscillations [92,112,113], the mechanical and enzymatic stress during preparation of isolated islets [114,115], as well as culture conditions and duration [99,116] have been put forward as possible sources of these differences. More specifically, alpha cells have been suggested as a potential source of local proglucagon peptides [117]. They are primarily situated in the mantle of pancreatic islets in mice, and this outer region is particularly susceptible to damage during the islet isolation process, potentially resulting in the loss of alpha cells during islet preparation. Given that both glucagon and GLP-1 have been shown to elevate the frequency of oscillations in beta cells, the diminished intra-islet alpha cell signalling could be a contributing factor to the observed decrease in beta cell oscillatory frequency in isolated islets [71,91,118,119]. Further, there may be a run-down of certain ion channels and changes in the expression [120, 121] with time, which obviously impact the identity and physiology of beta cells in the cultured isolated islets more than in the immediately used islets in slices [122,123]. This theory is at least partly supported by the finding that oscillations in mouse islets cultured for less than one day closely resemble oscillations in non-cultured islets [103,124] and in islets studied in vivo [125,126] or rapidly after the death of the animal [92], as well as the oscillations in tissue slices [27]. Until the influence of the above factors is fully understood, we can provide at least two practical suggestions. First, studies on isolated islets and tissue slices should always exactly state what the composition of the extracellular fluid was, and which type of oscillations were used for the network analyses, as well as provide details about the basic characteristics of these oscillations, i.e., their frequency and duration. Second, freshly microdissected islets or islets cultured for shorter time periods may yield results that are more closely comparable with results from tissue slices. Finally, the above advice also applies for studies utilizing yet other experimental preparations, such as in vivo imaging of isolated and transplanted mouse islets in the anterior chamber of the eye [127] and islets from other species, as mentioned in the preceding paragraph. In the present study, we used a range of different stimulatory concentrations. They are not intended to illustrate possible glucose-dependencies of different physiological and network metrics as these are covered elsewhere [13,27,33], but to demonstrate that the analytical tools work robustly across a range of frequently used stimulatory conditions. Given that the slow oscillations are rather glucose-insensitive in terms of their frequency in both slices [96] and islets [98] and that fast oscillations have comparable dose-response relationships in slices [27,33] and isolated islets [128130], we believe that the different concentrations we used did not introduce any critical bias and that most of our findings are applicable to concentrations beyond the range used here.

In conclusion, we would like to stress that the scope of network analyses has, in recent years, been extended to investigate intercellular interactions and functional connectivity patterns in different types of tissues. These encompass various kinds of neural assemblies [131], pituitary endocrine cells [132,133], astrocytes [134], yeast cells [135], distinct epithelial cell types [136,137], acinar cells [138], and hepatocytes [139]. As such, the insights we present herein hold relevance for comprehending the intricacies of collective cellular activity across diverse contexts, where the assessment of multicellular dynamics can be achieved through suitable imaging techniques. Moreover, in tandem with advancements in imaging methods, which are expected to soon enable the simultaneous high-resolution assessment of multiple variables defining multicellular activity, potentially even in three dimensions, it is imperative to stay attuned to progress on the computational front. Over recent years, a plethora of sophisticated methods has emerged for evaluating dynamic interactions within complex systems, such as multilayer networks [140,141], detection of higher-order interactions [142,143], information-theoretic metrics describing causal relationships [144,145], and deep learning-based methods [146,147]. These approaches hold substantial potential for further and more profound research, extending even into the realm of multicellular systems, as already demonstrated by some recent studies [13,56,148150]. We strongly believe that future progress in this field will rely on such interdisciplinary endeavors that combine cutting-edge experiments with innovative computational procedures. Along these lines, we anticipate a deeper understanding of how heterogeneous populations of interacting cells, placed within a dynamic and noisy environment, operate to ensure proper functionality, and how the regulatory mechanisms are altered in disease.

Materials and methods

Ethics statement

We conducted the study in strict accordance with all national and European recommendations on care and handling experimental animals, and all efforts were made to minimize the suffering of animals. Mice were used under protocols approved by the University of Colorado Institutional Animal Care and Use Committee (IACUC Protocol number: 00024) and The Administration of the Republic of Slovenia for Food Safety, Veterinary and Plant Protection (permit numbers: U34401-35/2018-2).

Animals and [Ca2+]i imaging in tissue slices

Slice preparation.

C57Bl6J and NMRI male and female mice were held in a temperature-controlled environment with a 12 h light/dark cycle and given continuous access to food and water. Preparation of mouse-derived acute pancreas tissue slices was executed as described previously in full [122]. In brief, after sacrifice with CO2 and cervical dislocation, the abdominal cavity is accessed via laparotomy and the papilla Vateri is clamped. 1.9% Low melting agarose dissolved in ECS containing (in mM) 125 NaCl, 26 NaHCO3, 6 glucose, 6 lactic acid, 3 myo-inositol, 2.5 KCl, 2 Na-pyruvate, 2 CaCl2, 1.25 NaH2PO4, 1 MgCl2, 0.5 ascorbic acid is heated to 40°C and injected through the bile duct. The pancreas is cooled with ice-cold ECS, extracted, and cut into tissue blocks, which are embedded in low melting point agarose and cut with a vibratome (VT 1000 S, Leica) to yield 140 μm slices. The slices are kept in HEPES-buffered saline (HBS) consisting of (in mM) 150 NaCl, 10 HEPES, 6 glucose, 5 KCl, 2 CaCl2, 1 MgCl2; titrated to pH = 7.4 with 1 M NaOH at room temperature and stained with a HBS staining solution containing 7 μM Calbryte 520 AM (AAT Bioquest), 0.03% Pluronic F-127 (w/v), and 0.12% dimethyl sulfoxide (v/v) for 50 min at room temperature. All chemicals were obtained from Sigma-Aldrich (St. Louis, Missouri, USA) unless stated otherwise. Individual tissue slices were placed into the recording chamber and used for one stimulation protocol. The recording chamber was continuously perifused with carbogenated ECS containing 6 mM glucose heated to 37°C at basal conditions. At 20–40 minutes, the perifusion was manually changed to stimulatory (8–12) mM glucose before it was returned to the basal glucose concentration.

Imaging.

Beta cell calcium dynamics were imaged using an upright confocal microscope system Leica TCS SP5 AOBS Tandem II with a 20X HCX APO L water immersion objective, NA 1.0, and an inverted confocal system Leica TCS SP5 DMI6000 CS with a 20X HC PL APO water/oil immersion objective, NA 0.7 (all from Leica Microsystems, Germany). A 488 nm argon laser was used to excite the fluorescent dye, and a Leica HyD hybrid detector operating in the 500–700 nm range was used to detect the fluorescence that was released (all from Leica Microsystems, Germany), as previously described [27,122]. The resolution used for time series acquisition was 512 X 512 pixels with a frequency of 2–10 Hz.

[Ca2+]i imaging in isolated islets

Islet isolation and culture.

Islets were isolated from mice under ketamine/xylazine anaesthesia (80 and 16 mg/kg) by collagenase delivery into the pancreas via injection into the bile duct. The collagenase-inflated pancreas was surgically removed and digested. Islets were handpicked and planted into the glass-bottom dishes (MatTek) using CellTak cell tissue adhesive (Sigma-Aldrich). Islets were cultured in RPMI medium (Corning, Tewksbury, MA) containing 10% fetal bovine serum, 100 U/mL penicillin, and 100 mg/mL streptomycin. Islets were incubated at 37C, 5% CO2 for 24–72 h before imaging.

Imaging.

An hour prior to imaging nutrition media from the isolated islets was replaced by an imaging solution (125 mM NaCl, 5.7 mM KCl, 2.5 mM CaCl2, 1.2 mM MgCl2, 10 mM HEPES, and 0.1% BSA, pH 7.4) containing 2 mM glucose and fluo4 AM [Ca2+]i sensitive dye (4 mM). After one hour the solution was replaced by dye-free imaging solution. During imaging the glucose level was raised from 2 mM to 11 mM. Islets were imaged using either a LSM780 system (Carl Zeiss, Oberkochen, Germany) with a 40x 1.2 NA objective or with an LSM800 system (Carl Zeiss) with 20x 0.8 NA PlanApochromat objective or a 40x 1.2 NA objective, with samples held at 37°C. The resolution was 512x512 pixels and time series were recorded with frequencies 1–2 Hz.

Pre-processing of recorded [Ca2+]i time series

Fluorescence signals of Calbryte 520 AM or Fluo-4 representing time series for manually selected regions of interest (ROIs), i.e., individual beta cells, were exported along with their corresponding coordinates using a custom software called ImageFiltering (copyright Denis Špelič) or ImageJ [151]. As both dyes can detect both fast- and slow-component in beta [61,152], data obtained by either dye was pre-processed equally. Time series that exhibited large artifacts, low signal-to-noise ratio, or dynamics inconsistent with beta cells were excluded after visual inspection. The recordings from tissue slices underwent band-pass filtering using a zero-lag filter to extract either the fast-activity component (with typical cut-off frequencies of 0.05 and 2.0 Hz) or the slow-activity component (with cut-off frequencies of 0.001 and 0.07 Hz). Similarly, the recordings from isolated islets underwent band-pass filtering to eliminate baseline drifts and capture the oscillatory component (with typical cut-off frequencies of 0.005 and 0.25 Hz). Fast-component signals from slices and oscillatory signals from isolated islets were further smoothed using an adjacency averaging procedure and then binarized by setting values to 1 (active state) for periods of increased [Ca2+]i signals or 0 (inactivity) for periods of low-amplitude signals. All subsequent analyses were performed either on the raw, filtered (fast or slow oscillatory component), or binarized cellular signals. The binarized signals were also used to calculate the relative active time. This metric represents the ratio of the time a given cell is in an active state, indicating thereby the overall cellular activity.

Evaluating synchronicity between [Ca2+]i traces

We use three different methods to quantify the similarity between recorded [Ca2+]i traces: the Pearson correlation, the coactivity, and the mutual information. The latter two are computed on the basis of binarized traces, whereas the Pearson correlation is computed on the basis of filtered traces. The Pearson correlation coefficient (PCi,j) between the i-th and j-th cell quantifies the linear relationship between the corresponding [Ca2+]i traces and is computed as: (1) where xi and xj are recorded time series for cells i and j, respectively, and σ represents their corresponding standard deviations. Eq (1) generates a value between -1 and 1, where -1 indicates completely anticorrelated time series and 1 implies identical time series. Although this is a popular similarity measure, it has some limitations. It cannot capture nonlinear relationships and can result in very high values in coupled oscillator systems like beta cell networks, which can limit its usefulness. Another popular method for measuring time series similarity is coactivity [132,153]. This method relies on binarized time series and measures the degree of simultaneous activity between cell pairs. It corresponds to the dot product of two normalized binary time series vectors, where each data point is one dimension of the vector. The coactivity coefficient between cell pairs i and j is calculated as: (2) where xb,i and xb,j are the binarized activity time series of cells i and j, respectively. This measure gives a value of 0 if there is no overlap in activity and 1 if the cells are completely coactive. The method quantifies the overlap of [Ca2+]i oscillations in a very straightforward manner but does not capture nonlinear relationships between cells.

Mutual information (MI) is another commonly used method to measure the similarity of time series [154]. It is a statistical measure that quantifies the amount of information shared between two random variables and the degree of dependence between them. In essence, it measures how knowledge about one variable can help predict the other. Mutual information is calculated based on the Shannon entropy of both time series (H(xb,i) and H(xb,j)) and their joint entropy (H(xb,i,xb,j)). For binarized time series of the cell pair i and j, MI is expressed as follows [155]: (3)

The value obtained by Eq (3) ranges from 0 (independent time series) to 1 (completely co-dependent time series). The Shannon entropy (H) of binarized time series i (xb,i) is defined as [155]: (4) where p(si) is the probability distribution that xb,i takes one of its possible values (si∈[0,1]). Values for p(si) can be easily calculated based on the frequency count (F(si)) of individual values si and the length (L) of the time series (p(si) = F(si)/L). The joint entropy of two binarized time series xb,i and xb,j is then calculated as [155]: (5)

In the case of two binarized time series, there are four possible value combinations for (si,sj). To ensure comparability between values of different cell pairs, the normalized mutual information () is calculated as: (6)

By using Eqs (1), (2), and (6), we can construct similarity matrices of size (N, N), whereby N stands for the number of cells, that encode the correlation, coactivity, and normalized mutual information between all cell pairs in individual recordings, respectively. Notably, MI captures also non-linear relationships between the discretized time series.

Network construction and analysis

The functional networks are extracted from similarity matrices, whereby different approaches can be used to this purpose. The most elemental way is to set up a predetermined threshold, so that a connection between a cell pair i and j is established only if their similarity coefficient (SCi,j) exceeds the pre-set similarity threshold (SCth). Specifically, cells i and j are connected if: (7) where SCi,j is based on one of the aforementioned similarity measures.

Alternatively, a variable similarity threshold technique can be used instead of a fixed threshold, which can create a network with a pre-set target average node degree, so the threshold is varied until a network with the target average node degree is designed. In our analyses the variable threshold was determined so that the resulting network had an average degree 8. This value was used to mimic the connectivity of realistic beta cell network architectures [156] and to obtain adequately dens networks suitable for analyses. However, it should be noted that previous studies have demonstrated that, within reasonable limits, the conclusions drawn from network analyses are not significantly influenced by the somewhat arbitrary choice of the average degree [13,96].

The third option is to utilize the multilayer minimum spanning tree (multilayer MST). In graph theory, a minimum spanning tree is defined as a graph that connects all the nodes with the minimum possible total edge weight, without forming any cycles. To construct the MST, the similarity coefficients of cell pairs (SCi,j) must be recomputed to an abstract distance measure (Di,j) using the following eq: (8)

Based on the computed abstract distances, an MST can be constructed with so-called greedy algorithms such as Kruskal’s [157] or Prim’s [158] algorithm. These algorithms create graphs with N-1 edges (N–number of nodes) which contain the lowest possible sum of edge weights (ΣDi,j) without creating any cycles. We expand this idea for the generation of a multilayer MST, where a single MST is computed sequentially for the same network, but already existing edges (i.e., cell pairs) are excluded from the calculation of the next MST layer. In our analyses we calculated four layers of MST’s, which yielded an average node degree of 8 (the average degree of the original MST is 2, and each of the three subsequent layers contributes an additional 2 degrees).

For each extracted network, we calculated several basic network parameters, such as average network node degree (kavg) and degree distribution, average clustering coefficient (Cavg) and clustering coefficient distribution, modularity (Q), global efficiency (E), relative largest component (Smax), and edge length distribution, and small-world coefficient (SW). See Ref. [159] for technical details and Ref. [56] for a physiological meaning of these specific network parameters.

Quantifying inter-network similarity

To quantify the similarity between networks constructed with different approaches, we computed the network similarity index (NSI) based on the sets of edges (A) for each individual network (α). For network pairs α and α′ NSI was calculated using the Jaccard similarity coefficient as follows: (9)

In other words, inter-network similarity is defined as the ratio between the cardinality, i.e., the total number of edges of the intersection of edges in networks α and α′, and the cardinality of the corresponding union. The resulting value of NSI ranges from 0 to 1, where 0 indicates no common edges between the networks and 1 indicates identical networks. This method was used to assess the similarity between functional networks derived from various oscillatory components and constructed with the above-described construction techniques. We additionally quantified the similarity of these networks with the postulated structural networks of islet cells, which we constructed as geometric networks by appropriate intercellular distance thresholding.

Methods for the time series processing, analyses of cellular signals, and network analyses were designed with Python programming language version 3.11.1, using the following packages: Numpy (https://numpy.org/), Matplotlib (https://matplotlib.org/), and NetworkX (https://networkx.org/). All code is available on the GitHub repository: https://github.com/MarkoSterk/beta_cell_analysis_suite

Supporting information

S1 Fig. Collective beta cell activity under the protocol of a glucose ramp.

A) Ca2+ traces of all responding beta cells in the slice (upper panel) and the corresponding raster plot of binarized fast Ca2+ oscillations. The glucose concentration was ramped from 6 mM to 12 mM, as indicated at the top. B) Functional beta cell networks extracted in different glucose concentrations and with different thresholding techniques. The fixed threshold approach (Rth = 0.8) leads to very different network structures under different stimulation levels. Under lower glucose, when the degree of correlated beta cell dynamics is low, the networks are sparse and segregated. With increasing stimulation, the networks become progressively more integrated and dense (i.e., average node degree kavg is increasing), highlighting the heightened intercellular coordination. Conversely, the fixed avg. degree and multilayer MST approaches fail to capture this behavior, as they enforce a fixed number of connections, irrespective of the level of coordinated intercellular activity. Furthermore, utilizing a fixed average degree under conditions of low multicellular activity results in exceedingly low thresholds (Rth < 0.5), thereby promoting the establishment of functional connections by chance, which introduces unpredictability into the network analysis. Consequently, techniques that enforce a fixed number of connections are unsuitable for experiments where the level of activity changes significantly.

https://doi.org/10.1371/journal.pcbi.1012130.s001

(TIF)

S2 Fig. Exploring the Impact of Oscillatory Components and Calcium Signal Processing on Functional Network Structure.

The figure presents four types of networks derived from analysis of the five different islets examined in Fig 5: i) A structural network modelled as a geometric network, wherein nearby cells are deemed connected. ii) A functional network derived from unprocessed signals. iii) A functional network extracted from the fast oscillatory component. iv) A functional network constructed based on the slow oscillatory component. All four networks were designed with a fixed average degree kavg = 8. Remarkably, across all five islets, the functional network based on the fast oscillatory component exhibits the fewest long-range connections and shows the highest similarity to the hypothesized structural network. In contrast, networks derived from unprocessed or slow-component signals display a greater proportion of long-range connections, exhibit similar characteristics to each other, and diverge significantly from the structural network.

https://doi.org/10.1371/journal.pcbi.1012130.s002

(TIF)

S3 Fig. Investigating the influence of oscillatory component on functional network structure in an isolated islet.

A) The average unprocessed (black) and extracted slow-component Ca2+ signal (blue) from a Gcamp mouse islet are depicted. The inset shows the corresponding derived fast-component signal (red). B) Different types of beta cell networks: structural (modelled as a geometric network) and three functional networks derived from the unprocessed, slow-component, and fast-component Ca2+ dynamics. Hub cells are highlighted in red. C) Inter-network similarity matrix quantifying the degree of overlap between the four networks. Evidently, the networks extracted from the unprocessed and slow-component traces are very similar, while the fast component network exhibits the highest degree of similarity with the structural network. In contrast, the similarity between the networks derived from unprocessed and slow-component signals and the structural network is notably lower, mirroring observations in tissue slices (see Figs 6 and S2).

https://doi.org/10.1371/journal.pcbi.1012130.s003

(TIF)

S4 Fig. Comparative analysis of functional intercellular network design methods.

A) Three representative beta cell signals (red line) and three alpha cell signals (blue line) subjected to the indicated stimulation protocol: 9 mM -> 10 mM -> 11 mM -> 11 mM glucose + μM epinephrine. This protocol was used to functionally discriminate alpha and beta cells, as the addition of 1 μM epinephrine activates alpha cells and inhibits beta cells. B) Functional networks were extracted using two methods: the fixed average degree method (left) and the four-layered multilayer minimum spanning tree (MST) method (right). The multilayer MST method enforced connections to all cells, including those with asynchronous dynamics, such as alpha cells. Consequently, alpha cells were integrated into the functional network despite their lack of correlation with the rest of the syncytium. This highlights the unsuitability of the MST method for network analyses involving elements with diverse dynamics. Alpha cells are indicated with blue circles and beta cells with red circles.

https://doi.org/10.1371/journal.pcbi.1012130.s004

(TIF)

Acknowledgments

We thank Jasmina Jakopiček, Nika Polšak, Rudi Mlakar, and Maruša Plesnik Rošer for their excellent technical assistance.

References

  1. 1. Haythorne E, Ashcroft FM. Metabolic regulation of insulin secretion in health and disease. The Biochemist. 2021;43: 4–8.
  2. 2. MacDonald PE, Joseph JW, Rorsman P. Glucose-sensing mechanisms in pancreatic beta-cells. Philos Trans R Soc Lond B Biol Sci. 2005;360: 2211–2225. pmid:16321791
  3. 3. Seino S, Sugawara K, Yokoi N, Takahashi H. β-Cell signalling and insulin secretagogues: A path for improved diabetes therapy. Diabetes, Obesity and Metabolism. 2017;19: 22–29. pmid:28880474
  4. 4. Ashcroft FM, Proks P, Smith PA, Ammälä C, Bokvist K, Rorsman P. Stimulus-secretion coupling in pancreatic beta cells. J Cell Biochem. 1994;55 Suppl: 54–65. pmid:7929618
  5. 5. Islam MS. Stimulus-Secretion Coupling in Beta-Cells: From Basic to Bedside. Adv Exp Med Biol. 2020;1131: 943–963. pmid:31646540
  6. 6. Rorsman P, Ashcroft FM. Pancreatic β-Cell Electrical Activity and Insulin Secretion: Of Mice and Men. Physiological Reviews. 2018;98: 117–214. pmid:29212789
  7. 7. Arrojo E Drigo R, Roy B, MacDonald PE. Molecular and functional profiling of human islets: from heterogeneity to human phenotypes. Diabetologia. 2020;63: 2095–2101. pmid:32894320
  8. 8. Benninger RKP, Kravets V. The physiological role of β-cell heterogeneity in pancreatic islet function. Nat Rev Endocrinol. 2022;18: 9–22. pmid:34667280
  9. 9. Skelin Klemen M, Dolenšek J, Slak Rupnik M, Stožer A. The triggering pathway to insulin secretion: Functional similarities and differences between the human and the mouse β cells and their translational relevance. Islets. 2017;9: 109–139. pmid:28662366
  10. 10. Fletcher PA, Marinelli I, Bertram R, Satin LS, Sherman AS. Pulsatile Basal Insulin Secretion Is Driven by Glycolytic Oscillations. Physiology. 2022;37: 216–223. pmid:35378996
  11. 11. Huising MO. Paracrine regulation of insulin secretion. Diabetologia. 2020;63: 2057–2063. pmid:32894316
  12. 12. Benninger RKP, Head WS, Zhang M, Satin LS, Piston DW. Gap junctions and other mechanisms of cell-cell communication regulate basal insulin secretion in the pancreatic islet. The Journal of physiology. 2011;589: 5453–66. pmid:21930600
  13. 13. Šterk M, Dolenšek J, Skelin Klemen M, Križančić Bombek L, Paradiž Leitgeb E, Kerčmar J, et al. Functional characteristics of hub and wave-initiator cells in β cell networks. Biophys J. 2023;122: 784–801. pmid:36738106
  14. 14. Benninger RKP, Hutchens T, Head WS, McCaughey MJ, Zhang M, Le Marchand SJ, et al. Intrinsic Islet Heterogeneity and Gap Junction Coupling Determine Spatiotemporal Ca2+ Wave Dynamics. Biophysical Journal. 2014;107: 2723–2733. pmid:25468351
  15. 15. Moullé VS. Autonomic control of pancreatic beta cells: What is known on the regulation of insulin secretion and beta-cell proliferation in rodents and humans. Peptides. 2022;148: 170709. pmid:34896576
  16. 16. Stožer A, Dolenšek J, Rupnik MS. Glucose-Stimulated Calcium Dynamics in Islets of Langerhans in Acute Mouse Pancreas Tissue Slices. Maedler K, editor. PLoS ONE. 2013;8: e54638. pmid:23358454
  17. 17. Ahrén B. Islet nerves in focus—defining their neurobiological and clinical role. Diabetologia. 2012;55: 3152–3154. pmid:23001378
  18. 18. Henquin J-C. Misunderstandings and controversies about the insulin-secreting properties of antidiabetic sulfonylureas. Biochimie. 2017;143: 3–9. pmid:28711685
  19. 19. Farnsworth NL, Walter R, Piscopio RA, Schleicher WE, Benninger RKP. Exendin-4 overcomes cytokine-induced decreases in gap junction coupling via protein kinase A and Epac2 in mouse and human islets. J Physiol. 2019;597: 431–447. pmid:30412665
  20. 20. Ahrén B, Yamada Y, Seino Y. The Insulin Response to Oral Glucose in GIP and GLP-1 Receptor Knockout Mice: Review of the Literature and Stepwise Glucose Dose Response Studies in Female Mice. Front Endocrinol. 2021;12: 665537. pmid:34122340
  21. 21. Miranda MA, Macias-Velasco JF, Lawson HA. Pancreatic β-cell heterogeneity in health and diabetes: classes, sources, and subtypes. Am J Physiol Endocrinol Metab. 2021;320: E716–E731. pmid:33586491
  22. 22. Nasteska D, Hodson DJ. The role of beta cell heterogeneity in islet function and insulin release. J Mol Endocrinol. 2018;61: R43–R60. pmid:29661799
  23. 23. Dominguez-Gutierrez G, Xin Y, Gromada J. Heterogeneity of human pancreatic β-cells. Mol Metab. 2019;27S: S7–S14. pmid:31500834
  24. 24. Pipeleers DG. Heterogeneity in pancreatic beta-cell population. Diabetes. 1992;41: 777–781. pmid:1612191
  25. 25. Katsuta H, Aguayo-Mazzucato C, Katsuta R, Akashi T, Hollister-Lock J, Sharma AJ, et al. Subpopulations of GFP-marked mouse pancreatic β-cells differ in size, granularity, and insulin secretion. Endocrinology. 2012;153: 5180–5187. pmid:22919061
  26. 26. Wojtusciszyn A, Armanet M, Morel P, Berney T, Bosco D. Insulin secretion from human beta cells is heterogeneous and dependent on cell-to-cell contacts. Diabetologia. 2008;51: 1843–1852. pmid:18665347
  27. 27. Stožer A, Skelin Klemen M, Gosak M, Križančić Bombek L, Pohorec V, Slak Rupnik M, et al. Glucose-dependent activation, activity, and deactivation of beta cell networks in acute mouse pancreas tissue slices. American Journal of Physiology-Endocrinology and Metabolism. 2021;321: E305–E323. pmid:34280052
  28. 28. Kravets V, Dwulet JM, Schleicher WE, Hodson DJ, Davis AM, Pyle L, et al. Functional architecture of the pancreatic islets reveals first responder cells which drive the first-phase [Ca2+] response. bioRxiv. 2021; 2020.12.22.424082. https://doi.org/10.1101/2020.12.22.424082
  29. 29. Gosak M, Stožer A, Markovič R, Dolenšek J, Perc M, Rupnik MS, et al. Critical and Supercritical Spatiotemporal Calcium Dynamics in Beta Cells. Frontiers in Physiology. 2017;8: 1106. pmid:29312008
  30. 30. Pedersen MG, Tagliavini A, Henquin J-C. Calcium signaling and secretory granule pool dynamics underlie biphasic insulin secretion and its amplification by glucose: experiments and modeling. Am J Physiol Endocrinol Metab. 2019;316: E475–E486. pmid:30620637
  31. 31. Scarl RT, Corbin KL, Vann NW, Smith HM, Satin LS, Sherman A, et al. Intact pancreatic islets and dispersed beta-cells both generate intracellular calcium oscillations but differ in their responsiveness to glucose. Cell Calcium. 2019;83: 102081. pmid:31563790
  32. 32. Curry DL, Bennett LL, Grodsky GM. Dynamics of insulin secretion by the perfused rat pancreas. Endocrinology. 1968;83: 572–84. pmid:4877098
  33. 33. Pohorec V, Križančić Bombek L, Skelin Klemen M, Dolenšek J, Stožer A. Glucose-Stimulated Calcium Dynamics in Beta Cells From Male C57BL/6J, C57BL/6N, and NMRI Mice: A Comparison of Activation, Activity, and Deactivation Properties in Tissue Slices. Frontiers in Endocrinology. 2022;13: 1–16. pmid:35399951
  34. 34. Benninger RKP, Zhang M, Head WS, Satin LS, Piston DW, Steven Head W, et al. Gap Junction Coupling and Calcium Waves in the Pancreatic Islet. Biophysical Journal. 2008;95: 5048–5061. pmid:18805925
  35. 35. Šterk M, Dolenšek J, Bombek LK, Markovič R, Zakelšek D, Perc M, et al. Assessing the origin and velocity of Ca2+ waves in three-dimensional tissue: Insights from a mathematical model and confocal imaging in mouse pancreas tissue slices. Communications in Nonlinear Science and Numerical Simulation. 2021;93: 105495.
  36. 36. Dolenšek J, Špelič D, Klemen M, Žalik B, Gosak M, Rupnik M, et al. Membrane Potential and Calcium Dynamics in Beta Cells from Mouse Pancreas Tissue Slices: Theory, Experimentation, and Analysis. Sensors. 2015;15: 27393–27419. pmid:26516866
  37. 37. Speier S, Gjinovci A, Charollais A, Meda P, Rupnik M. Cx36-mediated coupling reduces beta-cell heterogeneity, confines the stimulating glucose concentration range, and affects insulin release kinetics. Diabetes. 2007;56: 1078–86. pmid:17395748
  38. 38. Meda P. Gap junction proteins are key drivers of endocrine function. Biochim Biophys Acta Biomembr. 2018;1860: 124–140. pmid:28284720
  39. 39. Cigliola V, Chellakudam V, Arabieter W, Meda P. Connexins and β-cell functions. Diabetes Res Clin Pract. 2013;99: 250–259. pmid:23176806
  40. 40. Head WS, Orseth ML, Nunemaker CS, Satin LS, Piston DW, Benninger RKP. Connexin-36 gap junctions regulate in vivo first- and second-phase insulin secretion dynamics and glucose tolerance in the conscious mouse. Diabetes. 2012;61: 1700–1707. pmid:22511206
  41. 41. Westacott MJ, Ludin NWF, Benninger RKP. Spatially Organized β-Cell Subpopulations Control Electrical Dynamics across Islets of Langerhans. Biophysical Journal. 2017;113: 1093–1108. pmid:28877492
  42. 42. Yang YHC, Briant LJ, Raab CA, Mullapudi ST, Maischein H-M, Kawakami K, et al. Innervation modulates the functional connectivity between pancreatic endocrine cells. Sussel L, Bronner ME, Knox SM, Meyer D, editors. eLife. 2022;11: e64526. pmid:35373736
  43. 43. Hartig SM, Cox AR. Paracrine signaling in islet function and survival. J Mol Med (Berl). 2020;98: 451–467. pmid:32067063
  44. 44. Holter MM, Saikia M, Cummings BP. Alpha-cell paracrine signaling in the regulation of beta-cell insulin secretion. Frontiers in Endocrinology. 2022;13. Available from: https://www.frontiersin.org/articles/10.3389/fendo.2022.934775. pmid:35957816
  45. 45. Rutter GA, Hodson DJ. Beta cell connectivity in pancreatic islets: a type 2 diabetes target? Cell Mol Life Sci. 2015;72: 453–467. pmid:25323131
  46. 46. Pizarro-Delgado J, Fasciani I, Temperan A, Romero M, Gonzalez-Nieto D, Alonso-Magdalena P, et al. Inhibition of connexin 36 hemichannels by glucose contributes to the stimulation of insulin secretion. AJP: Endocrinology and Metabolism. 2014;306: E1354–E1366. pmid:24735890
  47. 47. Martinez C, Maschio DA, de Fontes CC, Vanzela EC, Benfato ID, Gazarini ML, et al. Early decrease in Cx36 is associated with increased cell adhesion molecules (CAMs) junctional content in mouse pancreatic islets after short-term high-fat diet feeding. Annals of Anatomy—Anatomischer Anzeiger. 2022;241: 151891. pmid:35114378
  48. 48. Carvalho CPF, Oliveira RB, Britan A, Santos-Silva JC, Boschero AC, Meda P, et al. Impaired β-cell-β-cell coupling mediated by Cx36 gap junctions in prediabetic mice. American Journal of Physiology-Endocrinology and Metabolism. 2012;303: E144–E151. pmid:22569071
  49. 49. Corezola do Amaral ME, Kravets V, Dwulet JM, Farnsworth NL, Piscopio R, Schleicher WE, et al. Caloric restriction recovers impaired β-cell-β-cell gap junction coupling, calcium oscillation coordination, and insulin secretion in prediabetic mice. American Journal of Physiology-Endocrinology and Metabolism. 2020;319: E709–E720. pmid:32830549
  50. 50. Clair JRS, Westacott MJ, Miranda J, Farnsworth NL, Kravets V, Schleicher WE, et al. Restoring Connexin-36 Function in Diabetogenic Environments Precludes Mouse and Human Islet Dysfunction. bioRxiv; 2023. p. 2020.11.03.366179. https://doi.org/10.1101/2020.11.03.366179
  51. 51. Saadati M, Jamali Y. The effects of beta-cell mass and function, intercellular coupling, and islet synchrony on [Formula: see text] dynamics. Sci Rep. 2021;11: 10268. pmid:33986325
  52. 52. Ravier MA, Güldenagel M, Charollais A, Gjinovci A, Caille D, Söhl G, et al. Loss of Connexin36 Channels Alters -Cell Coupling, Islet Synchronization of Glucose-Induced Ca2+ and Insulin Oscillations, and Basal Insulin Release. Diabetes. 2005;54: 1798–1807. pmid:15919802
  53. 53. Peercy BE, Sherman AS. Do oscillations in pancreatic islets require pacemaker cells? J Biosci. 2022;47: 14. pmid:35212266
  54. 54. Johnston NR, Mitchell RK, Haythorne E, Pessoa MP, Semplici F, Ferrer J, et al. Beta Cell Hubs Dictate Pancreatic Islet Responses to Glucose. Cell Metabolism. 2016;24: 389–401. pmid:27452146
  55. 55. Salem V, Silva LD, Suba K, Georgiadou E, Neda Mousavy Gharavy S, Akhtar N, et al. Leader β-cells coordinate Ca2+ dynamics across pancreatic islets in vivo. Nature Metabolism. 2019;1: 615–629. pmid:32694805
  56. 56. Stožer A, Šterk M, Paradiž Leitgeb E, Markovič R, Skelin Klemen M, Ellis CE, et al. From Isles of Königsberg to Islets of Langerhans: Examining the Function of the Endocrine Pancreas Through Network Science. Front Endocrinol (Lausanne). 2022;13: 922640. pmid:35784543
  57. 57. Kravets V, Dwulet JM, Schleicher WE, Hodson DJ, Davis AM, Pyle L, et al. Functional architecture of pancreatic islets identifies a population of first responder cells that drive the first-phase calcium response. PLoS Biol. 2022;20: e3001761. pmid:36099294
  58. 58. Félix-Martínez GJ, Godínez-Fernández JR. A primer on modelling pancreatic islets: from models of coupled β-cells to multicellular islet models. Islets. 2023;15: 2231609. pmid:37415423
  59. 59. Lei C-L, Kellard JA, Hara M, Johnson JD, Rodriguez B, Briant LJB. Beta-cell hubs maintain Ca 2+ oscillations in human and mouse islet simulations. Islets. 2018;10: 151–167. pmid:30142036
  60. 60. Loppini A, Chiodo L. Biophysical modeling of β-cells networks: Realistic architectures and heterogeneity effects. Biophysical Chemistry. 2019;254: 106247. pmid:31472460
  61. 61. Gosak M, Yan-Do R, Lin H, MacDonald PE, Stožer A. Ca2+ Oscillations, Waves, and Networks in Islets From Human Donors With and Without Type 2 Diabetes. Diabetes. 2022;71: 2584–2596. pmid:36084321
  62. 62. Aslanidi OV, Mornev OA, Skyggebjerg O, Arkhammar P, Thastrup O, Sørensen MP, et al. Excitation Wave Propagation as a Possible Mechanism for Signal Transmission in Pancreatic Islets of Langerhans. Biophysical Journal. 2001;80: 1195–1209. pmid:11222284
  63. 63. Chabosseau P, Yong F, Delgadillo-Silva LF, Lee EY, Melhem R, Li S, et al. Molecular phenotyping of single pancreatic islet leader beta cells by “Flash-Seq.” Life Sci. 2023;316: 121436. pmid:36706832
  64. 64. Zhang Q, Huising MO, Da Silva Xavier G, Hauge-Evans AC. Editorial: The pancreatic islet–a multifaceted hub of inter-cellular communication. Front Endocrinol. 2023;14. pmid:37265699
  65. 65. Satin LS, Zhang Q, Rorsman P. “Take Me To Your Leader”: An Electrophysiological Appraisal of the Role of Hub Cells in Pancreatic Islets. Diabetes. 2020;69: 830–836. pmid:32312899
  66. 66. Rutter GA, Ninov N, Salem V, Hodson DJ. Comment on Satin et al. “Take Me To Your Leader”: An Electrophysiological Appraisal of the Role of Hub Cells in Pancreatic Islets. Diabetes 2020;69:830–836. Diabetes. 2020;69: e10–e11. pmid:32312899
  67. 67. Satin LS, Rorsman P. Response to Comment on Satin et al. “Take Me To Your Leader”: An Electrophysiological Appraisal of the Role of Hub Cells in Pancreatic Islets. Diabetes 2020;69:830–836. Diabetes. 2020;69: e12–e13. pmid:32820057
  68. 68. Cappon G, Pedersen MG. Heterogeneity and nearest-neighbor coupling can explain small-worldness and wave properties in pancreatic islets. Chaos: An Interdisciplinary Journal of Nonlinear Science. 2016;26: 053103. pmid:27249943
  69. 69. Gosak M, Markovič R, Dolenšek J, Slak Rupnik M, Marhl M, Stožer A, et al. Network science of biological systems at different scales: A review. Physics of Life Reviews. 2018;24: 118–135. pmid:29150402
  70. 70. Briggs JK, Gresch A, Marinelli I, Dwulet JM, Albers DJ, Kravets V, et al. Beta-cell intrinsic dynamics rather than gap junction structure dictates subpopulations in the islet functional network. Satin LS, editor. eLife. 2023;12: e83147. pmid:38018905
  71. 71. Paradiž Leitgeb E, Kerčmar J, Križančić Bombek L, Pohorec V, Skelin Klemen M, Slak Rupnik M, et al. Exendin-4 affects calcium signalling predominantly during activation and activity of beta cell networks in acute mouse pancreas tissue slices. Front Endocrinol. 2024;14. pmid:38292770
  72. 72. Markovič R, Stožer A, Gosak M, Dolenšek J, Marhl M, Rupnik MS. Progressive glucose stimulation of islet beta cells reveals a transition from segregated to integrated modular functional connectivity patterns. Scientific Reports. 2015;5: 7845. pmid:25598507
  73. 73. Berglund ED, Li CY, Poffenberger G, Ayala JE, Fueger PT, Willis SE, et al. Glucose Metabolism In Vivo in Four Commonly Used Inbred Mouse Strains. Diabetes. 2008;57: 1790–1799. pmid:18398139
  74. 74. Hull RL, Willard JR, Struck MD, Barrow BM, Brar GS, Andrikopoulos S, et al. High fat feeding unmasks variable insulin responses in male C57BL/6 mouse substrains. Journal of Endocrinology. 2017;233: 53–64. pmid:28138002
  75. 75. Fontaine DA, Davis DB. Attention to Background Strain Is Essential for Metabolic Research: C57BL/6 and the International Knockout Mouse Consortium. Diabetes. 2016;65: 25–33. pmid:26696638
  76. 76. Emfinger CH, Clark LE, Yandell B, Schueler KL, Simonett SP, Stapleton DS, et al. Novel regulators of islet function identified from genetic variation in mouse islet Ca2+ oscillations. eLife. 2023;12. pmid:37787501
  77. 77. Marinelli I, Fletcher PA, Sherman AS, Satin LS, Bertram R. Symbiosis of Electrical and Metabolic Oscillations in Pancreatic β-Cells. Front Physiol. 2021;12: 781581. pmid:34925070
  78. 78. Šterk M, Barać U, Stožer A, Gosak M. Both electrical and metabolic coupling shape the collective multimodal activity and functional connectivity patterns in beta cell collectives: A computational model perspective. Phys Rev E. 2023;108: 054409. pmid:38115462
  79. 79. Luchetti N, Filippi S, Loppini A. Multilevel synchronization of human β-cells networks. Front Netw Physiol. 2023;3: 1264395. pmid:37808419
  80. 80. Marinelli I, Vo T, Gerardo-Giorda L, Bertram R. Transitions between bursting modes in the integrated oscillator model for pancreatic β-cells. J Theor Biol. 2018;454: 310–319. pmid:29935201
  81. 81. Bertram R, Satin LS, Sherman AS. Closing in on the Mechanisms of Pulsatile Insulin Secretion. Diabetes. 2018;67: 351–359. pmid:29463575
  82. 82. Wijk BCM van Stam CJ, Daffertshofer A. Comparing Brain Networks of Different Size and Connectivity Density Using Graph Theory. PLOS ONE. 2010;5: e13701. pmid:21060892
  83. 83. De Vico Fallani F, Richiardi J, Chavez M, Achard S. Graph analysis of functional brain networks: practical issues in translational neuroscience. Philosophical Transactions of the Royal Society B: Biological Sciences. 2014;369: 20130521. pmid:25180301
  84. 84. Theis N, Rubin J, Cape J, Iyengar S, Prasad KM. Threshold Selection for Brain Connectomes. Brain Connect. 2023;13: 383–393. pmid:37166374
  85. 85. Langer N, Pedroni A, Jäncke L. The Problem of Thresholding in Small-World Network Analysis. PLOS ONE. 2013;8: e53199. pmid:23301043
  86. 86. van den Heuvel MP, de Lange SC, Zalesky A, Seguin C, Yeo BTT, Schmidt R. Proportional thresholding in resting-state fMRI functional connectivity networks and consequences for patient-control connectome studies: Issues and recommendations. Neuroimage. 2017;152: 437–449. pmid:28167349
  87. 87. Joudaki A, Salehi N, Jalili M, Knyazeva MG. EEG-based functional brain networks: does the network size matter? PLoS One. 2012;7: e35673. pmid:22558196
  88. 88. Hassan M, Dufor O, Merlet I, Berrou C, Wendling F. EEG source connectivity analysis: from dense array recordings to brain networks. PLoS One. 2014;9: e105041. pmid:25115932
  89. 89. De Domenico M. Multilayer modeling and analysis of human brain networks. Gigascience. 2017;6: 1–8. pmid:28327916
  90. 90. Kajimura S, Margulies D, Smallwood J. Frequency-specific brain network architecture in resting-state fMRI. Sci Rep. 2023;13: 2964. pmid:36806195
  91. 91. Skelin Klemen M, Dolenšek J, Križančić Bombek L, Pohorec V, Gosak M, Slak Rupnik M, et al. The effect of forskolin and the role of Epac2A during activation, activity, and deactivation of beta cell networks. Frontiers in Endocrinology. 2023;14. pmid:37701894
  92. 92. Gilon P, Henquin JC. Influence of membrane potential changes on cytoplasmic Ca2+ concentration in an electrically excitable cell, the insulin-secreting pancreatic B-cell. Journal of Biological Chemistry. 1992;267: 20713–20720. pmid:1400388
  93. 93. Gilon P, Shepherd RM, Henquin JC. Oscillations of secretion driven by oscillations of cytoplasmic Ca2+ as evidences in single pancreatic islets. The Journal of biological chemistry. 1993;268: 22265–8. pmid:8226733
  94. 94. Satin LS, Butler PC, Ha J, Sherman AS. Pulsatile insulin secretion, impaired glucose tolerance and type 2 diabetes. Molecular Aspects of Medicine. 2015;42: 61–77. pmid:25637831
  95. 95. Henquin JC, Meissner HP, Schmeer W. Cyclic variations of glucose-induced electrical activity in pancreatic B cells. Pflugers Arch. 1982;393: 322–327. pmid:6750552
  96. 96. Zmazek J, Klemen MS, Markovič R, Dolenšek J, Marhl M, Stožer A, et al. Assessing Different Temporal Scales of Calcium Dynamics in Networks of Beta Cell Populations. Frontiers in Physiology. 2021;12: 337. pmid:33833686
  97. 97. Bergsten P. Role of oscillations in membrane potential, cytoplasmic Ca2+, and metabolism for plasma insulin oscillations. Diabetes. American Diabetes Association Inc.; 2002. p. 171. pmid:11815477
  98. 98. Nunemaker CS, Bertram R, Sherman A, Tsaneva-Atanasova K, Daniel CR, Satin LS. Glucose modulates [Ca2+]i oscillations in pancreatic islets via ionic and glycolytic mechanisms. Biophysical Journal. 2006;91: 2082–2096. pmid:16815907
  99. 99. Gilon P, Jonas JC, Henquin JC. Culture duration and conditions affect the oscillations of cytoplasmic calcium concentration induced by glucose in mouse pancreatic islets. Diabetologia. 1994;37: 1007–1014. pmid:7851679
  100. 100. Zhang M, Goforth P, Bertram R, Sherman A, Satin L. The Ca2+ Dynamics of Isolated Mouse β-Cells and Islets: Implications for Mathematical Models. Biophysical Journal. 2003;84: 2852–2870. pmid:12719219
  101. 101. Jonkers FC, Jonas J-C, Gilon P, Henquin J-C. Influence of cell number on the characteristics and synchrony of Ca 2+ oscillations in clusters of mouse pancreatic islet cells. The Journal of Physiology. 1999;520: 839–849. pmid:10545148
  102. 102. Emfinger CH, Lőrincz R, Wang Y, York NW, Singareddy SS, Ikle JM, et al. Beta-cell excitability and excitability-driven diabetes in adult Zebrafish islets. Physiol Rep. 2019;7: e14101. pmid:31161721
  103. 103. Valdeolmillos M, Santos RM, Contreras D, Soria B, Rosario LM. Glucose-induced oscillations of intracellular Ca2+ concentration resembling bursting electrical activity in single mouse islets of Langerhans. FEBS Letters. 1989;259: 19–23. pmid:2689228
  104. 104. Manning Fox JE, Gyulkhandanyan AV, Satin LS, Wheeler MB. Oscillatory Membrane Potential Response to Glucose in Islet β-Cells: A Comparison of Islet-Cell Electrical Activity in Mouse and Rat. Endocrinology. 2006;147: 4655–4663. pmid:16857746
  105. 105. Antunes CM, Salgado AP, Rosário LM, Santos RM. Differential patterns of glucose-induced electrical activity and intracellular calcium responses in single mouse and rat pancreatic islets. Diabetes. 2000;49: 2028–2038. pmid:11118004
  106. 106. Zimliki CL, Chenault VM, Mears D. Glucose-dependent and -independent electrical activity in islets of Langerhans of Psammomys obesus, an animal model of nutritionally induced obesity and diabetes. Gen Comp Endocrinol. 2009;161: 193–201. pmid:19167400
  107. 107. Pertusa JAG, Nesher R, Kaiser N, Cerasi E, Henquin J-C, Jonas J-C. Increased glucose sensitivity of stimulus-secretion coupling in islets from Psammomys obesus after diet induction of diabetes. Diabetes. 2002;51: 2552–2560. pmid:12145170
  108. 108. Bertuzzi F, Zacchetti D, Berra C, Socci C, Pozza G, Pontiroli AE, et al. Intercellular Ca2+ waves sustain coordinate insulin secretion in pig islets of Langerhans. FEBS Lett. 1996;379: 21–25. pmid:8566222
  109. 109. Bertuzzi F, Garancini P, Socci TC, Nano R, Taglietti MV, Santopinto M, et al. Lessons from in vitro perifusion of pancreatic islets isolated from 80 human pancreases. Cell Transplant. 1999;8: 709–712. pmid:10701499
  110. 110. Martin F, Soria B. [Ca*+], oscillations pancreatic islets. Cell Calcium. 1996; 20: 409.
  111. 111. Dolenšek J, Stožer A, Skelin Klemen M, Miller EW, Slak Rupnik M. The Relationship between Membrane Potential and Calcium Dynamics in Glucose-Stimulated Beta Cell Syncytium in Acute Mouse Pancreas Tissue Slices. Rakonczay Z, editor. PLoS ONE. 2013;8: e82374. pmid:24324777
  112. 112. Bolea S, Pertusa JAG, Martín F, Sanchez-Andrés JV, Soria B. Regulation of pancreatic β-cell electrical activity and insulin release by physiological amino acid concentrations. Pflügers Arch. 1997;433: 699–704. pmid:9049159
  113. 113. Jaffredo M, Bertin E, Pirog A, Puginier E, Gaitan J, Oucherif S, et al. Dynamic Uni- and Multicellular Patterns Encode Biphasic Activity in Pancreatic Islets. Diabetes. 2021;70: 878–888. pmid:33468514
  114. 114. Speier S, Rupnik M. A novel approach to in situ characterization of pancreatic β-cells. Pflügers Archiv. 2003;446: 553–558. pmid:12774232
  115. 115. Irving-Rodgers HF, Choong FJ, Hummitzsch K, Parish CR, Rodgers RJ, Simeonovic CJ. Pancreatic islet basement membrane loss and remodeling after mouse islet isolation and transplantation: impact for allograft rejection. Cell Transplant. 2014;23: 59–72. pmid:23211522
  116. 116. Roe MW, Spencer B, Lancaster ME, Mertz RJ, Worley JF, Dukes ID. Absence of effect of culture duration on glucose-activated alterations in intracellular calcium concentration in mouse pancreatic islets. Diabetologia. 1995;38: 876–877. pmid:7556996
  117. 117. Campbell SA, Johnson J, Light PE. Evidence for the existence and potential roles of intra-islet glucagon-like peptide-1. Islets. 2021;13: 32–50. pmid:33724156
  118. 118. Liu Y-J, Tengholm A, Grapengiesser E, Hellman B, Gylfe E. Origin of slow and fast oscillations of Ca 2+ in mouse pancreatic islets. The Journal of Physiology. 1998;508: 471–481. pmid:9508810
  119. 119. Zaborska KE, Jordan KL, Thorson AS, Dadi PK, Schaub CM, Nakhe AY, et al. Liraglutide increases islet Ca2+ oscillation frequency and insulin secretion by activating hyperpolarization-activated cyclic nucleotide-gated channels. Diabetes Obes Metab. 2022;24: 1741–1752. pmid:35546791
  120. 120. Ahn YB, Xu G, Marselli L, Toschi E, Sharma A, Bonner-Weir S, et al. Changes in gene expression in beta cells after islet isolation and transplantation using laser-capture microdissection. Diabetologia. 2007;50: 334–342. pmid:17180350
  121. 121. Negi S, Jetha A, Aikin R, Hasilo C, Sladek R, Paraskevas S. Analysis of Beta-Cell Gene Expression Reveals Inflammatory Signaling and Evidence of Dedifferentiation following Human Islet Isolation and Culture. PLOS ONE. 2012;7: e30415. pmid:22299040
  122. 122. Stožer A, Dolenšek J, Križančić Bombek L, Pohorec V, Slak Rupnik M, Klemen MS. Confocal Laser Scanning Microscopy of Calcium Dynamics in Acute Mouse Pancreatic Tissue Slices. Journal of Visualized Experiments. 2021; e62293. pmid:33938876
  123. 123. Marciniak A, Cohrs CM, Tsata V, Chouinard JA, Selck C, Stertmann J, et al. Using pancreas tissue slices for in situ studies of islet of Langerhans and acinar cell biology. Nat Protoc. 2014;9: 2809–2822. pmid:25393778
  124. 124. Santos RM, Rosario LM, Nadal A, Garcia-Sancho J, Soria B, Valdeolmillos M. Widespread synchronous [Ca2+]i oscillations due to bursting electrical activity in single pancreatic islets. Pflügers Archiv—European Journal of Physiology. 1991;418: 417–22. pmid:1876486
  125. 125. Sánchez-Andrés JV, Gomis A, Valdeolmillos M. The electrical activity of mouse pancreatic beta-cells recorded in vivo shows glucose-dependent oscillations. The Journal of Physiology. 1995;486: 223–228. pmid:7562637
  126. 126. Fernandez J, Valdeolmillos M. Synchronous glucose-dependent [Ca 2+] i oscillations in mouse pancreatic islets of Langerhans recorded in vivo. FEBS Letters. 2000;477: 33–36. pmid:10899306
  127. 127. Chen C, Chmelova H, Cohrs CM, Chouinard JA, Jahn SR, Stertmann J, et al. Alterations in β-Cell Calcium Dynamics and Efficacy Outweigh Islet Mass Adaptation in Compensation of Insulin Resistance and Prediabetes Onset. Diabetes. 2016;65: 2676–2685. pmid:27207518
  128. 128. Cook DL, Ikeuchi M. Tolbutamide as Mimic of Glucose on β-Cell Electrical Activity: ATP-Sensitive K+ Channels as Common Pathway for Both Stimuli. Diabetes. 1989;38: 416–421. pmid:2647550
  129. 129. Meissner HP, Schmelz H. Membrane potential of beta-cells in pancreatic islets. Pflügers Archiv—European Journal of Physiology. 1974;351: 195–206. pmid:4608967
  130. 130. Henquin JC. Adenosine triphosphate-sensitive K+ channels may not be the sole regulators of glucose-induced electrical activity in pancreatic B-cells. Endocrinology. 1992;131: 127–131. pmid:1611991
  131. 131. Schröter M, Paulsen O, Bullmore ET. Micro-connectomics: probing the organization of neuronal networks at the cellular scale. Nat Rev Neurosci. 2017;18: 131–146. pmid:28148956
  132. 132. Hodson DJ, Schaeffer M, Romanò N, Fontanaud P, Lafont C, Birkenstock J, et al. Existence of long-lasting experience-dependent plasticity in endocrine cell networks. Nat Commun. 2012;3: 605. pmid:22215080
  133. 133. Fazli M, Bertram R. Network Properties of Electrically Coupled Bursting Pituitary Cells. Front Endocrinol (Lausanne). 2022;13: 936160. pmid:35872987
  134. 134. Pires M, Raischel F, Vaz SH, Cruz-Silva A, Sebastião AM, Lind PG. Modeling the functional network of primary intercellular Ca2+ wave propagation in astrocytes and its application to study drug effects. J Theor Biol. 2014;356: 201–212. pmid:24813072
  135. 135. Mojica-Benavides M, van Niekerk DD, Mijalkov M, Snoep JL, Mehlig B, Volpe G, et al. Intercellular communication induces glycolytic synchronization waves between individually oscillating cells. Proc Natl Acad Sci U S A. 2021;118: e2010075118. pmid:33526662
  136. 136. Gosak M, Markovič R, Fajmut A, Marhl M, Hawlina M, Andjelić S. The Analysis of Intracellular and Intercellular Calcium Signaling in Human Anterior Lens Capsule Epithelial Cells with Regard to Different Types and Stages of the Cataract. PLoS One. 2015;10: e0143781. pmid:26636768
  137. 137. Stevenson AJ, Vanwalleghem G, Stewart TA, Condon ND, Lloyd-Lewis B, Marino N, et al. Multiscale imaging of basal cell dynamics in the functionally mature mammary gland. Proc Natl Acad Sci U S A. 2020;117: 26822–26832. pmid:33033227
  138. 138. Marolt U, Paradiž Leitgeb E, Pohorec V, Lipovšek S, Venglovecz V, Gál E, et al. Calcium imaging in intact mouse acinar cells in acute pancreas tissue slices. PLoS One. 2022;17: e0268644. pmid:35657915
  139. 139. Verma A, Antony AN, Ogunnaike BA, Hoek JB, Vadigepalli R. Causality Analysis and Cell Network Modeling of Spatial Calcium Signaling Patterns in Liver Lobules. Front Physiol. 2018;9: 1377. pmid:30337879
  140. 140. Boccaletti S, Bianconi G, Criado R, Del Genio CI, Gómez-Gardeñes J, Romance M, et al. The structure and dynamics of multilayer networks. Phys Rep. 2014;544: 1–122. pmid:32834429
  141. 141. Kivelä M, Arenas A, Barthelemy M, Gleeson JP, Moreno Y, Porter MA. Multilayer networks. Journal of Complex Networks. 2014;2: 203–271.
  142. 142. Battiston F, Cencetti G, Iacopini I, Latora V, Lucas M, Patania A, et al. Networks beyond pairwise interactions: Structure and dynamics. Physics Reports. 2020;874: 1–92.
  143. 143. Majhi S, Perc M, Ghosh D. Dynamics on higher-order networks: a review. Journal of The Royal Society Interface. 2022;19: 20220043. pmid:35317647
  144. 144. Porta A, Faes L. Wiener–Granger Causality in Network Physiology With Applications to Cardiovascular Control and Neuroscience. Proceedings of the IEEE. 2016;104: 282–309.
  145. 145. Reid AT, Headley DB, Mill RD, Sanchez-Romero R, Uddin LQ, Marinazzo D, et al. Advancing functional connectivity research from association to causation. Nat Neurosci. 2019;22: 1751–1760. pmid:31611705
  146. 146. Riaz A, Asad M, Alonso E, Slabaugh G. DeepFMRI: End-to-end deep learning for functional connectivity and classification of ADHD using fMRI. J Neurosci Methods. 2020;335: 108506. pmid:32001294
  147. 147. Ha S, Jeong H. Unraveling hidden interactions in complex systems with deep learning. Sci Rep. 2021;11: 12804. pmid:34140551
  148. 148. Choi HJ, Wang C, Pan X, Jang J, Cao M, Brazzo JA, et al. Emerging machine learning approaches to phenotyping cellular motility and morphodynamics. Phys Biol. 2021;18: 041001. pmid:33971636
  149. 149. Printz Y, Patil P, Mahn M, Benjamin A, Litvin A, Levy R, et al. Determinants of functional synaptic connectivity among amygdala-projecting prefrontal cortical neurons in male mice. Nat Commun. 2023; 14: 1667.
  150. 150. Li G, LeFebre R, Starman A, Chappell P, Mugler A, Sun B. Temporal signals drive the emergence of multicellular information networks. Proceedings of the National Academy of Sciences. 2022;119: e2202204119. pmid:36067282
  151. 151. Schindelin J, Arganda-Carreras I, Frise E, Kaynig V, Longair M, Pietzsch T, et al. Fiji: an open-source platform for biological-image analysis. Nat Methods. 2012;9: 676–682. pmid:22743772
  152. 152. Fletcher PA, Thompson B, Liu C, Bertram R, Satin LS, Sherman AS. Ca2+ release or Ca2+ entry, that is the question: what governs Ca2+ oscillations in pancreatic β cells? Am J Physiol Endocrinol Metab. 2023;324: E477–E487. pmid:37074988
  153. 153. Gava GP, McHugh SB, Lefèvre L, Lopes-dos-Santos V, Trouche S, El-Gaby M, et al. Integrating new memories into the hippocampal network activity space. Nat Neurosci. 2021;24: 326–330. pmid:33603228
  154. 154. Mijatovic G, Loncar-Turukalo T, Bozanic N, Milosavljevic N, Storchi R, Faes L. A Measure of Concurrent Neural Firing Activity Based on Mutual Information. Neuroinformatics. 2021;19: 719–735. pmid:33852134
  155. 155. Cover TM, Thomas JA. Elements of Information Theory. 1st ed. Wiley; 2001.
  156. 156. Zhang Q, Galvanovskis J, Abdulkader F, Partridge CJ, Göpel SO, Eliasson L, et al. Cell coupling in mouse pancreatic -Cells measured in intact islets of Langerhans. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences. 2008;366: 3503–3523. pmid:18632454
  157. 157. Kruskal JB. On the shortest spanning subtree of a graph and the traveling salesman problem. Proceedings of the American Mathematical Society. 1956. pp. 48–50.
  158. 158. Prim RC. Shortest connection networks and some generalizations. The Bell System Technical Journal. 1957;36: 1389–1401.
  159. 159. Boccaletti S, Latora V, Moreno Y, Chavez M, Hwang D. Complex networks: Structure and dynamics. Physics Reports. 2006;424: 175–308.