Spatial distribution of heterogeneity as a modulator of collective dynamics in pancreatic beta-cell networks and beyond

1 Introduction

Many non-linear systems exhibit excitable behaviour, whereby they exhibit large-amplitude oscillations in response to small-amplitude, transient perturbations. One prominent example in biology is electrical excitability in cells, which underlies the function of neurons (Izhikevich, 2000; De Maesschalck and Wechselberger, 2015; Wedgwood et al., 2021), cardiac cells (Majumder et al., 2018; Barrio et al., 2020), pituitary cells (Sanchez-Cardenas et al., 2010; Hodson et al., 2012) and pancreatic beta cells (Bertram et al., 2007; McKenna et al., 2016). Moreover, the study of these systems is broadly applicable, and such excitable dynamics are also observed in semiconductor lasers (Terrien et al., 2020; Terrien et al., 2021), social media networks (Mathiesen et al., 2013), epidemiology (Vannucchi and Boccaletti, 2004), and wildfires (Punckt et al., 2015). When excitable units are combined into networks, they can generate complex rhythms (Bittihn et al., 2017; Fretter et al., 2017; Hörning et al., 2017). Interestingly, such networks may also generate dynamics that occur over low-dimensional manifolds of the full system (Ashwin and Swift, 1992; Watanabe and Strogatz, 1993; Ott and Antonsen, 2009; Bick et al., 2020). For example, neurons in the pre-Bötzinger complex fire synchronously to induce the inspiratory and expiratory phases during breathing (Wittmeier et al., 2008; Gaiteri and Rubin, 2011).

Heterogeneity is ubiquitous in natural systems. Whilst often portrayed as a undesirable attribute, it can play an important role in governing network dynamics (Manchanda et al., 2017; Delgado et al., 2018; Lambert and Vanni, 2018). For example, neurons may coarsely be stratified into excitatory and inhibitory groups, with the former promoting firing behaviour in other neurons and the latter suppressing it. When coupled, these neuronal subtypes give rise to a variety of behaviours, including synchronisation, and enable the network to respond differentially to incoming inputs (Börgers and Kopell, 2003; Börgers et al., 2005; Kopell et al., 2010). The classification of neuronal subtypes is becoming ever finer (Gouwens et al., 2019; Lipovsek et al., 2021) and it remains an open question as to how this heterogeneity governs overall brain dynamics. Even when networks comprise only a single unit type, heterogeneity may still impact the global dynamics. For example, if the natural frequencies of nodes in a coupled oscillator network are too far apart, the network will be unable to synchronise and will instead display more complex rhythms (Ottino-Löffler and Strogatz, 2016). In pancreatic islets, heterogeneity has increasingly been acknowledged as crucial for healthy glucose metabolism (Nasteska et al., 2021), and classification of beta-cell populations in particular is an area of active investigation. Heterogeneity in cell-intrinsic properties (i.e., excitability, metabolic activity, and genetic profiles), network properties, and functional properties have all been demonstrated (Johnston et al., 2016; Westacott et al., 2017; Nasteska et al., 2021; Kravets et al., 2022; Šterk et al., 2023).

Here, we explore transitions between quiescent states and collective oscillations in coupled networks of heterogeneous, excitable nodes. For the node dynamics, we treat two types of system, the first being the FitzHugh–Nagumo (FHN) model, which is a prototypical model of electrical excitability in biological cells (Fitzhugh, 1961; Nagumo, 1962), is often used to investigate collective network dynamics, and has been used as a model for beta cells (Scialla et al., 2021; Pedersen et al., 2022). For the second, we consider a conductance-based model of pancreatic beta cells designed to more closely reproduce the signalling dynamics in these cells (Sherman et al., 1988). These cells remain at rest until they receive a significantly large electrical impulse or, in the case of beta cells, the extracellular concentration of glucose surpasses a threshold value (Ashcroft Frances M and Rorsman, 1989; Braun et al., 2008).

For the past 40 years, based on empirical evidence from rodent islets, it has generally been assumed that beta cells form a syncticium, in which activity of an entire islet of Langerhans can be described by considering the dynamics of a single cell (Dolenšek et al., 2013; Podobnik et al., 2020; Satin et al., 2020). A number of studies have challenged this perspective in recent years, in particular highlighting that some ‘leader cells’ disproportionately influence the activity of the entire network which is made up primarily of ‘follower cells’ (Johnston et al., 2016; Westacott et al., 2017; Salem et al., 2019; Benninger and Kravets, 2021). In particular, one hypothesis suggests that beta cells can be grouped into two overarching classes based on their degree of excitability. This then suggests that islets are composed of a small number (∼10%) of highly excitable cells, with the remainder being less excitable (Benninger and Hodson, 2018).

To this end, we diffusively couple cells having two different degrees of excitability according to two different network architecture classes. The first of these comes from the Watts–Strogatz model, which allows for the generation of connected graphs, with consistent mean degree, and a parameter for adjusting the balance between local and long-range connections. This form of coupling allows one to explore architectures ranging from regular to small world graphs and, as such, is well-suited to investigations of the general types of graphs that occur in the natural world (Watts and Strogatz, 1998). While these graphs are not directly relevant to the electrical coupling architecture of pancreatic beta-cell networks, there is potential that the dynamics of beta-cell networks share common features with small-world networks. For example, small-worldness has be observed in the functional connectivity of pancreatic beta-cell Ca²⁺ signals (Stožer et al., 2013), and long-range connections have been proposed in mathematical modelling studies (Barua and Goel, 2016). The second model architecture was designed to mimic the spatial configuration of electrical coupling (i.e., connexin36 gap junctions) within beta-cell networks. Cells are locally coupled via gap junctions (Benninger et al., 2011; Rorsman and Braun, 2013) and spatially arranged in a three-dimensonal lattice embedded within a sphere.

With these network structures in mind, we explore how the spatial organisation of the two subpopulations affects the propensity of the whole network to oscillate in a synchronous fashion. The work of Westacott et al. (2017) recently demonstrated that β-cells with similar electrical excitability are spatially correlated within pancreatic islets. Motivated by this finding, we develop a metric for evaluating the degree of this spatial non-uniformity in β-cell networks and an algorithm for generating non-uniform networks with a specified value of this metric. We demonstrate that spatial correlation in heterogeneity, or network sortedness, plays a critical role in modulating collective dynamics in these networks. Moreover, we show that these metrics and algorithms can be applied to networks that do not have an explicit notion of space (e.g., Watts–Strogatz networks) and that network sortedness is similarly important in defining collective dynamics on these general network architectures. This work aims to demonstrate the importance of cellular organisation when studying heterogeneity in pancreatic β-cell networks (and in networks more generally) through a case study in heterogeneous excitability. Furthermore, it provides a set of tools for introducing correlated heterogeneities in arbitrary networks and studying their impact on dynamics. This spatial organisation is likely to be of particular relevance to the study of human islets, since these have been shown possess a highly non-trivial architecture, particularly when compared with rodent islets, on which most of our knowledge of beta cell function is based (Cabrera et al., 2006).

The remainder of the manuscript is arranged as follows: In Section 2, we describe the single cell dynamics and graph structures, introduce a metric that captures how sorted a network is with respect to its heterogeneity, and present an algorithm that can generate networks with arbitrary sortedness. In Section 3, we investigate how dynamic transitions to synchronous activity depend on the degree of sortedness in the network and end in Section 4 with concluding remarks.

2 Methods

2.1 Mathematical model

2.1.1 Fitzhugh-Nagumo model (FHN)

The Fitzhugh-Nagumo model is a minimal model of an excitable unit. It exhibits relaxation oscillations when the external stimulus parameter, I, exceeds a critical value, where a Hopf bifurcation occurs. Here, we consider a network of N diffusively-coupled excitable nodes, which may, for example, be cells, each of which is described by the following model.

$\frac{d{V}_i}{dt}=V_i-\frac{V_i^3}{3}-W_i+GI-I_{\mathit{coup,i}},i=1,\dots N,(1)$

$\tau \frac{d{W}_i}{dt}=V_i+a-b{W}_i.(2)$

Here, V models the fast upstroke characteristic of excitability whilst W models the slow recovery (i.e., negative feedback). We replace the regular stimulus current term, I, with GI, where I is now the maximum stimulus current to a node, and thus defines the excitability of the node, whereas G ∈ [0, 1] is the degree of drive to the network. This allows us to consider a heterogeneous network model based on the excitability of nodes, but whilst including a network-wide drive parameter. A list of parameters is included in Supplementary Table S1.

2.1.2 Sherman-Rinzel-Keizer model (SRK)

The electrical activity of pancreatic beta cells is proportional to the extracellular concentration of glucose. For sufficiently high extracellular glucose, the cells exhibit bursting dynamics, in which their voltage periodically switches between high frequency oscillations and quiescence. The high frequency oscillations in voltage are correlated with the secretion of insulin from these cells, so that these bursting dynamics are tightly coupled to the cells’ functional role. Here, we consider a network of N diffusively-coupled excitable cells, each of which is described by the three-variable model.

$C_m\frac{d{V}_i}{dt}=-I_K\left(V_i,n_i\right)-I_{Ca}\left(V_i\right)-I_{K-Ca}\left(V_i,c_i\right)-I_L\left(V_i\right)-I_{\mathit{coup,i}},i=1,\dots N,(3)$

$\frac{d{n}_i}{dt}=\frac{n_{\infty }\left(V_i\right)-n_i}{{\tau }_n\left(V_i\right)},(4)$

$\frac{d{c}_i}{dt}=-f\left(\alpha I_{Ca}\left(V_i\right)+k_{Ca}{c}_i\right).(5)$

