GABA Regulates Ca²⁺ Oscillations and Synchronization in Pancreatic Beta Cells

1. Introduction

Pancreatic beta cells play a central role in glucose homeostasis by coupling metabolic stimulation to insulin secretion, and impairment of this function is a major determinant of diabetes mellitus. In beta cells, nutrient metabolism, ATP production, electrical activity, Ca²⁺ influx, and intracellular Ca²⁺ oscillations are tightly interconnected components of stimulus–secretion coupling. Because pulsatile and coordinated beta-cell activity is essential for efficient insulin release, identifying endogenous regulators of beta-cell Ca²⁺ dynamics remains important for understanding both normal islet physiology and beta-cell dysfunction in diabetes [1].

Among such regulators, gamma-aminobutyric acid (GABA) has attracted attention for several decades. Early work demonstrated that pancreatic beta cells contain unusually high concentrations of GABA, establishing the endocrine pancreas as a major non-neuronal site of GABA accumulation [2]. This concept was extended by the demonstration that glutamic acid decarboxylase (GAD) and GABA colocalize with synaptic-like microvesicles in beta cells, suggesting specialized machinery for GABA storage and secretion [3]. Subsequently, regulated exocytosis of GABA-containing synaptic-like microvesicles was directly demonstrated in pancreatic beta cells [4], and GABA was later shown to function as an autocrine excitatory transmitter in human beta cells through GABA_A receptor-dependent signaling [5]. Together, these studies established that GABA is not merely a by-product of amino acid metabolism, but an endogenous signaling molecule with functional relevance in pancreatic islets.

In parallel with this signaling perspective, a substantial body of work emphasized the metabolic role of GABA in beta-cell physiology. Broader analyses of amino acid metabolism highlighted the importance of mitochondrial coupling and non-glucose metabolic pathways in insulin secretion [6]. This view was strengthened by the finding that glucose suppresses GABA release from pancreatic beta cells through increased GABA shunt activity [7]. Shortly thereafter, it was shown that glucose-promoted GABA metabolism contributes to insulin secretion and is associated with changes in ATP content and ATP/ADP ratio in beta cells [8]. Additional mechanistic work on glutamate dehydrogenase and amino acid transporters further supported the idea that glutamate/GABA-related pathways are functionally integrated with beta-cell metabolism and secretory activity [9,10]. This metabolic interpretation was explicitly developed in later reviews that argued for a specific role of GABA metabolism and the GABA shunt in sustained insulin secretion and beta-cell bioenergetics [11,12,13].

More recent work brought GABA into the context of oscillatory and collective beta-cell behavior. A major advance came with the demonstration that human beta cells release GABA from cytosolic pools in a pulsatile manner and that this release imposes a synchronizing rhythm on pulsatile insulin secretion, with a critical contribution of the volume-regulated anion channel (VRAC) [14]. More recent reviews have synthesized growing evidence that islet GABA signaling combines receptor-mediated auto/paracrine actions with metabolic coupling through the GABA shunt, particularly in human islets [15,16]. Most importantly for the present study, Ferreira et al. [17] showed that beta-cell-specific loss of GAD65 and GAD67 abolishes endogenous islet GABA and results in abnormal Ca²⁺ oscillations characterized by altered active-phase duration and reduced amplitudes, together with impaired insulin secretion. These findings strongly support the view that endogenous beta-cell-derived GABA is an important determinant of islet Ca²⁺ rhythm generation and coordinated activity.

Despite these advances, the mechanistic basis by which GABA regulates beta-cell Ca²⁺ oscillations remains incompletely understood. In particular, it is still unresolved how GABA-dependent metabolic effects and GABA-dependent modulation of intercellular signaling are translated into quantitative changes in oscillatory frequency, amplitude, entrainment, and the eventual synchronization of beta-cell populations. This unresolved issue is especially important because coordinated oscillatory activity at the islet level is essential for robust pulsatile insulin secretion. In our recent work, we proposed the Dual Anaplerotic Model (DAM), in which the GABA shunt was incorporated as an anaplerotic component of beta-cell metabolism, functionally linked to oscillatory redistribution of mitochondrial and cytosolic metabolites and to ATP generation [18]. That framework now provides a basis for extending the model toward explicit Ca²⁺ dynamics and for interpreting recent experimental observations on GABA-dependent beta-cell behavior.

In the present study, we use an extended mathematical model to investigate how GABA links intracellular metabolism to Ca²⁺ oscillations and coordinated beta-cell dynamics. Building on the DAM framework, we incorporate explicit dynamics of cytosolic Ca²⁺, endoplasmic reticulum Ca²⁺, ATP, and a regulatory variable governing Ca²⁺ influx, while representing extracellular GABA as a delayed interstitial signal that feeds back onto cellular excitability. This formulation enables us to address two interconnected questions. First, we ask how altered GABA production affects the amplitude, frequency, and temporal organization of Ca²⁺ oscillations in a single beta cell, and which underlying mechanisms dominate this response. Second, we ask whether the same GABA-dependent processes can provide an effective coupling pathway between non-identical beta cells and thereby promote entrainment and synchronization.

By combining explicit ATP- Ca²⁺ dynamics with metabolically grounded GABA dynamics and delayed intercellular GABA signaling, the model goes beyond previous conceptual descriptions of GABA action and provides a unified mechanistic framework for interpreting both single-cell and collective beta-cell behavior. In this way, the study not only offers a mechanistic explanation for the experimentally observed effects of GABA deficiency on Ca²⁺ oscillations [17], but also provides a plausible theoretical basis for the synchronizing role of pulsatile beta-cell-derived GABA within the islet [14,15]. Overall, our findings support the view that GABA is an important regulator of beta-cell oscillatory dynamics and collective activity, acting through a functional link between mitochondrial metabolism, ATP production, Ca²⁺ influx, and intercellular phase coordination.

2. Model

The mathematical model developed in this study is based on established frameworks of intracellular Ca²⁺ dynamics and incorporates key mechanisms described in previous models, such as the Integrated Oscillator Model (IOM) [19]. These mechanisms are reduced to a minimal formulation that preserves the essential oscillatory behavior arising from Ca²⁺ exchange between the cytosol and the endoplasmic reticulum, in line with classical minimal models of calcium oscillations [20,21,22]. In addition, the model explicitly incorporates coupling between Ca²⁺ and metabolic dynamics, thereby reproducing the characteristic phase relationship between Ca²⁺ and metabolic/ATP oscillations reported in previous studies [19,22,23]. Within this framework, the primary aim of the model is to introduce GABA-dependent processes in order to investigate how GABA modulates Ca²⁺ oscillations, their phase relationships, and synchronization in pancreatic beta cells.

We first introduce a Single-Cell Model, which describes intracellular Ca²⁺ dynamics within an individual beta cell. The model integrates metabolic and electrophysiological components, with particular emphasis on the role of GABA in mitochondrial ATP production via the anaplerotic GABA shunt. In addition, GABA release is incorporated as a dynamic process that introduces delayed negative feedback on Ca²⁺ influx, thereby modulating the temporal characteristics of intracellular Ca²⁺ oscillations.

The framework is subsequently extended to a Two-Cell Coupled Model, in which individual cellular oscillators are interconnected through GABA-mediated coupling. This formulation enables us to examine how the coupling between oscillatory Ca²⁺ dynamics and oscillatory GABA signaling at the single-cell level can facilitate phase alignment and promote synchronization of Ca²⁺ oscillations across cells.

2.1. Single-Cell Model

The cellular processes governing intracellular Ca²⁺ oscillations in a single β-cell are schematically illustrated in Figure 1. The model includes Ca²⁺ fluxes between the cytosol and the endoplasmic reticulum (ER), as well as Ca²⁺ influx across the plasma membrane through voltage-gated Ca²⁺ channels (VGCCs). VGCC opening is regulated by the membrane potential, such that depolarization promotes Ca²⁺ entry. Membrane depolarization is, in turn, favored by closure of ATP-sensitive K⁺ channels (K_ATP channels), driven by an increase in local ATP concentration near the plasma membrane. This local ATP signal is associated with PEP-cycle activity and ATP delivery to the microdomain surrounding K_ATP channels, whereas mitochondrial oxidative phosphorylation provides the principal energetic background for this process [23,24]. Within this framework, ATP production is additionally influenced by the GABA shunt, which acts as an anaplerotic pathway [18] feeding succinate into the TCA cycle and thereby reinforcing mitochondrial oxidative metabolism. In this manner, GABA metabolism contributes to ATP generation and functionally couples mitochondrial metabolism to the regulation of Ca²⁺ dynamics. The model further includes net cytosolic GABA production (${GABA}_{\mathrm{c}\mathrm{y}\mathrm{t}}$), as described in detail in our previous study [18], as well as its release into the interstitial space through VRAC, thereby generating the interstitial GABA signal $GAB{A}_{\mathrm{i}\mathrm{s}}$. Interstitial GABA is assumed to inhibit Ca²⁺ influx through mechanisms involving both ${\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}_{\mathrm{A}}$ and ${\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}_{\mathrm{B}}$ receptors [17], as will be described in more detail below. In effective form, these pathways are represented in the model as a delayed negative-feedback loop through which GABA release suppresses Ca²⁺ entry and modulates the temporal characteristics of intracellular Ca²⁺ oscillations.

