We built our biophysical model of human β-cell networks based on a comprehensive electrophysiological description of human β-cells activity, also accounting for metabolites oscillations (Pedersen, 2010; Riz et al., 2014; Loppini et al., 2015), and added gap-junction coupling between adjacent cells, defined by the islet structure as detailed in the next subsection.

The equation for the membrane potential dynamics of the ith cell is d V i d t = − I ion, i + ∑ j ∈ Ω i I coup, j i , (1)

where I _ion,i is sum of the membrane ionic currents, Ω _i is the cell neighborhood and I _coup,ji = g _c(V _i − V _j ) is the gap-junction current between i ^th and j ^th cells.

Among the ionic currents included in the electrical compartment, the ATP-sensitive potassium (KATP) current is modulated by a variable mimicking the slow periodic ATP increasing and decreasing downstream of the glycolytic oscillations, triggering action potential firing and bursting. The glycolytic compartment, models the dynamics of intermediate metabolites whose production follows the glucose uptake by the cell, i.e., glucose 6-phosphate (G6P), fructose 6-phosphate (F6P), fructose 1,6-bisphosphate (FBP), dihydroxyacetone-phosphate (DHAP) and glyceraldehyde 3-phosphate (G3P). In our model, the metabolic coupling between nearby cells is realized imposing G6P and F6P diffusion through gap junctions. In particular, the dynamical equation regulating the total concentration of G6P and F6P (G6PF6P) for the i ^th cell is d [ G 6 P F 6 P ] i d t = V GK , i − V PFK − P G 6 PF 6 P ∑ j ∈ Ω i G 6 P F 6 P i − G 6 P F 6 P j , (2) where V _GK and V _PFK are the glucokinase and phosphofructokinase reaction rates, and P _G6PF6P denotes the gap-junction permeability to metabolites, i.e., the strength of metabolic coupling. As reference case, we set the electrical conductance and metabolites permeability of gap junctions to g _c = 0.01 nS/pF and P _G6PF6p = 0.01 ms⁻¹, in line with previous estimations based on electrophysiological data and experimental measures of metabolites diffusion among adjacent β-cells (Kohen et al., 1979; Loppini et al., 2015). These values are also used in other modeling studies on β-cells networks (Loppini et al., 2015; Loppini et al., 2017; Loppini and Chiodo, 2019; Saadati and Jamali, 2021). Biological heterogeneity among cells was reproduced by varying the conductance of the voltage-sensitive potassium channels (g _Kv) for the electrical compartment and V _GK for the metabolic one, keeping the other parameters fixed at their default values (Loppini et al., 2015). In particular, g _Kv can impose a very fast busting or spiking response with variable frequency, and V _GK sets the intrinsic metabolic oscillation frequency. In order to simulate a fast spiking activity within the slow bursting active period, we set mean value and standard deviation of g _Kv to 0.215 nS/pF and 2%, respectively, while for V _GK to 0.0556 mM/s and 10%. The effect of these parameters was tested in uncoupled populations and single-cell simulations. This allowed for the identification of a V _GK cut-off value setting the onset of metabolic oscillations. Low values of V _GK not inducing metabolic oscillations at the single-cell level were discarded.

A full list of model equations, their detailed description and the complete set of parameters can be found in refs. Riz et al. (2014), Loppini et al. (2015).

It is known that the spatial organization of β-cells in pancreatic islets and islets composition differ across species (Brissova et al., 2005; Cabrera et al., 2006; Steiner et al., 2010). In humans, β-cells appear to be intermingled with α-, γ-, and δ-cells, and are about 50% of the islet cells. We built human-like β-cells networks by stacking 2D layers of cells in a hexagonal-closest packing (HCP) structure (Nittala et al., 2007; Nittala and Wang, 2008; Stamper et al., 2014; Cappon and Pedersen, 2016), with 12 layers, each including 12 × 12 cells. A single cell was modeled as a sphere of radius r _c = 6.5 μm. To build the final islet structure, we extracted a spherical region with 65 μm radius within the HCP stacked layers and imposed a site percolation probability of 50% (p_site = 0.50) to match data on β-cell percentage. Nearest-neighbor couplings were defined based on structural contacts between cells, i.e., for cells whose center-to-center distance was equal to 2r _c. We further deleted from the structure isolated cells resulting from percolation (less than 5). The total number of cells included in the reconstructed islet was N = 378, with an average and maximum number of neighbors equal to ≃ 5 and 10, respectively. The structural adjacency matrix a _ij describing nearest-neighbors was then used in model equations to define the electrical and metabolic coupling terms (Eqs 1, 2).

Synchronization phenomena have been investigated by computing the correlation index between all pairs of β-cell binarized potential and metabolic signals (Rodgers and Nicewander, 1988; Stožer et al., 2013b): R i j = ⟨ x i t − ⟨ x i t ⟩ ( x j t − ⟨ x j t ⟩ ) ⟩ σ i σ j , (3)

