The effect of unhealthy β-cells on insulin secretion in pancreatic islets

  • Yang Pu1,

    Affiliated with

    • Saangho Lee2,

      Affiliated with

      • David C Samuels2,

        Affiliated with

        • Layne T Watson1, 3 and

          Affiliated with

          • Yang Cao1Email author

            Affiliated with

            BMC Medical Genomics20136(Suppl 3):S6

            DOI: 10.1186/1755-8794-6-S3-S6

            Published: 11 November 2013



            Insulin secreted by pancreatic islet β-cells is the principal regulating hormone of glucose metabolism and plays a key role in controlling glucose level in blood. Impairment of the pancreatic islet function may cause glucose to accumulate in blood, and result in diabetes mellitus. Recent studies have shown that mitochondrial dysfunction has a strong negative effect on insulin secretion.


            In order to study the cause of dysfunction of pancreatic islets, a multiple cell model containing healthy and unhealthy cells is proposed based on an existing single cell model. A parameter that represents the function of mitochondria is modified for unhealthy cells. A 3-D hexagonal lattice structure is used to model the spatial differences among β-cells in a pancreatic islet. The β-cells in the model are connected through direct electrical connections between neighboring β-cells.


            The simulation results show that the low ratio of total mitochondrial volume over cytoplasm volume per β-cell is a main reason that causes some mitochondria to lose their function. The results also show that the overall insulin secretion will be seriously disrupted when more than 15% of the β-cells in pancreatic islets become unhealthy.


            Analysis of the model shows that the insulin secretion can be reinstated by increasing the glucokinase level. This new discovery sheds light on antidiabetic medication.


            With 25-30 percent of adults in the developed world having a high risk of diabetes, and 24 million people in the United States with diabetes, diabetes mellitus ranks as a principal cause of death [1]. The incidence of diabetes mellitus is expected to double within the coming 20 years. After a meal, healthy individuals secrete insulin into the bloodstream signaling the consumption of glucose produced from the digested food. Thus insulin has an important signaling role in the maintenance of low blood glucose levels via glucose consumption. Type II diabetes may result from disruption of this signaling process [2, 3]. Hence understanding the insulin secretion mechanism is vitally important.

            The main pathways involved in glucose-stimulated insulin secretion are shown in Figure 1. Glucose is imported into β-cells and undergoes the glycolysis process. Three primary products, pyruvate, malate, and nicotinamide adenine dinucleotide plus hydrogen (NADH), are produced in glycolysis [4] and transported into mitochondria. In mitochondria pyruvate and malate are consumed in the tricarboxylic acid (TCA) cycle to produce fumarate. Adenosine triphosphate (ATP), the most important energy resource in the cell, is produced from fumarate and NADH through oxidative phosphorylation and released to the cytoplasm [5]. The increasing ATP level in the cytoplasm blocks the ATP-sensitive K+ channel [68], resulting in a higher intracellular K + level and depolarization of the cell membrane. This results in the opening of voltage-dependent Ca2+ channels. Cytoplasmic Ca2+ is transported into the β-cell from the extracellular fluid through the voltage-gated Ca2+ channels [911]. When both the cytoplasmic ATP and Ca2+ levels rise high enough, the vesicles containing insulin fuse with the cell plasma membrane, releasing insulin from the cell. The increasing formation and export of ATP to the entire cell will raise the Ca2+ levels in the mitochondria as well. High mitochondrial Ca2+ levels eventually collapse the mitochondrial membrane potential, shutting off ATP production, possibly by opening a low-conductance state of the permeability transition pore (PTP) in the mitochondrial membrane [1214]. The rise and fall in the mitochondrial membrane potential and intramitochondrial Ca2+ concentration is the mitochondrial oscillation. The frequency of this oscillation depends on how rapidly the mitochondria can build up its membrane potential and transport Ca2+. Besides the metabolic pathway in a single β-cell, cells in each pancreatic islet are also connected through direct electrical connections between neighboring β-cells. This connection is sufficient to synchronize both electrical bursting activity and metabolic oscillations [15].
            Figure 1

            Metabolic pathway in the β -cell.

            The cells, in fact, have differing capabilities in consuming glucose, since their differing physical properties cause different rates of electron transport chain function, ATP production, Ca2+ uptake and release, and cell membrane potential. Furthermore, some β-cells do not function normally as just described. The term unhealthy is used to describe cells that cannot consume glucose normally. Unhealthy cells interact strongly with adjacent healthy cells, because of their electrical connections. The total insulin secretion is the aggregation of the insulin secreted by all pancreatic islet cells. Its dynamics may be destroyed by a small fraction of unhealthy cells even though a large majority of the cells are still healthy. Consequently, understanding how a cohort of unhealthy cells can affect the insulin secretion process in other healthy cells is very important, and can lead to a better understanding of the cause of diabetes.

            This present work proposes a mathematical model to study scenarios with unhealthy cells leading to insulin oscillation failure. The new model is based on a mathematical insulin secretion model of the oscillation and metabolic pathway of insulin secretion in a single β-cell, due to Bertram et al. [1618]. Bertram's original model is extended to model spatial and coupling effects [14] due to the electrical connections between neighboring β-cells. Mitochondrial dysfunction is believed to have a strong negative effect on insulin secretion [1921], and also mitochondrial dysfunction is believed to be a main cause of unhealthy cells, based on the recent data of [2224]. The ratio of the total mitochondrial volume to the cytoplasm volume per β-cell, a parameter corresponding to mitochondria function, is used to define unhealthy cells, as distinguished from healthy cells. The effect on the islet insulin oscillation of unhealthy cells coupled with healthy cells is then studied in detail. Multiple cell simulations demonstrate insulin oscillation malfunction when the fraction of unhealthy cells exceeds approximately 15%. To verify this result, a simplified coupling topology is used to study the effect of one unhealthy cell on neighboring healthy cells. The latter study confirmed the previous more complicated simulation.

            In order to understand the cause of insulin secretion failure resulted from unhealthy cells, an eight-cell model is built and studied in details. We discovered a critical difference of the dynamics between healthy cells and unhealthy cells. Our simulation results demonstrate that stimulating glucokinase can make unhealthy pancreatic islets function normally. Based on the discovery, a possible strategy for antidiabetic medicine is proposed. Our strategy is consistent with recent antidiabetic medicine development [2530] that identifies glucokinase as a major drug target.


            Single β-cell modeling

            The mathematical model of a single pancreatic β-cell is based on the deterministic model introduced by Bertram et al. [1618]. The model has four components: glycolytic, mitochondrial metabolism, cytoplasmic intermediate, and plasma membrane. Each component of the model is connected by variables shared between components as illustrated in Figure 2.
            Figure 2

            Model components and interconnections. The whole model is divided into four components: glycolytic, mitochondrial metabolism, cytoplasmic intermediate, and plasma membrane.

            The glycolytic component of the model is a simplified model of the glycolysis pathway [31] that converts glucose into pyruvate. This pathway generates oscillations via a critical enzyme phosphofructokinase (PFK), the rate limiting variable in glycolysis and an important control. PFK is allosterically inhibited by ATP [32]. As a result high ATP levels will cause glycolysis to slow down. PFK phosphorylates fructose-6-phosphate (F6P) to fructose-1,6-bisphosphate (FBP), which is a substrate for glyceradehyde 3-P dehydrogenase (GPDH). The GPDH reaction is in a rapid equilibrium, therefore used as the input to the mitochondrial component, replacing pyruvate dehydrogenase (PDH) in the model. The glycolysis model consists of equations for the concentrations [G6P] and [FBP]. Here [G6P] is the glucose 6-phosphate concentration, which is assumed to be in rapid equilibrium with [F6P], and the Js are fluxes.
            The model for the mitochondrial component is based on the detailed Magnus-Keizer model of mitochondrial metabolism [33]. The tricarboxylic acid cycle (TCA cycle) and oxidative phosphorylation occur in mitochondria. The TCA cycle is a series of enzyme-catalyzed chemical reactions whose main role is to convert pyruvate to NADH and reduced flavin adenine dinucleotide (FADH), and then feed them into oxidative phosphorylation, which uses energy released by the oxidation of nutrients to produce ATP [34]. The model includes an expression for Ca2+-dependent dehydrogenases of the citric acid cycle, yielding NADH. The NADH supplies electrons for the electron transport chain, which produces a membrane potential Δψ across the mitochondrial inner membrane. The flux of protons down the electrical gradient through the F1F0-ATP synthase converts ADP to ATP. Finally, Ca2+ enters the mitochondria through Ca2+ uniporters and is transported out through Na+/Ca2+ exchangers [35]. The mitochondrial compartment has variables for the NADH, ADP, and mitochondrial Ca2+ concentrations, and for Δψ The ATP concentration is calculated from the ADP concentration by assuming that the sum of the two is conserved. Differential equations for the mitochondrial variables are
            There are five variables in the cytoplasmic intermediate and plasma membrane components, which describe the plasma membrane potential V and the cytosolic ADP concentration, the fraction n of the opened delayed rectifying K+ channels and cytosolic Ca2+ concentration, as well as the ER Ca2+ concentration. The differential equations are

            where I K , I Ca , I K(Ca), and I K(ATP) are the ionic current on the membrane, C is the membrane capacitance, f c is the fraction of free Ca2+ in the cytosol, J mem is the flux of Ca2+ across the plasma membrane, J m is the flux of Ca2+ out of the mitochondria, which is scaled by the mitochondria/cytosol volume ratio κ, J er is the flux out of the ER, f er is the fraction of free Ca2+ in the ER, V c and V er are the volumes of the cytosolic and ER compartments, respectively, n is the fraction of the open delayed rectifying K+ channels, τ is the relaxation time for the open and close reactions of the delayed rectifying K+ channels to reach equilibrium, and the steady state function for n is n (V) = (1 + exp(−(V + 16)/5))−1. For more details on these terms in the differential equations of the model, refer to Bertram et al.

            Multiple β-cells modeling

            There are about 1, 000 β-cells in each pancreatic islet. They interact with neighboring β-cells through direct electrical connections. These interactions are modeled with a change in (3a):
            where I j is the current for cell j, gc is the electrical coupling conductance, and N (j) is the index set for all neighbor cells of cell j. The neighbors surrounding a cell j are detected by considering the position of a cell j within a 3-D hexagonal lattice (different from the 3-D von Neumann neighborhoods of cellular automata theory). As shown in Figure 3, each internal cell is connected with six cells in the same layer, three cells in the upper layer and three other cells in the lower layer. For example, cell 23 in Figure 3 is connected with cells 18, 19, 22, 24, 27, 28 in its own layer, and connected with cells 7, 11, 12 in the upper layer, and cells 34, 35, 39 in the lower layer. Face, edge, and corner cells have between three and nine neighbors, depending on their location. Table 1 is used to calculate the coordinates of the neighbors of cell (a, b, c) in the 3-D hexagonal lattice.
            Figure 3

            Three layers of β -cells in 3-D hexagonal lattice.

            Table 1

            The coordinates of the neighbors of the cell with coordinates (a, b, c)


            Coordinates (x, y, z)


            (a − 1, b, c)

            Upper Layer

            (a − 1, b + 1, c + 1)


            (a − 1, b, c + 1)


            (a, b, c − 1)


            (a, b, c + 1)

            Same Layer

            (a, b − 1, c − 1)


            (a, b − 1, c)


            (a, b + 1, c)


            (a, b + 1, c + 1)


            (a + 1, b − 1, c − 1)

            Lower Layer

            (a + 1, b, c − 1)


            (a + 1, b, c)

            Besides the spatial differences, each cell has its own parameter set. Among all the parameters of the single cell model, κ, the ratio of mitochondria volume to cytosol volume, represents the function of mitochondria. The single cell study has shown that the oscillation behavior of a cell is very sensitive to κ, and thus κ is used to describe the cell differences. The value of κ is different for each cell in the 3-D structure, while all other cell parameters are set the same.

            Results and discussion

            The differential equations for the model, implemented in Fortran 90, were solved using LSODE in ODEPACK [36]. Generated data were plotted using Matlab. The simulations are performed on a Mac OS 2.4 GHz Intel(R) Core 2 Duo CPU with 4GB memory.

            Single β-cell simulation

            The simulation model for multiple cells has at its core the single cell model. Since the total amount of ATP and ADP (in the single cell model) is conserved, only ADP appears in the differential equation

            where J hyd is the hydrolysis rate of ATP to ADP, κ (κ = 0.0733 in the standard model [37]) is the ratio of the total functional mitochondrial volume to the cytoplasm volume, and J ANT is the flux through the adenine nucleotide translocator, which exchanges ADP and ATP molecules between the cytoplasm and the mitochondria. With κ = 0.0733, both the cell membrane potential and insulin secretion show periodic bursts of rapid oscillation, as illustrated in Figure 1 of Additional file 1. The slowest component of the compound bursting is due to oscillatory glycolysis, reflected by an oscillatory FBP concentration. The oscillatory glycolysis causes slow oscillations in ADP, which superimpose with the faster ADP oscillations driven by Ca2+. The multimodal ADP rhythm leads to oscillations in the conductance of the ATP-dependent potassium channel, which drives the burst episodes of the membrane potential V, which then results in compound bursting of intracellular calcium, leading to pulsatile insulin secretion.

            In truth, though, κ varies between cells; specifically, cells with some mitochondria weakened due to aging will effectively have a κ smaller than the default value. It is therefore paramount to study the sensitivity of the cell membrane potential V and insulin secretion I patterns to the parameter κ. For this work, Bertram's single cell model is modified by varying the parameter κ to represent different levels of functional mitochondria within different cells. In terms of changes in the patterns of V, Ca2+, and I, the bursting behavior of the model is classified into four categories as a function of the volume of functioning mitochondria in the β-cell, in order of occurrence: burst formation, periodic burst, burst loss, and decoupling.

            The bursting pattern starts to form at around 85% of the regular mitochondrial volume, which is when the total volume of the mitochondria is about 6% of the cytoplasmic volume. If the volume of functioning mitochondria is lower than this limit, insufficient ATP is produced to allow the insulin burst. In the burst formation category (Additional file 1: Figure 2) each burst consists of a small number of spikes in plasma membrane voltage, cytoplasmic calcium ion concentration, and the insulin secretion rate. If the mitochondrial volume fraction is increased beyond the value used in the standard model, the behavior of the burst begins to change in pattern, entering the burst loss category (Additional file 1: Figure 3). In the burst loss category every other burst decreases in duration until it is lost completely, while the remaining bursts lengthen their duration. The net effect is a decrease in the fraction of time spent in each burst.

            When the volume fraction of the mitochondria rises over about 0.095 the oscillation patterns undergo a radical alteration (Additional file 1: Figure 4). The frequency of the cell membrane voltage oscillation within the bursts becomes very large. While the membrane voltage continues to undergo bursts of rapid oscillations separated by periods of no oscillation, this no longer drives the calcium ion concentration, which then normally drives the insulin secretion. Thus this behavior will be called the decoupling category. While insulin release does undergo a long and slow periodic variation, the rapid spikes that characterize the normal burst no longer occur.
            Figure 4

            Bursting fraction for the cell membrane potential V , the calcium ion concentration Ca 2+ , and the insulin secretion rate I . The X-axis represents the ratio of functional mitochondrial volume and the cytoplasmic volume. The Y-axis represents the bursting fraction.

            In order to quantify the development and eventual loss of the bursting behavior, define the "bursting fraction" as the portion of the total time taken up by the bursts, where a burst duration is defined as the time from the first to the last peak in a burst. Bertram's single cell original model has a bursting fraction of approximately 0.55 (meaning that a burst in the model variables' oscillations is occurring for 55% of the time). When the cell mitochondrial volume fraction is changed, this bursting fraction changes markedly (cf. Figure 4), and the patterns of the calcium ion concentration Ca2+, the cell membrane potential V, and the insulin secretion rate I all alter noticeably. Figure 4 shows the ranges of the mitochondrial volume ratios corresponding to the four classifications: burst formation, periodic bursts, burst loss, and decoupling. Since the classifications are qualitative, the divisions between the classifications are blurred and the ranges are drawn overlapping. The membrane voltage, calcium ion concentration, and insulin release all have nearly identical bursting fractions for the first three categories (formation, periodic busts, and loss). In the decoupling category the bursting fractions in the calcium ion concentration and insulin release rapidly fall to zero, while the voltage bursting fraction holds relatively constant. The normal bursting behavior of the single β-cell model occurs only within a narrow range of mitochondrial volumes, from roughly 7% to 8% of the cellular volume.

            The reason for the absence of the bursting function for low mitochondrial volume is straightforward. With a lower volume of functioning mitochondria the cytoplasmic ATP levels are lower and are insufficient to trigger the remaining stages of the signaling pathway for insulin release (Figure 1). The change in behavior at higher mitochondrial volumes, where cytoplasmic ATP levels are slightly higher, has a more complex mechanism. This model of the insulin secretion is built around a fast central oscillator [10, 38, 39] consisting of the plasma membrane potential V and the fraction n of open delayed rectifying potassium ion channels, denoted by K+ in Figure 1. As the symbol implies, the opening and closing of these potassium channels is controlled by the cytoplasmic ATP concentration. During the bursts these two variables execute a fast oscillation (Figure 5), which drives the other oscillations of this model system. Between bursts, this oscillation stops and the V and n values move over to the "tail of the Q" in Figure 4, until the next burst begins. The period of this central oscillator driving the whole model is controlled by the ATP dependence on the opening of the potassium channels. Higher ATP levels lead to a faster oscillation of V and n, which (in this model) becomes too fast for the dynamics of the cytoplasmic calcium ion concentration to follow, decoupling the rest of the system from the central oscillator. Therefore, altering the mitochondrial/cytosol volume ratio in the model, which corresponds to the amount of functional mitochondria in the cell, alters the amount of ATP production and through this the level and pattern of insulin secretion. In order to reflect cell differences, the value of κ is different for each cell in the 3-D structure, with all other cell parameters being the same.
            Figure 5

            Phase plane showing relationship between plasma membrane potential V an activation variable for the delayed rectifying K + channels n . The X-axis represents the fraction of open delayed rectifying channels. The Y-axis represents the membrane potential with a unit of mV .

            The insulin bursts are totally absent when the volume fraction of the mitochondria is either less than 0.06 or greater than 0.095, which suggests an alternative definition for "healthy cell" and "unhealthy cell". Define a healthy cell as having a mitochondria volume fraction κ between 0.06 and 0.095, and a cell with κ out of this range as an unhealthy cell. It is plausible that some mitochondria lose their function, effectively resulting in a smaller κ. The ensuing numerical experiments use κ = 0.05 for unhealthy cells and the standard value κ = 0.0733 for healthy cells.

            Multiple β-cells simulation

            Although there are about 1,000 cells in each pancreatic islet, for multiple β-cells simulations, consider first the case with 125 cells coupled spatially in a 3-D hexagonal lattice. The justification for using 125 rather than 1000 cells is a pragmatic one--the CPU time for simulating 1,000 cells is rather long. In a 1,000 cell heterogeneous model, each single cell model has ten variables, yielding a 10,000-dimensional ODE, which required 26 hours to solve till time t = 2 × 106. Furthermore, the oscillation patterns observed from 125 cells are qualitatively very similar to those observed from 1,000 cells.

            If there are no unhealthy cells at all, the membrane potentials of all the cells synchronize after the coupling is turned on at time 400,000 milliseconds (ms) by changing gc from 0 to 150. This simulation of 125 cells without unhealthy cells is shown in Figures 1 (membrane potential) and 2 (total insulin secretion) in Additional file 2. The curves in the membrane potential plot are out of phase at time t = 0, but soon after 400,000 ms, these curves coalesce (see Additional file 2: Figure 1). Because the insulin levels of some cells are high while those of other cells are low, the total insulin is relatively flat before synchronization. Immediately after the coupling is turned on, the total insulin secretion shows bursts and its value rises to a hundred times that of a single cell, because there are more than a hundred cells synchronized and releasing insulin in phase.

            Focus on total insulin secretion to see how unhealthy cells, through the 3-D coupling in the hexagon structure, affect the total insulin secretion. To save computational time the coupling is turned on at the beginning of the simulations, t = 0. Figure 3 in Additional file 2 shows the resulting total insulin behavior with 10% of the cells being unhealthy spread uniformly in the 3-D hexagonal structure. The total insulin, as in the case of 100% healthy cells, shows periodic oscillations and maintains a reasonable level. When the percent of cells being unhealthy increases to 15%, the oscillations of total insulin still look normal (see Additional file 2: Figure 4), but now some bursts have fewer spikes. As the percentage of cells being unhealthy increases to 20% and 30% from 10% and 15% of cells being unhealthy, the spikes within each burst become much less numerous (shown in Additional file 2: Figures 5 and 6). These bursts are also much more irregular, and even more significantly, totally disappear after 2.25 × 106 ms (Additional file 2: Figure 6). In summary, the cohort of unhealthy cells dominates the global behavior, resulting in a level of total insulin too low to maintain proper pancreatic islet function. The conclusion is that if there is more than approximately 15% of cells being unhealthy, the function of the pancreatic islet will be severely affected.

            To verify the conclusions from simulations with 125 cells, simulations were also performed for the model with 1,000 cells. This 1000 cell model cannot be run to as long a time as the 125 cell model, but observe that when, for a certain time, insulin secretion cannot generate enough spikes, the pancreatic islet can be considered to be malfunctioning. Hence the simulation monitors the numbers of spikes in the last five bursts, and if the mean number of spikes is below three, deems that the overall system is malfunctioning and halts. The simulations for 1,000 cells with 30%, 20%, 15% of cells being unhealthy, all considered malfunctioning systems, are shown in Figures 7, 8, and 9 of Additional file 2 respectively. A simulation for 10% of cells being unhealthy found no malfunction in an extremely long time (comparable to the time for the 125 cell model runs). In summary, the 125 cell and 1000 cell simulations yield similar conclusions: If the percentage of cells being unhealthy is larger than approximately 15%, the system will malfunction; at 15% of cells being unhealthy, the system still functions but is close to malfunctioning; below 10% of cells being unhealthy, the system can definitely function very well.

            Another ratio, the number of links between an unhealthy cell and a healthy cell to the total number of links between all cells, is also calculated for each simulation; the value of this ratio is quite close to the number of unhealthy cells to total number of cells ratio, see Table 2. Therefore, for further analysis, using 125 cells and the ratio of number of unhealthy cells to total number of cells should suffice.
            Table 2

            Comparison of Two Different Ratios

            Ratio of Links

            Ratio of Cells









            Simplified multiple β-cells simulation

            In order to understand further how unhealthy cells affect the healthy cells, consider a simplified case with only one unhealthy cell in the 3-D structure, as shown in Figure 6. The red dot represents the unhealthy cell, the blue dots healthy cells. The unhealthy cell is connected with all the healthy cells, while the healthy cells themselves are connected through three rings corresponding to three layers in the 3-D structure. The experiment starts with two cells, one healthy and one unhealthy. Then the number of healthy cells is increased one by one. If the total number of cells is smaller than 13, the corresponding subgraph of Figure 6 is used. Figure 1 in Additional file 3 shows the total insulin of one unhealthy cell and one healthy cell (coupling starts at 3 × 105 ms). According to the "malfunction criterion" (the average number of spikes in the last five bursts is smaller than three), the system malfunctions. Figure 2 in Additional file 3 shows the total insulin for one unhealthy cell and three healthy cells. The simulations show that when the number of healthy cells reaches three (marginally) or four, the total insulin shows functioning bursts. This suggests that one unhealthy cell needs at least four healthy cells to make the whole system function well based on the topology in Figure 6. To study the case of two unhealthy cells with exactly the same topology and connected to the same number of healthy cells, the topology in Figure 7 is used. The number of healthy cells needed to rescue two unhealthy cells is eight (marginally) or nine, which is consistent with the result of one unhealthy cell in Figure 6. Therefore one unhealthy cell needs at least four healthy cells to make the overall system function well, i.e., the system probably can tolerate at most 20% of cells being unhealthy. This result matches with the previous simulations with 125 and 1,000 cells in the 3-D topology.
            Figure 6

            Topology for one unhealthy cell. A simplified case with only one unhealthy cell in the 3-D structure. The red dot (center dot) represents the unhealthy cell, the blue dots are healthy cells.

            Figure 7

            Topology for two unhealthy cells. The simplified topology for two adjacent unhealthy cells in two layers.

            The cause of insulin secretion failure

            To study possible reasons for the failure of insulin secretion caused by unhealthy cells, consider a 2 × 2 × 2 model of eight cells, in which three unhealthy cells are uniformly distributed. Figures in Additional file 4 show the behaviors of all the variables in this model for each of the eight cells in each figure. The insulin secretion starts to fail at around time 2 × 105 ms as shown in Figure 8. Except for the plot of [G6P] (see Additional file 4: Figure 10), it is hard to distinguish healthy cells from unhealthy cells in all the plots, such as that of membrane potential in Figure 5 of Additional file 4. In Figure 10 of Additional file 4 eight curves are partitioned into two groups. The upper group consists of healthy cells, while the lower one consists of unhealthy cells. These two groups of curves do not intersect with each other after the failure, but they are mixed together before the failure. This observation suggests that G6P might be closely related to the insulin secretion failure. In order to test this hypothesis, the two groups of curves were manually reset to see whether insulin secretion can be reinstated: if the insulin secretion fails for 5 × 105 ms, the values of [G6P] in all cells are reset to 279, the initial value of [G6P] used in the simulation. The oscillations of insulin secretion are resumed (see Additional file 5: Figures 1 and 2). The black vertical lines in these figures are the times when the values of [G6P] are equalized. The results show that the divergence of the two groups of [G6P] leads to the disappearance of insulin secretion oscillations.
            Figure 8

            Insulin secretion failure of the eight cells system. There are three unhealthy cells uniformly distributed in the 2 X 2 X 2 topology

            Since an unhealthy cell has a smaller mitochondria/cytosol volume ratio than that of a healthy cell, unhealthy cells usually have lower ATP (from oxidative phosphorylation) levels. ATP has a negative feedback to phosphofructokinase (PFK) in glycolysis. When the ATP level is lower, the PFK reaction becomes faster, which directly consumes G6P faster leading to lower G6P levels in unhealthy cells than in healthy cells. As shown in Figure 10 of Additional file 4 the levels of G6P in healthy cells are higher than those in unhealthy cells. At the same time, because of the low production rate of ATP in unhealthy cells, the ratio of ATP to ADP gets lower. ATP-sensitive K+ channels in the plasma membrane are activated by ADP and inactivated by ATP, so the ratio of these nucleotides determines the fraction of open ATP-sensitive K+ channels. When the ATP/ADP ratio is low there is an increase in the number of open ATP-sensitive K+ channels, which results in the difficulty of membrane depolarization. Voltage-dependent Ca2+ channels are blocked. Since the insulin is secreted when Ca2+ exceeds a certain level, the blocked Ca2+ channels will reduce insulin secretion. Therefore, raising the levels of G6P in unhealthy cells back to normal will bring back normal insulin secretion. This analysis suggests the hypothesis that increasing the value of any substance ahead of the ATP synthesis in the pathway, such as the substances FBP and NADH in mitochondria, reinstates the oscillations. Manually resetting the values of either [FBP] or [NADH] in all the cells to the same initial value confirms the conjecture: not only can G6P restart the insulin oscillations, but also FBP and NADHm can restart the oscillations.

            Since G6P and FBP are both substances in the glycolysis pathway, if the rate of glycolysis were faster, the oscillations of insulin might be resumed as well. In order to test this conjecture, the glucokinase level is raised fourfold after the insulin failure is detected. Glucokinase is the enzyme that phosphorylates glucose to glucose-6-phosphate (G6P). Figures 3 and 4 in Additional file 5 show that the oscillations of insulin are resumed when the glucokinase level increased. Two oscillation periods are shown in Figure 5 of Additional file 5 (from the earlier 125 cell model): The longer period (slow oscillations) is caused by the glycolytic oscillations, while Ca2+ feedback is responsible for the fast oscillations. Only the short period remains in Figure 6 of Additional file 5. This observation suggests that if the longer period could be reestablished in Figure 6 of Additional file 5 the insulin oscillations might be reinstated. Since the glycolysis pathway is responsible for the slow oscillations, that pathway may need some stimulus. The simulation results in Figures 4 and 5 of Additional file 5 demonstrate that by increasing the glucokinase level in the glycolysis pathway, the insulin oscillations can be resumed after the failure.

            Figure 7 in Additional file 5 shows the simulation result when the three unhealthy cells are placed together in the 2 × 2 × 2 grid. With a high ratio of unhealthy cells to total cells (more than .20), the insulin secretion fails. If the glucokinase level in the unhealthy cells is raised fourfold and then kept unchanged during the simulation, the total insulin secretion becomes normal (see Additional file 5: Figure 8). From Figure 9 in Additional file 5 all the unhealthy cells look like the healthy cells. Even though those unhealthy cells are still unhealthy, if their glucokinase levels are high enough, the overall behavior of the system will be dominated by the healthy cells. If there is no healthy cell, only unhealthy cells with high glucokinase levels, the system can't be repaired, see Figure 9. Therefore, the unhealthy cells are rendered harmless by high glucokinase level and other healthy cells in the network.
            Figure 9

            Total insulin secretion of the eight unhealthy cells system with the stimulation on glucokinase of the cells. All the cells are placed in the 2 × 2 × 2 topology


            In conclusion, the insulin secretion of pancreatic islets usually will be functionally destroyed when there are more than 20% of cells being unhealthy among all cells. The more unhealthy cells there are, the more irregular insulin secretion will be. Increasing the level of glucokinase can make the pancreatic islet function normally when there is a high fraction of unhealthy cells by increasing the glucose absorption of the glycolysis pathway. This has implications for the clinical treatment of type II diabetes. Currently there are three classes of medications used to treat type II diabetes. The first treatment is to increase the amount of insulin secreted by the pancreas by inhibiting the opened delayed rectifying K+ channels. The second treatment is to increase the sensitivity of target organs to insulin. The third treatment is to decrease the rate at which glucose is absorbed from the gastrointestinal tact, which is a method to reduce the glucose uptake from food. It appears that all the available treatments are insufficient to stem the tide. Therefore, new treatments are currently under investigation including the development of therapeutic agents with novel action mechanisms. Recently, researchers have identified glucokinase as an outstanding drug target for developing antidiabetic medicines [2530]. Assuming the mathematical models are valid, the simulation results demonstrate that stimulating glucokinase can make unhealthy pancreatic islets function normally, consistent with new antidiabetic medicines. Such multiple cell models are good candidates for guiding the development of the next generation of antidiabetic medicines.



            adenosine diphosphate


            adenosine triphosphate


            endoplasmic reticulum




            flavin adenine dinucleotide




            glucose 6-phosphate


            glyceradehyde 3-P dehydrogenase


            nicotinamide adenine dinucleotide plus hydrogen


            pyruvate dehydrogenase




            permeability transition pore


            tricarboxylic acid.



            This work is based on the conference paper "The effect of unhealthy beta-cells in synchronized insulin secretion", which appeared in the 2012 IEEE International Conference on Bioinformatics and Biomedicine (BIBM 2012) [14]. This work was partially supported by the National Science Foundation under awards CCF-0726763 and CCF-0953590, and the National Institutes of Health under award GM078989.


            The publication costs for this article were funded by the corresponding author with the above sources of support.

            This article has been published as part of BMC Medical Genomics Volume 6 Supplement 3, 2013: Selected articles from the IEEE International Conference on Bioinformatics and Biomedicine 2012: Medical Genomics. The full contents of the supplement are available online at http://​www.​biomedcentral.​com/​bmcmedgenomics/​supplements/​6/​S3.

            Authors’ Affiliations

            Department of Computer Science, Virginia Polytechnic Institute and State University
            Center for Human Genetics Research, Vanderbilt University
            Department of Mathematics, Virginia Polytechnic Institute and State University


            1. Ford ES, Giles WH, Mokdad AH: Increasing prevalence of the metabolic syndrome among us adults. Diabetes Care 2004, 27:2444–2449.PubMedView Article
            2. Cerasi E: Aetiology of type II diabetes. In Insulin: Molecular Biology to Pathology. Oxford University Press 1992, 347–392.
            3. Nesher R: Modeling phasic insulin release: immediate and time-dependent effects of glucose. Diabetes 2002, 51:S53-S59.PubMedView Article
            4. Nielsen K, Sorensen PG, Hynne F, Buss HG: Sustained oscillations in glycolysis: An experimental and theoretical study of chaotic and complex periodic behavior and of quenching of simple oscillations. Biophys Chem 1998, 72:49–62.PubMedView Article
            5. Wiederkehr A, Wolheim CB: Minireview: implication of mitochondria in insulin secretion and action. Endocrinology 2006, 147:2643–2649.PubMedView Article
            6. Rorsman P: The pancreatic beta-cell as a fuel sensor: an electrophysiologist's viewpoint. Diabetologia 1997, 40:487–495.PubMedView Article
            7. Lang J: Molecular mechanisms and regulation of insulin exocytosis as a paradigm of endocrine secretion. European journal of biochemistry 1999, FEBS 259: 3–17.View Article
            8. Tsuura Y, Ishida H, Okamoto Y, Kato S, Sakamoto K, Horie H, Okada Y, Seino Y: Glucose sensitivity of ATP-sensitive K+ channels is impaired in beta-cells of the GK rat. A new genetic model of NIDDM. Diabetes 1993, 42:1446–1453.PubMedView Article
            9. Wolheim CB, Maechler P: Beta cell mitochondria and insulin secretion: messenger role of nucleotides and metabolites. Diabetes 2002, 51:37–42.View Article
            10. Bergsten P, Grapengiesser E, Gylfe E, Tengholm A, Gylfe E, Tengholm A, Hellman B: Synchronous oscillations of cytoplasmic Ca 2+ and insulin release in glucose-stimulated pancreatic islets. J Biol Chem 1994, 269:8749–8753.PubMed
            11. Jung SK, Kauri LM, Qian WJ, Kennedy RT: Correlated oscillations in glucose consumption, oxygen consumption, and intracellular free Ca 2+ in single islets of Langerhans. J Biol Chem 2000, 275:6642–6650.PubMedView Article
            12. Kindmark H, Kohler M, et al.: Glucose induced oscillation in cytoplasmic free Ca 2+ concentration precede oscillations in mitochondrial membrane potential in the pancreatic beta cell. J Biol Chem 2001, 276:34530–34536.PubMedView Article
            13. Yang Pu, Saangho Lee, David Samuels, Layne Watson, Yang Cao: Hybrid modeling and simulation of insulin secretion pathway in pancreatic islets. In Proceeding of BIBE 2010. IEEE Computer Society; 2010:156–161.
            14. Yang Pu, Saangho Lee, David Samuels, Layne Watson, Yang Cao: The effect of unhealthy ß-cells in synchronized insulin secretion. Bioinformatics and Biomedicine (BIBM), 2012 IEEE International Conference on: 4–7 October 2012 2012, 1–4.
            15. Zimliki CL, Mears D, Sherman A: Three roads to islet bursting: emergent oscillations in coupled phantom bursters. Biophys J 2004, 87:193–206.PubMedView Article
            16. Richard Bertram, Leslie SatinS, et al.: Interaction of glycolysis and mitochondrial respiration in metabolic oscillations of pancreatic islets. Biophysical Journal 2007, 92:1544–1555.View Article
            17. Richard Bertram, Rudy ArceoC II: A mathematical study of the differential effects of two SERCA isoforms on Ca 2+ oscillations in pancreatic islets. Bulletin of Mathematical Biology 2008, 70:1251–1271.View Article
            18. Bertram R, Gram Pedersen M, Luciani DS, Sherman A: A simplified model for mitochondrial ATP production. Journal of Theoretical Biology 2006, 243:575–586.PubMedView Article
            19. Soejima A, Inoue K, Takai D, Kaneko M, Ishihara H, Oka Y, Hayashi JI: Mitochondrial DNA Is Required for Regulation of Glucose-stimulated Insulin Secretion in a Mouse Pancreatic Beta Cell Line, MIN6. J Biol Chem 1996, 271:26194–26199.PubMedView Article
            20. Malaisse WJ, Hutton JC, Kawazu S, Herchuelz A, Valverde I, Sener A: The stimulus-secretion coupling of glucose-induced insulin release. XXXV. The links between metabolic and cationic events. Diabetologia 1979, 16:331–341.PubMedView Article
            21. Kennedy E, Maechler P, Wollheim CB: Effects of depletion of mitochondrial DNA in metabolism secretion coupling in INS-1 cells. Diabetes 1998, 47:374–380.PubMedView Article
            22. Maechler P, Wollheim CB: Mitochondrial function in normal and diabetic beta-cell. Nature 2001, 414:807–812.PubMedView Article
            23. Wollheim CB: Beta-cell mitochondria in the regulation of insulin secretion: a new culprit in Type II diabetes. Diabetologia 2000, 43:265–277.PubMedView Article
            24. Mizukami H, Wada R, Koyama M, et al.: Augmented beta cell loss and mitochondrial abnormalities in sucrose-fed GK rats. Virchows Arch 2008, 452:383–392.PubMedView Article
            25. Grimsby J, Sarabu R, Corbett WL, et al.: Allosteric activators of glucokinase: potential role in diabetes therapy. Science 2003, 301:370–373.PubMedView Article
            26. Sarabu R, Grimsby J: Targetting glucokinase activation for the treatment of type 2 diabetes: a status review. Current Opinion in Drug Discovery and Development 2005, 8:631–637.PubMed
            27. Guertin KR, Grimsby J: Small molecule glucokinase activators as glucose lowering agents: a new paradigm for diabetes therapy. Current Medicinal Chemistry 2006, 13:1839–1843.PubMedView Article
            28. Matschinsky FM, Magnuson MA, Zelent B, et al.: The network of glucokinase-expressing cells in glucose homeostasis and the potential of glucokinase activators for diabetes therapy. Diabetes 2006, 55:1–12.PubMedView Article
            29. Pal M: Recent advances in glucokinase activators for the treatment of type 2 diabetes. Drug Discovery Today 2009, 14:784–792.PubMedView Article
            30. Matschinsky Fm, Zelent B, Doliba N, Li CH, et al.: Glucokinase activators for diabetes therapy. Diabetes Care 2011, 34:S236-S243.PubMedView Article
            31. Nan Jiang, Roger Cox, John HancockM: A kinetic core model of the glucose-stimulated insulin secretion network of pancreatic beta cells. Mamm Genome 2007, 18:508–520.View Article
            32. Uyeda K, Furuya E, Luby L: The effect of natural and synthetic D-fructose 2,6-biphosphate on the regulatory kinetic-properities of liver and muscle phosphofructokinases. J Biol Chem 1981, 256:8394–8399.PubMed
            33. Magnus G, Keizer J: Model of beta-cell mitochondrial calcium handling and electrical activity. I. Cytoplasmic variables. Am J Physiol 1998, 274:C1158-C1173.PubMed
            34. Cortassa S, Aon MA, et al.: An integrated model of cardiac mitochondrial energy metabolism and calcium dynamics. Biophys J 2003,84(4):2734–2755.PubMedView Article
            35. Wu F, Yang F, Vinnakota KC, Beard DA: Computer modeling of mitochondrial tricarboxylic acid cycle, oxidative phosphorylation, metabolite transport, and electrophysiology. J Biol Chem 2007, 282:24525–24537.PubMedView Article
            36. Alan Hindmarsh: ODEPACK, A Systematized Collection of ODE Solvers. Scientific Computing 1983, 55:C64.
            37. Dean P: Ultrastructural Morphometry of the Pancreatic Beta-Cell. Diabetologia 1973, 9:115–119.PubMedView Article
            38. Liu YJ, Tengholm A, et al.: Origin of slow and fast oscillations of Ca2 + in mouse pancreatic islets. J Physiol 1998, 508:471–481.PubMedView Article
            39. Valdeolmillos M, Santos R, Contreras D, Soria B, Rosario L: Glucose-induced oscillations of intracellular Ca 2+ concentration resembling bursting electrical-activity in single mouse islets of langerhans. FEBS LETTERS 1989, 259:19–23.PubMedView Article


            © Pu et al; licensee BioMed Central Ltd. 2013

            This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://​creativecommons.​org/​licenses/​by/​2.​0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The Creative Commons Public Domain Dedication waiver (http://​creativecommons.​org/​publicdomain/​zero/​1.​0/​) applies to the data made available in this article, unless otherwise stated.