## Abstract

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], [Ca^{2+}], 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 (V_{m}), 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 (V_{L}) analysis (Cha et al., 2009), we could quantify contributions of individual membrane currents to the slow changes in V_{m} 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 V_{L} 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: V_{m}, seven gating variables of ion channels, three state variables of I_{NaCa}, four ionic concentrations in the cytosol ([Na^{+}]_{i}, [K^{+}]_{i}, and [Ca^{2+}]_{i}) or in ER ([Ca^{2+}]_{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 Ca^{2+}) through ion channels and exchangers should be included in calculating the time derivatives of both the membrane potential (dV_{m}/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}, [Ca^{2+}]_{i}, [K^{+}]_{i}, and [Ca^{2+}]_{ER}) and V_{m} was avoided by applying the charge conservation law:^{2+}]_{ER}/dt from the model, leaving 17 independent variables. The initial value of each variable (for example, V_{m}(0)) are presented in Table S1 in our companion paper (Cha et al., 2011).

### 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 P_{CaV} 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 V_{m} 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.

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]_{i}s) and [Ca^{2+}]_{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 I_{CaV} (f_{us}) (Fig. S2 A). When f_{us} 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 I_{NaCa} (I_{1}) was also classified as a slow variable because it has a comparable time constant to the f_{us}. [Ca^{2+}]_{i} oscillated rapidly over the time span of a single action potential (fast Ca^{2+} ripple) superimposed on a slow drift of the Ca^{2+} plateau during the burst (inset of Fig. 4 A or Fig. 2 in Cha et al., 2011). In this study, [Ca^{2+}]_{i} was treated as a slow variable because it took ∼4 s for [Ca^{2+}]_{i} to reach a new steady state when other slow variables including [Ca^{2+}]_{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 (V_{m}, d_{CaV}, U_{CaV}, r_{KDr}, q_{KDr}, m_{KCa(BK)}, h_{KCa(BK)}, E_{i_total}, and I_{2}), after fixing nine slow variables ([ATP], [MgADP], [Re], [Na^{+}]_{i}, [K^{+}]_{i}, [Ca^{2+}]_{i}, and [Ca^{2+}]_{ER}, in addition to f_{us} and I_{1}) 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 ([S_{i}]s) are fixed, or when f_{us} was also fixed together with [S_{i}]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 at http://www.jgp.org/cgi/content/full/jgp.201110612/DC1.

## 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 V_{m}, [ATP], [Ca^{2+}]_{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 V_{m} 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 V_{m} 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], [Ca^{2+}]_{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 V_{m}, 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 V_{m} 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 I_{KATP} conductance continuously decreases, and the resulting membrane depolarization (Fig. 1 A) increases the equilibrium level of [Ca^{2+}]_{i} via a fractional activation of I_{CaV}, even though I_{SOC} is deactivating (Fig. 1 C). Elevated [Ca^{2+}]_{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^{+}/Ca^{2+} exchange, which is accelerated by increased [Ca^{2+}]_{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 (P_{CaV}) of I_{CaV}, a dominant current in generating action potentials. The bifurcation diagram showed a typical S-shaped curve of EP in the P_{CaV}–V_{m} 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 P_{CaV} (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 P_{CaV} 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 P_{CaV} 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.

In the P_{CaV}–V_{m} plane in Fig. 2 A, the membrane excitability at the standard P_{CaV} 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, EP_{1}, EP_{2}, and EP_{3} on the I-V curve, correspond to the intersections of the EP curve with the black vertical lines at P_{CaV} of 48.9 pA mM^{−1} in Fig. 2 A. The resting V_{m} is at EP_{1}, and action potentials, if triggered by an appropriate stimulus, oscillate around EP_{3}. 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 (V_{m} vs. P_{CaV}) were obtained using the values of [S]_{i}s, [Ca^{2+}]_{ER}, f_{us}, and I_{1} 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 P_{CaV} (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 [Ca^{2+}]_{i} fluctuated rapidly during the burst, as shown in the inset in Fig. 4 A. Because of these rapid changes in [Ca^{2+}]_{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 [Ca^{2+}]_{i}. The mode changes from D to C or from C to B occurred slightly earlier at the local minimums of the Ca^{2+} transient than at the maximums.

Fig. 4 B indicates that a sequence of bifurcation occurred during one cycle of bursting rhythm. At 0 s, V_{m} remains stable at EP_{1}. During the initial 7.6 s, the diagram shifts leftward with reference to the standard P_{CaV}, and EP_{1} and EP_{2} approach each other. At 7.6 s, EP_{1} coalesces with EP_{2} 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 V_{m} to move from EP_{1} 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, V_{m} returns to the stable EP_{1}, 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 I_{KATP}. In Fig. 5 B, the bifurcation diagram was calculated with respect to the ultraslow inactivation gate of I_{CaV} (f_{us}). The value of f_{us} also decreased toward the LP bifurcation point, indicating that f_{us} was also a key factor in the burst termination. Namely, the decrease of f_{us} reduces the inward I_{CaV} and leads to gradual hyperpolarization. The concurrent reduction of Ca^{2+} influx through I_{CaV} lowers the plateau level of [Ca^{2+}]_{i}, overcoming the opposite effect of the accumulation of Ca^{2+} in the ER. After [Ca^{2+}]_{i} reached a peak at 12 s, it facilitates the burst termination (Fig. 5 C) by decreasing inward I_{NaCa} and I_{TRPM}. The accumulation of intracellular Na^{+} also had a significant effect on the termination of the burst through I_{NaK} at 8 mM [G] (Fig. 5 D). In contrast, [K^{+}]_{i} and I_{1} (inactivation state of I_{NaCa}) 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].

## 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 [Ca^{2+}]_{i} and Ca^{2+}-dependent activation of I_{KCa} (Chay and Keizer, 1983; Sherman et al., 1988) or Ca^{2+}-dependent inactivation of I_{CaV} (Chay and Kang, 1988). Subsequently, other hypothetical mechanisms include voltage-dependent slow inactivation of I_{Ca} (Chay, 1990; Keizer and Smolen, 1991), [ATP], or [ADP] to modulate the conductance of I_{KATP} (Smolen and Keizer, 1992), or [Ca^{2+}]_{ER} to activate I_{SOC} (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 (f_{us}) via I_{CaV}, intracellular ATP metabolism via I_{KATP}, [Ca^{2+}]_{i} via I_{NaCa} and I_{TRPM}, and [Na^{+}]_{i} via I_{NaK}. These mechanisms make a positive contribution to termination of the burst against the minor and opposite effects of [K^{+}]_{i} or I_{1} (inactivation state of I_{NaCa}). 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 I_{KATP} but, at the same time, retards it by depressing I_{NaK}; [Ca^{2+}]_{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 [Ca^{2+}]_{i} as the bifurcation parameter, which was the sole slow variable in their model. The bifurcation diagram was superimposed with a Ca^{2+} nullcline in V_{m}–[Ca^{2+}]_{i} space, and thereby, the overall behavior of the system was easily tracked along the V_{m}–[Ca^{2+}]_{i} diagram guided by both the Ca^{2+} nullcline and bifurcation points. In the study of Bertram and Sherman (2004), their model was composed of three slow variables, [Ca^{2+}]_{i}, [Ca^{2+}]_{ER}, and [ADP], and bifurcation diagrams were presented with respect to [Ca^{2+}]_{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 [Ca^{2+}]_{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 V_{m} 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 Ca^{2+} sequestration (Chay and Keizer, 1983) or the Ca^{2+}-binding affinity to Ca^{2+} channels (Chay, 1993). In more complex models, a hypothetical parameter proportional to the proton motive force (Keizer and Magnus, 1989) or the conductance of I_{KATP} (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 [Ca^{2+}]_{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 (s_{1}) and the other with a much larger time constant (s_{2}). Three slow variables were assigned to [Ca^{2+}]_{i} and [ATP] or [Ca^{2+}]_{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 I_{CaV} and I_{NaCa}. 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 f_{us} might correspond to s_{1}, and [Na^{+}]_{i} or [Ca^{2+}]_{ER} to s_{2}, 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 f_{us} were predominant, whereas relatively large variations in [Na^{+}]_{i} or [Ca^{2+}]_{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 Ca^{2+}-dependent consumption of ATP (k_{ATP,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

## 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.

## Footnotes

- Abbreviations used in this paper:
- EP
- equilibrium point
- LC
- limit cycle
- LP
- limit point

- Submitted: 7 February 2011
- Accepted: 8 June 2011

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/).