The variable V represents the membrane voltage of a β-cell, n is the activation variable for K⁺ ion channels, and c is the cytosolic Ca²⁺. This system was adapted from the Sherman–Rinzel–Keizer model, which describes the dynamics of electrical activity in pancreatic β-cells in the presence of glucose (Sherman et al., 1988). The intrinsic dynamics of the voltage, V given by Eq. 3 are driven by K⁺ (I_K), Ca²⁺ (I_Ca), and Ca²⁺-activated K⁺ (I_K−Ca) ionic currents, with a rate governed by the whole cell capacitance given by C_m. These currents are described via.

$I_K\left(V,n\right)={\bar{g}}_{K}n\left(V-V_k\right),(6)$

$I_{Ca}\left(V\right)={\bar{g}}_{Ca}{m}_{\infty }\left(V\right)h_{\infty }\left(V\right)\left(V-V_{Ca}\right),(7)$

$I_{K-Ca}\left(V,c\right)={\bar{g}}_{K-Ca}\frac{c}{K_d+c}\left(V-V_k\right),(8)$

$I_L\left(V\right)={\bar{g}}_L\left(1-G\right)\left(V-V_K\right).(9)$

In Eqs 6–9, ${\bar{g}}_X$ denotes the maximal conductance of the channel X where X ∈ {K, Ca, K − Ca, L} where L signifies a leak channel; V_X are the reversal potentials of the respective channels, m and n are the proportion of open activating gates for the Ca²⁺ and K⁺ channels, respectively; h is the proportion of open inactivating Ca²⁺ channels; c is the cytosolic concentration of Ca²⁺; and G is the extracellular concentration of glucose, which provides a global drive to promote activity and is taken to be homogeneous across the network. The activation of I_K−Ca is a function of free intracellular Ca²⁺ concentration and is defined by a Hill-type function with disassociation constant K_d. The current I_coup,i captures the influence of the coupling between cells and will be discussed in Section 2.2.4.

The dynamics for n, m, and h follow exponential decay to their state values given by

$x_{\infty }\left(V\right)=\frac{1}{1+\exp \left[\left(V_x-V\right)/S_x\right]},x\in \left\{h,m,n\right\},(10)$

at a rate given by the voltage-dependent time constant.

${\tau }_n\left(V\right)=\frac{\bar{\tau }}{\exp \left[\left(V-\bar{V}\right)/{\kappa }_1\right]+\exp \left[-\left(V-\bar{V}\right)/{\kappa }_2\right]}(11)$

for n. In Eq. 10, V_x represents the activation (inactivation) thresholds for m and n (h) and S_x represents the sensitivity of the channels around this point. Finally, Eq. 5 describes the evolution of the concentration of cytosolic Ca²⁺, which decays and is pumped out of the cell following a combined linear process with rate k_Ca and enters the cell via the Ca²⁺ ion channel at a rate given by the scale factor α. The parameter f specifies the fraction of free to bound Ca²⁺ in the cell, where bound Ca²⁺ plays no role in the relevant dynamics in our model.

To expose the dependence of our system on glucose, we introduced a hyperpolarising leak current given by Eq. 9 that explicitly depends on the glucose concentration G. For an isolated cell (i.e., without coupling) with the parameters specified in Supplementary Table S2, the system describing each node exhibits steady state behaviour for low G and passes through a bifurcation as G ∈ [0, 1] is increased, as shown in Figure 1.

FIGURE 1

FIGURE 1. Excitability of single SRK cells. The voltage traces from 3 cells with varying levels of intrinsic excitability $({\bar{g}}_L)$, but the same level of drive (G = 0.7). The red, blue, and black traces show decreasing levels of excitability with values of ${\bar{g}}_L=60$, 120, and 180, respectively. More excitable cells have a shorter interburst interval. The solid black line represents a Hopf bifurcation as a function of both G and g_L. At the lowest drive (G = 0), the Hopf bifurcation occurs for ${\bar{g}}_L=45.21$ pS. The dotted lines represent “level sets” of the $({\bar{g}}_L,G)$ parameter space, along which the behaviour of the single cell is identical. Data for the bifurcation diagram was computed using XPP 8.0 (Ermentrout, 2002).

The bursting dynamics in our model are of the fold-homoclinic type under the classification specified in (Izhikevich, 2000). This classification is based on separation of the full system into a fast subsystem (Eqs 3, 4) and a slow subsystem (Eq. 5), treating the slow subsystem variables (in this case, c) as parameters in the fast subsystem. During each bursting cycle, the slow evolution of c pushes the fast subsystem through bifurcations that initiate and terminate oscillatory behaviour. In particular, when c decreases to a small enough value, the fast subsystem passes through a fold bifurcation in which a stable steady state and a saddle steady state collide and annihilate one another. Following this, the system exhibits stable periodic activity, during which c increases according to Eq. 5. When c increases to a sufficiently large value, the fast subsystem passes through a homoclinic bifurcation that destroys the periodic orbit and the system returns to the original stable steady state. Following this, c decreases until it once again reaches the fold point and the cycle repeats.

2.2 Model simulations