2.1.1. Ca²⁺ Dynamics

Intracellular Ca²⁺ dynamics is described by the following two differential equations:

$${\displaystyle \frac{d{Ca}_{\mathrm{c}\mathrm{y}\mathrm{t}}}{dt}}=J_{\mathrm{i}\mathrm{n}}-J_{\mathrm{o}\mathrm{u}\mathrm{t}}+J_{\mathrm{E}\mathrm{R},\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}-J_{\mathrm{E}\mathrm{R},\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}+J_{\mathrm{E}\mathrm{R},\mathrm{c}\mathrm{h}},$$

(1)

$${\displaystyle \frac{d{Ca}_{\mathrm{E}\mathrm{R}}}{dt}}=J_{\mathrm{E}\mathrm{R},\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}-J_{\mathrm{E}\mathrm{R},\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}-J_{\mathrm{E}\mathrm{R},\mathrm{c}\mathrm{h}} ,$$

(2)

which represent the temporal evolution of Ca²⁺ concentration in the cytosol (${Ca}_{\mathrm{c}\mathrm{y}\mathrm{t}}$) and in the endoplasmic reticulum (${Ca}_{\mathrm{E}\mathrm{R}}$), respectively. The model accounts for Ca²⁺ exchange between the cytosol and the ER, as well as Ca²⁺ fluxes across the plasma membrane.

In the following, the individual flux terms are described in more detail. The values of the model constants given below represent reference values in arbitrary units.

Within this framework, the flux $J_{\mathrm{E}\mathrm{R},\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}$ represents the active uptake of Ca²⁺ from the cytosol into the ER, mediated by the sarco/endoplasmic reticulum Ca²⁺-ATPase (SERCA). In the present model, this flux is described by

$$J_{\mathrm{E}\mathrm{R},\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}=k_{\mathrm{E}\mathrm{R},\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}\cdot C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}} ,$$

(3)

where $k_{\mathrm{E}\mathrm{R},\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}=4$, assuming a linear dependence of pump activity on cytosolic Ca²⁺ concentration. In this formulation, the dependence of the flux on ATP concentration is neglected, as its variations in the cytosol are relatively small. This simplification is justified by experimental evidence showing that cytosolic ATP levels exhibit relatively small fluctuations compared to those in the submembrane compartment [25] and in the microdomain associated with ${\mathrm{K}}_{\mathrm{A}\mathrm{T}\mathrm{P}}$ channels, where local ATP dynamics are sufficient to regulate channel permeability [23]. Instead, the model explicitly incorporates Ca²⁺ dependence, as SERCA activity is tightly regulated by cytosolic Ca²⁺ levels. An increase in $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ enhances Ca²⁺ binding to the pump, thereby increasing its turnover rate and facilitating Ca²⁺ sequestration into the ER. Capturing this dependence is essential for accurately reproducing intracellular Ca²⁺ dynamics, particularly the feedback mechanisms governing Ca²⁺ oscillations and ER loading.

The flux $J_{\mathrm{E}\mathrm{R},\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}$ describes passive Ca²⁺ leak from the ER back into the cytosol, reflecting the basal permeability of the ER membrane. In the present model, this flux is described by

$$J_{\mathrm{E}\mathrm{R},\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}=k_{\mathrm{E}\mathrm{R},\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}\cdot C{a}_{\mathrm{E}\mathrm{R}} ,$$

(4)

where $k_{\mathrm{E}\mathrm{R},\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{k}}=0.5$, and the leak rate is assumed to increase linearly with the ER luminal Ca²⁺ concentration. This assumption reflects the fact that passive Ca²⁺ efflux from the ER depends on the ER Ca²⁺ load, with higher $C{a}_{\mathrm{E}\mathrm{R}}$ generating a larger concentration gradient and thereby a greater leak flux into the cytosol. Notably, the dependence of the flux on the concentration difference between the ER and the cytosol is not modeled explicitly. This simplification is justified by the fact that Ca²⁺ concentration in the ER is several orders of magnitude higher than in the cytosol. Consequently, the concentration gradient effectively follows variations in $C{a}_{\mathrm{E}\mathrm{R}}$, which supports the approximation of the leak flux as being proportional to $C{a}_{\mathrm{E}\mathrm{R}}$.

The term $J_{\mathrm{E}\mathrm{R},\mathrm{c}\mathrm{h}}$ represents Ca²⁺ release from the ER through Ca²⁺-permeable ER channels and accounts for channel-mediated mobilization of stored Ca²⁺ into the cytosol. In the present model, this flux is described by

$$J_{\mathrm{E}\mathrm{R},\mathrm{c}\mathrm{h}}=k_{\mathrm{E}\mathrm{R},\mathrm{c}\mathrm{h}}\cdot {\displaystyle \frac{C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}^{n_{\mathrm{c}\mathrm{h}}}}{K_{\mathrm{c}\mathrm{h}}^{n_{\mathrm{c}\mathrm{h}}}+C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}^{n_{\mathrm{c}\mathrm{h}}}}}\cdot C{a}_{\mathrm{E}\mathrm{R}} ,$$

(5)

with $k_{\mathrm{E}\mathrm{R},\mathrm{c}\mathrm{h}}=1$, $K_{\mathrm{c}\mathrm{h}}=0.4$, and $n_{\mathrm{c}\mathrm{h}}=4$. This formulation assumes a sigmoidal dependence of channel activation on cytosolic Ca²⁺ concentration, reflecting the cooperative nature of Ca²⁺-induced Ca²⁺ release (CICR). The multiplicative dependence on ${Ca}_{\mathrm{E}\mathrm{R}}$ accounts for the fact that the magnitude of Ca²⁺ release is also governed by the ER Ca²⁺ load, with higher luminal Ca²⁺ levels providing a stronger driving force for Ca²⁺ efflux.

In addition to intracellular Ca²⁺ exchange, the model also includes Ca²⁺ transport across the plasma membrane, represented by the fluxes $J_{\mathrm{i}\mathrm{n}}$ and $J_{\mathrm{o}\mathrm{u}\mathrm{t}}$. The efflux term $J_{\mathrm{o}\mathrm{u}\mathrm{t}}$ describes Ca²⁺ extrusion from the cytosol to the extracellular space via plasma membrane Ca²⁺ clearance mechanisms, predominantly Ca²⁺-ATPases and related outward transport processes. In the present model, this flux is described by

$$J_{\mathrm{o}\mathrm{u}\mathrm{t}}=k_{\mathrm{o}\mathrm{u}\mathrm{t}}\cdot C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}} ,$$

(6)

with $k_{\mathrm{o}\mathrm{u}\mathrm{t}}=1$, assuming a linear dependence of Ca²⁺ extrusion on cytosolic Ca²⁺ concentration. The proportional dependence on $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ captures the first-order approximation of Ca²⁺ extrusion kinetics, where the rate of transport increases with substrate availability. This is consistent with the behavior of plasma membrane Ca²⁺ pumps and exchangers, which respond to elevated cytosolic Ca²⁺ by enhancing Ca²⁺ clearance, thereby contributing to the restoration of basal intracellular Ca²⁺ levels.

The influx term $J_{\mathrm{i}\mathrm{n}}$ denotes Ca²⁺ entry from the extracellular space into the cytosol, primarily through voltage-gated Ca²⁺ channels (VGCCs). In the present model, this flux is described by

$$J_{\mathrm{i}\mathrm{n}}=g_{\mathrm{i}\mathrm{n}} x {\displaystyle \frac{K_{\mathrm{i}\mathrm{n}}}{K_{\mathrm{i}\mathrm{n}}+C_{\mathrm{c}\mathrm{y}\mathrm{t}}}},$$

(7)

where $g_{\mathrm{i}\mathrm{n}}$ represents the maximum influx rate modulated by intracellular signaling, and $x$ denotes the channel activation variable. In the model, we additionally include an inhibitory dependence of the influx term $J_{\mathrm{i}\mathrm{n}}$ on $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$, thereby phenomenologically representing calcium-dependent inactivation (CDI) of voltage-gated Ca²⁺ channels. In this way, an increase in intracellular Ca²⁺ reduces channel availability and thereby limits further Ca²⁺ influx [26,27]. Here, $K_{\mathrm{i}\mathrm{n}}=0.2$ denotes the half-saturation constant controlling the dependence of the influx on cytosolic Ca²⁺.

