Phase modulation of insulin pulses enhances glucose regulation and enables inter-islet synchronization

Insulin is secreted in a pulsatile manner from multiple micro-organs called the islets of Langerhans. The amplitude and phase (shape) of insulin secretion are modulated by numerous factors including glucose. The role of phase modulation in glucose homeostasis is not well understood compared to the obvious contribution of amplitude modulation. In the present study, we measured Ca2+ oscillations in islets as a proxy for insulin pulses, and we observed their frequency and shape changes under constant/alternating glucose stimuli. Here we asked how the phase modulation of insulin pulses contributes to glucose regulation. To directly answer this question, we developed a phenomenological oscillator model that drastically simplifies insulin secretion, but precisely incorporates the observed phase modulation of insulin pulses in response to glucose stimuli. Then, we mathematically modeled how insulin pulses regulate the glucose concentration in the body. The model of insulin oscillation and glucose regulation describes the glucose-insulin feedback loop. The data-based model demonstrates that the existence of phase modulation narrows the range within which the glucose concentration is maintained through the suppression/enhancement of insulin secretion in conjunction with the amplitude modulation of this secretion. The phase modulation is the response of islets to glucose perturbations. When multiple islets are exposed to the same glucose stimuli, they can be entrained to generate synchronous insulin pulses. Thus, we conclude that the phase modulation of insulin pulses is essential for glucose regulation and inter-islet synchronization.