Simulations were conducted using Matlab 2019B. The dynamical systems were solved using ode15s, the relative tolerance set to 10^–5, and explicit Jacobians were provided. The code was run on the University of Birmingham BlueBEAR HPC running RedHat 8.3 (x86_64) (see http://www.birmingham.ac.uk/bear for more details). Each set of simulations ran over 16 cores using a maximum of 128 GB RAM (32 GB was sufficient in most cases). All code used in the project is freely available for download from: github.com/dgalvis/network_spatial.

2.2.1 Initial conditions

For the FHN model, the initial conditions $y_i(0)=\left(V_i(0),W_i(0)\right)$ for node i = 1, … , N were sampled independently from the distributions.

$V_i\left(0\right)\sim \mathcal{N}\left(-1,{\left(1/6\right)}^2\right),W_i\left(0\right)=0.(12)$

For the SRK model, the initial conditions $y_i(0)=\left(V_i(0),n_i(0),c_i(0)\right)$ for node i = 1, … , N were sampled independently from the distributions

$V_i\left(0\right)\sim \mathcal{N}\left(-68,{\left(68/6\right)}^2\right),n_i\left(0\right)=0,c_i\left(0\right)\sim \mathcal{N}\left(0.57,{\left(0.57/6\right)}^2\right).(13)$

The term $\mathcal{N}(\mu ,{\sigma }^2)$ represents a normal distribution with mean μ and variance σ². Throughout, we use Y (0) to denote the set of initial conditions across the whole network model, i.e., $Y(0)=\left(y_1(0),\dots ,y_N(0)\right)$.

2.2.2 Excitability and drive in the single-cell FHN model

The term GI is used to represent node excitability (I) and network drive (G). A Hopf bifurcation occurs at GI = 0.33. During the study, we partition the nodes into two sub-populations based on their excitability, with one population being highly excitable and the other being significantly less so.

2.2.3 Excitability and drive in the single-cell SRK model

The ionic current I_L (Eq. 9) is a hyperpolarising current that can be used to adjust the excitability of each cell and to determine the activation level of the network model. In particular, the maximum conductance ${\bar{g}}_L$ determines the excitability of a cell. As this value increases, the cell becomes less excitable, that is, for a given value of G, cells with higher ${\bar{g}}_L$ are less likely to burst. This behaviour is summarised in Figure 1, which shows a two parameter bifurcation diagram showing the transition from quiescent to bursting behaviour under simultaneous variation of $({\bar{g}}_L,G)$, which occurs via a Hopf bifurcation of the full system (Eqs 3–5). For G = 0, this Hopf bifurcation occurs at ${\bar{g}}_L=45.21$ pS. For non-zero values of G, the bifurcation curve is defined via $(1-G){\bar{g}}_L=45.21$ pS, as can be seen by examining the form of the Eqs 3, 9. Note that when G = 1, system (Eqs 3–5) matches that of (Sherman et al., 1988). We use the observations about the link between ${\bar{g}}_L$ and excitability to partition the network into two sub-populations, one being highly excitable, the other being significantly less excitable.

2.2.4 Graph structure and coupling

In this work, we consider two types of graph structure. The first is the Watts-Strogatz model of graph generation to create small-world graphs (WS graphs). We use N = 1000 nodes, a mean degree of D = 12, and unless stated, we use a rewiring probability of β = 0.2. The second type is motivated by pancreatic β-cell networks, which are arranged into roughly spherical clusters called islets of Langerhans (which also encompass other cell types that are disregarded in our model), which each contain ∼ 1,000 β-cells. To capture this, we arrange N = 1, 018 nodes on a hexagonal close packed (hcp) lattice embedded within a sphere (βC graph). Each node is connected to all of its nearest-neighbours so that the number of connections of nodes away from the boundary of the sphere is equal to the coordination number 12 whilst nodes on the boundary have fewer connections. This type of connectivity scheme is commonly employed in β-cell network modelling studies (Benninger et al., 2014; Westacott et al., 2017; Dwulet et al., 2019, 2021), as it qualitatively captures the network architecture. In reality, pancreatic islet architecture is more complex in terms of islet shape, cell shape, and the layout of cells in the tissue (Stožer et al., 2013; Šterk et al., 2023); all of which could be considered as additional sources of heterogeneity. Moreover, there may be heterogeneity in the number of connections that β-cells form amongst each other (Cabrera et al., 2006). We choose to employ this homogeneous lattice structure so that we could consider heterogeneity in cell-intrinsic excitability alone without introducing additional sources of heterogeneity due to the graph structure.

The dominant form of coupling between beta cells in the islets is through gap junctions, which allow small molecules, including charged ions to pass directly from a cell to its adjacent neighbours. Mathematically, this is represented through the inclusion of the diffusive term I_coup,i in Eq. 1 or Eq. 3 (we use the same coupling for both FHN and SRK models) given by

$I_{\mathit{coup,i}}={\bar{g}}_{\mathit{coup}}\sum _{j\in J_i}\left(V_i-V_j\right),(14)$

where J_i is the set of all nodes to which node i is coupled.

Together, the two models for node dynamics (FHN and SRK models) and the two types of graph structure (WS and βC graphs) give rise to four combinations of network models, which we denote by WS-FHN, WS-SRK, βC-FHN, and βC-SRK. Note that throughout this work, we use the term graph to refer specifically to the set of links between nodes. We use the term network or network model to refer to the entirety of the model including choice of dynamics on the node (FHN or SRK), connectivity structure, and choice of parameters.

2.2.5 Heterogeneity

For the FHN model, we consider networks consisting of two sub-populations of nodes distinguished by their excitability (i.e., by their I values). Population 1 is highly excitable (I = 2) and population 2 is less excitable (I = 1). Similarly for the SRK model, we consider networks consisting of two sub-populations of nodes distinguished by their excitability (i.e., by their ${\bar{g}}_L$ values). Population 1 is highly excitable (${\bar{g}}_L=60$ pS) and population 2 is less excitable (${\bar{g}}_L=100$ pS). We then consider the range over G for which population 1 nodes are intrinsically active (i.e., G such that population 1 is active when ${\bar{g}}_{\mathit{coup}}=0$). We then consider the effects of degree of sortedness between the two subpopulations (see Section 2.3.1), graph structure (WS or βC), global network drive (G), and global coupling strength $({\bar{g}}_{\mathit{coup}})$ on the collective dynamics of the network.

2.3 Measuring sortedness

2.3.1 Definition

To track the degree of sortedness in the network, we define a node sortedness measure that, for a given node, measures the proportion of neighbours that are of the same population type. For a general network with nodes attributed to $K\in \mathbb{N}$ populations, the node sortedness, A_i, is defined as

$A_i=\frac{1}{|J_i|}\sum _{j\in J_i}{\chi }_{ij},{\chi }_{ij}=\sum _{k=1}^K{\mu }_i^{\left(k\right)}{\mu }_j^{\left(k\right)},i=1,\dots ,N,{\mu }_i^{\left(k\right)}=\left\{\begin{array}{cc}1,\hfill & i\in P_k\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}\right.,(15)$

where the population sets P_k contain the indices of the nodes within population k = 1, … , K and form a partition over the node indices {1, 2, … , N}, J_i is the set of indices of nodes that are coupled to node i, ${\mu }_i^{(k)}$ is an indicator function that takes value 1 if i belongs to population k and value 0 otherwise, and χ_ij is an indicator function that takes value 1 when node i and j belong to the same population and value 0 otherwise. For each population, the average node sortedness is defined via.

${\bar{A}}_k=\frac{1}{|P_k|}\sum _{n\in P_k}{A}_n,k=\mathrm{1,2},\dots K.(16)$

Finally, the network sortedness is defined as

$\mathcal{A}=\frac{1}{K-1}\left(-1+\sum _{k=1}^K{\bar{A}}_k\right).(17)$

where $\mathcal{A}\in [-1/(K-1),1]$ and, for the present case with K = 2, $\mathcal{A}\in [-\mathrm{1,1}]$. For a network in which populations are assigned to nodes following a uniformly random distribution, $\mathcal{A}\approx 0$ since ${\bar{A}}_k$ is approximately equal to N_k/N where N_k, k = 1, 2 is the number of nodes in population k. And therefore ${\sum }_k{\bar{A}}_k\approx 1$. An illustration of the computation of the sortedness metrics (Eqs 15–17) is shown in Supplementary Figure S1A.

2.3.2 Modified network sortedness for beta cell networks

For a βC lattice network, ${\bar{A}}_k$ (as defined in Eq. 16) is maximised when population 1 nodes are arranged in a single cluster with a minimal number of connections to population 2 nodes. This naturally occurs at the edges of the domain, since any cluster of population 1 nodes in the domain interior must be surrounded by population 2 nodes. We are interested in the dynamics that arise as the small population of highly excitable nodes forms clusters within the lattice, hence, we wish to remove this tendency for clusters to form at the domain boundary. To overcome this, we use a modified definition of the node sortedness (Eq. 15).

${\tilde{A}}_i=\frac{1}{J}\sum _{j\in J_i}{\chi }_{ij}+\frac{{\mu }_i^{\left(2\right)}\left(J-|J_i|\right)}{J},i=1,\dots N,(18)$

where J = 12 is the number of connections that interior lattice nodes possess. For nodes with |J_i| < J (i.e., nodes on the domain boundary) the additional term in Eq. 18 compared to Eq. 15 incorporates a further J − |J_i| connections to population 2 nodes for the purposes of calculating node sortedness values. This procedure is equivalent to assuming that the lattice defining our domain is embedded within a larger lattice of population 2 nodes.

2.4 Modifying network sortedness

2.4.1 The sorting algorithm

Here, we describe our approach for generating networks with different network sortedness. The algorithm works by exchanging the population type of nodes from different populations randomly to increase (or decrease) $\mathcal{A}$. The algorithm begins by randomly permuting the order of the N indices. The first N₁ indices of the permuted sequence are attributed to P₁, with the remaining N₂ indices attributed to P₂, yielding a distribution of population 1 nodes that is uniformly random in space.

On each iteration, a, of the algorithm, pairs of nodes (from different populations) are sampled without replacement from a joint probability density function (pdf) P(X = i, Y = j) = f (i, j), i ∈ P₁, j ∈ P₂, where X and Y are random integer variables indicating the node selected from population 1 and 2, respectively. The population types of these nodes are then exchanged, that is, if i ∈ P₁ and j ∈ P₂, then i is added to P₂ and removed from P₁ and vice versa for j. The network sortedness (Eq. 17) is then recomputed for the adjusted population sets. If the exchange leads to an increase (decrease) in $\mathcal{A}$, the exchange is accepted and the algorithm proceeds to iteration a + 1. If the exchange does not lead to an increase (decrease) in $\mathcal{A}$, the exchange is rejected and indices i and j are placed back in P₁ and P₂, respectively. In this case, a new pair of nodes is drawn from f and the process is repeated until either: a pair whose exchange leads to an increase (decrease) in $\mathcal{A}$ is found and the algorithm proceeds to the next iteration; or it is determined that no such pair exists, at which point the algorithm terminates. An example of one iteration of this algorithm is depicted in Supplementary Figure S1B. We refer to the algorithm in which swaps are accepted only if they lead to an increase (decrease) in $\mathcal{A}$ as the forward (backward) algorithm. We define ${\mathcal{A}}_a$ to be the evaluation of $\mathcal{A}$ of the network after a iterations, where a ∈ {0 … a_final}. Running the algorithm to convergence (a_final) produces the sets ${\mathcal{P}}_k={\{P_k^a\}}_{a=0}^{a_{\text{final}}}$ containing the population sets after each iteration.

Figure 2 depicts a WS network (Figure 2A, top row) and a βC network (Figure 2B, top row) at five iterations of the sorting algorithm. The WS network can be described by a set of nodes equispaced along a ring, and when the rewiring probability is small, the majority of connections are between nodes that are close to one another on the ring. For our choice of D = 12 (and β = 0.2), the majority of connections will occur between an arbitrary node and the D/2 = 6 closest nodes in either direction. Figure 2A shows that as $\mathcal{A}$ increases, the population 1 nodes (blue points) begin to form clusters along the ring due to the high interconnectivity amongst near neighbours. They do not necessarily form a single cluster at convergence (a = 314 in this example), but due to the rewirings, many connections between clusters exist (blue lines). In the case of a βC network (Figure 2B), we see the emergence of localised population 1 clusters culminating in a single cluster at convergence (a = 239 in this example).

FIGURE 2

FIGURE 2. Dynamics on networks defined by the sorting algorithm. (A) Top Row: Five iterations of the sorting algorithm on a WS graph (with β = 0.2, D = 12). The majority of connections in this graph $(\approx 80\%)$ are between a node and its D/2 = 6 closest neighbours in each direction along the ring. The other $\approx 20\%$ of connections are randomly selected for rewiring. In each panel, only the connections within population 1 (blue lines) are shown, and only the rewiring connections are visible. As $\mathcal{A}$ increases, the population 1 nodes (blue points) begin to form clusters with many connections between clusters. The black boxes are shown zoomed in on Supplementary Figure S1C. Bottom Row: The dynamics of V (WS-FHN) averaged across population 1 (blue) and 2 (black) are shown for a fixed value of G = 0.253, for strong coupling ${\bar{g}}_{\mathit{coup}}=0.1$, and on each of the five iterations of the sorting algorithm. A phase transition from quiescence to global, synchronised activity occurs with respect to $\mathcal{A}$. (B) Top Row: Five iterations of the sorting algorithm on the β-cell graph. Initially, the population 1 cells (blue) are spread uniform randomly throughout the lattice. As $\mathcal{A}$ increases, they begin to form localised clusters culminating in a single cluster at convergence (a = 239). Bottom Row: The dynamics of V (βC-FHN) averaged across population 1 (blue) and 2 (black) are shown for a fixed value of G = 0.253, for strong coupling ${\bar{g}}_{\mathit{coup}}=0.1$, and on each of the five iterations of the sorting algorithm. A phase transition from quiescence to global, synchronised activity again occurs with respect to $\mathcal{A}$. (C) The average number of peaks $\bar{P}$ as a function of $(G,\mathcal{A})$ for one run of the sorting algorithm on a WS-FHN network (left) and one run on the βC-FHN network (right) and for 20 equispaced values of G ∈ [0.15, 0.345]. As sortedness increases, the drive necessary for global, synchronised activation of the network decreases.

2.4.2 Node selection probabilities

In this section, we formulate the node selection pdf used in the network sortedness adjustment algorithm. We assume that the selection of node from P₁ is independent of the selection of node from P₂:

$f\left(i,j\right)=f_{P_1}\left(i\right)f_{P_2}\left(j\right),i\in P_1,j\in P_2.(19)$

One choice would set f₁ and f₂ to be uniform over P₁ and P₂, respectively. This is the choice we take for the WS networks.

Empirical observations of the algorithm applied to the βC networks, however, demonstrate that clusters of population 1 nodes tend to form at the edge of the domain. As discussed in Section 2.3.2, we wish to avoid this scenario. The tendency for clusters to form near the edge occurs because of the spherical nature of our lattice domain. In particular, a uniform choice for f₁ and f₂ means that nodes at the centre of the domain are less likely to be selected under a uniformly random sampling of indices than those at the edge because the number of nodes in the network increases superlinearly with respect to the domain radius. Therefore, we derive choices for $f_{P_k}$ that equalise the probability of a node being selected on the basis of its radial coordinate. The heuristic for generating $f_{P_1}$ will be the same as that for generating $f_{P_2}$ up to the population identity.

Denote the radial distance from the origin of node $i\in {\mathbb{N}}_N$ by $r_i={(x_i^2+y_i^2+z_i^2)}^{1/2}\in {\mathbb{R}}_{\ge 0}$ where $(x_i,y_i,z_i)\in {\mathbb{R}}^3$ are the Cartesian coordinates of the location of the node. We define a sequence of intervals, ${\mathcal{I}}_n=[(n-1)\delta r,n\delta r]$, for n = 1, … 8 with δr = r_max/8 where r_max = max_i{r_i} so that each node is assigned to exactly one interval. The set of nodes from P_k belonging to a given interval ${\mathcal{I}}_n$ is given by $R_{n,P_k}=\{i\in {\mathbb{N}}_N\mid r_i\in {\mathcal{I}}_n,i\in P_k\}$. Using these set definitions, the pdf $f_{P_k}$ may be defined as

$f_{P_k}\left(i\right)=\frac{1}{Q|R_{n_i,P_k}|},i\in {\mathbb{N}}_N.(20)$

where $R_{n_i,P_k}$ is such that $r_i\in {\mathcal{I}}_{n_i}$ and Q is a normalisation factor ensuring that ${\sum }_{i\in P_k}{f}_{P_k}(i)=1$. This choice for $f_{P_k}$ reweights the probability of a given node being selected by a factor proportional to the number of nodes from the same population within a spherical annulus with inner and outer radii specified by the boundaries of the intervals ${\mathcal{I}}_n$. This reweighting favours selecting nodes closer to the centre of the domain over those more distal.

Pseudocode for the βC network generation and network sortedness manipulation algorithm is provided in the Supplementary Section S2.

2.5 Evaluation of collective dynamics

To characterise the network dynamics, we consider two features based on the Ca²⁺ trajectories across all nodes, namely, the mean number of peaks $(\bar{P})$ and the time-averaged degree of phase synchronisation $(\bar{R})$ calculated as the average magnitude of the Kuramoto order parameter (see Supplementary Section S3). The mean number of Ca²⁺ peaks across all nodes is proportional to the network participation, that is, the fraction of nodes that undergo oscillation. The value of $\bar{R}$ captures the network coordination, tracking how closely the phases of the Ca²⁺ trajectories stay to one another across the simulation duration. We additionally define ${\bar{P}}_k$ and ${\bar{R}}_k$ where k ∈ {1, 2} to be the mean number of peaks in Ca²⁺ and the time-averaged degree of phase synchronisation across nodes in population k, respectively.

3 Results

3.1 For strong coupling, increasing sortedness decreases the necessary drive for a phase transition from quiescence into global and synchronised activity

We ran the sorting algorithm to convergence once on a typical WS network and again on the βC network to produce the sets ${\mathcal{P}}_k$ for both network types (denote them WS-${\mathcal{P}}_k$ and βC-${\mathcal{P}}_k$ for k ∈ {1, 2}). For the WS network, this means specifically that we ran the Watts-Strogatz model of graph generation once to create a single small-world graph defining the connectivity structure of the network. We then uniform-randomly selected 10% of the nodes to be in population 1, thus defining the population sets $P_k^0$. These were used as the initial configuration for the sorting algorithm, which was run to convergence (Figure 2A, top row) producing WS-${\mathcal{P}}_k$. For the βC network, the connectivity structure is fixed. We defined $P_k^0$ as above and ran the sorting algorithm to convergence (Figure 2B, top row) producing βC-${\mathcal{P}}_k$. For both network structures (WS and βC), we then simulated the FHN dynamics on the network models defined by $P_k^a$ for all a ∈ [0, a_final] over 20 equispaced values of G ∈ [0.15, 0.345] with strong coupling $({\bar{g}}_{\mathit{coup}}=0.1)$. In this case, we used the same set of initial conditions across runs.

Figure 2C shows the average number of peaks $\bar{P}$ as a function of (G, $\mathcal{A}$). We found that as sortedness increases, the drive necessary to activate the networks decreases. Moreover, we found that for strong coupling, this phase transition results in a transition from quiescence to global, synchronised activity ($\bar{R}$ not shown but takes a value $\bar{R}\approx 1$ when $\bar{P}\ne 0$). This was the case for both the WS-FHN (left panel) and βC-FHN networks (right panel). Figures 2A, B (bottom rows) shows the dynamics of V (for WS-FHN and βC-FHN resp.) averaged over population 1 and 2 using a fixed value of G. This illustrates the strong synchronisation between populations (for strong coupling) as well as the phase transition as a function of $\mathcal{A}$ alone.

We next verified that this observed relationship is robust to changes in ${\mathcal{P}}_k$ and Y (0). To do this, we first calculated a Latin hypercube of pairs (G, a). For each point, we used a different seed network $P_k^0$ and ran the sorting algorithm to find $P_k^a$. We then ran dynamics on the network model using a different set of initial conditions Y (0). Note that for the WS models, each seed network was based both on a different WS graph and a different initial choice of indices for the population 1 nodes. For the SRK model, we used G ∈ [0.2, 0.6], ${\bar{g}}_{\mathit{coup}}=10$ (strong coupling), and 5000 (G, a) pairs. For the FHN model, we used G ∈ [0.15, 0.34], ${\bar{g}}_{\mathit{coup}}=0.1$, and 10,000 (G, a) pairs. We then calculated the features ($\bar{P}$ and $\bar{R}$) for all the points $(G,\mathcal{A})$ (where $\mathcal{A}$ is calculated from $P_k^a$ for each point) and used them to train a Gaussian process regression (GP) model (using the fitrgp function in MATLAB). We then estimated the features on a regular grid of points $(G,\mathcal{A})$ and used that to generate heatmaps for the features. Supplementary Figure S2 shows the process of generating these heatmaps for $\bar{P}$ in the case of the WS-FHN network models.

We found that for each of the four combinations of network models, increasing $\mathcal{A}$ results in a decreased G required for a phase transition from quiescence to global, synchronised activity. Figure 3A shows $\bar{P}$ for the WS networks and Figure 3B shows the same for the βC networks. The synchronisation parameter $\bar{R}$ is again not shown but takes values $\bar{R}\approx 1$ when $\bar{P}\ne 0$. For the heatmaps in Figures 3A, B, the $\bar{P}$ half-maximum value is shown as a black curve.

FIGURE 3

FIGURE 3. Phase transitions with respect to $\mathcal{A}$ for strong coupling. (A) The phase transition from quiescence to global, synchronised activity requires lower drive as sortedness increases in WS networks and for both models (FHN and SRK) when coupling is strong (${\bar{g}}_{\mathit{coup}}=0.1$ for FHN and ${\bar{g}}_{\mathit{coup}}=10$ for SRK). (B) This is likewise the case for the βC-FHN and βC-SRK model for strong coupling. (C) As β for the Watts-Strogatz model increases, the model transitions from a regular model with no rewiring to a random network. Here, we show the phase transition for three values of the rewiring probability. We found that the more regular a network, the lower the threshold G for the phase transition across the range of $\mathcal{A}$.

For the WS networks, the rewiring parameter β is related to the small-worldness of the network and inversely related to the regularity of the graph. Moreover, the maximum possible value of $\mathcal{A}$ will occur for β = 0, when population 1 nodes form a single cluster. We found that, in general, the sorting algorithm converges to larger values of $\mathcal{A}$ for smaller values of β. Figure 3C shows the $\bar{P}$ half-maximum value (using the Latin hypercube protocol described above for WS-FHN) for β ∈ {0.1, 0.2, 0.4}. We found that for a given value of $\mathcal{A}$, the G value at the phase transition is smaller for smaller values of β. This supports the hypothesis that the networks with more localised connectivity, such as the WS graphs when β ≈ 0 and the βC graph, tend to have a lower activation threshold than networks with increased small-worldness (also compare Figures 3A, B). This is perhaps unsurprising because in the networks with localised connectivity, clusters of excitable nodes only interact with the less excitable nodes along some boundary. For example, in the βC network, the nodes in the centre of a population 1 cluster will be completely separate from population two nodes. Whereas in the networks with small-world properties (WS, β ≠ 0), very few population 1 nodes (if any) will have a node sortedness of 1.

3.2 Weak coupling and high sortedness leads to wave propagation in the FHN networks

Next, we studied the influence of coupling strength on the relationship between $\mathcal{A}$, G, and the features of collective dynamics ($\bar{P}$ and $\bar{R}$). To do this, we calculated a Latin hypercube of triplets $({\bar{g}}_{\mathit{coup}},G,a)$. As before, we used different seed networks to identify $P_k^a$ and ran the dynamics using different initial conditions Y (0). For the FHN model, we used G ∈ [0.15, 0.34], ${\bar{g}}_{\mathit{coup}}\in [0.02,0.1]$, and 20,000 $({\bar{g}}_{\mathit{coup}},G,a)$ triplets. For the SRK model, we used G ∈ [0.2, 0.6], ${\bar{g}}_{\mathit{coup}}\in [\mathrm{2,10}]$, and 10,000 triplets. We again calculated the features ($\bar{P}$ and $\bar{R}$) for all the points $({\bar{g}}_{\mathit{coup}},G,\mathcal{A})$ (where $\mathcal{A}$ is calculated from $P_k^a$ for each point) and used them train a GP model. We used these to draw isosurfaces (level sets) of $\bar{P}$ and $\bar{R}$ as functions of $({\bar{g}}_{\mathit{coup}},G,\mathcal{A})$.

For the WS-FHN networks, we found that the intra-population isosurfaces (e.g., ${\bar{P}}_1$ and ${\bar{P}}_2$) corresponding to the phase transition between quiescence and activity separate for weak coupling and high $\mathcal{A}$. Figure 4A (left) shows the isosurfaces $\{({\bar{g}}_{\mathit{coup}},G,a)\mid {\bar{P}}_k=0.5\ast \max {\bar{P}}_k\}$ for k = 1 (blue, population 1) and k = 2 (black, population 2). For population 1, the phase transition threshold (G) decreases as $\mathcal{A}$ increases for all values of ${\bar{g}}_{\mathit{coup}}\in [0.02,0.1]$. However, for population 2, this relationship was only observed for strong coupling. For weak coupling and high sortedness, ${\bar{P}}_2$ is an increasing function of $\mathcal{A}$. Figure 4B (left) shows the heatmap of $\bar{P}$ for ${\bar{g}}_{\mathit{coup}}=0.02$ as well as the contour curves $\{({\bar{g}}_{\mathit{coup}}=0.02,G,a)\mid {\bar{P}}_k=0.5\ast \max {\bar{P}}_k\}$ for k = 1 (blue, population 1) and k = 2 (black, population 2). In the region between these contours, we found a region where all of population 1 activates, but many population 2 nodes do not. Figure 4B (diamond) shows a raster plot of all nodes (population 1 in blue) arranged so that neighbours on the raster plot correspond to neighbours on the ring (see Figure 2A). Population 1 clusters become active and activate the population 2 nodes neighbouring the cluster boundaries. The population 2 nodes then display wave-like propagation (i.e., time-delayed) that terminates before all the nodes have become active. The waves occur because in the region, sortedness is sufficiently high that the primary signal activating many population 2 nodes comes via local connections from previously activated population 2 nodes. Moreover, the long-range connections from population 1 nodes are weak and insufficient to activate population 2 in a synchronous manner [as in Figure 4B (square)]. To show this, we reran the experiment using β = 0 (no rewirings), which is shown in Supplementary Figure S3. We found that regardless of choice of coupling strength and sortedness, if the network activates, it produces waves emanating from one or more clusters of population 1 nodes. Moreover, we did not observe wave termination when β = 0, implying that it is the rewiring of connections that leads to the wave termination. Since the activation of a population 2 node is dependent on the strength of local connections in the chain and rewiring decreases the number of local connections, there exist some population 2 nodes (or clusters) that do not activate, thus terminating the wave.

FIGURE 4

FIGURE 4. Weak coupling in the FHN networks. (A) Left panel: The isosurfaces for the ${\bar{P}}_1$ (blue) and ${\bar{P}}_2$ (black) half-maxima in the WS-FHN networks. Above these surfaces, population 1 and 2, respectively transition to population-wide activity. These surfaces separate for weak coupling and high sortedness, bounding a regime where population 1 is active but population 2 is only partially active. Right panel: The isosurfaces for $\bar{R}$ illustrate the boundary between low and high synchrony in the WS-FHN networks. For weak coupling and high sortedness, a region of low synchrony emerges. This region encompasses the region where population 1 is active and population 2 is partially active. (B) Heatmaps for $\bar{P}$ (left panel) and $\bar{R}$ (right panel) when ${\bar{g}}_{\mathit{coup}}=0.02$ (weak). The intersection of each isosurface in A with the plane ${\bar{g}}_{\mathit{coup}}=0.02$ is drawn as a contour line. Raster plots illustrate typical activity in the different regions bounded by these contours and nodes are arranged according to location along the ring. Square shows the synchronised activity when sortedness is low. Diamond shows the terminating waves that occur when sortedness is high but drive is low. Circle shows the non-terminating waves that occur for high sortedness and intermediate drive. (C) Left panel: The isosurfaces for the ${\bar{P}}_1$ (blue) and ${\bar{P}}_2$ (black) half-maxima in the βC-FHN networks. Right panel: The isosurfaces for $\bar{R}$ illustrate the boundary between low and high synchrony in the βC-FHN networks. Low synchrony emerges across sortedness and drive when coupling is weak. For low sortedness and when drive is very close to the transition between quiescence and activity [see (D) Diamond], there exists an region of low synchrony that persists for higher values of ${\bar{g}}_{\mathit{coup}}$ (i.e., yellow bump). (D) Heatmaps for $\bar{P}$ (left panel) and $\bar{R}$ (right panel) when ${\bar{g}}_{\mathit{coup}}=0.02$ (weak). For the raster plots, nodes are ordered according to radial distance from the centre of the lattice. Square shows the fast but wave-like propagation of activity through the lattice for intermediate drive and sortedness. Diamond shows the asynchronous activity for low sortedness and low drive. Circle shows fast wave-like propagation when sortedness and drive are high.

Figure 4A (right) shows the isosurfaces $\{({\bar{g}}_{\mathit{coup}},G,a)\mid \bar{R}=\bar{r}\}$ for values between $\bar{r}=0.5$ (red) and $\bar{r}=0.7$ (yellow). These isosurfaces illustrate the boundary between low and high phase synchrony for the network as a whole. We found that for weak coupling and high $\mathcal{A}$, a region of low synchrony emerges. Figure 4B (right) shows the heatmap of $\bar{R}$ for ${\bar{g}}_{\mathit{coup}}=0.02$ and the contour curves corresponding to Figure 4A isosurfaces. This low synchrony is a result of the wave propagation that occurs in the region. Note that the region of low synchrony encapsulates the population 2 phase transition. There exists a region where population 2 is fully activated, but the synchrony is still low. Figure 4B (circle) shows an example where waves emit from the population 1 clusters, but the drive G to population 2 nodes is high enough to prevent rewiring-based termination of the waves. If G is increased further, synchrony re-emerges due to the fact that even weak (and few) long-range connections from population 1 can synchronise the primed population 2 nodes.

For the βC-FHN networks, we repeated the experiment as described above. In this case, however, the raster plots are arranged by radial distance from the centre of the lattice (Figure 4D rasters). Unlike in the case of WS-FHN networks, ${\bar{P}}_1$ and ${\bar{P}}_2$ did not separate for weak coupling. We found that across the range of ${\bar{g}}_{\mathit{coup}}$, the phase transition threshold is a decreasing function of $\mathcal{A}$ and that it represents a phase transition of the entire network [Figures 4C, D (left)]. When coupling is weak (Figure 4D shows ${\bar{g}}_{\mathit{coup}}=0.02$), the synchrony of the network decreases within the entirety of the active region. This is due to the propagation of waves from population 1 clusters, similar to the case of WS-FHN with no rewirings, however, these waves propagate quickly due to the supralinearity in the number of nodes being recruited over time. Figure 4D (square and circle) show two examples where a cluster of population 1 nodes near the centre of the lattice emit such waves. Finally, we found a region for low $\mathcal{A}$ and low G (in the active region) where asynchronous activity was observed [Figure 4D (diamond)]. This is the only region where we observed asynchrony, and it corresponds to the region where a maximal value of ${\bar{g}}_{\mathit{coup}}$ is necessary to synchronise the network (Figure 4C, see the bump in the yellow isosurface).

3.3 Weak coupling and high sortedness leads to resonance in the SRK networks

For the WS-SRK networks, we found that the intra-population isosurfaces corresponding to the phase transition between quiescence and activity separate for weak coupling and high $\mathcal{A}$ (Figure 5A left panel). The population 1 phase transition threshold (G) decreases as $\mathcal{A}$ increases for all values of ${\bar{g}}_{\mathit{coup}}\in [\mathrm{2,10}]$. For population 2, this relationship is observed for strong coupling, but for weak coupling and high sortedness, ${\bar{P}}_2$ is an increasing function of $\mathcal{A}$. Figure 5B (left) shows the heatmap of $\bar{P}$ for ${\bar{g}}_{\mathit{coup}}=2$ as well as the contour curves $\{({\bar{g}}_{\mathit{coup}}=2,G,a)\mid {\bar{P}}_k=0.5\ast \max {\bar{P}}_k\}$ for k = 1 (blue, population 1) and k = 2 (black, population 2). Within this region, we observe a region of 2:1 resonance, wherein population 1 is active and synchronised, whilst population 2 is active and synchronised but with half the frequency of population 1 (Figure 5B circle). Nodes in the SRK model produce bursts of electrical activity resulting in the slow build-up of Ca²⁺ (i.e., c). As Ca²⁺ concentration increases, hyperpolarising K⁺ channels open and terminate the burst. During the inactive period, Ca²⁺ concentration decreases and the cell becomes increasingly excitable over time. At this level of coupling, the population 2 nodes are still capable of being driven to activity by nodes in population 1. However, they also require a longer recovery time (interburst interval) leading to the observed resonance. It is worth noting that wave-like activity was not generally observed, due to the fact that our variable of interest, c, is very slow with respect to the coupled variable, V (see Figure 1). Thus, the time delays in activation between nodes occur on the order of oscillations in V (i.e., action potentials), which are very short relative to the oscillation frequency of c (i.e., bursts). This is illustrated in Supplementary Figure S4, which shows the behaviour of the WS-SRK networks when β = 0. For strong coupling, the activity of the network is still nearly synchronous despite only having local connections. For weak coupling, we did observe waves (and lowered synchrony), however, they are fast relative to the period of c dynamics (compare Supplementary Figure S4B triangle to Supplementary Figure S3B circle). Interestingly, as in the case of WS-FHN networks, we did not observe the separation in ${\bar{P}}_1$ and ${\bar{P}}_2$ when β = 0 and there was no 2:1 resonance region. We conclude, therefore, that it is the rewiring that leads to the emergence of 2:1 resonance. The lowered number of local connections in the chain at the boundaries between population 1 and population 2 clusters do have the effect of terminating a “wave” (where synchronised activity is the limit as wave speed tends to infinity), as in the case of the WS-FHN networks, however, it only occurs half the time because of the slow recovery dynamics of the SRK model. The tradeoff is that the increased number of rewiring connections allows synchronisation in the 2:1 resonance region (compare Figure 5B circle to Supplementary Figure S4B triangle).

FIGURE 5

FIGURE 5. Weak coupling in the SRK networks. (A) Left panel: The isosurfaces for the ${\bar{P}}_1$ (blue) and ${\bar{P}}_2$ (black) half-maxima in the WS-SRK networks. Above the surfaces, population 1 and 2, respectively, transition to population-wide activity. These surfaces separate for weak coupling and high sortedness, revealing a 2:1 resonance regime. Right panel: The isosurfaces for $\bar{R}$ illustrate the boundary between low and high synchrony in the WS-SRK networks. For weak coupling and high sortedness, low synchrony emerges at the transition where population 2 becomes fully active. (B) Heatmaps for $\bar{P}$ (left panel) and $\bar{R}$ (right panel) when ${\bar{g}}_{\mathit{coup}}=2$ (weak). The intersection of each isosurface in A with the plane ${\bar{g}}_{\mathit{coup}}=2$ is drawn as a contour line. Raster plots illustrate the activity and nodes are arranged according to location along the ring. Diamond shows fully synchronous activity. Circle shows the 2:1 resonance. Square shows the irregular activity in the low synchrony region. (C) Left panel: The isosurfaces for the ${\bar{P}}_1$ (blue) and ${\bar{P}}_2$ (black) half-maxima in the βC-SRK networks. Right panel: The isosurfaces for $\bar{R}$ show that synchrony is high everywhere the network is active in the βC-SRK networks. (D) Heatmaps for $\bar{P}$ (left panel) and $\bar{R}$ (right panel) when ${\bar{g}}_{\mathit{coup}}=2$ (weak). For the raster pots, nodes are ordered according to radial distance from the centre of the lattice. Square shows fully synchronous activity for low sortedness. Diamond shows fully synchronous activity for intermediate sortedness. Circle shows the 2:1 resonance.

Figure 5A (right) shows the isosurfaces $\{({\bar{g}}_{\mathit{coup}},G,a)|\bar{R}=\bar{r}\}$ for values between $\bar{r}=0.5$ (red) and $\bar{r}=0.7$ (yellow). These isosurfaces show that the region of lowered phase synchrony does not fully encompass the 2:1 resonance regime. Rather, it surrounds the ${\bar{P}}_2$ isosurface for weak coupling and high sortedness. Figure 5B (right) shows the heatmap of $\bar{R}$ for ${\bar{g}}_{\mathit{coup}}=2$ and the contour curves corresponding to Figure 5A isosurfaces. This region of lowered synchrony corresponds to the transition between the 2:1 resonance region (Figure 5B circle) and the fully synchronised region (Figure 5B diamond). Within this region, we found a region with irregular activity (Figure 5B square) with a mixture of network-wide, synchronous bursts and wave-like activity. This is due to the variability in local and non-local connections among the population 2 nodes interacting with the slow recovery of Ca²⁺ very near the population 2 phase transition. Consider the raster plot shown in Figure 5B (square) and note that the population 1 nodes fire regular, synchronous bursts (beats). With a 3:1 resonance, the local and global connections together are strong enough to elicit a (roughly) synchronous burst. On the other beats, the small number of global connections are too weak to induce synchrony (as the population 2 Ca²⁺ concentration is still recovering) and instead wave propagation or termination occurs.

For the βC-SRK networks, we repeated the experiment and found that a 2:1 resonance regime emerged in this case as well. Figures 5C, D (left panel) shows the separation in ${\bar{P}}_1$ and ${\bar{P}}_2$ for weak coupling and high $\mathcal{A}$, and Figure 5D (circle) shows an example raster plot. Unlike in the case of the WS-SRK networks, synchrony is strong everywhere, including the 2:1 resonance region and the transition to the fully active region (Figures 5C, D right panel). We did find some examples of the irregular behaviour described above (notice the slight dark patch in Figure 5D right panel at the transition boundary), but this is less prominent in the βC-SRK networks, likely due to the regularity of the graph.

4 Discussion

In this manuscript, we demonstrated how transitions to globally-coordinated activity are dependent on the degree of sortedness in population excitability. We used two prototypical models of cellular excitability where a small population was highly excitable, whilst a larger population was less excitable. This scenario has recently been hypothesised to be relevant in beta cell networks, in particular in the description of “hubs,” populations displaying high functional connectivity (Johnston et al., 2016), ‘wave-initiators’ (also known as “leaders”) (Šterk et al., 2023), and “first-responders” (Kravets et al., 2022). As the global drive to the network was increased, activity across the network transitioned from a globally inactive state to one in which subsets of nodes became active and synchronised their activity. By perturbing the spatial distribution of the highly excitable population, we showed that the drive strength at which such transitions occur is dependent on the sortedness of the network. Moreover, we highlighted the presence of other types of network solution, including partially synchronised states, slow waves, and 2:1 resonant states and explored how these depended on sortedness. These results have implications for insulin secretion in the pancreatic islets of Langerhans, and more general implications regarding transitions to synchrony and other forms of collective dynamics in networks of coupled, excitable units.

Pulsatile insulin secretion is a key component of healthy islet function and is driven by the coordinated electrical activity of β-cell populations (Bertram et al., 2007). In turn, the coordination of electrical activity is facilitated by local cell-cell communication through the gap junctional network (Ravier et al., 2005). Recently, many studies have demonstrated the existence of β-cell heterogeneity and its importance in determining network dynamics (Stožer et al., 2013; Benninger et al., 2014; Johnston et al., 2016; Westacott et al., 2017; Salem et al., 2019; Benninger and Kravets, 2021; Nasteska et al., 2021; Kravets et al., 2022; Šterk et al., 2023). Thus, a key question in this work has been to study the potential impact of the cytoarchitecture (e.g., spatial layout of heterogeneity) in modulating these dynamics. Our work suggests that it may be key to fully understanding the impact of heterogeneity on collective dynamics. We chose to focus on one form of heterogeneity, the cell-intrinsic excitability, as a case study; however, it is likely that cytoarchitecture of β-cell networks will modulate the impact of heterogeneity more generally.

We focussed on heterogeneity in cell-intrinsic excitability because of its critical role in delineating subpopulations of pancreatic beta cell networks. High excitability is a feature of both the first-responder cells in the transient phase of glucose-stimulated insulin secretion (GSIS) (Salem et al., 2019; Benninger and Kravets, 2021; Kravets et al., 2022) and the wave-initiator population during the pulsatile phase of GSIS (Benninger et al., 2014). The existence of these populations has been heavily implicated as a feature of healthy pancreatic islets, as has their critical roles in driving the collective dynamics required for GSIS. However, defining these subpopulations is complex because they likely arise as an emergent property of the network, relying on an interplay of both cell-intrinsic properties (i.e., excitability and metabolism) and network properties (i.e., coupling strength and non-uniformity of heterogeneity) (Benninger et al., 2014; Westacott et al., 2017; Salem et al., 2019; Kravets et al., 2022; Šterk et al., 2023). For example, it appears that both first-responder cells (Kravets et al., 2022) and wave-initiator cells (Šterk et al., 2023) may have weaker local connectivity, a property that would isolate them from the hyperpolarising influence of less-excitable cells and could improve network responsiveness. Although our model did not include heterogeneous coupling; we did observe that weaker coupling lowers the activation threshold of the network (Figures 4A, 5A). Likewise, we demonstrated that the spatial aggregation of highly-excitable nodes, a network property previously demonstrated by Westacott et al. (2017), can also have this effect. On the other hand, our work also illustrates potential tradeoffs in these network features. When coupling is too weak, the network can fail to synchronise and waves can fail to propagate. When sortedness is too high, excitable clusters can fail to recruit the whole network.

A variety of simplifications have been used in modelling studies to provide insight into the collective dynamics of heterogeneous coupled oscillators. For example, in specific cases of heterogeneity under the assumption of isotropic coupling and in the limit of infinitely many nodes, the network state can be mapped exactly to a low dimensional description, in which boundaries representing dynamical transitions can be found in closed form (Watanabe and Strogatz, 1993; Ott and Antonsen, 2009; Montbrió et al., 2015). For networks with a large, but not infinite number of nodes, and with general types of structure and heterogeneity, such a transformation may not be available, meaning that exact prediction of transitions between dynamical states is largely intractable. Our approach was to perform Monte-Carlo sampling across networks and initial conditions to elucidate common behaviours and trends regarding such transitions, taking advantage of the robustness of results across trials. In spite of the analytical intractability of the system, our results are consistent with the notion that the dynamics of the full network could be projected to a low dimensional manifold (e.g., the manifold representing full synchronisation of the subpopulations). In this regard, there are a number of approaches under active development, such as polynomial chaos expansion (Ghanem and Red-Horse, 2017), dynamic mode decomposition (Kutz et al., 2016), and proper orthogonal decomposition (Berkooz et al., 1993), that may facilitate construction of an approximate low dimensional system that captures the core dynamical features of the full heterogeneous network, thus allowing for a deeper exploration of network state transitions in the future.

To perform our study, we developed an algorithm that perturbs the sortedness of the network in a directed manner. The algorithm can be adapted to a range of domain geometries and network architectures, as showcased in our studies of Watts–Strogatz networks. Although our study focussed on conditions in which there are only two different populations, Section 2 discusses how our metrics can be extended to networks with more population types. Moreover, while we focussed on the case where the cells were homogeneous within a subpopulation, small changes to the algorithm could allow for the study of multimodal parameter distributions. We hope that this work will facilitate further study into the interplay of node-intrinsic heterogeneity and network features in driving the function of excitable networks, in particular beta-cell networks. To this end, future work will focus on adapting the sortedness algorithm to multiple parameters, which would allow us to further interrogate the importance of coupling strength, metabolic activity, electrical excitability, and neighbourhood on the emergence of functional roles in heterogeneous, excitable networks.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/dgalvis/network_spatial.

Author contributions

DG, DH, and KW conceived and designed the research; DG and KW developed the algorithm, model, and codebase; DG, DH, and KW interpreted the results; DG and KW prepared the figures; DG and KW drafted the manuscript; DG, DH, and KW edited and revised the manuscript.

Funding

DG acknowledges funding from the University of Birmingham Dynamic Investment Fund and the EPSRC Centre grant EP/N014391/2. DH was supported by MRC (MR/N00275X/1 and MR/S025618/1) and Diabetes UK (17/0005681) Project Grants, as well as a UKRI ERC Frontier Research Guarantee Grant (EP/X026833/1). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Starting Grant 715884 to DH). KW acknowledges funding from the MRC Fellowship MR/P01478X/1 and the Hub for Quantitative Modelling in Healthcare EP/T017856/1.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnetp.2023.1170930/full#supplementary-material