The parameter $g_{\mathrm{i}\mathrm{n}}$ is further defined as

$$g_{\mathrm{i}\mathrm{n}}=g_{\mathrm{i}\mathrm{n},0}-k_{\mathrm{i}\mathrm{n},\mathrm{G}}\cdot GAB{A}_{\mathrm{i}\mathrm{s}} ,$$

(8)

where $g_{\mathrm{i}\mathrm{n},0}=8$ is the baseline maximal ${\mathrm{C}\mathrm{a}}^{2+}$ influx rate in the absence of inhibitory GABAergic input, and $k_{\mathrm{i}\mathrm{n},\mathrm{G}}=500$ quantifies the strength by which interstitial GABA suppresses ${\mathrm{C}\mathrm{a}}^{2+}$ influx. Here, $GAB{A}_{\mathrm{i}\mathrm{s}}$ denotes the interstitial GABA concentration, which follows the intracellular GABA concentration with a time delay, as described in more detail below. In this way, we assume that the delayed interstitial GABA signal modulates Ca²⁺ influx via a rapid inhibitory feedback loop. This term provides a phenomenological description of the fast ionotropic action of extracellular GABA on beta-cell excitability and membrane-potential-dependent Ca²⁺ entry, consistent with the current view that GABA released from beta cells acts through both GABA_A and GABA_B receptors and that exogenous GABA suppresses islet Ca²⁺ oscillations through a mechanism involving both receptor classes [15,17].

In modeling the influx $J_{\mathrm{i}\mathrm{n}}$ (Eq. 7), we additionally introduce a slower inhibitory feedback loop via the regulatory variable $x$. The gating variable $x$follows first-order kinetics,

$${\displaystyle \frac{dx}{dt}}={\displaystyle \frac{x_{\mathrm{\infty },\mathrm{G}}-x}{{\tau }_x}},$$

(9)

where ${\tau }_x=1$ determines the timescale of adaptation of the gating process. The steady-state value ${\mathrm{x}}_{\mathrm{\infty },\mathrm{G}}$ depends on both ATP and GABA levels,

$$x_{\infty ,\mathrm{G}}=\left\{\begin{array}{cc}{x}_{\infty ,0}-k_{x,\mathrm{G}} GAB{A}_{\mathrm{i}\mathrm{s}};& ATP>AT{P}_{\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{n}}\\ 0;& ATP<AT{P}_{\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{n}}\end{array} \right.,$$

(10)

where $x_{\infty ,0}=1$ is the baseline steady-state activation level in the absence of GABA, and $k_{x,\mathrm{G}}=100$ determines the strength of GABA-dependent inhibition of channel activation. Since $x$ positively regulates the influx $J_{\mathrm{i}\mathrm{n}}$ (Eq. 7), a reduction in $x_{\mathrm{\infty },\mathrm{G}}$ mediated by $GAB{A}_{\mathrm{i}\mathrm{s}}$ (Eq. 10) implies that delayed interstitial GABA gradually reduces the effective drive for Ca²⁺ influx on a slower timescale. Physiologically, this slow inhibitory component is intended to phenomenologically represent a metabotropic ${\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}_{\mathrm{B}}$-like pathway acting on slower excitability and adaptation processes, consistent with inhibitory G-protein (Gi/o)-dependent inhibition of adenylyl cyclase and the resulting reduction in Ca²⁺ influx [15]. Together, these terms (Eqs. 7–10) define a physiologically plausible delayed negative feedback loop through which interstitial GABA ($GAB{A}_{\mathrm{i}\mathrm{s}}$) can suppress excessive Ca²⁺ entry [17]. The parameter $AT{P}_{\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{n}}=0.65$ represents the threshold ATP concentration above which Ca²⁺ channels can become active (open), whereas below this threshold the channels remain closed.

2.1.2. ATP Dynamics and Coupling to the DAM

In the present model, ATP dynamics is described at the level of the submembrane ATP pool, rather than bulk cytosolic ATP, because ATP oscillations in the submembrane region are substantially more pronounced and are more directly relevant for the regulation of Ca²⁺ influx and K_ATP-dependent excitability than ATP in the bulk cytosol [23]. The temporal evolution of ATP is therefore described by:

$${\displaystyle \frac{dATP}{dt}}=J_{\mathrm{A}\mathrm{T}\mathrm{P},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}}-J_{\mathrm{A}\mathrm{T}\mathrm{P},\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}-J_{\mathrm{A}\mathrm{T}\mathrm{P},\mathrm{u}\mathrm{s}\mathrm{e}},$$

(11)

where $J_{\mathrm{A}\mathrm{T}\mathrm{P},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}}$ denotes ATP production, $J_{\mathrm{A}\mathrm{T}\mathrm{P},\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}$ represents ATP consumption by Ca²⁺-ATPases, and $J_{\mathrm{A}\mathrm{T}\mathrm{P},\mathrm{u}\mathrm{s}\mathrm{e}}$ accounts for other ATP-consuming cellular processes.

The ATP production term is written as

$$J_{\mathrm{A}\mathrm{T}\mathrm{P},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}}={\displaystyle \frac{v_{\mathrm{P}\mathrm{E}\mathrm{P}}}{1+{\left(\frac{{Ca}_{\mathrm{c}\mathrm{y}\mathrm{t}}}{K_{\mathrm{C}\mathrm{a}}}\right)}^n}}+v_{\mathrm{T}\mathrm{C}\mathrm{A}}·{Ca}_{\mathrm{c}\mathrm{y}\mathrm{t}},$$

(12)

where $K_{\mathrm{C}\mathrm{a}}=0.2$ is the half-saturation constant controlling the inhibitory effect of cytosolic Ca²⁺ on the PEP-dependent component of ATP production, and $n=4$ is the corresponding Hill coefficient determining the steepness of this dependence. Equation (12) separates ATP production into two physiologically distinct components. The first term describes ATP production associated with the cataplerotic phase and the associated PEP-cycle activity, whereas the second term represents the oxidative mitochondrial contribution associated with TCA-cycle activity.

The inverse dependence of the first term on $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ is an essential feature of the model. It reflects the idea that ATP production via the PEP cycle is most effective when cytosolic Ca²⁺ is low, i.e., during the cataplerotic phase, when carbon is preferentially routed through pyruvate carboxylase (PC) rather than pyruvate dehydrogenase (PDH), and when PEP cycling becomes an important source of ATP delivery to the submembrane region [18,23]. In this regime, the availability of carbon skeletons for PEP cycling depends strongly on the anaplerotic input provided by the GABA shunt. Within the DAM framework [18], this cataplerotic contribution is associated with the flux $J_{13}$, which describes carbon redistribution toward the GABA-related pool during the cataplerotic phase of the oscillation cycle. In the present formulation, this flux is identified with the effective GABA production flux, $J_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}}\equiv J_{13}$. Accordingly, the effective strength of the PEP-dependent ATP-producing branch is written as

$$v_{\mathrm{P}\mathrm{E}\mathrm{P}}=v_{\mathrm{P}\mathrm{E}\mathrm{P},0}+k_{\mathrm{P}\mathrm{E}\mathrm{P},\mathrm{G}}\cdot J_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}},$$

(13)

where $v_{\mathrm{P}\mathrm{E}\mathrm{P},0}=0.1$ is the basal contribution of the PEP-related ATP-producing component, and $k_{\mathrm{P}\mathrm{E}\mathrm{P},\mathrm{G}}=1$ determines how strongly this component is enhanced by the cataplerotic flux $J_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}}=J_{13}$. In this way, the model links GABA-dependent cataplerotic carbon redistribution directly to the PEP-cycle-associated ATP production that is most effective during the low- Ca²⁺ phase of the oscillation.

The second term in Eq. (12), $v_{\mathrm{T}\mathrm{C}\mathrm{A}}\cdot C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$, represents the oxidative component of ATP production. This term captures the fact that mitochondrial ATP synthesis in the submembrane region is also stimulated during the oxidative phase, when elevated cytosolic Ca²⁺ promotes mitochondrial metabolism, in part through activation of PDH and enhanced TCA-cycle turnover. In the DAM framework, this oxidative contribution is linked to the GABA shunt flux $J_{32}$, which describes the return of carbon from the GABA-related pool back into the left part of the TCA cycle, thereby reinforcing oxidative ATP production during the phase in which fresh carbon enters the cycle through GABA-shunt-dependent anaplerosis [18]. In the present formulation, this flux is identified with the effective GABA-shunt contribution, $J_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{s}\mathrm{h}\mathrm{u}\mathrm{n}\mathrm{t}}\equiv J_{32}$. We therefore write

