In our companion paper, the physiological functions of pancreatic β cells were analyzed with a new β-cell model by time-based integration of a set of differential equations that describe individual reaction steps or functional components based on experimental studies. In this study, we calculate steady-state solutions of these differential equations to obtain the limit cycles (LCs) as well as the equilibrium points (EPs) to make all of the time derivatives equal to zero. The sequential transitions from quiescence to burst–interburst oscillations and then to continuous firing with an increasing glucose concentration were defined objectively by the EPs or LCs for the whole set of equations. We also demonstrated that membrane excitability changed between the extremes of a single action potential mode and a stable firing mode during one cycle of bursting rhythm. Membrane excitability was determined by the EPs or LCs of the membrane subsystem, with the slow variables fixed at each time point. Details of the mode changes were expressed as functions of slowly changing variables, such as intracellular [ATP], [Ca2+], and [Na+]. In conclusion, using our model, we could suggest quantitatively the mutual interactions among multiple membrane and cytosolic factors occurring in pancreatic β cells.

## INTRODUCTION

In our companion paper (see Cha et al. in this issue), we constructed a detailed model of a pancreatic β cell based on published electrophysiological measurements of ion channels or exchangers. The model successfully reconstructed three representative electrical activities in response to a varying glucose concentration ([G]): the quiescent states of the membrane potential (Vm), bursting activity with alternate burst–interburst events, and continuous firing of action potentials. It was suggested that the burst–interburst cycle is generated by the interactions of channels or transporters with intracellular ions and/or metabolic intermediates. By applying lead potential (VL) analysis (Cha et al., 2009), we could quantify contributions of individual membrane currents to the slow changes in Vm during the interburst period, and suggested distinct ionic mechanisms underlying the bursting rhythm at different [G].

In our companion paper (Cha et al., 2011), the dynamic behavior of the model was calculated by time-based integration of ordinary differential equations. For example, [G]-dependent activity in a β cell was examined by integrating 18 differential equations until a constant pattern was observed at each [G]. However, this time-based simulation will never be more than an approximation because of the uncertainty that always remains in defining the steady states even after a long integration time, and in discriminating different patterns of membrane excitability explicitly. These problems may be overcome if steady-state solutions of the differential equations are obtained with respect to continuous variation in [G] or slowly varying cytosolic substances. To achieve this aim, we have applied bifurcation analysis to the comprehensive β-cell model developed in our companion paper.

Based on the steady-state solutions, the cellular responses in a β cell will be explicitly described in terms of “mode changes” of system behavior. We focused on three main questions: what kind of modes contribute to membrane excitability in β cells; when does each mode change occur during the burst–interburst rhythm; and, finally, which are the important intracellular factors underlying the mode changes? Collectively, with the results of time-based simulations and VL analysis in our companion paper (Cha et al., 2011), the results of the bifurcation analysis clarified the physiological roles of several intracellular factors in promoting modal changes in β-cell function. The results will be compared with those of the bifurcation analysis to a series of simple β-cell models reported by Chay, Keizer, and colleagues in the past few decades.

## MATERIALS AND METHODS

### The new β-cell model as a nondrift system with unique solutions

The structure and individual components of our β-cell model were fully described in our companion paper (Cha et al., 2011). In brief, the model is composed of 18 variables: Vm, seven gating variables of ion channels, three state variables of INaCa, four ionic concentrations in the cytosol ([Na+]i, [K+]i, and [Ca2+]i) or in ER ([Ca2+]ER), and three metabolic substrate concentrations ([ATP], [MgADP], and [Re]). Time-dependent changes of these variables are described by ordinary differential equations (refer to the supplemental material in Cha et al., 2011).

The new β-cell model is a nondrift system; that is, the full 18-variable system has steady-state solutions to make all the time derivatives in the model zero simultaneously. A necessary condition for the existence of these steady-state solutions is that all of the cation fluxes (Na+, K+, and Ca2+) through ion channels and exchangers should be included in calculating the time derivatives of both the membrane potential (dVm/dt) and the intracellular ion concentrations (d[X+]/dt). No intracellular ion concentration was fixed arbitrarily in the model. Additionally, to avoid unlimited changes of the metabolic compounds, the sum of NADH and NAD was set to be constant, like that of ATP and ADP. Using these approaches, the simulated responses of our β-cell model were all completely reversible.

To obtain the unique solutions in the full system (Fig. 1), the redundancy in calculating four ion concentrations ([Na+]i, [Ca2+]i, [K+]i, and [Ca2+]ER) and Vm was avoided by applying the charge conservation law:

$[Ca2+]ER=voli fER2 volER{Cmvoli F(Vm−Vm(0))−([Na+]i−[Na+]i(0))−([K+]i−[K+]i(0))−2fi([Ca2+]i−[Ca2+]i(0))}+[Ca2+]ER(0),$
(A8)

as derived in the Appendix. This procedure removed d[Ca2+]ER/dt from the model, leaving 17 independent variables. The initial value of each variable (for example, Vm(0)) are presented in Table S1 in our companion paper (Cha et al., 2011).

Figure 1.

Changes in EPs and LCs in the whole β-cell model by varying [G]. Four bifurcation diagrams showing continuous changes in EPs or LCs for Vm (A), [ATP] (B), [Ca2+]i (C), and [Na+]i (D) with respect to [G]. Stable EPs, unstable EPs, stable LCs, and unstable LCs are indicated by black, red, blue, and yellow lines, respectively. For LCs, the maximum and minimum values in oscillations were plotted. The amplitude of oscillation in [ATP] and [Na+]i was small, and the maximum and minimum curves of LCs fused with each other. AUTO failed to find unstable LCs at [G] < 9 mM. Open circles indicate bifurcation points; HB, Hopf bifurcation from a stable EP to an unstable EP (at 6.9 mM [G]); TR, Torus bifurcation from a stable LC to an unstable LC (at 18.84 mM [G]).