where x _i and σ _i are the binarized signals and the corresponding standard deviation for the i ^th cell. The binarized signals were constructed by evaluating at first the onset time in cells activation, both on membrane potential (V) and FBP, by means of a signal derivative thresholding, and at second, by imposing an active time duration of 15 ms for V, and 1.2 s for FBP. Therefore the binarized signals were equal to 0 and 1 outside and inside the active time, respectively. This choice allowed us to match the action potential duration in V and the initial rise in FBP.

The correlation indices were further thresholded to obtain electrical and metabolic functional networks. In particular, we investigated the behavior of the networks by varying the correlation thresholds, R_V,thr and R_FBP,thr, in the range of [0,1]. We finally performed a motif analysis by considering, at the three network levels (structural, electrical, and metabolic), the pairwise multiplex connection motifs, for a total of 2³ = 8 possible motifs. In this context, a multiplex pair-wise connection motif represents how β-cell pairs are connected through the structural, electrical and metabolic layers, and serves to compare cells connectivity at different levels. We investigated connection motifs at increasing distances over the islet, quantifying synchronization spatial features in comparison to the underlying structural architecture. To choose the optimal value of the correlation thresholds to be used in the motif analysis, we studied four network parameters (average node degree, average local clustering coefficient, average local efficiency, and number of connected components) as functions of the correlation thresholds for both electrical and metabolic dynamics. Specifically, the average node degree is evaluated as the mean of the nodes degree k _i = ∑ _j∈N a _ij . The average local clustering is the mean of the local clustering coefficient of the nodes c _i = 2e _i /(k _i (k _i − 1)), where e _i is the number of connections in the local subgraph of node i, i.e., G _i . Similarly, the average local efficiency is the mean of the local efficiency of the nodes E i = ( 1 / ( k i ( k i − 1 ) ) ) ∑ j ≠ j ′ ∈ G i 1 / d j j ′ , where d _jj′ is the distance between nodes j and j′ (Latora and Marchiori, 2001; Boccaletti et al., 2006). In relation to intercellular networks, the node degree quantifies the number cells to which a cell is connected, the local clustering coefficient measures the tendency of cells to develop dense local connection patterns, the local efficiency is linked the clustering coefficient and quantifies how efficient is the local intercellular communication, while the number of connected components measures the amount of disjoint connected cells groups in the network. In the context of functional networks, high degree, clustering and efficiency, and a low number of connected components denote strong correlations in cells activity encompassing large groups of cells.

β-cells architecture was reproduced with MATLAB R2021a (The MathWorks, Inc.), together with the data analysis. The numerical integration of the biophysical model is implemented in a C++ algorithm using a fourth-order Runge–Kutta scheme for the electrical subsystem, and an Euler scheme for the metabolic subsystem in parallel, due to the non-restrictive slow dynamics of the glycolytic oscillator. In each case, a fixed time step of 0.02 ms was used. Numerical accuracy was verified with a MATLAB implementation for both the electrical and the glycolytic compartments, using a stiff integrator.