References

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

CrossRef Full Text | Google Scholar

Ashwin, P., and Swift, J. W. (1992). The dynamics of n weakly coupled identical oscillators. J. Nonlinear Sci. 2, 69–108. doi:10.1007/bf02429852

CrossRef Full Text | Google Scholar

Barrio, R., Coombes, S., Desroches, M., Fenton, F., Luther, S., and Pueyo, E. (2020). Excitable dynamics in neural and cardiac systems. Commun. Nonlinear Sci. Numer. Simul. 86, 105275. doi:10.1016/j.cnsns.2020.105275

PubMed Abstract | CrossRef Full Text | Google Scholar

Barua, A. K., and Goel, P. (2016). Isles within islets: The lattice origin of small-world networks in pancreatic tissues. Phys. D. Nonlinear Phenom. 315, 49–57. doi:10.1016/J.PHYSD.2015.07.009

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Benninger, R. K., Hutchens, T., Head, W. S., McCaughey, M. J., Zhang, M., Marchand, S. J. L., et al. (2014). Intrinsic islet heterogeneity and gap junction coupling determine spatiotemporal Ca²⁺ wave dynamics. Biophysical J. 107, 2723–2733. doi:10.1016/J.BPJ.2014.10.048

CrossRef Full Text | Google Scholar