Figure 1.

Changes in EPs and LCs in the whole β-cell model by varying [G]. Four bifurcation diagrams showing continuous changes in EPs or LCs for Vm (A), [ATP] (B), [Ca2+]i (C), and [Na+]i (D) with respect to [G]. Stable EPs, unstable EPs, stable LCs, and unstable LCs are indicated by black, red, blue, and yellow lines, respectively. For LCs, the maximum and minimum values in oscillations were plotted. The amplitude of oscillation in [ATP] and [Na+]i was small, and the maximum and minimum curves of LCs fused with each other. AUTO failed to find unstable LCs at [G] < 9 mM. Open circles indicate bifurcation points; HB, Hopf bifurcation from a stable EP to an unstable EP (at 6.9 mM [G]); TR, Torus bifurcation from a stable LC to an unstable LC (at 18.84 mM [G]).

### Bifurcation analysis

In experimental studies using real β cells, a steady-state response to a new experimental condition is observed after a certain delay, for example, after applying a new level of [G]. The time-based simulation is straightforward to reconstruct this experimental protocol by repeating the integration of model differential equations with a tiny time step. In contrast to this time-based simulation, the bifurcation analysis directly solves multiple differential equations, which are responsible for the kinetic behavior of the model. Namely, the steady-state solutions of the model are directly obtained by setting all of the differential equations to zero and solving them simultaneously. The solutions are usually represented in a bifurcation diagram, which shows changes in equilibrium points (EPs) and limit cycles (LCs) when one parameter is systematically varied on the x axis, for example, [G] in Fig. 1 or PCaV in Fig. 2. An EP corresponds to a steady-state point at which the system permanently stays unless any perturbations are applied, and an LC is a steady-state periodic solution, for example, a sustained oscillation in Vm and substrate concentrations. In this paper, LCs are represented with a pair of lines that indicate the maximum and minimum values during an oscillation. The stability of an EP can be explicitly determined by an eigenvalue of a Jacobian matrix. A system eventually converges to a stable EP, such as the resting membrane potential and accompanying stable intracellular substrate concentrations. A system can also stay at a hypothetical unstable EP, but any perturbation will cause the system to leave from the EP, like a ball perfectly balanced on the peak of a hill. Similarly, an oscillation (LC) can also be stable or unstable.

Figure 2.

Mode changes of membrane excitability by varying PCaV. (A) Bifurcation diagrams showing continuous changes in Vm of EPs or LCs as a function of PCaV. Stable EPs, unstable EPs, stable LCs, and unstable LCs are indicated by black, red, blue, and yellow lines, respectively. The black line for a stable EP corresponds to the resting membrane potential, and the two blue lines for a stable LC correspond to the amplitude of the action potentials. The slow variables were fixed to the following values: [ATP] = 2.64 mM, [MgADP] = 0.0591 mM, [Re] = 0.61 mM, [Na+]i = 5.87 mM, [K+]i = 127 mM, [Ca2+]i = 0.124 µM, [Ca2+]ER = 0.0247 mM, fus = 0.843, and I1 = 0.152. Open circles indicate bifurcation points; LP, LP bifurcation; HB, Hopf bifurcation; PD, period doubling bifurcation. Black vertical lines with the blue italic numeral indicate control values of PCaV (48.9 pA mM−1) in the β-cell model. EP1, EP2, and EP3 (gray dots) are the intersections of the EP curve with the black vertical line. Gray vertical lines passing through the corresponding bifurcation points separate individual modes, as indicated at the top. (B) Steady-state I-V relationship. Zero current potentials correspond to EP1, EP2, and EP3 in A.

Figure 2.

Mode changes of membrane excitability by varying PCaV. (A) Bifurcation diagrams showing continuous changes in Vm of EPs or LCs as a function of PCaV. Stable EPs, unstable EPs, stable LCs, and unstable LCs are indicated by black, red, blue, and yellow lines, respectively. The black line for a stable EP corresponds to the resting membrane potential, and the two blue lines for a stable LC correspond to the amplitude of the action potentials. The slow variables were fixed to the following values: [ATP] = 2.64 mM, [MgADP] = 0.0591 mM, [Re] = 0.61 mM, [Na+]i = 5.87 mM, [K+]i = 127 mM, [Ca2+]i = 0.124 µM, [Ca2+]ER = 0.0247 mM, fus = 0.843, and I1 = 0.152. Open circles indicate bifurcation points; LP, LP bifurcation; HB, Hopf bifurcation; PD, period doubling bifurcation. Black vertical lines with the blue italic numeral indicate control values of PCaV (48.9 pA mM−1) in the β-cell model. EP1, EP2, and EP3 (gray dots) are the intersections of the EP curve with the black vertical line. Gray vertical lines passing through the corresponding bifurcation points separate individual modes, as indicated at the top. (B) Steady-state I-V relationship. Zero current potentials correspond to EP1, EP2, and EP3 in A.

Bifurcation refers to a behavioral change in a system, such as the change in the cellular electrical activity from the resting potential to the spontaneous action potential generation when the extracellular [G] is increased. A bifurcation is generally indicated with an emergence or disappearance of EPs or LCs, or with a change in their stability. The following kinds of bifurcation were observed in this study: limit point (LP) bifurcation: two EPs or LCs approach and annihilate each other; Hopf bifurcation: an EP loses stability and a new LC appears; period doubling bifurcation: the system switches to a new behavior with twice the period of the original system; and Torus bifurcation: an LC becomes unstable, and the system oscillates around the unstable LC. Mathematical details of classifying bifurcations could be skipped. In this study, the steady-state solutions were numerically obtained using AUTO implemented in XPPAUT (Ermentrout, 2002), a specific computational tool for bifurcation analysis. For a more extensive explanation about the bifurcation analysis to the cellular models, see Fall et al. (2002).

### Fast and slow decomposition of the system to determine membrane excitability