At first, we investigated the islet dynamics with electric and metabolic coupling strengths as previously used in other modeling studies, i.e., g _c = 0.01 nS/pF and P _G6PF6P = 0.01 ms⁻¹. At the single-cell level, the pattern is a slow bursting oscillation driven by glucose metabolism with a bursting period of ≃ 160 s (frequency of ≃ 0.37 min⁻¹) and active phase duration of ≃ 45 s, resembling slow electrical activity observed in human β-cells (Braun et al., 2012; Rorsman and Braun, 2013; Riz et al., 2014). The membrane potential is characterized by a fast action potential firing intermingled to silent periods, while the intracellular calcium shows a compound oscillation with fast and slow components following the membrane potential oscillations and metabolites dynamics, respectively (Figure 1A). As a representative variable of the glycolytic component, we focused on FBP, whose oscillations are in phase with membrane potential active phases and with the slow calcium component. To analyze cells coordination at the population level, we investigated the activation times of cells on both membrane potential and FBP, as extracted from the peaks of signals derivative. We observed a coherent and synchronized activity all over the islet, characterized by in-phase action potentials within the active periods both between near and distant β-cells. Dephasing in cells activity was instead observed at the onset and termination of the active phases. Indeed, cells activation at the beginning of the active phases followed a wave pattern propagating through the islet. This pattern was observed on membrane potential, on calcium, and on FBP and was progressively lost when moving within the central part of the active phases, where cells were more synchronized and electrically activating in small groups. The wave pattern and the synchronized regime just after are shown in Figure 1B (top) by means of raster plots computed on the FBP and membrane potential activation times. To further dissect the spatial organization of the activation wave, we evaluated, for each burst, the FBP activation time for every cell with respect to the first responding cell in the islet (ΔT onset) and analyzed its spatial spreading over the islet structure (Figure 1C, left). At g _c = 0.01 nS/pF and P _G6PF6P = 0.01 ms⁻¹, ΔT onset was in the range 0–0.45 s with an organized spatial distribution indicating an activation wave triggered by peripheral cells and spreading toward the opposite area of the islet. To get deeper insights in wave dynamics, we analyzed the correlation of ΔT onset versus the values of parameters we varied in the β-cells population (V _GK and g _Kv) and found a negative correlation with the glucokinase reaction rate V _GK (Figure 1D, top), suggesting that the wave initiators were also among the cells with increased glucokinase reaction rate. Given that the diameter of our simulated islet was ≃ 125 μm, we estimated a wave velocity of about 280 μm/s at this coupling conditions, a value significantly higher than the ones experimentally measured in real human islets, close to ≃ 10 μm/s (Gosak et al., 2022). Based on this aspect, we further tested the model at a lower metabolites permeability, i.e., at P _G6PF6P = 0.001 ms⁻¹, keeping unaltered the electrical coupling. In this conditions, same bursting dynamics was recovered at the single-cell level, but at the population level, cells dephasing in slow oscillations was more pronounced, especially at bursting onset (Figure 1B, bottom). We observed an increasing in ΔT onset, ranging in this case between 0 and 3.5 s. Also, its spatial distribution followed a similar organization as the one computed at higher metabolites permeability, still indicating a wave spreading phenomena triggering the islet activation (Figure 1C, right). Same anti-correlation between ΔT onset versus V _GK was computed, as for the high permeability case (Figure 1D, bottom). Interestingly, in this condition we could obtain a wave velocity of ≃ 36 μm/s, i.e., on the same order of the one measured experimentally. In the following, the functional network analysis is performed on both models with high and low permeability to G6PF6P. We refer to the high and low permeability islet models as I1 and I2.

We extracted functional networks from the simulated dynamics by a thresholding procedure performed on the pairwise correlation indices between binarized signals of membrane potential and FBP (see Methods). The correlation matrices evaluated from simulated signals are shown in Figure 2, in comparison to the structural adjacency matrix. In particular, the functional network is extracted by filtering only the pairwise correlations higher than a selected threshold, resulting in a binary adjacency matrix describing the functional topology (functional connectivity matrix). In order to identify the correct threshold for both the electrical layer (V) and the metabolic layer (FBP), we focused on the behavior of four metrics commonly used in complex networks, namely, the average degree, average local clustering, average local efficiency and number of components. Their variations versus the threshold for both I1 and I2 are shown in Figure 3. The network efficiency is not reported, having a similar variation and slightly higher values with respect to the clustering. As expected, network metrics were observed to decrease at increasing thresholds, except for the average clustering in the metabolic layer, which showed a minimal increase at thresholds close to 1. This behavior was mainly due to the fragmentation of the functional network into small connected components. When moving from 0 to 1 in terms of threshold, the functional networks range from fully connected networks (or almost fully) to fully disconnected nodes. Given its higher correlations in cells activity, the I1 model showed higher network metrics with respect to I2. To select specific thresholds for both the models and layers, we compared the metrics, minimizing the number of components, limiting the average degree to ≃ 20% of the total number of cells, and guarantying an average local clustering and local efficiency in the range 0.7–0.8 and 0.8–0.9, respectively. The aim was to avoid too sparse functional networks analyzing differences both in short- and long-range correlations between electrical and metabolic layers, therefore selecting thresholds were 0.9 and 0.95 for V and FBP for I1, respectively, and 0.3 (V) and 0.75 (FBP) for I2. Such a choice led to a relative abundance of strongly connected nodes. We also tested higher values of thresholds giving rise to more sparse networks and a reduced number of hubs, obtaining results similar to the ones here presented. The extracted functional connectivity matrices are shown in Figure 4 (panels A, B, D, and E). To analyze the link between the extracted functional networks and wave activation propagation, we checked on the relation of nodes degree with respect to ΔT onset for both I1 and I2 and for both FBP and V activations. For the membrane potential, we focused only on the first action potential wave at the active phase onset. Interestingly, for every condition, the more connected nodes were the ones activating at intermediate times (see Figure 5A). A similar bell-shaped curve was also obtained for the I2 model at higher thresholds (0.85 and 0.9 on the metabolic layer), proving that a similar relation holds in sparser functional networks (Figure 5B). In particular, activation waves in our simulations were triggered by a small peripheral population of cells with a relatively low number of functional connections, recruiting nearby connected cells, further activating inner cells, and eventually propagating the activation to the last responding cells with functional connectivity comparable to wave initiators. The occurrence of higher functional connections in intermediate activating cells can be ascribed to the fact that their activity is correlated to both first and last responders, as well as showing a high correlation among themselves. In other words, at intermediate times, a subpopulation of cells synchronously activates and delivers the activation signal from first to last responders.