Introduction
Hormones are long-distance messengers that regulate numerous body functions, including metabolism. Many hormones are secreted in an oscillatory manner [1]. Insulin, one of the most important hormones in human physiology, is also secreted in a pulsatile manner [2]. In a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 amount of Ca 2+ oscillation data (60-100 samples for each condition), which were sufficient to probe the glucose-dependent modulations. Here we excluded two samples (from total 66 samples) for 15 mM that showed prolonged active phases without regular oscillations. A recent study has reported that 15 mM is a bistable glucose concentration that can generate both oscillatory and stationary Ca 2+ dynamics [30]. First, we examined the duration of the active and silent phases of the Ca 2+ oscillations (See Materials and methods). We defined the ratio of active to silent phases as the "plateau fraction" [13]. The islets generated higher plateau fractions at higher glucose concentrations ( Fig 1B). Second, we quantified the period of Ca 2+ oscillations using a Fourier transformation of their time traces. Dominant periods showed a Ushaped glucose dependence (Fig 1C). Compared to 8 mM (4.1±1.1 min) and 15 mM glucose (3.8±0.9 min), 10 mM and 12 mM glucose showed shorter periods (3.3±0.6 and 2.9±0.7 min, respectively; P<0.01). We also noticed longer dominant periods accompanied by larger variations. The dominant periods were not dependent on islet size (Fig 1D).
We observed that these nontrivial features of the plateau factions and dominant periods could be simply described by the "active rotator model," which has been widely used to understand oscillatory phenomena in physics [31][32][33]. Suppose that the spontaneous Ca 2+ oscillations can be approximated as a periodic function, where r and θ are the amplitude and phase at time t, respectively, and unity is added to prevent a negative Ca 2+ concentration. For the moment, we ignore amplitude by setting r = 1 to focus on the phase of the Ca 2+ oscillations. The phases θ = 0 and π represent the peak and valley of the Ca 2+ oscillations, respectively. If the phase linearly increases with a constant angular velocity (dθ / dt = ω), the Ca 2+ oscillations become a pure sinusoidal wave with equal durations for the active and silent phases. However, the islet Ca 2+ concentrations oscillate in phase-modulated waves that are dependent on the glucose concentration. The modulation can be captured by the active rotator model, where the strength μ of the pinning force controls the degree of phase modulation by breaking the rotational symmetry and forcing oscillators to remain in the active phase (θ = 0) for a longer time and in the silent phase (θ = π) for a shorter time for positive μ and vice versa for negative μ. In the presence of the pinning force, the effective period of the active rotator is 2p= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi o 2 À m 2 p . Thus, the pinning force always makes the rotator slower regardless of its sign. These features of the active rotator model are consistent with the glucose-dependent plateau fraction and dominant period of the islet Ca 2+ oscillations (Fig 2A). Thus, we hypothesize that islets have an intrinsic period in their Ca 2+ oscillations, but the period is modified by phase modulations that are dependent on the glucose concentration. Specifically, the phase modulation was absent at 12 mM glucose because in this condition, the durations of the active and silent phases were equal ( Fig 1B). Furthermore, 8 mM glucose could impose a negative μ that increased the silent phase and the period of Ca 2+ oscillations, whereas 15 mM glucose could impose a positive μ that increased the active phase and the period of Ca 2+ oscillations. We quantified the μ values to fit the Ca 2+ oscillations under different glucose concentrations. One practical issue with the estimation is that Ca 2+ oscillations have heterogeneous periods between islets, even for the same glucose concentration ( Fig 2B). Therefore, we considered a distribution of intrinsic periods. For simplicity, we assumed that the periods follow a Gaussian distribution with a mean ( " o) and standard deviation (δω) Thus, different glucose concentrations have different μ values that modulate the Gaussian distribution. With the maximum likelihood estimation of the parameters, our hypothesis could fit the period distributions of the measured Ca 2+ oscillations ( Fig 2B). The estimated μ values for glucose concentration (G) could be fit with a nonlinear function, with parameters " m ¼ 0:016 min -1 , G m 0 ¼ 11:6 mM and δG μ = 1.2 mM (Fig 2C). The phase modulation was gradual at~12 mM glucose. In conclusion, the phase modulation of the simple rotator model was sufficient to explain the glucose-dependent changes in the plateau fraction and the U-shaped dominant period of the islet Ca 2+ oscillations.

Optimal glucose oscillations can entrain islet Ca 2+ oscillations
The phase modulation of islet Ca 2+ oscillations by μ(G) is the response to glucose stimuli. In the present study, we simulated Ca(t), which represents the response of the active rotators to glucose concentrations that alternate between G 1 and G 2 with a period of T ( Fig 3A). For the simulation, we used active rotators for which the intrinsic period and pinning force followed the previously determined distributions, ω and μ(G). Because the rotators experienced common glucose stimuli, the optimal driving signals of the glucose alternation could entrain the rotators to generate synchronous oscillations. For G 1 = 8 mM, we examined the degree of synchronization between the rotators for various amplitudes, G 2 -G 1 , and periods, T ( Fig 3B). The glucose alternation could more effectively provide synchronization when the amplitudes were larger and had an optimal period (3-4 min). We then repeated the same simulation with G 1 = 10 mM, and μ(G) was less sensitive to the glucose changes than in the simulation using G 1 = 8 mM ( Fig 3C). The optimal driving periods did not change, but larger amplitudes were required to induce synchronization.
To confirm these predictions, we measured islet Ca 2+ oscillations stimulated by alternating glucose concentrations ( Fig 3A). Glucose concentrations that alternated between G 1 = 8 mM and G 2 = 10 mM with a period of T = 4 min could entrain islets to generate synchronous Ca 2+ oscillations (Table 1). However, glucose alternation with a smaller amplitude (0.5 mM) or with shorter (2.7 min) or longer (5.3 min) periods was not able to entrain islets significantly. We also repeated the same experiment with G 1 = 10 mM and confirmed that G 1 = 10 mM was less effective than G 1 = 8 mM for entraining islets, as predicted in the simulations.

Phase modulation of insulin pulses has a functional role in glucose homeostasis
Thus far, we have examined how islet Ca 2+ oscillations respond to constant and alternating glucose stimuli. Our ultimate question is whether the phase modulation of the islet Ca 2+ oscillations has a functional role in glucose regulation. To answer this question, one needs to consider insulin secretion by glucose and glucose regulation by insulin. Given that the intracellular Ca 2+ concentration is highly correlated with insulin secretion [25,34], we assumed that insulin pulses have the same phase, θ(t), as Ca 2+ oscillations. Therefore, insulin secretion can also be described by a periodic function, where r(G) and θ(t) are the glucose-dependent amplitude and time-varying phase, respectively, of the insulin pulses. Note that r(G) is the amplitude of the insulin pulses, not the Ca 2+ oscillations. Thus, the amplitude and the modulated phase of the insulin pulses contribute to insulin secretion. In particular, we defined the latter contribution as a shape factor, s(G), which is the average insulin secretion with a fixed amplitude, r = 1 ( Fig 4A). Thus, the average insulin secretion becomes hI(G)i = r(G) Á s(G). Given the information of hI(G)i and s(G), we could infer r (G) (Fig 4B). Notably, hI(G)i is most sensitive at G = 15 mM, but r(G) is most sensitive at G = 8 mM, which is closer to the normal glucose range.
To complete the glucose-insulin regulation loop (Fig 5A), one also needs to understand glucose regulation by insulin and insulin secretion by glucose. Thus, we model the glucose regulation by insulin, where G in is the external glucose infusion rate, and v is the glucose clearance rate. The second term in Eq (5) represents the glucose clearance, which is proportional to the insulin level (I) and the present glucose concentration (G). For the simulation with multiple islets, I should be the average insulin secretion of the islets (See Materials and methods). Now, we have a complete model of glucose-insulin regulation to examine the functional role of the phase modulation of insulin pulses. We compared two models: Model I considers a phase modulation (μ 6 ¼ 0), whereas Model II ignores it (μ = 0). Thus, Model II always generates pure sinusoidal waves of insulin with a constant shape factor, s(G) = 1. Given a large glucose infusion (G in = 8 mM/min), Model I maintained a lower glucose concentration and showed synchronous insulin pulses and an oscillatory glucose concentration, unlike Model II (Fig 5B). Stationary glucose concentrations for G in could also be mathematically estimated. The stationary condition(dG / dt = 0) implies that the glucose infusion and clearance are balanced: In addition to the balance equation, we have another equation, glucose-dependent insulin secretion: hI(G)i = r(G)Ás(G) for Model I and hI(G)i = r(G) for model II due to s(G) = 1. These two equations uniquely determine the stationary glucose concentrations (Fig 5C). The analysis demonstrated that Model I can maintain the stationary glucose concentration within a narrower range than Model II, given various glucose infusion rates. This result occurs because the shape factor s(G) effectively modifies insulin secretion by suppressing/enhancing insulin secretion for small/large glucose infusions. We also checked this conclusion with simulations ( Fig  5D). In addition to the enhanced glucose regulation, the phase modulation led islets to generate synchronous insulin pulses (Fig 5E). However, when the glucose concentration was raised to a higher level (G = 10−12mM) with a large glucose infusion, the synchronization of the insulin pulses decreased because the phase modulation became insensitive to the glucose change occurring at the higher glucose concentration (Fig 2C). This result may explain the diminished insulin pulsatility in diabetic patients with prolonged hyperglycemia.

Discussion
In this study, we showed that the shape (phase) of islet Ca 2+ oscillations was modulated depending on glucose concentration. Using the Ca 2+ oscillation data, we developed a mathematical model that incorporated (i) insulin secretion by glucose and (ii) glucose regulation by insulin. The model demonstrated that the glucose-dependent shape modulation of insulin pulses contributed to tighter regulation of the normal glucose level, compared to their amplitude modulation alone. The additional shape modulation could suppress unnecessary insulin secretion at low glucose, while it could enhance demanding insulin secretion at high glucose. Furthermore, the shape modulation provided a common cue for physically-separate islets to generate synchronous insulin pulses. Otherwise, incoherent insulin pulses from numerous islets would lead to diminished pulsatility of circulating insulin levels. Therefore, we concluded that the shape modulation of insulin pulses was crucial for glucose regulation and inter-islet synchronization.
The diminished pulsatility of circulating insulin levels has been observed in diabetic conditions (refer to a review of [11]). We have found that the shape modulation of insulin pulses was relatively insensitive at high glucose levels (10-14 mM) (Fig 2C). Thus when the external glucose stimulus was alternating between G 1 and G 2 , starting from G 1 = 10 mM was less effective to synchronize the islet Ca 2+ oscillations than starting from G 1 = 8 mM (Fig 3B  and 3C). In the glucose-insulin loop simulation, when glucose levels stayed high with large glucose infusion rates (G in > 2 mM/min), the inter-islet synchronization was poor (Fig 5E). The asynchronous insulin pulses at hyperglycemic conditions could result in the diminished pulsatility of circulating insulin levels. The incoherent insulin pulses at prolonged hyperglycemia may suggest another reasoning for the loss of insulin pulsatility under diabetes, in addition to the deteriorated pulse generation of individual islets. This mechanism, however, cannot explain the loss of insulin pulsatility for non-diabetic subjects without hyperglycemia [12].
We modeled the phase dynamics of the Ca 2+ oscillations by adopting the active rotator model instead of describing the detailed molecular processes inside of the cells. Rather than using the constant pinning force (μ) of the original active rotator model, we coupled the pinning force μ(G) to an external stimulus, glucose (G). Thus, the simple rotator model could successfully capture the dependence of the plateau fractions and dominant periods of the Ca 2+ oscillations on the glucose concentration. The modified rotator model is closely related to the theta model (also known as the Ermentrout-Kopell canonical model), which is the simplest model for describing neuronal firing without the details of ionic movement through membranes [35]. A subtle difference is that the external driver is independent of the internal phase in the theta model, whereas the external driver (glucose) is governed by the internal phase (insulin secretion) of the islets in the modified rotator model. Our model is a phenomenological model that drastically simplifies insulin secretion, but precisely incorporates the observed phase modulation of insulin pulses. Therefore, our conclusion would likely be robust, independent of the detailed molecular mechanism of phase modulation. The simple model clearly revealed a link between phase modulation and inter-islet synchronization, and this link might not have been obviously captured in the complex models. The plateau fraction of the Ca 2+ oscillations was linked to their U-shaped dominant periods of the glucose concentration. A more surprising finding was that the plateau fraction and the U-shaped feature of single islets were linked to synchronization between islets.
To probe target molecules, e.g., voltage-gated Ca 2+ channel β 3 subunit in β cells [36], for phase modulation, however, more realistic models are required. An electrophysiological model has been developed to describe the islet Ca 2+ oscillations and explain inter-islet synchronization [17,18]. A realistic model will contribute to the formulation of testable hypotheses that apply insights from the simple rotator model. Currently, it is not feasible to noninvasively monitor insulin secretion from multiple islets that experience in vivo glucose regulation. Nevertheless, in an in vitro experiment, a glucose-insulin feedback loop was constructed using microfluidics, and the experiment demonstrated that the feedback loop could entrain islets to generate synchronous Ca 2+ oscillations [22].
These wave-like messengers, or hormones, can encode information in their amplitude, frequency, and phase. Hormonal frequency/phase modulation has received less attention than amplitude modulation. The frequency modulation of the luteinizing hormone for reproduction is one rare example [37,38]. As more hormone measurements at high temporal resolution become available, it is expected that frequency/phase modulation will be more ubiquitously observed in the oscillations of hormones other than insulin.

Experimental animals and islet preparation
All experimental procedures were approved by the Pohang University of Science and Technology Institutional Animal Care and Use Committee (POSTECH IACUC, Korea). C57BL/6J male mice (Jackson Laboratory, Bar Harbor, ME, USA) at 8-12 weeks of age were used. The animals were maintained with a 12 h light/dark cycle with free access to water. Mice were sacrificed by cervical dislocation under anesthesia with CO 2 . Islets were isolated from the pancreas of male C57BL/6J mice (25-30 g) after collagenase digestion. Briefly, 3 mL of 1 mg/mL collagenase P was injected into the pancreas via the common bile duct or by direct injection to expand the pancreas. The pancreas was excised and cut into small pieces, which were then digested with collagenase P to obtain free islets. Islets were then hand-picked under a dissecting microscope and placed in RPMI 1640 containing 11 mM glucose, 10% fetal bovine serum, and 1% penicillin-streptomycin. Islets were cultured at 37˚C in a 95% O 2 and 5% CO 2 mixture for one day before the experiment.

Measurements of Ca 2+ oscillations in islets
To monitor the intracellular concentration of the Ca 2+ ions, each batch of islets was incubated in a buffer solution containing 3 mM glucose with 2 μM fura-2/AM at 37˚C and 5% CO 2 for 30 min. After incubation, the islets were attached to coverslips using Puramatrix Hydrogel (BD Biosciences, Bedford, MA). The islets were placed in an open perifusion chamber (Live cell instrument, Seoul, South Korea). The chamber was mounted on an inverted fluorescence microscope (IX71, Olympus, Tokyo, Japan) and maintained at 37˚C. All of the buffer solutions with various glucose concentrations were delivered to the chambers containing the islets using a fast-flow solution-switching system (VC6; Warner Instruments, Hamden, CT, USA). The flow rate of this system was 3 mL/min, and all of the solutions were maintained at 37˚C. The intracellular Ca 2+ was measured using the 340/380 nm fluorescence ratio with a fluorescence microscope. The microscope was equipped with a 300 W xenon arc lamp (Lambda DG-4, Sutter Instruments, Novato, CA) containing the appropriate filters for excitation of fura-2 at 340 nm and 380 nm. Fluorescence images were acquired with a 100 ms exposure every 1 sec by a CCD camera (ImagEM, Hamamatsu Photonics, Hamamatsu, Japan). The 340/380 nm fluorescence ratios for all islets were obtained and analyzed using MetaFluor software (MDS Analytical Technologies, Sunnyvale, USA). Each islet was defined as a region of interest (ROI). Signals from the ROI were corrected by subtracting the background signal. Islets were perfused within the chamber with 3 mM glucose for 5 min prior to recording the fura-2/AM fluorescence. We used two protocols to measure the islet Ca 2+ oscillations under constant and alternating glucose conditions. The first protocol started with the glucose concentration at a basal level (3 mM) for 100 sec, followed by an increase to a higher glucose concentration (8, 10, 12, or 15 mM) for 3,000 sec and then a return to the basal glucose concentration for 120 sec. At the end, we depolarized the islets with 25 mM KCl for 250 sec to check their vitality. The second protocol started with the glucose concentration at the basal level for 100 sec, followed by a stimulation with a higher glucose concentration (8 or 10 mM) for 2,000 sec and then alternations of the glucose concentrations between 8 and 8/8.5/10 mM or 10 and 10/10.5/12 mM with a period of 160, 240, or 320 sec. The alternations were repeated 8 times, the glucose concentration returned to the basal level for 120 sec, and finally, the islets were depolarized with 25 mM KCl for 250 sec.

Data analysis
We removed a linear trend in the time traces of the Ca 2+ oscillations to focus on their phase by using the function "detrend" in MATLAB (MathWorks, Natick, MA). This modification set the mean of the time traces to zero. We computed the duration of continuous periods of positive values in the detrended time traces. The time fraction of the positive events was defined as the plateau fraction. We then applied the Fourier transformation to the detrended time traces and obtained their power spectra, and the maximum values defined their peak frequencies. We used the inverse value of the peak frequency as the dominant period. We used the time trace data only if their dominant periods were shorter than 10 min. To measure the degree of synchronization between islets, we adopted a stroboscopic approach that calculates the closeness of the phases between two oscillators [20,39]. Because the synchronization index (λ) was developed to calculate the coherence between two oscillators, we computed λ for every pair of islet Ca 2+ oscillations and defined their average value as our synchronization index [22]. Thus, λ = 0 and 1 represent complete desynchronization and synchronization, respectively.

Parameter estimation
We hypothesized that the islets had intrinsic Ca 2+ oscillations with a predetermined frequency (ω), but they generated the oscillations with different frequencies under different glucose conditions due to glucose-dependent phase modulation by μ(G). The modified frequencies at G = 8, 10, 12, and 15 mM glucose became o G ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Assuming that the intrinsic frequency (ω) followed a Gaussian distribution with a mean ( " o) and standard deviation (δω), the modified frequency (ω G ) should follow a certain distribution, Q(ω G ). We estimated the likelihood values of x ¼ ð " o; do; m 8 ; m 10 ; m 12 ; m 15 Þ that could make Q(ω G ) closest to the measured distribution, P(ω G ). To quantify the distance between the two distributions, we used the Kullback-Leibler divergence [40]: D(P || Q) = S ω P(ω)log P(ω)/Q(ω). We defined a total cost, C = S G D G (P || Q), for the four glucose conditions. Then, we conducted an importance sampling of x with a Monte-Carlo simulation and obtained 10 maximum likelihood values of x from~400,000 samples that produced a lower cost (C). Because the parameters were not entirely independent, we set μ 12 = 0 to fix the scale of ω using the evidence that phase modulation was absent in the Ca 2+ oscillations at 12 mM glucose (plateau fraction = 0.50±0.06). The final estimations obtained from the top 10 likelihood values of x were " o ¼ 0:305 AE 0:004 min -1 , δω = 0.037±0.003 min -1 , μ 8 = −0.163±0.011 min -1 , μ 10 = −0.028 ±0.020 min -1 , μ 12 = 0min -1 , and μ 15 = −0.136±0.018 min -1 . From these four values of μ G , we obtained the nonlinear function mðGÞ ¼ " msinh½ðG À G m 0 Þ=dG m , where " m ¼ 0:016 min -1 , G m 0 ¼ 11:6 mM, δG μ = 1.2 mM, and the hyperbolic function sinh(Z) = (e Z + e −Z )/2. For the fitting of the nonlinear function, we used the nonlinear least-squares Marquardt-Levenberg algorithm in GNUPLOT. For insulin secretion, we estimated its amplitude from the following published data: Bergsten [41], Henquin [42], and Gylfe [43]. We normalized the three data sets by establishing that insulin secretion became maximal (1 as the dimensionless value) at 20 mM. We hypothesized that the insulin secretion is determined by both the amplitude, r(G), and the shape factor, s(G). The phase dynamics determined the shape factor, CaðtÞdt. For the amplitude, we assumed a sigmoidal function, rðGÞ ¼ 0:25½1 þ tanh½ðG À G r 0 Þ=dG r , where tanh(Z) = (e Z −e −Z )/(e Z +e −Z ). We then estimated the maximum likelihood values of G r 0 ¼ 7:6 mM and δG r = 8.0 mM that minimized the mismatch, χ 2 = S G [hI(G)i−r(G)Ás(G)] 2 .

Model simulation
The Ca 2+ oscillation of the nth islet was described by the active rotator model, dy n dt ¼ o n À mðGÞcosy n with a Gaussian random number, o n ¼ " o AE do, and the pinning function, μ(G). Then, we used Ca n (t) = 1+cos θ n (t) as the normalized Ca 2+ concentration to examine its response to constant and alternating glucose concentrations. For the glucose-insulin feedback simulation, we defined the average insulin, I ¼ N À 1 X N n¼1 rðGÞð1 þ cosy n Þ, with an amplitude of r(G) for insulin secretion. In the simulation, we used N = 20 islets. Then, the glucose regulation followed dG dt ¼ G in À n I Á G: The above N + 1 differential equations were numerically integrated using the 2 nd order Runge-Kutta method with a sufficiently small time step, Δt = 0.01 [44].

Statistical analysis
Differences between unpaired observations were evaluated with Student's t-test. The data in the graphs and tables are presented as the means ± standard deviations. Findings were considered statistically significant at P<0.01.

Author Contributions
Conceptualization: BL KL JK SH SR JJ.