Our companion paper (Cha et al., 2011) demonstrated that the ionic mechanism for membrane excitation is continuously modified by the time-dependent variations in cytosolic substrates during burst–interburst rhythm. If the rate of these time-dependent changes in cytosolic substrate concentrations is slow enough to scarcely affect the configuration of individual action potentials, but could determine the burst–interburst cycle, the membrane excitability at a given time point might be examined by fixing the substrate concentrations. We applied bifurcation analysis by fixing the intracellular substrate concentrations ([S]is) and [Ca2+]ER at a specific time point of the burst–interburst rhythm to determine the modes of the membrane excitability. After fixing these concentrations, however, we occasionally found that a slow change in the membrane excitability still occurred solely because of the time-dependent changes of the ultraslow inactivation gate of ICaV (fus) (Fig. S2 A). When fus was fixed in addition to the substrate concentrations, the membrane excitability stayed in the same state as observed at the fixing moment. That is, a steady rhythm of repetitive action potentials or a steady resting potential was established depending on the fixed time (Fig. S2 B). Thus, fixing these variables largely satisfied the requirement for fast and slow decomposition to define the membrane excitability at a given time point. Based on the range of time constants during bursting rhythm (Fig. S2 C), the inactivation state of INaCa (I1) was also classified as a slow variable because it has a comparable time constant to the fus. [Ca2+]i oscillated rapidly over the time span of a single action potential (fast Ca2+ ripple) superimposed on a slow drift of the Ca2+ plateau during the burst (inset of Fig. 4 A or Fig. 2 in Cha et al., 2011). In this study, [Ca2+]i was treated as a slow variable because it took ∼4 s for [Ca2+]i to reach a new steady state when other slow variables including [Ca2+]ER were fixed. As a consequence, to determine the mode of membrane excitability, we calculated steady-state solutions for the membrane system consisting of nine fast variables (Vm, dCaV, UCaV, rKDr, qKDr, mKCa(BK), hKCa(BK), Ei_total, and I2), after fixing nine slow variables ([ATP], [MgADP], [Re], [Na+]i, [K+]i, [Ca2+]i, and [Ca2+]ER, in addition to fus and I1) at a given time of observation. Refer to Fig. S1 or the supplemental material in our companion paper (Cha et al., 2011) for definitions of individual variables.

### Online supplemental material

Fig. S1 shows the interactions between model variables in the cytosolic and membrane systems. Fig. S2 shows model behaviors when all the intracellular concentrations ([Si]s) are fixed, or when fus was also fixed together with [Si]s. Fig. S2 C demonstrates the time constants of the model variables during bursting. Fig. S3 shows a simulation that is started from an unstable EP at 16 mM [G]. The source code of this model is provided as a PDF file for XPPAUT (Ermentrout, 2002). The supplemental material is available.

## RESULTS

### Mode changes in cellular activity induced by varying [G]

Extracellular glucose affects [ATP] and [MgADP] through ATP production pathways; the equilibrium values of the remaining variables are readjusted through reciprocal interactions among the metabolic compounds, membrane excitation, and intracellular ion concentrations (Fig. S1). Accordingly, different patterns (or modes) of cellular activity are established by varying [G], as shown by the time-based simulation in our companion paper (Cha et al., 2011). In the bifurcation analysis, underlying mode changes are revealed explicitly by the steady-state solutions (EPs and LCs) of the differential equations for all 17 variables in the model. Fig. 1 demonstrates changes in Vm, [ATP], [Ca2+]i, and [Na+]i of EPs or LCs when [G] is varied. At [G] < 6.9 mM, one stable EP exists (Fig. 1, black lines), predicting that the cellular state will always return exactly to this EP regardless of the initial condition. Namely, for [G] < 6.9 mM, EP defines the resting membrane potential and the corresponding steady-state values of substrate concentrations. For 6.9 < [G] < 30 mM, EP becomes unstable (Fig. 1, red lines) indicating that spontaneous bursts of action potentials start to take place. Around the unstable EP, an LC is present represented by the yellow and blue lines (Fig. 1), which give the peak and the most negative potentials of individual action potential in the case of stable LCs. For 6.9 < [G] < 18.84 mM, the LC is unstable and the action potential configuration slowly varies during the intermittent spike burst period. The LC converts from unstable to stable at [G] = 18.84 mM, predicting the continuous spike burst in the time-based simulation. (refer to Fig. 2 of Cha et al., 2011). Thus, the sequential transitions from quiescence to burst–interburst oscillations and then to continuous firing with increasing [G] were defined objectively by the EPs or LCs obtained by solving the whole set of equations.

The unstable EP for 6.9 < [G] < 30 mM (Fig. 1, red lines) indicates a hypothetical condition, under which the membrane potential and the intracellular composition of ions and substrates can remain constant unless any perturbation is applied. The unstable EP (Fig. 1, red lines) is usually difficult to observe in experiments using real cells, or in time-based simulations, but its existence can be convinced by starting the simulation from that particular set of values of all variables defined by the EP in the bifurcation analysis. An example is shown in Fig. S3 A, where Vm remained near the initial values for ∼1 s before spontaneous activity started. This is typical behavior for any unstable EP. The “spontaneous” deviation from the EP is most probably a result of numerical errors in integration. An alternative way to put the system on the unstable EP in real experiments or in the time-based simulation might be provided by the voltage-clamp experiments. Fig. S3 B shows a protocol of clamping Vm at the equilibrium potential given by the unstable EP for [G] = 16 mM. Under the voltage clamp, values of the remaining 16 variables, including [ATP], [Ca2+]i, and [Na+]i, spontaneously relaxed to their equilibrium levels defined by the unstable EP (Fig. S3 C). This response is expected for a stable EP but not for an unstable EP. We ran the bifurcation analysis again and found that EP became stable when Vm, the pivotal factor in the membrane system, was fixed (not depicted). When the membrane was released from the voltage clamp (Fig. S3 B), the EP became unstable again, and Vm escaped after a similar latency as in Fig. S3 A. Obviously, voltage-clamp experiments are awaited to confirm this prediction in real cells.