Based on the extracted functional networks, the islet synchronization properties were studied by adopting a multiplex description composed of three layers: structural, electrical, and metabolic. In particular, we looked at the pairwise connection motifs in the network, analyzing their relative frequency at increasing distances from a reference cell. The analyzed motifs included all the possible connections combinations on the three layers: structural/electrical/metabolic (M1), electrical/metabolic (M2), structural/metabolic (M3), structural/electrical (M4), structural (M5), electrical (M6), metabolic (M7), disconnected (M8) (see Figure 4). Specifically, we analyzed M(r) as the relative frequency of motifs within a distance in the range [0, r] (Figures 4C, F). As expected, all the motifs involving a structural connection (M1, M3, M4, and M5) decreased monotonically and strongly at increasing distances for both I1 and I2. For such motifs, and consistently with the reduced value of permeability, I2 showed a reduced fraction of M1 and an increased fraction of M4, indicating a slight loss of functional metabolic coupling. Interestingly, a combined electrical and functional connection without a structural link (M2) was significantly present at distances within the range 20–50 μm for both I1 and I2, indicating a significant spatial extension of combined electric and metabolic coordination. Concerning the frequency of exclusive electrical and metabolic functional connections (M6 and M7), they were found to be more abundant at r > 40–50 μm, with long-range metabolic connections outnumbering long-range electrical ones.

Given the importance of biological noise in β-cells dynamics due to stochastic ion channel gating on the cell membrane, we investigated the effects of an additional stochastic term in the membrane potential dynamics on synchronization and functional multiplex motifs. In particular, we added a noise term in the membrane potential equation with shape σW(t), where σ represents the noise intensity, i.e., the standard deviation, and W(t) is a Gaussian white-noise process with zero mean and covariance ⟨W(t)W(t′)⟩ = δ(t − t′). In our simulations, we used σ = 0.2, and we tested the effect of such a stochastic perturbation on the I2 model, closer than I1 to the experimental observations in terms of activation waves and globally less synchronized. Interestingly, the network metrics evaluated on the binarized membrane potential revealed an increased electrical synchronization across the β-cells population at low thresholds (Figures 6A–C). This consideration is evident from the behavior of the average degree and number of components, indicating that the functional network exhibits a fully coupled network at very low R_V,thr values in these conditions. When moving to higher thresholds (≃ 0.3), the average number of connections per node is slightly higher with respect to the deterministic case (63 vs. 57), while the reduced value of the local clustering is mainly due to the appearance of locally larger coupled functional aggregates having a clustering value lower than the smaller assemblies extracted in the absence of noise. The frequency of the multiplex connection motifs (Figure 6D) nicely summarizes this increase in membrane potential synchronization, showing higher frequencies for every motif including a functional electrical connection. In this regard, combined electrical and metabolic coordination between cells pairs not structurally connected was reinforced within small groups when a moderate noise was added as input.

Recent findings on β-cells network activity pointed out the existence of different sub-populations of β-cells. To investigate the impact of multiple sub-populations, we tested the I2 model response imposing a bimodal random sorting of the glucokinase reaction rate and defined subsets of normally responding (subset A) and highly-active cells (subset B), similarly to other modeling investigations (Dwulet et al., 2019). Specifically, we extracted the V _GK with mean and standard deviation of 0.0756 mM/s and 3% for A, and 0.1556 mM/s and 3% for B. Also, subset B included 10% of the total number of cells, in line with experimental findings on hub cells percentage in pancreatic islets, i.e., a sub-population with a characteristic fingerprint, also including increased metabolic activity. At the cell scale, slow bursting was only minimally affected, with a slight increase in bursting frequency ( ≃ 0.44 min⁻¹). Also, the wave pattern activation was still present and characterized by a slight increase in time-lags between wave initiators and last-responding cells, indicating a mild decrease in the computed velocity ( ≃ 30 μm/s). To address the subset of cells originating the wave, we analyzed the V _GK parameter with respect to the ΔT onset, as defined previously, and found that subset B was the main pacemaker of wave activity (Figure 7). No correlation between V _GK and ΔT onset was instead found within group A. We further analyzed functional network properties at different thresholds (Figure 8) and found no substantial differences with respect to the uni-modal case for I2 (Figure 3). At the same thresholds used for the uni-modal I2, more functionally connected cells, i.e., functional hubs, still overlapped with cells activating at intermediate time within the activation wave front (Figure 7). Further, the multiplex functional connection motifs for the bimodal case revealed correlation patterns consistent with the uni-modal model (Figure 8).