$$v_{\mathrm{T}\mathrm{C}\mathrm{A}}=v_{\mathrm{T}\mathrm{C}\mathrm{A},0}+k_{\mathrm{T}\mathrm{C}\mathrm{A},\mathrm{G}}\cdot J_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{s}\mathrm{h}\mathrm{u}\mathrm{n}\mathrm{t}},$$

(14)

where $v_{\mathrm{T}\mathrm{C}\mathrm{A},0}=1$ is the basal oxidative contribution to ATP production, and $k_{\mathrm{T}\mathrm{C}\mathrm{A},\mathrm{G}}=2$ quantifies the extent to which this component is increased by the flux $J_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{s}\mathrm{h}\mathrm{u}\mathrm{n}\mathrm{t}}=J_{32}$. Thus, the model explicitly links GABA-shunt-mediated return of carbon into the TCA cycle to the oxidative ATP-producing branch that predominates during the high- Ca²⁺ phase of the oscillation.

A key point is that the fluxes $J_{13}$ and $J_{32}$, as well as the dynamics of cytosolic GABA, are not introduced here as ad hoc functions, but are taken directly from the previously published DAM formulation [18]. In the original DAM study, the time courses of $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ and ATP were prescribed from experimental fits, whereas the remaining metabolic pools and inter-pool fluxes were modeled explicitly. In the present study, this logic is extended one step further: the same DAM subsystem is retained for the calculation of $J_{13}$, $J_{32}$, and $GAB{A}_{\mathrm{c}\mathrm{y}\mathrm{t}}$, but $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ and ATP are now computed self-consistently from the differential equations of the present model. Thus, the DAM provides the internal metabolic structure required to determine how GABA-shunt-dependent carbon redistribution feeds back onto ATP production and, indirectly, onto Ca²⁺ dynamics.

ATP consumption by Ca²⁺-transporting ATPases is described by

$$J_{\mathrm{A}\mathrm{T}\mathrm{P},\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}=k_{\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}\cdot {Ca}_{\mathrm{c}\mathrm{y}\mathrm{t}},$$

(15)

where $k_{\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}=1.5$ denotes the effective rate constant for ATP consumption by Ca²⁺-ATPases. This term represents the ATP cost of Ca²⁺ sequestration and extrusion. In this effective formulation, ATP consumption is assumed to increase linearly with cytosolic Ca²⁺, reflecting the enhanced activity of Ca²⁺-ATPases when intracellular Ca²⁺ levels are elevated.

Finally, we include a further ATP consumption term,

$$J_{\mathrm{A}\mathrm{T}\mathrm{P},\mathrm{u}\mathrm{s}\mathrm{e}}=k_{\mathrm{u}\mathrm{s}\mathrm{e}}\cdot ATP,$$

(16)

where $k_{\mathrm{u}\mathrm{s}\mathrm{e}}=0.1$ denotes the effective rate constant for ATP consumption by other ATP-dependent cellular processes that are not modeled explicitly. These include, for example, ATP utilization in biosynthetic reactions, maintenance processes, cyclic nucleotide metabolism, and broader inhibitory effects associated with high ATP availability, including feedback on metabolic enzymes. This term therefore represents a generic ATP sink that prevents unbounded ATP accumulation and contributes to shaping the oscillatory ATP profile.

Together, Eqs. (11)–(16) define the ATP dynamics of the present single-cell model. Relative to the original DAM formulation, the essential novelty here is that ATP is no longer prescribed externally, but emerges dynamically from the interplay between Ca²⁺-dependent ATP consumption, PEP-cycle-associated ATP production, oxidative ATP production, and the metabolic redistribution of carbon through the GABA shunt. This explicit treatment of ATP dynamics is crucial, because ATP acts as the central link between mitochondrial metabolism, GABA-shunt activity, and Ca²⁺-dependent excitability in the present model.

2.1.3. DAM-Based GABA Dynamics

To describe GABA dynamics in a metabolically consistent way, we retained the DAM subsystem introduced in our previous study [18], including the same pool variables $P_0$, $P_1$, $P_2$, and $P_3$, as well as the same inter-pool fluxes $J_0$, $J_{01}(=J_{21})$, $J_{12}$, $J_{13}$, and $J_{32}$. In the original DAM formulation, these equations were used together with experimentally fitted time courses of $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ and ATP. In the present study, however, $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ and ATP are no longer prescribed externally, but are computed self-consistently from the differential equations introduced above. Thus, the DAM subsystem is retained in its original form, whereas its coupling to the rest of the model is extended by replacing the previously imposed $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ and ATP dynamics with explicitly calculated variables.

Within this framework, the variable $P_3$, representing the GABA pool in the original DAM model, directly provides the cytosolic GABA dynamics used in the present study. Accordingly, cytosolic GABA is not introduced here through an additional phenomenological equation, but follows directly from the original DAM equations. This preserves the metabolic interpretation of GABA as an integral part of the oscillatory redistribution of carbon between the glycolytic pool, the two TCA-cycle-related pools, and the GABA pool, as described previously [18].

As in the DAM model, the flux $J_{13}$ represents the cataplerotic branch associated with transfer from the right half of the TCA cycle toward the GABA pool, whereas $J_{32}$ represents the oxidative return flux from the GABA pool back toward the left half of the TCA cycle. These two fluxes are of particular importance in the present model because they provide the metabolic link between the DAM subsystem and ATP production. Specifically, $J_{13}$ modulates the PEP-cycle-related contribution to ATP production, whereas $J_{32}$ modulates the oxidative TCA-related contribution, as introduced in Eqs. (13) and (14). In this way, the previously developed DAM structure is preserved, while its metabolic output is now dynamically coupled to the explicitly modeled ATP and Ca²⁺ oscillations.

The remaining DAM equations and parameter values are the same as in the original publication [18] and are therefore not repeated here in full. This approach avoids unnecessary duplication while maintaining full consistency with the previously established metabolic framework. The essential extension introduced in the present study is that the DAM subsystem now operates within a closed dynamical model in which $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$, $C{a}_{\mathrm{E}\mathrm{R}}$, ATP, and GABA mutually interact. Consequently, GABA dynamics is no longer driven by externally imposed oscillatory inputs, but emerges from the coupled metabolic-calcium system itself.

To account for the fact that extracellular GABA does not follow cytosolic GABA instantaneously, the interstitial GABA concentration $GAB{A}_{\mathrm{i}\mathrm{s}}$ is modeled as a delayed function of cytosolic GABA,

$$GAB{A}_{\mathrm{i}\mathrm{s}}\left(t\right)=r_{\mathrm{G}}{GABA}_{\mathrm{c}\mathrm{y}\mathrm{t}}(t-{\tau }_{\mathrm{G}}),$$

(17)

where $r_{\mathrm{G}}=0.01$ denotes the proportionality factor relating cytosolic and interstitial GABA levels, and ${\tau }_{\mathrm{G}}=1$ represents a short effective delay associated with GABA release, diffusion in the interstitial space, and uptake/clearance. The introduction of the factor $r_G<1$ is physiologically motivated by the fact that beta cells maintain a high intracellular GABA pool, estimated to be in the millimolar range, whereas interstitial GABA in the islet is thought to lie in the nanomolar-to-low-micromolar range [14,15,16]. At present, however, a quantitative phase delay between cytosolic and interstitial GABA oscillations in pancreatic beta cells has not been experimentally established. Available evidence indicates that human beta cells release GABA directly from a cytosolic pool via VRAC in a pulsatile manner, with secretion periods typically in the 4–10 min range, while taurine transporter (TauT)-mediated uptake contributes to maintaining low interstitial GABA levels; together, these findings support a rapid release-clearance cycle but do not provide a direct numerical estimate of the lag between intracellular and extracellular GABA dynamics [14,15,16]. Recent work further showed that endogenous beta-cell-derived GABA is required for proper islet Ca²⁺ oscillation dynamics and suggested that GABA signaling may operate as a delayed negative feedback mechanism, yet the magnitude of this delay has likewise not been measured directly [17]. Accordingly, in the absence of direct experimental constraints, interstitial GABA is represented here by a small effective delay.

2.2. Two-Cell Coupled Model

To demonstrate how GABA promotes coordinated activity between pancreatic beta cells, we extended the single-cell model to a system of two coupled cells (Figure 2). Each cell retained the same intracellular dynamics as in the single-cell case, including Ca²⁺ handling, ATP production, and GABA oscillations, while coupling was introduced through a shared interstitial GABA signal. In this way, intracellular metabolic oscillations are converted into an intercellular coupling signal capable of entraining Ca²⁺ dynamics between neighboring beta cells.

The common interstitial GABA signal is defined as the average value of the interstitial GABA level associated with each individual cell. The arithmetic mean of the interstitial GABA level of the two-cell system is given by