The [G]–EP relationships (Fig. 1, black and red lines) give important clues on the principal mechanisms of the glucose-induced signal transductions, free from the cyclic oscillation in time-based simulations. With increasing [G], the accelerated ATP production raises the equilibrium level of [ATP] in a dose-dependent manner (Fig. 1 B). Then, the IKATP conductance continuously decreases, and the resulting membrane depolarization (Fig. 1 A) increases the equilibrium level of [Ca2+]i via a fractional activation of ICaV, even though ISOC is deactivating (Fig. 1 C). Elevated [Ca2+]i accelerates ATP consumption, resulting in a small inflection in the equilibrium [ATP]–[G] curve (Fig. 1 B). The biphasic relationship between [Na+]i and [G] is determined by two factors (Fig. 1 D). The initial falling phase in [Na+]i is a result of activation of the Na+/K+ pump through an increase in [ATP], and the rising phase is attributed to the increase in Na+ influx via the forward mode of Na+/Ca2+ exchange, which is accelerated by increased [Ca2+]i. With the saturation of ATP production >20 mM [G], all relationships become flat.

### Mode changes in membrane excitability induced by slow changes in the intracellular substrates

Time-base simulations, as well as experimental studies, have led to the hypothesis that membrane excitability is cyclically modulated by slow changes in intracellular substrates during the burst–interburst rhythm. To prove this hypothesis mathematically, we separated the membrane system from the cytosolic factors by adopting the same strategy of the fast and slow decomposition of the system variables, as has been established in previous bifurcation analyses of simple β-cell models (see Materials and methods). The nine slow variables were fixed at each time point of the bursting activity when solving the differential equations of the nine fast variables to obtain EPs or LCs in the fast membrane system. Based on the EPs and LCs, we can define the membrane excitability, as described below.

#### Definition of modes of membrane excitability.

In discriminating modes of membrane excitability in our new β-cell model, we found it convenient to calculate EPs and LCs by varying the amplitude factor (PCaV) of ICaV, a dominant current in generating action potentials. The bifurcation diagram showed a typical S-shaped curve of EP in the PCaV–Vm plane (Fig. 2 A). The black and red lines indicate stable and unstable EPs, as in Fig. 1, and LCs are shown by blue (stable) or yellow (unstable) lines. The black vertical line in Fig. 2 A indicates the standard PCaV (48.9 pA mM−1), and the intersections with the bifurcation diagram correspond to the EPs or LCs in the control system.

Six distinct modes of membrane excitability (Fig. 2 A, top of the panel, from A to F) were defined by finding the values of PCaV for each bifurcation point (Fig. 2 A, gray lines). Each mode has a different number of stable or unstable EPs and LCs, and shows specific electrophysiological characteristics in response to current injections into the cell (Fig. 3 and summarized in Table I). Modes A and F with extremely small or large PCaV are nonexcitable. Mode B has two unstable EPs and one stable EP, and a single action potential is induced by an electrical stimulus. Mode C has two stable states with spontaneous action potentials (stable LC) and a resting potential (stable EP). Accordingly, the stable firing can be switched on and off by applying a depolarizing or a hyperpolarizing current pulse, respectively. Mode D with a stable LC shows a stable oscillation regardless of stimuli. Mode E is also bistable with stable oscillations and a depolarized resting potential.

Table I.

Membrane excitability

 Mode EPs or LCs Physiology A 1 sEP Nonexcitable mode B 1 sEP, 2 uEP Single action potential mode C 1 sEP, 2 uEP, sLC Bistable mode with a stable firing and a quiescent state D 1 uEP, sLC Stable firing mode E 1 sEP, sLC, uLC Bistable mode with a stable firing and a quiescent state in depolarized potential F 1 sEP Nonexcitable mode with damped oscillation
 Mode EPs or LCs Physiology A 1 sEP Nonexcitable mode B 1 sEP, 2 uEP Single action potential mode C 1 sEP, 2 uEP, sLC Bistable mode with a stable firing and a quiescent state D 1 uEP, sLC Stable firing mode E 1 sEP, sLC, uLC Bistable mode with a stable firing and a quiescent state in depolarized potential F 1 sEP Nonexcitable mode with damped oscillation

sEP, stable EP; uEP, unstable EP; sLC, stable LC; uLC, unstable LC.

Figure 3.

Representative membrane responses to current injections in different modes. The amplitudes of injected current pulses are indicated inside panels. Duration of the pulse was 5 ms, except 15 ms in Mode E. In each panel, a different value of PCaV was used for simulation (units in pA mM−1): Mode A, PCaV was 30; Mode B, 45; Mode C, 48.9; Mode D, 60; Mode E, 65; Mode F, 69. The same values were used for the slow variables as in Fig. 2. The time axes were identical for all panels.

Figure 3.

Representative membrane responses to current injections in different modes. The amplitudes of injected current pulses are indicated inside panels. Duration of the pulse was 5 ms, except 15 ms in Mode E. In each panel, a different value of PCaV was used for simulation (units in pA mM−1): Mode A, PCaV was 30; Mode B, 45; Mode C, 48.9; Mode D, 60; Mode E, 65; Mode F, 69. The same values were used for the slow variables as in Fig. 2. The time axes were identical for all panels.

In the PCaV–Vm plane in Fig. 2 A, the membrane excitability at the standard PCaV is determined to be at mode C under the given set of slow variables. The corresponding bistable membrane excitation is obvious from the N-shaped I-V diagram with three intersections, drawn with the same set of slow variables (Fig. 2 B). The zero current potentials, EP1, EP2, and EP3 on the I-V curve, correspond to the intersections of the EP curve with the black vertical lines at PCaV of 48.9 pA mM−1 in Fig. 2 A. The resting Vm is at EP1, and action potentials, if triggered by an appropriate stimulus, oscillate around EP3. The I-V curve, however, fails to determine if the action potential is repetitive or singular because of lack of information about LCs.