Benninger, R. K., and Kravets, V. (2021). The physiological role of β-cell heterogeneity in pancreatic islet function. Nat. Rev. Endocrinol. 18, 0123456789–123456822. doi:10.1038/s41574-021-00568-0

CrossRef Full Text | Google Scholar

Benninger, R. K. P., Head, W. S., Zhang, M., Satin, L. S., and Piston, D. W. (2011). Gap junctions and other mechanisms of cell-cell communication regulate basal insulin secretion in the pancreatic islet. J. Physiology 589, 5453–5466. doi:10.1113/jphysiol.2011.218909

CrossRef Full Text | Google Scholar

Berkooz, G., Holmes, P., and Lumley, J. L. (1993). The proper orthogonal decomposition in the analysis of turbulent flows. Annu. Rev. Fluid Mech. 25, 539–575. doi:10.1146/annurev.fl.25.010193.002543

CrossRef Full Text | Google Scholar

Bertram, R., Satin, L. S., Pedersen, M. G., Luciani, D. S., and Sherman, A. (2007). Interaction of glycolysis and mitochondrial respiration in metabolic oscillations of pancreatic islets. Biophysical J. 92, 1544–1555. doi:10.1529/biophysj.106.097154

PubMed Abstract | CrossRef Full Text | Google Scholar