$$GAB{A}_{\mathrm{i}\mathrm{s},\mathrm{a}\mathrm{v}\mathrm{g}}={\displaystyle \frac{1}{N}}{\sum }_{i=1}^{N}GAB{A}_{\mathrm{i}\mathrm{s},i}，N=2$$

(18)

where $GAB{A}_{\mathrm{i}\mathrm{s},i}=r_{\mathrm{G}}{GABA}_{\mathrm{c}\mathrm{y}\mathrm{t},i}(t - {\tau }_{\mathrm{G}})$ denotes the interstitial GABA signal associated with the $i$-th cell (see Eq. 17).

As illustrated schematically in Figure 2, the shared interstitial GABA signal, $GAB{A}_{\mathrm{i}\mathrm{s},\mathrm{a}\mathrm{v}\mathrm{g}}$, acts on both neighboring cells through ${\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}_{\mathrm{A}}$ and ${\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}_{\mathrm{B}}$ receptor-mediated pathways. Through these interactions, the common interstitial GABA pool modulates Ca²⁺ influx into each cell and thereby influences the timing and amplitude of intracellular Ca²⁺ oscillations. Because both cells are exposed to the same extracellular inhibitory GABA signal, this coupling provides a mechanism for phase alignment of their Ca²⁺ dynamics. In this way, GABA-mediated intercellular coupling can promote synchronization of Ca²⁺ signals in two-cell as well as multicellular systems, thereby supporting coordinated pulsatile insulin secretion at the islet level.

3. Results

The mathematical model defined by Eqs. (1)–(18) was used to investigate how GABA regulates Ca²⁺ oscillations and coordinated beta-cell dynamics. We begin by analyzing the consequences of altered GABA production, since this question is directly motivated by recent experimental observations showing that reduced endogenous GABA availability leads to decreased Ca²⁺ oscillation amplitude and prolonged oscillatory periods in pancreatic beta cells [17]. The first part of the Results section therefore examines how changes in GABA production affect the amplitude, frequency, and temporal organization of cytosolic Ca²⁺ oscillations in a single beta cell. We then dissect the underlying mechanisms by separating the metabolic contribution of GABA to ATP production from its extracellular signaling effect on Ca²⁺ influx, thereby identifying the dominant processes responsible for the observed oscillatory phenotype. Finally, we extend the analysis to coupled cells and show how delayed interstitial GABA signaling provides a physiologically plausible pathway for phase adjustment, entrainment, and synchronization of beta-cell oscillations. In this way, the Results section proceeds from experimentally motivated single-cell effects of altered GABA production to the emergence of coordinated multicellular dynamics.

3.1. GABA Controls Amplitude and Frequency of Ca²⁺ Oscillations in a Single Beta Cell

Within the proposed single-cell framework, we first examine how GABA production modulates intracellular Ca²⁺ oscillations. This question is directly motivated by the experimental findings of Ferreira et al. [17], who showed that endogenous beta-cell-derived GABA is essential for maintaining properly shaped Ca²⁺ oscillations, whereas loss of GABA signaling leads to a reduced oscillation amplitude and a prolonged oscillation period.

To investigate this effect, we varied the relative strength of the flux representing GABA production. Specifically, we introduced the dimensionless parameter $p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}}\in \left[\mathrm{0,1}\right]$, which rescales the DAM flux $J_{13}$ [18]. In the present model, $J_{13}$ is identified with the effective GABA production flux, $J_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}}$, such that

$$J_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}} :=p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}} J_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}}.$$

(19)

Here, $p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}}=0$ corresponds to complete suppression of GABA production, whereas $p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}}=1$ retains the original DAM flux unchanged. In this way, $p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}}$ provides a direct control parameter for intracellular GABA generation. Because cytosolic GABA in turn determines both the metabolic contribution of GABA to ATP production and the delayed interstitial GABA signal, varying $p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}}$ alters the coupled dynamics of GABA, ATP, and Ca²⁺ oscillations, as illustrated in Figure 3.

Figure 3A shows the model dynamics for the reference condition, $p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}}=1.0$, for which GABA production remains unchanged relative to the original DAM formulation. The upper panel displays oscillations of intracellular GABA, $GAB{A}_{\mathrm{c}\mathrm{y}\mathrm{t}}$, together with the corresponding interstitial GABA signal, $GAB{A}_{\mathrm{i}\mathrm{s}}$. As expected from Eq. (17), $GAB{A}_{\mathrm{i}\mathrm{s}}$ follows the intracellular oscillations with a small delay and reaches concentrations that are approximately two orders of magnitude lower than those of $GAB{A}_{\mathrm{c}\mathrm{y}\mathrm{t}}$. The lower panel shows the corresponding oscillations of cytosolic calcium, $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$, and local ATP concentration near $K_{\mathrm{A}\mathrm{T}\mathrm{P}}$ channels, illustrating the coupled metabolic and signaling dynamics under reference conditions.

Figure 3B shows the corresponding model dynamics under strongly reduced GABA production, $p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}}=0.1$. Under these conditions, both intracellular and interstitial GABA concentrations are substantially reduced relative to the reference case. This reduction is accompanied by pronounced changes in Ca²⁺ dynamics. In particular, the amplitudes of $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ oscillations are reduced and the oscillation frequency decreases, corresponding to a prolongation of the oscillation period. Thus, diminished GABA production leads to both weaker and slower cytosolic Ca²⁺ oscillations. These predictions are in good qualitative agreement with the experimental observations of Ferreira et al. [17], who reported that impaired endogenous GABA signaling is associated with reduced Ca²⁺ oscillation amplitude and prolonged oscillatory periods.

Figure 3C summarizes the dependence of the model output on $p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}}$. The upper panel shows that the average intracellular GABA concentration, $GAB{A}_{\mathrm{c}\mathrm{y}\mathrm{t},\mathrm{a}\mathrm{v}\mathrm{g}}$, increases monotonically with increasing GABA production. The lower panel shows the corresponding changes in Ca²⁺ oscillatory behavior. Both the minimum and maximum values of $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ increase with $p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}}$, and the oscillation frequency rises as well. The model therefore predicts that GABA production affects not only the amplitude of cytosolic Ca²⁺ oscillations, but also their temporal organization. Higher GABA production is associated with larger-amplitude and higher-frequency calcium oscillations, whereas reduced GABA production leads to smaller amplitudes and lower oscillation frequencies.

Because GABA influences the system through more than one pathway, including intracellular metabolic effects on ATP production and extracellular signaling effects on Ca²⁺ influx, we next analyze these mechanisms separately. This allows us to distinguish the relative contributions of the metabolic and signaling actions of GABA to the modulation of Ca²⁺ oscillatory dynamics.

3.1.1. Impact of GABA-Dependent ATP Production on Ca²⁺ Oscillations

In the analysis above, variation of the overall GABA production rate altered both intracellular GABA metabolism and extracellular GABA signaling. We next isolate the metabolic contribution of GABA and examine how GABA-dependent ATP production affects Ca²⁺ oscillations (Figure 4). For this purpose, GABA production itself is kept unchanged, whereas its contribution to ATP generation is selectively reduced through the parameters $k_{\mathrm{P}\mathrm{E}\mathrm{P},\mathrm{G}}$ and $k_{\mathrm{T}\mathrm{C}\mathrm{A},\mathrm{G}}$. As defined in Eqs. (13) and (14), these parameters determine how strongly GABA production enhances ATP generation through the PEP-related and TCA-related components, respectively. This approach allows us to assess how the intracellular metabolic action of GABA, independently of its receptor-mediated extracellular effects, contributes to the modulation of Ca²⁺ oscillatory dynamics.

Figure 4A shows the effect of the parameter $k_{\mathrm{P}\mathrm{E}\mathrm{P},\mathrm{G}}$, which determines the GABA-dependent contribution to ATP production through the PEP-related branch, on the properties of cytosolic Ca²⁺ oscillations. As $k_{\mathrm{P}\mathrm{E}\mathrm{P},\mathrm{G}}$ decreases, corresponding to a progressive reduction of the GABA effect on this component of ATP production, both the minimum and maximum values of $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ oscillations decrease, accompanied by a clear reduction in oscillation amplitude. The oscillation frequency also decreases, with the decline being most pronounced at lower values of $k_{\mathrm{P}\mathrm{E}\mathrm{P},\mathrm{G}}$.

Figure 4B shows the corresponding effect of the parameter $k_{\mathrm{T}\mathrm{C}\mathrm{A},\mathrm{G}}$, which determines the GABA-dependent contribution to ATP production through the TCA-related branch. In this case, decreasing $k_{\mathrm{T}\mathrm{C}\mathrm{A},\mathrm{G}}$ likewise lowers both the minimum and maximum values of $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ oscillations, whereas the oscillation amplitude remains comparatively stable over most of the parameter range. The oscillation frequency nevertheless decreases, with the effect becoming somewhat more pronounced at higher values of $k_{\mathrm{T}\mathrm{C}\mathrm{A},\mathrm{G}}$.