#### Time-dependent mode changes at 8 mM [G].

To disclose specific modes of membrane excitability during the burst–interburst rhythm at 8 mM [G], bifurcation diagrams (Vm vs. PCaV) were obtained using the values of [S]is, [Ca2+]ER, fus, and I1 at successive time points, measured from the time-based simulation shown in Fig. 4 A (gray line). Fig. 4 B shows bifurcation diagrams at six representative time points during spontaneous bursting activity. EPs or LCs were measured at the standard value of PCaV (48.9 pA mM−1, black vertical lines) in individual bifurcation diagrams and were plotted together with the result of time-based simulation (Fig. 4 A). Finally, we determined the mode of membrane excitability at each time point on the basis of the definitions summarized in Table I and showed sequential mode changes at the top of Fig. 4. Mode C (bistable mode) observed at the beginning of the record is followed by mode D (stable firing mode) at 7.6 s, and action potentials are initiated after a delay. The burst was terminated by extinction of the stable LC through a mode change from D to B (single action potential mode) via the temporal appearance of mode C. After cessation of the burst, the membrane returns again to mode C. It should be noted that [Ca2+]i fluctuated rapidly during the burst, as shown in the inset in Fig. 4 A. Because of these rapid changes in [Ca2+]i, the position of the EP or LC curves fluctuated slightly along the abscissa in Fig. 4 B, and we obtained bifurcation diagrams at the extremes of [Ca2+]i. The mode changes from D to C or from C to B occurred slightly earlier at the local minimums of the Ca2+ transient than at the maximums.

Figure 4.

Time-dependent mode changes in membrane excitability during one burst–interburst cycle. (A) Time-based simulation of Vm at 8 mM [G] (gray continuous line). Colored dots are EP1, EP2, EP3, and LC (min and max) in the membrane system, which were measured from the bifurcation diagrams calculated with fixed values of [S]is, [Ca2+]ER, fus, and I1 at corresponding time points. Stable EPs, unstable EPs, stable LCs, and unstable LCs are indicated by black, red, blue, and yellow dots, respectively. The unstable LC (yellow) was only observed at the moment of mode change from Mode B to C (19.5 s) during <0.1 s, whereas it was not observed at the switch from Mode C to B during the burst (15.5 s). During the burst period, two sets of EPs and LCs were demonstrated at the sequential maximum or minimum of [Ca2+]i. These two EP3s are almost superimposed in the figure. The mode of membrane excitability at each moment is indicated at the top. (Inset) Trace of [Ca2+]i during the same burst–interburst cycle. (B) Bifurcation diagrams at six representative time points in A, as indicated in the top left part of the figure. During the burst, three time points were selected from those of sequential minimums of [Ca2+]i (11, 14.6, and 15.5 s). The same color codes were used for the dots as those in A. Black vertical lines were drawn at PCaV = 48.9 pA mM−1, the standard value in the β-cell model.

Figure 4.

Time-dependent mode changes in membrane excitability during one burst–interburst cycle. (A) Time-based simulation of Vm at 8 mM [G] (gray continuous line). Colored dots are EP1, EP2, EP3, and LC (min and max) in the membrane system, which were measured from the bifurcation diagrams calculated with fixed values of [S]is, [Ca2+]ER, fus, and I1 at corresponding time points. Stable EPs, unstable EPs, stable LCs, and unstable LCs are indicated by black, red, blue, and yellow dots, respectively. The unstable LC (yellow) was only observed at the moment of mode change from Mode B to C (19.5 s) during <0.1 s, whereas it was not observed at the switch from Mode C to B during the burst (15.5 s). During the burst period, two sets of EPs and LCs were demonstrated at the sequential maximum or minimum of [Ca2+]i. These two EP3s are almost superimposed in the figure. The mode of membrane excitability at each moment is indicated at the top. (Inset) Trace of [Ca2+]i during the same burst–interburst cycle. (B) Bifurcation diagrams at six representative time points in A, as indicated in the top left part of the figure. During the burst, three time points were selected from those of sequential minimums of [Ca2+]i (11, 14.6, and 15.5 s). The same color codes were used for the dots as those in A. Black vertical lines were drawn at PCaV = 48.9 pA mM−1, the standard value in the β-cell model.

Fig. 4 B indicates that a sequence of bifurcation occurred during one cycle of bursting rhythm. At 0 s, Vm remains stable at EP1. During the initial 7.6 s, the diagram shifts leftward with reference to the standard PCaV, and EP1 and EP2 approach each other. At 7.6 s, EP1 coalesces with EP2 and disappears (LP bifurcation of EPs), and thus the membrane system shifts to the stable LC, corresponding with spontaneous action potentials. It takes about another 3 s for Vm to move from EP1 to LC because of the very small inward current charging the membrane capacitance. Once spontaneous oscillations are initiated, the diagram shifts rightward until the stable LC disappears by an LP bifurcation at 15.5 s. Finally, Vm returns to the stable EP1, and the whole cycle of events repeats again. Here, it should be noted that the movement of the bifurcation diagram is entirely the result of time-dependent changes in the slow intracellular factors.

### Which slow variables are responsible for termination of the burst?