Bick, C., Goodfellow, M., Laing, C. R., and Martens, E. A. (2020). Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: A review. J. Math. Neurosci. 10, 9. doi:10.1186/s13408-020-00086-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Bittihn, P., Berg, S., Parlitz, U., and Luther, S. (2017). Emergent dynamics of spatio-temporal chaos in a heterogeneous excitable medium. Chaos 27, 093931. doi:10.1063/1.4999604

PubMed Abstract | CrossRef Full Text | Google Scholar

Börgers, C., Epstein, S., and Kopell, N. J. (2005). Background gamma rhythmicity and attention in cortical local circuits: A computational study. Proc. Natl. Acad. Sci. U. S. A. 102, 7002–7007. doi:10.1073/pnas.0502366102

PubMed Abstract | CrossRef Full Text | Google Scholar

Börgers, C., and Kopell, N. (2003). Synchronization in networks of excitatory and inhibitory neurons with sparse, random connectivity. Neural Comput. 15, 509–538. doi:10.1162/089976603321192059

PubMed Abstract | CrossRef Full Text | Google Scholar

Braun, M., Ramracheya, R., Bengtsson, M., Zhang, Q., Karanauskaite, J., Partridge, C., et al. (2008). Voltage-gated ion channels in human pancreatic beta-cells: Electrophysiological characterization and role in insulin secretion. Diabetes 57, 1618–1628. doi:10.2337/db07-0991

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabrera, O., Berman, D. M., Kenyon, N. S., Ricordi, C., Berggren, P.-O., and Caicedo, A. (2006). The unique cytoarchitecture of human pancreatic islets has implications for islet cell function. Proc. Natl. Acad. Sci. U. S. A. 103, 2334–2339. doi:10.1073/pnas.0510790103