Taken together, these results indicate that GABA-dependent ATP production in both metabolic branches contributes to the maintenance of higher cytosolic Ca²⁺ levels and higher oscillation frequencies. However, the PEP-related branch has a more prominent effect on oscillation amplitude, whereas the TCA-related branch primarily affects oscillation timing while leaving the amplitude relatively preserved.

To further explain the changes in Ca²⁺ oscillation amplitude and frequency shown in Figure 4, we next examine how selective blockade of the GABA-dependent contribution to ATP production reshapes the coupled dynamics of local ATP near $K_{\mathrm{A}\mathrm{T}\mathrm{P}}$ channels and cytosolic Ca²⁺. These two variables are tightly interdependent, and their temporal interplay provides the mechanistic basis for the changes in Ca²⁺ oscillatory behavior observed in the model.

Figure 5 shows the dynamics of local ATP near $K_{\mathrm{A}\mathrm{T}\mathrm{P}}$ channels and cytosolic calcium, $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$, under reference conditions (Fig. 5A), after blocking the GABA-dependent contribution to ATP production through the PEP-related branch (Fig. 5B), and after blocking the corresponding contribution through the TCA-related branch (Fig. 5C). Here, $AT{P}_{\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{n}}$ denotes the threshold ATP level required to activate the Ca²⁺ influx pathway in the model (Eq. 10). When ATP rises above this threshold, Ca²⁺ influx is initiated and $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ begins to increase. The subsequent rise in cytosolic Ca²⁺ activates Ca²⁺-removal processes, which in turn reduce $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ and terminate the calcium pulse. Thus, the timing and magnitude of ATP excursions relative to $AT{P}_{\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{n}}$ determine both the onset and the strength of the corresponding Ca²⁺ response.

When the GABA-dependent contribution to ATP production through the PEP-related branch is blocked ($k_{\mathrm{P}\mathrm{E}\mathrm{P},\mathrm{G}}=0$), ATP generation during the cataplerotic phase becomes slower. This branch corresponds to the GABA-dependent support of the PEP cycle, which is most effective when $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ is low. As a result, the local ATP concentration near $K_{\mathrm{A}\mathrm{T}\mathrm{P}}$ channels rises more gradually, as reflected by the shallower ATP upstroke in Figure 5B. Consequently, ATP exceeds the threshold value $AT{P}_{\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{n}}$ only weakly, leading to a smaller Ca²⁺ influx and, therefore, to a lower peak of $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$. Because the Ca²⁺ rise remains reduced, the subsequent oxidative support of ATP production also remains limited. Together, these effects lead to both a lower amplitude and a lower frequency of Ca²⁺ oscillations.

When the GABA-dependent contribution to ATP production through the TCA-related branch is blocked ($k_{\mathrm{T}\mathrm{C}\mathrm{A},\mathrm{G}}=0$), ATP production during the oxidative phase is reduced because the contribution of the GABA shunt to oxidative ATP production, mediated by the return of carbon from the GABA pool into the TCA cycle, is abolished. Under these conditions, the local ATP concentration declines more strongly during the active phase of the oscillation. During the following cataplerotic phase, ATP production increases again; however, because ATP starts from a lower level, more time is required to reach the threshold value $AT{P}_{\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{n}}$. This primarily prolongs the oscillation period and thus lowers the oscillation frequency. In contrast, the amplitude of the Ca²⁺ oscillations remains relatively preserved, because the late ATP rise still produces a sufficient overshoot above $AT{P}_{\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{n}}$ to trigger a pronounced Ca²⁺ influx, as shown in Figure 5C.

Taken together, the results in Figure 5 show that the GABA-dependent ATP-producing contribution through the PEP cycle is particularly important for generating a rapid and sufficiently large ATP rise during the cataplerotic phase, which in turn determines both the amplitude and the frequency of the subsequent Ca²⁺ oscillations. By contrast, the GABA-dependent contribution associated with the GABA shunt and oxidative TCA-cycle metabolism primarily supports ATP during the oxidative phase and therefore affects oscillation timing more strongly than oscillation amplitude.

3.1.2. Impact of GABA-Dependent Ca²⁺ Influx on Ca²⁺ Oscillations

We next examine how interstitial GABA modulates Ca²⁺ entry into the cell and, consequently, the dynamics of cytosolic Ca²⁺ oscillations. In the model, this regulation acts through the influx term $J_{\mathrm{i}\mathrm{n}}$, which describes Ca²⁺ entry from the extracellular space into the cytosol. As summarized in Eqs. (7)–(10), $GAB{A}_{\mathrm{i}\mathrm{s}}$ affects $J_{\mathrm{i}\mathrm{n}}$ through two inhibitory mechanisms. First, it directly reduces the effective maximal Ca²⁺ influx rate by decreasing $g_{\mathrm{i}\mathrm{n}}$. Second, it gradually lowers the effective channel activation variable $x$, thereby introducing a slower inhibitory feedback on Ca²⁺ entry. Together, these two actions represent the combined fast and delayed effects of interstitial GABA on Ca²⁺ influx and provide the basis for analyzing how extracellular GABA signaling reshapes Ca²⁺ oscillatory dynamics.

To quantify this effect, we introduce the dimensionless parameter $p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{i}\mathrm{s}}$, which rescales the interstitial GABA signal according to

$$GAB{A}_{\mathrm{i}\mathrm{s}} :=p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{i}\mathrm{s}}\cdot GAB{A}_{\mathrm{i}\mathrm{s}}.$$

(20)

The parameter $p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{i}\mathrm{s}}$ ranges from 0 to 1 and therefore controls the strength of the $GAB{A}_{\mathrm{i}\mathrm{s}}$-dependent inhibition in Eqs. (8) and (10). A value of $p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{i}\mathrm{s}}=0$ corresponds to a complete blockade of the interstitial GABA effect on Ca²⁺ influx, whereas $p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{i}\mathrm{s}}=1$ preserves the full inhibitory contribution of $GAB{A}_{\mathrm{i}\mathrm{s}}$ used in the reference model.

The resulting changes in Ca²⁺ dynamics are shown in Figure 6. Figure 6A compares representative time courses of $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ for the reference case ($p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{i}\mathrm{s}}=1$) and for complete blockade of the interstitial GABA effect ($p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{i}\mathrm{s}}=0$). Blocking the effect of interstitial GABA on $J_{\mathrm{i}\mathrm{n}}$ leads to a slight increase in the amplitude of $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ oscillations, consistent with the inhibitory action of $GAB{A}_{\mathrm{i}\mathrm{s}}$ on Ca²⁺ entry. In contrast, the oscillation frequency changes only modestly. This behavior is also evident in Figure 6A, where the two calcium traces remain very similar in shape but gradually develop a phase shift over time. The shaded region highlights the resulting slow beating interval, which arises from the small difference in oscillation frequency between the two conditions.

This behavior is quantified in Figure 6B, which shows the dependence of $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t},\mathrm{m}\mathrm{i}\mathrm{n}}$, $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}$, and oscillation frequency on $p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{i}\mathrm{s}}$. As the strength of interstitial GABA signaling is reduced, both the minimum and maximum values of $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ increase slightly, indicating a modest increase in oscillation amplitude. By contrast, the oscillation frequency changes only weakly over the entire parameter range. Thus, Figure 6B confirms that the primary effect of interstitial GABA on the single-cell oscillator is not a strong shift in intrinsic oscillation frequency, but rather a comparatively small modulation of Ca²⁺ influx that manifests mainly through subtle changes in oscillation amplitude and timing.

Although the isolated effect of interstitial GABA on Ca²⁺ influx produces only a modest change in oscillation frequency, it still plays an important mechanistic role by adjusting oscillatory phase. In the present model, the intrinsic rhythm of the single-cell oscillator is determined primarily by ATP dynamics, whereas GABA-dependent modulation of $J_{\mathrm{i}\mathrm{n}}$ acts mainly as a phase-adjusting mechanism. Accordingly, when GABA acts only through Ca²⁺ influx, the resulting $C{a}_{\mathrm{c}\mathrm{y}\mathrm{t}}$ traces remain similar in overall amplitude and period but progressively drift in phase, as illustrated by the beating pattern in Figure 6A. Thus, the GABA-dependent regulation of Ca²⁺ influx does not constitute the dominant mechanism setting the intrinsic oscillation frequency; rather, it fine-tunes oscillatory timing by shifting phase. This distinction becomes especially important in the coupled-cell setting, where even a moderate GABA-dependent phase adjustment of Ca²⁺ entry can facilitate entrainment between neighboring cells and thereby promote efficient synchronization, as examined in the next section.