Termination of the burst occurs from the disappearance of the stable LC (LP bifurcation), accompanied by the change in the membrane excitability from Mode C to B at 8 mM [G]. To examine the role of the slow variables in terminating the burst, we constructed a bifurcation diagram by treating a single slow factor as a bifurcation parameter. The other slow variables were fixed at the values obtained at 14.6 s, just before the LP bifurcation. In Fig. 5 A, the bifurcation parameter is [ATP]. A stable LC (Fig. 5 A, blue lines) appeared over a higher range of [ATP]. During the burst period, [ATP] gradually decreased, as indicated by gray vertical lines sampled every 1 s. This monotonic decrease in [ATP] clearly promotes the burst termination through the membrane system approaching the LP bifurcation point. As demonstrated in our companion paper (Cha et al., 2011), the terminating effect of [ATP] was mediated through gradual activation of outward IKATP. In Fig. 5 B, the bifurcation diagram was calculated with respect to the ultraslow inactivation gate of ICaV (fus). The value of fus also decreased toward the LP bifurcation point, indicating that fus was also a key factor in the burst termination. Namely, the decrease of fus reduces the inward ICaV and leads to gradual hyperpolarization. The concurrent reduction of Ca2+ influx through ICaV lowers the plateau level of [Ca2+]i, overcoming the opposite effect of the accumulation of Ca2+ in the ER. After [Ca2+]i reached a peak at 12 s, it facilitates the burst termination (Fig. 5 C) by decreasing inward INaCa and ITRPM. The accumulation of intracellular Na+ also had a significant effect on the termination of the burst through INaK at 8 mM [G] (Fig. 5 D). In contrast, [K+]i and I1 (inactivation state of INaCa) have minor and opposite effects on the mode change (Fig. 5, E and F). It should be noted that the relative importance of the slow factors on membrane excitability may be altered at a different [G].

Figure 5.

Effects of individual slow variables on the mode change toward termination of the burst. Each bifurcation diagram was obtained by varying single slow variable, [ATP] (A), fus (B), [Ca2+]i (C), [Na+]i (D), [K+]i (E), and I1 (F) as a bifurcation parameter, with the remaining slow variables fixed at their values at 0.9 s before the LP bifurcation (open circles). Stable EPs, unstable EPs, stable LCs, and unstable LCs are indicated by black, red, blue, and yellow lines, respectively. Six gray vertical lines in A–D indicate the values of the corresponding slow variables sampled at an interval of 1 s (from 11 to 16 s in Fig. 4 A). In C, the sequential values of the minimum [Ca2+]i were sampled. In E and F, two gray lines indicate the corresponding values at 11 and 16 s. Arrows represent the sampling sequence. We confirmed that the bifurcation diagrams were qualitatively the same when slow variables were fixed at different time points during the burst.

Figure 5.

Effects of individual slow variables on the mode change toward termination of the burst. Each bifurcation diagram was obtained by varying single slow variable, [ATP] (A), fus (B), [Ca2+]i (C), [Na+]i (D), [K+]i (E), and I1 (F) as a bifurcation parameter, with the remaining slow variables fixed at their values at 0.9 s before the LP bifurcation (open circles). Stable EPs, unstable EPs, stable LCs, and unstable LCs are indicated by black, red, blue, and yellow lines, respectively. Six gray vertical lines in A–D indicate the values of the corresponding slow variables sampled at an interval of 1 s (from 11 to 16 s in Fig. 4 A). In C, the sequential values of the minimum [Ca2+]i were sampled. In E and F, two gray lines indicate the corresponding values at 11 and 16 s. Arrows represent the sampling sequence. We confirmed that the bifurcation diagrams were qualitatively the same when slow variables were fixed at different time points during the burst.

## DISCUSSION

In the new model developed in our companion paper (Cha et al., 2011), the dynamics of a pancreatic β cell were described with 18 differential equations containing individual functional components. The present study focused on solving the differential equations directly, using bifurcation analysis. Examination was performed at two different levels: the entire system (Fig. 1) to investigate the whole cell response to varying [G], and at the level of the fast membrane subsystem (Figs. 2, 4, and 5) to explore the time-dependent changes in membrane excitability at a given [G]. The results of bifurcation analysis were supplemented by time-based simulations (Fig. 3 and Fig. 2 in Cha et al., 2011) to deduce physiologically relevant conclusions. Based on the number and stability of the steady-state solutions (EPs and LCs), we discriminated three different regions depending on [G]: quiescent, bursting, and continuous firing activity (Fig. 1). We also discriminated six kinds of modes of membrane excitability depending on slow cytosolic factors (Fig. 2). The exact timing of transition between modes of membrane excitability was also indicated on time-based simulation. We investigated the ionic mechanisms for the initiation of the burst with lead potential analysis in our companion paper (Cha et al., 2011), and those for the termination of the burst with bifurcation analysis in this paper. Thus, the two papers together provide a novel mathematical description about the β-cell bursting activity and the role of individual functional units or molecules in the response to glucose.

### Comparison with previous studies in respect to bifurcation analysis

In the past few decades, bifurcation analysis has been used to clarify the mathematical principles underlying bursting activity in pancreatic β cells. In most studies, extremely simplified models have been used because they were amenable to mathematical analysis. In this study, we applied bifurcation analysis to a complex β-cell model based on extensive experimental data, with the aim of clarifying important physiological mechanisms in reference to individual functional components. We demonstrated that the bursting activity in this complex system is generated with the same basic principles as before.

#### Interactions between fast and slow systems.

Previous modeling studies consistently indicated that the bursting rhythm is generated by the reciprocal interactions between fast and slow subsystems. These studies presented useful hypotheses for a negative feedback mechanism between an ionic current with one single slow variable. For example, early β-cell models hypothesized an interaction between [Ca2+]i and Ca2+-dependent activation of IKCa (Chay and Keizer, 1983; Sherman et al., 1988) or Ca2+-dependent inactivation of ICaV (Chay and Kang, 1988). Subsequently, other hypothetical mechanisms include voltage-dependent slow inactivation of ICa (Chay, 1990; Keizer and Smolen, 1991), [ATP], or [ADP] to modulate the conductance of IKATP (Smolen and Keizer, 1992), or [Ca2+]ER to activate ISOC (Chay, 1996, 1997). With our detailed cell model, we demonstrated that the following multiple slow variables work in concert to generate the burst–interburst rhythm: slow inactivation (fus) via ICaV, intracellular ATP metabolism via IKATP, [Ca2+]i via INaCa and ITRPM, and [Na+]i via INaK. These mechanisms make a positive contribution to termination of the burst against the minor and opposite effects of [K+]i or I1 (inactivation state of INaCa). Moreover, our model showed that even a single slow variable has complex influences on the membrane excitability. For example, [ATP] promotes the burst termination by activating IKATP but, at the same time, retards it by depressing INaK; [Ca2+]i stabilizes the oscillation during the initial phase of the burst but terminates the burst in the late phase.