PubMed Abstract | CrossRef Full Text | Google Scholar

De Maesschalck, P., and Wechselberger, M. (2015). Neural excitability and singular bifurcations. J. Math. Neurosci. 5, 29. doi:10.1186/s13408-015-0029-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Delgado, M. d. M., Miranda, M., Alvarez, S. J., Gurarie, E., Fagan, W. F., Penteriani, V., et al. (2018). The importance of individual variation in the dynamics of animal collective movements. Philosophical Trans. R. Soc. B Biol. Sci. 373, 20170008. doi:10.1098/rstb.2017.0008

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Dwulet, J. A. M., Ludin, N. W., Piscopio, R. A., Schleicher, W. E., Moua, O., Westacott, M. J., et al. (2019). How heterogeneity in glucokinase and gap-junction coupling determines the islet [ca2+] response. Biophysical J. 117, 2188–2203. doi:10.1016/J.BPJ.2019.10.037

CrossRef Full Text | Google Scholar

Dwulet, J. M., Briggs, J. K., and Benninger, R. K. (2021). Small subpopulations of β-cells do not drive islet oscillatory [Ca2+] dynamics via gap junction communication. PLoS Comput. Biol. 17, 10089488–e1009031. doi:10.1371/journal.pcbi.1008948

CrossRef Full Text | Google Scholar

Ermentrout, B. (2002). Simulating, analysing and animating dynamical systems: A guide to xppaut for researchers and students (siam).

Google Scholar

Fitzhugh, R. (1961). Impulses and physiological states in theoretical models of nerve membrane. Biophysical J. 1, 445–466. doi:10.1016/s0006-3495(61)86902-6

CrossRef Full Text | Google Scholar

Fretter, C., Lesne, A., Hilgetag, C. C., and Hütt, M. T. (2017). Topological determinants of self-sustained activity in a simple model of excitable dynamics on graphs. Sci. Rep. 7, 42340–42411. doi:10.1038/srep42340

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaiteri, C., and Rubin, J. E. (2011). The interaction of intrinsic dynamics and network topology in determining network burst synchrony. Front. Comput. Neurosci. 5, 10. doi:10.3389/fncom.2011.00010

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghanem, R., and Red-Horse, J. (2017). Polynomial chaos: Modeling, estimation, and approximation. Cham: Springer International Publishing, 521–551. doi:10.1007/978-3-319-12385-1_13

CrossRef Full Text | Google Scholar

Gouwens, N. W., Sorensen, S. A., Berg, J., Lee, C., Jarsky, T., Ting, J., et al. (2019). Classification of electrophysiological and morphological neuron types in the mouse visual cortex. Nat. Neurosci. 22, 1182–1195. doi:10.1038/s41593-019-0417-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Hodson, D. J., Schaeffer, M., Romanò, N., Fontanaud, P., Lafont, C., Birkenstock, J., et al. (2012). Existence of long-lasting experience-dependent plasticity in endocrine cell networks. Nat. Commun. 3, 605. doi:10.1038/ncomms1612

PubMed Abstract | CrossRef Full Text | Google Scholar

Hörning, M., Blanchard, F., Isomura, A., and Yoshikawa, K. (2017). Dynamics of spatiotemporal line defects and chaos control in complex excitable systems. Sci. Rep. 7, 7757–7810. doi:10.1038/s41598-017-08011-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Izhikevich, E. M. (2000). Neural excitability, spiking and bursting. Int. J. Bifurcation Chaos 10, 1171–1266. doi:10.1142/s0218127400000840

CrossRef Full Text | Google Scholar

Johnston, N. R., Mitchell, R. K., Haythorne, E., Trauner, D., Rutter, G. A., Hodson, D. J., et al. (2016). Beta cell hubs dictate pancreatic islet responses to Glucose. Cell Metab. 24, 389–401. doi:10.1016/j.cmet.2016.06.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Kopell, N., Kramer, M. A., Malerba, P., and Whittington, M. A. (2010). Are different rhythms good for different functions? Front. Hum. Neurosci. 4, 187–189. doi:10.3389/fnhum.2010.00187