3.2. Intercellular GABA Couples and Entrains Ca²⁺ Oscillations in Two Beta Cells

We next examined whether the GABA-dependent modulation of Ca²⁺ influx, which in the single-cell analysis acted primarily as a phase-adjusting mechanism, is sufficient to promote entrainment and synchronization between two coupled beta cells. For this purpose, we used the two-cell model introduced in Section 2.2, in which each cell retains its own intracellular metabolic and Ca²⁺ dynamics, whereas coupling is mediated through the shared interstitial GABA signal $GAB{A}_{\mathrm{i}\mathrm{s},\mathrm{a}\mathrm{v}\mathrm{g}}$. In this way, intracellular GABA oscillations generated by the metabolic subsystem are transformed into a common extracellular coupling signal that acts back on both cells and modulates Ca²⁺ influx in each of them.

The corresponding dynamics are shown in Figure 7. At the beginning of the simulation (Fig. 7A), the two cells exhibit a clear phase difference in their cytosolic Ca²⁺ oscillations, indicating that they are initially not entrained. Consistent with this, the two interstitial GABA signals are also phase-shifted during the initial regime. After the onset of GABA-mediated coupling, both cells begin to sense the shared interstitial GABA signal $GAB{A}_{\mathrm{i}\mathrm{s},\mathrm{a}\mathrm{v}\mathrm{g}}$, which progressively reduces the phase difference between the two oscillators by jointly modulating Ca²⁺ influx. To make the late stage of this process more visible, Figure 7 contains a break in the time axis. In the later regime shown in Fig. 7B, the two interstitial GABA signals overlap and form a common oscillatory signal, while the corresponding cytosolic Ca²⁺ traces also overlap almost completely, indicating effective entrainment and near-complete synchronization.

As an additional quantitative confirmation of this transition, synchronization was characterized by a time-dependent Pearson correlation coefficient, $corr$, calculated over an exponentially weighted moving window, following the general approach of Pozzi et al. [28]. In the present context, this coefficient serves as a compact measure of the transition from an initially phase-shifted regime to a synchronized state and does not require further methodological detail here. As shown in Fig. 7C, the correlation is initially low or negative, consistent with the phase-shifted oscillations in Fig. 7A, and then rises progressively toward values close to 1 after GABA-mediated coupling is switched on, indicating convergence of the two Ca²⁺ oscillators to a common rhythm.

These results show that even a relatively modest GABA-dependent modulation of Ca²⁺ influx is sufficient to coordinate oscillatory beta-cell activity when transmitted through the interstitial space. In mechanistic terms, ATP-dependent metabolic dynamics primarily determine the intrinsic oscillatory rhythm of each cell, whereas intercellular GABA provides the coupling signal that adjusts phase and gradually aligns the two oscillators. Thus, the model supports a division of roles in which metabolism sets the intrinsic pace of the oscillations, while extracellular GABA-mediated signaling enables entrainment and synchronization between neighboring cells. This interpretation is consistent with the experimental view that beta-cell-derived GABA can function as an intercellular coordinating signal within the islet [14].

To further evaluate the physiological relevance of this coupling mechanism, we next examined how strongly the parameters of the two cells may differ while synchronization is still maintained. This question is important because beta cells within the islet are intrinsically heterogeneous and therefore do not share identical metabolic, signaling, or excitability properties. In the model, such physiological variability is represented by differences in selected parameter values between the two cells. We therefore asked whether GABA-mediated coupling remains sufficiently robust to entrain cells that differ in their intrinsic oscillatory characteristics. Rather than analyzing all model parameters, we selected a representative set of nine parameters according to their mechanistic roles in the system. Specifically, the selected set includes parameters describing the metabolic component of the model, namely GABA production, the PEP-related and TCA-related contributions to ATP generation, and the two ATP-consumption terms; parameters representing extracellular GABA action on Ca²⁺ influx, including the overall strength of the interstitial GABA signal and the fast and slow inhibitory coupling pathways; and a parameter characterizing intrinsic cellular excitability through the ATP threshold for activation of Ca²⁺ influx. In this way, the selected set spans the principal determinants of the intrinsic oscillatory rhythm, the strength of intercellular GABA-mediated phase adjustment, and the baseline responsiveness of the Ca²⁺-influx pathway. Table 1 summarizes the corresponding synchronization ranges and thereby illustrates not only the robustness of the model itself, but also the capacity of GABA-mediated entrainment to synchronize metabolically and functionally non-identical beta cells.

In non-identical cells, GABA-mediated coupling does not necessarily induce complete in-phase synchronization of the two Ca²⁺ signals ($corr=1$), because the cells represent non-identical oscillators with slightly different intrinsic frequencies and intrinsic dynamics. Thus, coupling does not fully equalize their signals, but rather stabilizes their relative phase, allowing the cells to oscillate with a common rhythm while maintaining a finite phase difference ($corr<1$). Table 1 shows the ranges of selected parameters in the second cell relative to the first cell for which GABA-mediated coupling still supports stable phase locking with a moving Pearson correlation coefficient of $corr>0.7$, corresponding to a Ca²⁺ signal phase lag of up to approximately 10% of the oscillation period.

Column A of Table 1 shows the synchronization ranges obtained for the reference strength of the interstitial GABA signal ($p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{i}\mathrm{s}}=1.0$). To assess how stronger interstitial GABA action influences synchronization between the two cells, Column B shows the corresponding ranges for the increased value ($p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{i}\mathrm{s}}=1.2$). Comparison of Columns A and B shows that increasing the influence of interstitial GABA generally broadens the range of parameter mismatch over which stable phase locking is still achieved ($corr>0.7$). This indicates that a stronger shared GABA signal more effectively compensates for differences in the intrinsic properties of the two cells and thereby enhances the robustness of synchronization.

A notable exception is the parameter $k_{\mathrm{P}\mathrm{E}\mathrm{P},\mathrm{G}}$, for which increasing $p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{i}\mathrm{s}}$ reduces the upper limit of the synchronization range. A positive deviation of $k_{\mathrm{P}\mathrm{E}\mathrm{P},\mathrm{G}}$ in the second cell enhances the PEP-related contribution to ATP production, causing local ATP to rise more rapidly and to reach the threshold for Ca²⁺ influx earlier. As a result, the corresponding Ca²⁺ pulse is advanced in time relative to that of the first cell. This interpretation is consistent with Figure 5A and 5B, which show that the PEP-related branch strongly affects the rate of ATP rise and the timing of threshold crossing. Moreover, Figure 4A indicates that increasing $k_{\mathrm{P}\mathrm{E}\mathrm{P},\mathrm{G}}$ above its reference value ($k_{\mathrm{P}\mathrm{E}\mathrm{P},\mathrm{G}}=1$) changes oscillation frequency less strongly than decreasing it below the reference value, which explains the relatively broad positive synchronization range under reference coupling. When interstitial GABA signaling is strengthened ($p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{i}\mathrm{s}}=1.2$), however, the delayed GABA-dependent inhibition of Ca²⁺ influx is also enhanced. Under these conditions, larger positive deviations of $k_{\mathrm{P}\mathrm{E}\mathrm{P},\mathrm{G}}$ lead to a less favorable balance between earlier ATP-dependent triggering and stronger delayed inhibitory feedback, thereby narrowing the upper bound of the synchronization range. In other words, stronger interstitial GABA generally improves phase locking, but in the case of excessive enhancement of the PEP-related ATP-producing branch it can also accentuate the timing mismatch that limits stable entrainment.

To further extend the analysis, we asked whether synchronization could be reinforced not only by paracrine GABA-mediated coupling, but also by a weak effective electrical interaction between the cells. This question is physiologically relevant because beta cells within the islet are connected not only through extracellular signaling, but also through direct electrical coupling via gap junctions, predominantly formed by connexin-36 (Cx36) channels. Experimental studies have shown that Cx36-mediated gap-junction coupling is important for the synchronization of glucose-induced Ca²⁺ oscillations and pulsatile insulin secretion in beta cells and intact islets [29,30,31]. In addition, experimental and computational studies have demonstrated that gap-junction coupling, together with beta-cell heterogeneity, can shape coordinated Ca²⁺ wave propagation and population-level synchronization within the islet [32]. Because the present model does not include an explicit equation for membrane potential, we represented this effect phenomenologically by adding an additional term to Eq. (8), which defines the conductance of Ca²⁺ influx,

$$g_{\mathrm{i}\mathrm{n},\mathrm{i}}=g_{\mathrm{i}\mathrm{n},0}-k_{\mathrm{G}}GAB{A}_{\mathrm{i}\mathrm{s},\mathrm{a}\mathrm{v}\mathrm{g}}+k_{\mathrm{e}\mathrm{l}}(C{a}_{\mathrm{c}\mathrm{y}\mathrm{t},j}-C{a}_{\mathrm{c}\mathrm{y}\mathrm{t},i}),$$