#### Transitions between an EP and an LC.

Results in this study are consistent with the conclusions of previous studies, namely that burst–interburst rhythm at a given [G] can be explained by repetitive transitions between a stable EP and a stable LC, although our graphical analysis is largely different. In the most successful and straightforward presentation using a simple model (Sherman et al., 1988), a bifurcation diagram was obtained using [Ca2+]i as the bifurcation parameter, which was the sole slow variable in their model. The bifurcation diagram was superimposed with a Ca2+ nullcline in Vm–[Ca2+]i space, and thereby, the overall behavior of the system was easily tracked along the Vm–[Ca2+]i diagram guided by both the Ca2+ nullcline and bifurcation points. In the study of Bertram and Sherman (2004), their model was composed of three slow variables, [Ca2+]i, [Ca2+]ER, and [ADP], and bifurcation diagrams were presented with respect to [Ca2+]i, by fixing the other two variables. Therefore, the combined effects of multiple slow variables were only indirectly inferred by comparing bifurcation diagrams calculated with two different values of [ADP] or [Ca2+]ER. Our model is an even more complex system, with eight slow variables, so it was extremely difficult to prepare such a straightforward bifurcation diagram. To overcome this problem, we developed an alternative way to show the transition of Vm between an EP and an LC along the time axis. To show the net effects of the slow variables during the bursting rhythm, we used the values of all the slow variables at each time point to calculate EPs and LCs. In our presentation, time-dependent changes in the mode of membrane excitability were explicitly identified on the record of the time-based simulation (Fig. 4 A).

#### Modal changes in behavior of the whole system.

We found that the bursting activity ceased by decreasing [G], and the firing became uninterrupted by increasing [G] in our β-cell model. These transitions from quiescence to bursting or to continuous firing in the whole system were consistently observed in several previous studies, but only indirectly. Namely, changes in glucose were mimicked by increasing either the rate of cytosolic Ca2+ sequestration (Chay and Keizer, 1983) or the Ca2+-binding affinity to Ca2+ channels (Chay, 1993). In more complex models, a hypothetical parameter proportional to the proton motive force (Keizer and Magnus, 1989) or the conductance of IKATP (Smolen and Keizer, 1992; Bertram et al., 1995) was changed, instead of calculating the reaction pathways for the glucose signal. Therefore, these simple models did not directly address the key question of how changes in [G] modulate cellular activity. In this study, we simulated the whole spectrum of [G] dependency by implementing the changes in metabolic status after changes in [G] (Fig. 1). This enabled us to estimate the values of [G] at which the cell changes its electrophysiological characteristics. Although we successfully reproduced cellular response to a range of [G] relevant to experimental measurements in the mouse, improvements in the metabolic components are still awaited to get a deeper insight into the effects of glucose. In the future, we would especially like to examine the effects of [Ca2+]i on the production of ATP through the tricarboxylic acid cycle and oxidative phosphorylation.

#### A wide range of cycle length in bursting activity.

In a simple model possessing a single slow variable, the range of burst cycle length was rather limited if compared with the experimental observations. Bertram et al. (2000) demonstrated that burst cycle length can be varied over a wider range by including two different slow processes, one with a relatively small time constant (s1) and the other with a much larger time constant (s2). Three slow variables were assigned to [Ca2+]i and [ATP] or [Ca2+]ER in a subsequent study (Bertram and Sherman, 2004). In the present study, we used nine slow variables, including seven substrate concentrations as well as the slow gating of ICaV and INaCa. This resulted in a wide range of burst cycle length when external [G] was changed (Fig. 2 in Cha et al., 2011). Based on the theory of Bertram et al. (2000), comparing the time courses of the slow variables implies that [ATP] ([MgADP]) or fus might correspond to s1, and [Na+]i or [Ca2+]ER to s2, in our model. That is, the duration of one cycle of the bursting rhythm is short at a low [G] (8 mM) where the variations in [ATP] or fus were predominant, whereas relatively large variations in [Na+]i or [Ca2+]ER govern a much slower rhythm at a high [G] (16 mM). It should be noted that the rate of change in [ATP] was largely affected by the Ca2+-dependent consumption of ATP (kATP,Ca; Eq. S103 in Cha et al., 2011), and that of [Na+]i was determined by NaK or NaCa activity in our model. Thus, for more reliable reproduction of the experimental burst duration, the above parameters should be improved in the future, when more extensive experimental data are available.

### Conclusion

In conclusion, collectively with the time-based integrations, steady-state solutions of our differential equations explicitly proved the hypothesis about the roles of individual ion channels and transporters in generating β-cell electrical activity in relation to their complex interactions with slow variables. Furthermore, working hypotheses for new experiments can be obtained from a mathematical model with detailed membrane components and cytosolic mechanisms. Although quantitative predictions from any mathematical model are dependent on how correctly the individual model components are described, this study and our companion paper (Cha et al., 2011), using bifurcation analysis, lead potential analysis, and time-based simulations, provide a framework to improve an objective understanding of this complex system.

## Appendix