PubMed Abstract | CrossRef Full Text | Google Scholar

Kravets, V., Dwulet, J. A. M., Schleicher, W. E., Hodson, D. J., Davis, A. M., Pyle, L., et al. (2022). Functional architecture of pancreatic islets identifies a population of first responder cells that drive the first-phase calcium response. PLOS Biol. 20, e3001761. doi:10.1371/JOURNAL.PBIO.3001761

PubMed Abstract | CrossRef Full Text | Google Scholar

Kutz, J. N., Brunton, S. L., Brunton, B., and Proctor, J. L. (2016). Dynamic mode decomposition: Data-driven modelling of complex systems. SIAM.

Google Scholar

Lambert, D., and Vanni, F. (2018). Complexity and heterogeneity in a dynamic network. Chaos, Solit. Fractals 108, 94–103. doi:10.1016/j.chaos.2018.01.024

CrossRef Full Text | Google Scholar

Lipovsek, M., Bardy, C., Cadwell, C. R., Hadley, K., Kobak, D., and Tripathy, S. J. (2021). Patch-seq: Past, present, and future. J. Neurosci. 41, 937–946. doi:10.1523/JNEUROSCI.1653-20.2020

PubMed Abstract | CrossRef Full Text | Google Scholar

Majumder, R., Feola, I., Teplenin, A. S., de Vries, A. A., Panfilov, A. V., and Pijnappels, D. A. (2018). Optogenetics enables real-time spatiotemporal control over spiral wave dynamics in an excitable cardiac system. eLife 7, 410766–e41155. doi:10.7554/eLife.41076

CrossRef Full Text | Google Scholar

Manchanda, K., Bose, A., and Ramaswamy, R. (2017). Collective dynamics in heterogeneous networks of neuronal cellular automata. Phys. A Stat. Mech. its Appl. 487, 111–124. doi:10.1016/j.physa.2017.06.021

CrossRef Full Text | Google Scholar

Mathiesen, J., Angheluta, L., Ahlgren, P. T., and Jensen, M. H. (2013). Excitable human dynamics driven by extrinsic events in massive communities. Proc. Natl. Acad. Sci. U. S. A. 110, 17259–17262. doi:10.1073/pnas.1304179110

PubMed Abstract | CrossRef Full Text | Google Scholar

McKenna, J. P., Ha, J., Merrins, M. J., Satin, L. S., Sherman, A., and Bertram, R. (2016). Ca2+ effects on ATP production and consumption have regulatory roles on oscillatory islet activity. Biophysical J. 110, 733–742. doi:10.1016/j.bpj.2015.11.3526

CrossRef Full Text | Google Scholar

Montbrió, E., Pazó, D., and Roxin, A. (2015). Macroscopic description for networks of spiking neurons. Phys. Rev. X 5, 021028–21115. doi:10.1103/PhysRevX.5.021028

CrossRef Full Text | Google Scholar

Nagumo, J., Arimoto, S., and Yoshizawa, S. (1962). An active pulse transmission line simulating nerve axon. Proc. IEEE 117, 2061–2070. doi:10.1109/jrproc.1962.288235

CrossRef Full Text | Google Scholar

Nasteska, D., Fine, N. H., Ashford, F. B., Cuozzo, F., Viloria, K., Smith, G., et al. (2021). Pdx1low mafalow beta-cells contribute to islet function and insulin release. Nat. Commun. 12, 674–719. doi:10.1038/s41467-020-20632-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Ott, E., and Antonsen, T. M. (2009). Long time evolution of phase oscillator systems. Chaos 19, 023117. doi:10.1063/1.3136851

PubMed Abstract | CrossRef Full Text | Google Scholar

Ottino-Löffler, B., and Strogatz, S. H. (2016). Kuramoto model with uniformly spaced frequencies: Finite- N asymptotics of the locking threshold. Phys. Rev. E 93, 062220. doi:10.1103/PhysRevE.93.062220

PubMed Abstract | CrossRef Full Text | Google Scholar

Pedersen, M. G., Brøns, M., and Sørensen, M. P. (2022). Amplitude-modulated spiking as a novel route to bursting: Coupling-induced mixed-mode oscillations by symmetry breaking. Chaos Interdiscip. J. Nonlinear Sci. 32, 013121. doi:10.1063/5.0072497

CrossRef Full Text | Google Scholar

Podobnik, B., Korošak, D., Skelin Klemen, M., Stožer, A., Dolenšek, J., Slak Rupnik, M., et al. (2020). β cells operate collectively to help maintain glucose homeostasis. Biophysical J. 118, 2588–2595. doi:10.1016/j.bpj.2020.04.005

CrossRef Full Text | Google Scholar

Punckt, C., Bodega, P. S., Kaira, P., and Rotermund, H. H. (2015). Wildfires in the lab: Simple experiment and models for the exploration of excitable dynamics. J. Chem. Educ. 92, 1330–1337. doi:10.1021/ed500714f

CrossRef Full Text | Google Scholar

Ravier, M. A., ldenagel, M. G., Charollais, A., Gjinovci, A., Caille, D., hl, G. S., et al. (2005). Loss of connexin36 channels alters beta-cell coupling, islet synchronization of glucose-induced Ca2+ and insulin oscillations, and basal insulin release. DIABETES 54, 1798–1807. doi:10.2337/diabetes.54.6.1798

PubMed Abstract | CrossRef Full Text | Google Scholar

Rorsman, P., and Braun, M. (2013). Regulation of insulin secretion in human pancreatic islets. Annu. Rev. Physiology 75, 155–179. doi:10.1146/annurev-physiol-030212-183754

PubMed Abstract | CrossRef Full Text | Google Scholar

Salem, V., Silva, L., Suba, K., Georgiadou, E., Neda Mousavy Gharavy, S., Akhtar, N., et al. (2019). Leader β-cells coordinate ca2+ dynamics across pancreatic islets in vivo. Nat. Metab. 1, 615–629. doi:10.1038/s42255-019-0075-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanchez-Cardenas, C., Fontanaud, P., He, Z., Lafont, C., Meunier, A. C., Schaeffer, M., et al. (2010). Pituitary growth hormone network responses are sexually dimorphic and regulated by gonadal steroids in adulthood. Proc. Natl. Acad. Sci. U. S. A. 107, 21878–21883. doi:10.1073/pnas.1010849107

PubMed Abstract | CrossRef Full Text | Google Scholar

Satin, L. S., Zhang, Q., and Rorsman, P. (2020). Take Me to Your Leader”: An electrophysiological appraisal of the role of hub cells in pancreatic islets. Diabetes 69, 830–836. doi:10.2337/dbi19-0012

PubMed Abstract | CrossRef Full Text | Google Scholar

Scialla, S., Loppini, A., Patriarca, M., and Heinsalu, E. (2021). Hubs, diversity, and synchronization in fitzhugh-nagumo oscillator networks: Resonance effects and biophysical implications. Phys. Rev. E 103, 052211. doi:10.1103/PhysRevE.103.052211

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Šterk, M., Dolenšek, J., Klemen, M. S., Bombek, L. K., Leitgeb, E. P., Kerčmar, J., et al. (2023). Functional characteristics of hub and wave initiator cells in beta cell networks. Biophysical J. 122, 784–801. doi:10.1016/j.bpj.2023.01.039

CrossRef Full Text | Google Scholar

Stožer, A., Gosak, M., Dolenšek, J., Perc, M., Marhl, M., Rupnik, M. S., et al. (2013). Functional connectivity in islets of langerhans from mouse pancreas tissue slices. PLoS Comput. Biol. 9, 1002923. doi:10.1371/journal.pcbi.1002923

CrossRef Full Text | Google Scholar

Terrien, S., Pammi, V. A., Broderick, N. G. R., Braive, R., Beaudoin, G., Sagnes, I., et al. (2020). Equalization of pulse timings in an excitable microlaser system with delay. Phys. Rev. Res. 2, 023012. doi:10.1103/physrevresearch.2.023012

CrossRef Full Text | Google Scholar

Terrien, S., Pammi, V. A., Krauskopf, B., Broderick, N. G., and Barbay, S. (2021). Pulse-timing symmetry breaking in an excitable optical system with delay. Phys. Rev. E 103, 012210. doi:10.1103/PhysRevE.103.012210

PubMed Abstract | CrossRef Full Text | Google Scholar

Vannucchi, F., and Boccaletti, S. (2004). Chaotic spreading of epidemics in complex networks of excitable units. Math. Biosci. Eng. 1, 49–55. doi:10.3934/mbe.2004.1.49

PubMed Abstract | CrossRef Full Text | Google Scholar

Watanabe, S., and Strogatz, S. H. (1993). Integrability of a globally coupled oscillator array. Phys. Rev. Lett. 70, 2391–2394. doi:10.1103/PhysRevLett.70.2391

PubMed Abstract | CrossRef Full Text | Google Scholar

Watts, D. J., and Strogatz, S. H. (1998). Collective dynamics of ‘small-world’ networks. Nature 393, 440–442. doi:10.1038/30918

PubMed Abstract | CrossRef Full Text | Google Scholar

Wedgwood, K. C., Słowiński, P., Manson, J., Tsaneva-Atanasova, K., and Krauskopf, B. (2021). Robust spike timing in an excitable cell with delayed feedback. J. R. Soc. Interface 18, 20210029. doi:10.1098/rsif.2021.0029

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Wittmeier, S., Song, G., Duffin, J., and Poon, C. S. (2008). Pacemakers handshake synchronization mechanism of mammalian respiratory rhythmogenesis. Proc. Natl. Acad. Sci. U. S. A. 105, 18000–18005. doi:10.1073/pnas.0809377105

PubMed Abstract | CrossRef Full Text | Google Scholar