(21)

where $k_{\mathrm{e}\mathrm{l}}$ determines the strength of weak effective electrical coupling. This term does not represent direct diffusion of Ca²⁺ between cells; rather, it provides an effective description of how electrical coupling can influence cellular excitability and, consequently, voltage-dependent Ca²⁺ entry. The results of this case are summarized in Column C of Table 1 ($p_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{i}\mathrm{s}}=1.2$, $k_{\mathrm{e}\mathrm{l}}=1$). Comparison with Column B shows that the addition of weak electrical coupling further broadens the ranges of all selected parameters for which stable phase locking is maintained ($corr>0.7$). This indicates that weak electrical coupling can complement GABA-mediated synchronization and further increase the robustness of coordinated activity between non-identical cells.

Viewed together, the results summarized in Table 1 show that GABA-mediated intercellular coupling is sufficiently robust to entrain beta cells even in the presence of physiologically plausible heterogeneity, and that this robustness can be further enhanced by weak effective electrical coupling. The model therefore supports the view that extracellular GABA can act as an effective coordinating signal that aligns the phases of non-identical beta-cell oscillators and promotes synchronized Ca²⁺ activity at the islet level, while electrical coupling provides an additional stabilizing influence on this process. This interpretation is consistent with experimental evidence showing that pulsatile GABA release from beta cells can contribute to the coordination of islet activity [14], with the broader view that extracellular GABA participates in receptor-mediated auto/paracrine regulation of islet function [15], and with recent findings that endogenous beta-cell-derived GABA is required for proper Ca²⁺ oscillation dynamics and insulin secretion [17].

4. Discussion

The present study provides a mechanistic framework linking intracellular GABA metabolism to both single-cell Ca²⁺ oscillations and intercellular synchronization in pancreatic beta cells. By combining the previously developed DAM framework [18] with explicit equations for cytosolic Ca²⁺, ER Ca²⁺, ATP, and GABA-dependent regulation of Ca²⁺ influx, we were able to connect metabolic, signaling, and collective aspects of beta-cell function within a single dynamical model. In this way, the model extends earlier conceptual work on the metabolic role of the GABA shunt in beta cells [8,12,13] toward an explicit explanation of how GABA shapes Ca²⁺ oscillatory dynamics and coordinated beta-cell behavior.

A central result of the study is that reduced GABA production leads to lower-amplitude and lower-frequency Ca²⁺ oscillations. This prediction is in good qualitative agreement with the recent experiments of Ferreira et al. [17], who showed that endogenous beta-cell-derived GABA is required for properly shaped Ca²⁺ oscillations and normal insulin secretion. In the present model, this phenotype is not imposed phenomenologically, but emerges from the coupled effects of reduced intracellular GABA metabolism and weaker extracellular GABA signaling. Thus, the model provides a mechanistic explanation for how impaired endogenous GABA availability can alter both the magnitude and timing of beta-cell Ca²⁺ signals.

The analysis further shows that the metabolic action of GABA is best understood within the dual-anaplerotic logic of the DAM. In this framework, GABA contributes to ATP production in two related but functionally distinct ways. First, through the GABA shunt, it supports oxidative ATP production during the oxidative phase by returning carbon from the GABA pool into the TCA cycle. Second, by increasing the effective carbon volume available for subsequent cataplerosis, it indirectly promotes the PEP cycle and thereby strengthens local ATP production during the cataplerotic phase. The model indicates that this PEP-related contribution is especially important for generating a rapid and sufficiently large ATP rise near K_ATP channels, which in turn is crucial for both the amplitude and the frequency of the ensuing Ca²⁺ oscillations. By contrast, the oxidative TCA-related contribution primarily stabilizes ATP during the active phase and therefore affects oscillation timing more strongly than oscillation amplitude. This interpretation is consistent with the broader view that GABA metabolism is functionally integrated with beta-cell bioenergetics and insulin secretion [8,13,23,24].

A second important result is that extracellular GABA signaling, when considered in isolation, exerts only a modest effect on the intrinsic oscillation frequency of a single cell, but plays a key role in phase adjustment. In the model, this extracellular action is represented by a fast inhibitory effect on the effective Ca²⁺-influx term and a slower inhibitory effect mediated through the regulatory variable $x$. Physiologically, these two components correspond to a simplified description of ionotropic and metabotropic GABAergic actions, respectively, and are consistent with the current view that beta-cell-derived GABA acts through both auto/paracrine receptor-mediated pathways and metabolic coupling [15,16,17]. The present results therefore suggest a functional division of labor between the two major GABA actions in the model: ATP dynamics primarily determine the intrinsic oscillatory rhythm, whereas extracellular GABA-dependent modulation of Ca²⁺ influx primarily adjusts oscillatory phase.

This distinction becomes particularly important at the multicellular level. Once two cells are coupled through the shared interstitial GABA signal, even a relatively small GABA-dependent phase adjustment of Ca²⁺ influx is sufficient to entrain the two oscillators and drive them toward synchronized activity. Importantly, this synchronization persists over finite ranges of parameter mismatch between the two cells, which is physiologically relevant because beta cells within the islet are intrinsically heterogeneous in their metabolic activity, signaling properties, and excitability. The model therefore supports the view that extracellular GABA can act as an effective coordinating signal that aligns the phases of non-identical beta-cell oscillators. In addition, our analysis shows that weak effective electrical coupling can further broaden the parameter ranges over which stable phase locking is maintained, indicating that paracrine GABA signaling and electrical interactions may act in a complementary manner to stabilize coordinated beta-cell activity. This interpretation is consistent with experimental evidence that pulsatile GABA release from beta cells can contribute to the coordination of islet activity [14], with the broader view that extracellular GABA participates in receptor-mediated auto/paracrine regulation of islet function [15], with the observation that endogenous GABA is required for proper Ca²⁺ oscillation dynamics [17], and with studies showing that Cx36-mediated coupling supports synchronized beta-cell Ca²⁺ activity and pulsatile insulin secretion [29,30,31,32].

From a modeling perspective, the present framework does not aim to replace detailed electrophysiological models of beta-cell activity, but rather to complement them. Previous mathematical studies have provided important insight into bursting dynamics, the interaction between metabolic and electrical oscillations, oscillations in K_ATP conductance, phantom bursting, and the behavior of coupled and multicellular islet models [22,33,34,35,36,37]. Our model adds to this line of work by explicitly incorporating intracellular GABA dynamics, GABA-dependent ATP production, delayed interstitial GABA signaling, and a phenomenological weak electrical-coupling term into a minimal Ca²⁺–ATP oscillatory framework. In this sense, it can be viewed as a metabolically and paracrinely enriched module that captures how paracrine GABA signaling and electrical interactions may jointly contribute to beta-cell coordination, and that could in future work be integrated with more detailed membrane-potential-based models of beta-cell electrophysiology and islet network behavior.

The model also has several limitations. Membrane potential and individual ionic currents are not described explicitly, and the effects of GABA_A and GABA_B receptors are represented only phenomenologically through effective modulation of Ca²⁺ influx and the slow regulatory variable $x$. Likewise, the interstitial GABA signal is represented by a simple delayed relation to cytosolic GABA, and the present multicellular analysis is restricted to a pair of coupled cells. In addition, electrical coupling is not modeled through explicit gap-junction currents or membrane-potential equations, but only through a weak effective term acting on the Ca²⁺-influx conductance. These simplifications mean that the model is primarily suited to the slow oscillatory timescale considered here and does not yet address the faster and ultrafast electrical bursting dynamics that also contribute to beta-cell activity. Future extensions should therefore combine the present GABA-centered metabolic framework with more detailed electrophysiological descriptions, explicit gap-junction coupling, and larger multicellular islet architectures, thereby linking GABA metabolism, membrane excitability, electrical connectivity, and spatially distributed islet synchronization in a unified model [19,35,37].

In conclusion, the present study identifies GABA as a dual regulator of beta-cell dynamics: metabolically, it shapes the intrinsic amplitude and frequency of Ca²⁺ oscillations through its effects on ATP production, and paracrinely, it promotes entrainment and synchronization by adjusting the phase of Ca²⁺ influx between neighboring cells. The model further suggests that this GABA-mediated coordinating effect can be reinforced by weak electrical coupling, thereby increasing the robustness of synchronized activity in heterogeneous beta-cell populations. Together, these findings provide a coherent mechanistic framework that connects GABA metabolism, local ATP dynamics, Ca²⁺ oscillations, intercellular phase coordination, and coordinated collective activity in pancreatic islets.