$dVmdt=−INa,tot+ICa,tot+IK,tot+IinjectCm$
(A1)
$d[Na+]idt=−INa,totvoli F or −INa,tot=voli F d[Na+]idt$
(A2)
$d[K+]idt=−IK,tot−Iinjectvoli F or −IK,tot−Iinject=voli F d[K+]idt$
(A3)
$d[Ca2+]idt=fivoli(−ICa,tot2F−JSERCA+Jleak)$
(A4)
$d[Ca2+]ERdt=fERvolER(JSERCA−Jleak)$
(A5)

By combining Eqs. A4 and A5,

$d[Ca2+]idt=fivoli(−ICa,tot2 F−d[Ca2+]ERdtvolERfER),or−ICa,tot=2F(volifid[Ca2+]idt+volERfERd[Ca2+]ERdt).$
(A6)

By combining Eqs. A1, A2, A3, and A6,

$d[Ca2+]ERdt=voli fER2⋅volER(Cmvoli FdVmdt−d[Na+]idt−d[K+]idt−2fid[Ca2+]idt).$
(A7)

By integrating both terms with t from t = 0,

$[Ca2+]ER=voli fER2 volER{Cmvoli F(Vm−Vm(0))−([Na+]i−[Na+]i(0))−([K+]i−[K+]i(0))−2fi([Ca2+]i−[Ca2+]i(0))}+[Ca2+]ER(0).$
(A8)

## Acknowledgments

We thank Professor T. Powell for fruitful discussion and for improving the English of this paper.

This work was supported by the Biomedical Cluster Kansai project; a Grant-in-Aid (22590216 to C.Y. Cha and 22390039 to A. Noma) from the Ministry of Education, Culture, Sports, Science and Technology of Japan; and the Ritsumeikan-Global Innovation Research Organization at Ritsumeikan University.

Lawrence G. Palmer served as editor.

## References

References
Bertram
R.
,
Sherman
A.
.
2004
.
A calcium-based phantom bursting model for pancreatic islets
.
Bull. Math. Biol.
66
:
1313
1344
.
Bertram
R.
,
Butte
M.J.
,
Kiemel
T.
,
Sherman
A.
.
1995
.
Topological and phenomenological classification of bursting oscillations
.
Bull. Math. Biol.
57
:
413
439
.
Bertram
R.
,
Previte
J.
,
Sherman
A.
,
Kinard
T.A.
,
Satin
L.S.
.
2000
.
The phantom burster model for pancreatic beta-cells
.
Biophys. J.
79
:
2880
2892
.
Cha
C.Y.
,
Himeno
Y.
,
Shimayoshi
T.
,
Amano
A.
,
Noma
A.
.
2009
.
A novel method to quantify contribution of channels and transporters to membrane potential dynamics
.
Biophys. J.
97
:
3086
3094
.
Cha
C.Y.
,
Nakamura
Y.
,
Himeno
Y.
,
Wang
J.
,
Fujimoto
S.
,
Inagaki
N.
,
Earm
Y.E.
,
Noma
A.
.
2011
.
Ionic mechanisms and Ca2+ dynamics underlying the glucose sensing in pancreatic β cells: a simulation study
.
J. Gen. Physiol.
138
:
21
37
.
Chay
T.R.
1990
.
Effect of compartmentalized Ca2+ ions on electrical bursting activity of pancreatic beta-cells
.
Am. J. Physiol.
258
:
C955
C965
.
Chay
T.R.
1993
.
The mechanism of intracellular Ca2+ oscillation and electrical bursting in pancreatic beta-cells
.
Adv. Biophys.
29
:
75
103
.
Chay
T.R.
1996
.
Electrical bursting and luminal calcium oscillation in excitable cell models
.
Biol. Cybern.
75
:
419
431
.
Chay
T.R.
1997
.
Effects of extracellular calcium on electrical bursting and intracellular and luminal calcium oscillations in insulin secreting pancreatic beta-cells
.
Biophys. J.
73
:
1673
1688
.
Chay
T.R.
,
Kang
H.S.
.
1988
.
Role of single-channel stochastic noise on bursting clusters of pancreatic beta-cells
.
Biophys. J.
54
:
427
435
.
Chay
T.R.
,
Keizer
J.
.
1983
.
Minimal model for membrane oscillations in the pancreatic beta-cell
.
Biophys. J.
42
:
181
190
.
Ermentrout
B.
2002
.
Simulating, Analyzing, and Animating Dynamical System: A Guide to XPPAUT for Researchers and Students
.
SIAM Press
,
Philadelphia
.
290 pp
.
Fall
C.P.
,
Marland
E.S.
,
Wagner
J.M.
,
Tyson
J.J.
.
2002
.
Computational Cell Biology
.
Springer-Verlag
,
New York, Inc.
445 pp
.
Keizer
J.
,
Magnus
G.
.
1989
.
ATP-sensitive potassium channel and bursting in the pancreatic β cell. A theoretical study
.
Biophys. J.
56
:
229
242
.
Keizer
J.
,
Smolen
P.
.
1991
.
Bursting electrical activity in pancreatic beta cells caused by Ca2+- and voltage-inactivated Ca2+ channels
.
Proc. Natl. Acad. Sci. USA.
88
:
3897
3901
.
Sherman
A.
,
Rinzel
J.
,
Keizer
J.
.
1988
.
Emergence of organized bursting in clusters of pancreatic β-cells by channel sharing
.
Biophys. J.
54
:
411
425
.
Smolen
P.
,
Keizer
J.
.
1992
.
Slow voltage inactivation of Ca2+ currents and bursting mechanisms for the mouse pancreatic β-cell
.
J. Membr. Biol.
127
:
9
19
.

Abbreviations used in this paper:

• EP

equilibrium point

•
• LC

limit cycle

•
• LP

limit point

This article is distributed under the terms of an Attribution–Noncommercial–Share Alike–No Mirror Sites license for the first six months after the publication date (see http://www.rupress.org/terms). After six months it is available under a Creative Commons License (Attribution–Noncommercial–Share Alike 3.0 Unported license, as described at http://creativecommons.org/licenses/by-nc-sa/3.0/).