|Home | About | Journals | Submit | Contact Us | Français|
Pancreatic β-cells release insulin in response to increased glucose levels. Compared to isolated β-cells, β-cells embedded within the islets of Langerhans network exhibit a coordinated and greater insulin secretion response to glucose. This coordinated activity is considered to rely on gap-junctions. We investigated the β-cell electrophysiology and the calcium dynamics in islets in response to glucose gradients. While at constant glucose the network of β-cells fires in a correlated fashion, a glucose gradient induces a sharp division into an active and an inactive part. We hypothesized that this sharp transition is mediated by the specific properties of the gap-junctions. We used a mathematical model of the β-cell electrophysiology in islets to discuss possible origins of this sharp transition in electrical activity. In silico, gap-junctions were required for such a transition. However, the small width of transition was only found when a stochastic variability of the expression of key transmembrane proteins, such as the ATP-dependent potassium channel, was included. The agreement with experimental data was further improved by assuming a delay of gap-junction currents, which points to a role of spatial constraints in the β-cell. This result clearly demonstrates the power of mathematical modeling in disentangling causal relationships in complex systems.
The electrophysiology of β-cells is most relevant for understanding insulin secretion (Rorsman and Renstrom, 2003). Exocytosis of insulin loaded granules is dependent on intracellular calcium signaling. Glucose transport and metabolism induces the synthesis of adenosine triphosphate (ATP) in the mitochondria of β-cells. This increases the cellular ATP to adenosine diphosphate (ADP) ratio and inhibits ATP-sensitive potassium (K-ATP) channels. The reduced outflow of potassium depolarizes the β-cell and opens voltage-dependent calcium channels. Oscillations between calcium inwards and potassium outwards currents are induced leading to bursts of the membrane potential. A slow increase in the intracellular calcium is triggered that at some level interrupts the oscillation. The increased calcium level is at the same time the triggering event for an increased exocytosis frequency of insulin carrying granules (Gilon et al., 2002).
In experiments, a prominent difference in the stability of bursts in excitable β-cells was observed depending on whether the β-cells are isolated or in their natural environment within the islet of Langerhans (Zhang et al., 2003). While isolated β-cells exhibit unstable bursting, β-cells in an islet context show strongly correlated activity. This stability is attributed to the connection of β-cells with gap-junctions (Ravier et al., 2005; Rocheleau et al., 2006). However, the detailed role of gap-junctions for β-cell electrophysiology and calcium dynamics is not well understood and a matter of debate.
We are convinced that a mathematical model of β-cells coupled by gap-junctions will help to disentangle their impact on β-cell electrical activity in islets. So far, two β-cells have been coupled by gap-junction (Jo et al., 2005) based on a model of the β-cell electrophysiology (Sherman et al., 1988; Sherman, 1996). These are based on a subset of transmembrane currents and well reproduce the behavior of β-cells. A possible role of noise was considered to induce stabilization of electrical bursting (de Vries and Sherman, 2000; Jo et al., 2005). Bigger maps of coupled cells were considered in view of bursting (de Vries, 2001) and more recently with particular emphasis on pancreatic β-cells (Giugliano et al., 2000; Benninger et al., 2008).
A previously developed quantitative model of the β-cell electrophysiology (Meyer-Hermann, 2007) is fully based on measured single transmembrane properties and is, thus, considered as a good starting point for a map of coupled β-cells. Here, this model was combined with the gap-junction model in (Benninger et al., 2008) using β-cells arranged in a regular grid with defined neighborhood relations. The electrical activity throughout the modeled islet was considered for different applied glucose levels and for stable glucose gradients along the islet.
The modeling results were compared to experimental results (Rocheleau et al., 2004). It was found in silico that the gap-junctions induce a correlation of β-cell activity in the islet. The stability of the correlation can be perturbed with stochastic variations in some quantities (as found before in (Jo et al., 2005). A delay of the adaptation of the gap-junction current to changed intracellular properties turns out to induce nonlinear (chaotic) effects in the islet. When applying a glucose gradient along the islet a sharp division of the islet into an active and an inactive part is found in the experiment. This transition could be reproduced in the model by assuming a stochastic variation in K-ATP channel expression levels.
The goal of this work was to develop a mathematical model of the β-cell electrophysiology that allows an analysis of the impact of gap-junctions between β-cells coupled in islets of Langerhans on electrical activity and calcium dynamics. The resulting mathematical model and the analysis of gap-junctions under different conditions are presented. In particular, we applied different levels of glucose to the whole islet as well as glucose gradients along the islet.
We investigated the activity of β-cells in an islet that are connected by gap-junctions. A two-dimensional network of 27 β-cells arranged in a block of 3×9 cells was used. Each β-cell is assumed to exhibit exactly the same properties as in the model for isolated cells (Meyer-Hermann, 2007). The extension of this model to gap-junction coupling is described as follows.
Rate equations are used to describe the dynamics of intracellular ion concentrations. External ion concentrations are treated as a thermodynamic bath and are assumed to be constant. The ions treated explicitly are sodium (N), potassium (K), and calcium (C). The dynamics depend on the conductance of various membrane proteins, where each corresponds to a single term in the following differential equations.
Here, Ix denote the membrane currents, where x specifies the type of membrane protein, which are listed in Table Table1.1. Negative currents are defined to bring positive ions into the cell. The Ix are defined as single channel/protein conductance and corresponding currents can be derived from experiments [see Meyer-Hermann (2007)]. n in the term for the gap-junction stands for the number of connected neighbor cells. ρx denote the corresponding surface densities of the membrane proteins. ξ is the surface to volume factor that translates the transmembrane current into a space-averaged intracellular change in the ion concentration. The factors αNCX,Na,K account for the stoichiometry of exchangers. zCa=2 is the valence of calcium ions. The leakage currents JK,Na,Ca are defined to guarantee a stable resting state of the cell. In the resting state all other currents are balanced by the leakage currents.
Calmodulin is a calcium buffer. It is assumed to bind calcium faster than the typical time scale of the ion dynamics and the membrane potential. Thus, we use the rapid buffer approximation. This gives rise to the additional factor (1+xb) in the equation for calcium Eq. 1.
where b0 is the total buffer concentration and Kb is the corresponding dissociation constant.
The dynamics of these ions build up the membrane potential V.
where the same terms as in Eq. 1 appear. Note that all currents are defined to be electrical currents. Thus, the contributions enter without additional factor in Eq. 3 but have additional factors in Eq. 1 in the case that a single electrical charge does not correspond to a single ion passing the membrane. Cm is the membrane capacitance and replaces the surface to volume ratio ξ used in the case of ion concentrations. Typical values are Cm=0.01 pF/μm2 (Chay, 1997) or 0.0009 pF/μm2 [for neurons, Gentet et al. (2000)].
All currents IXX,YY (where the subscript is one of the acronyms in Table Table1)1) are derived from single transmembrane protein characteristics. This is described in detail in (Meyer-Hermann, 2007). In the following the additional gap-junction current is specified.
Model of gap-junction currents: The term entering all ionic differential equations of cell i and which was previously used in (Benninger et al., 2008) is in more general terms
where the sum goes over all n connected neighbors j and the currents are different for different ions. Assuming the density of gap-junctions to be identical for all cell-cell connections we have ρgap,ijρgap.
β-cells are connected with the gap-junction protein connexin-36 (Moreno et al., 2004). The whole cell conductance mediated by these junctions is in the order of 500–1,000 pS (Zhang et al., 2008). A single gap-junction has a conductance of 5.1±0.7 pS (Moreno et al., 2004). One cell-cell connection is assumed to be concentrated within 1/6 of the cell surface. Assuming 30 gap-junctions per contact infers
Indeed, then for six neighbors we have
using the physiological single gap-junction conductance of 5.1 pS, which is in agreement with the value 950±96 pS (Zhang et al., 2008).
Note that if a cell has only three, four, or five neighbors (such at the border) the resulting total conductance is reduced according to the real number of neighbors. This assumes that the connections between two cells are always of the same type. In 2D-simulations, the density is still calculated according to a 3D-cell with six connections because the density is related to the surface of the 3D-cell. However, in 2D the maximum number of connections is 2D=4 and the total possible gap-current is, thus, reduced by a factor of 2/3 (see the factor 1/6 in Eq. 4, which in principle, could also be replaced by 1/(2D) in a different strategy with D the dimension of the simulation space).
Connexin-36 has a very weak voltage gating (Moreno et al., 2004) so that a constant conductance can be safely assumed. The channels between β-cells are much bigger than the ions, and allow molecules of up to 1 kDa or 2 kDa to pass. Nevertheless, the resulting current depends on the potential difference in both connected cells and the ion-specific reversal potential.
The asymptotic potential difference is
with the latter term representing the Nernst-potential
x=K,Na,Ca and R=8.315 J/(K mol) is the Rydberg (molar) gas constant. The temperature T has to be adopted to the experimental situation or the in vivo case. xi,j are the intracellular ionic concentrations of the cells i and j, respectively. Note that the currents are symmetric by construction, i.e.,
In this approach it is assumed that the point of ionic currents between the extracellular medium and the cytosol is separated from the location of the gap-junction, such that changes in the ionic intracellular concentrations become effective on the level of the gap-junction with a time delay . The delay can be estimated from the diffusion constant and the distance that has to be overcome together with typical diffusion constants. The distance is assumed to be in the order of the cell radius around 6 μm. The diffusion constant is in the range of 223 μm2/s (Allbritton et al., 1992). This gives rise to
where the diffusion constant was simply assumed equal for all ions. One might argue whether the delay also applies to the membrane potential Vi−Vj or to the Nernst-term only. One might also argue that the gap-junctions themselves give rise to a delay in activity, which would be in analogy to other transmembrane proteins. In addition, on the scale of 161 ms the local changes in ions are already equilibrated in the sense that sodium and potassium exchange is achieved on the scale of milliseconds. Thus, a big part of the concentration changes might never reach the gap-junctions. Taken together these thoughts make the exact value of the delay rather uncertain and we, therefore, also evaluated delays in the range of a few to 500 ms in the simulations.
In the resting state Vi=Vj and xi=xj holds, implying a vanishing Nernst-potential and gap-junction current. This also implies that the leakage currents, which are determined to guarantee the resting state of the cells, are not altered by gap-junctions.
In the limit τ→0s the current equation for gap-junctions simplifies to
The glucose level is set to 10 mM throughout the islet. The cell activity in this in silico islet with β-cells connected by gap-junctions turns out to be correlated and fully synchronized. The synchronous behavior of all β-cells is the same as the one published for an isolated in silico β-cell at 10 mM stimulation (see Meyer-Hermann (2007); data not shown). The synchronization is well supported by experiment, where more than 90% of electrically active cells in an islet oscillated synchronously (Benninger et al., 2008).
Next the glucose concentration was varied from 2 mM to 12 mM at both ends of the long direction of the modeled islet. In experiment synchronous calcium oscillations are seen in the side of the islet stimulated with 12 mM glucose. There is a sharp transition of less than two cells width, where oscillations are halted. The oscillations in calcium are halted at a point in the islet that is stimulated with an average of 7.2±0.6 mM glucose, determined by NAD(P)H imaging (Rocheleau et al., 2004).
We asked whether and under which conditions this transition between an active and an inactive part of the islet is found in silico. Applying the same glucose gradient, the β-cells exhibited a degree of excitation that continuously followed the glucose level (Fig. (Fig.11 right panel). Thus, the sharp transition found in experiment was not reproduced. There are three striking observations: first, bursting at 10 mM glucose is not comparable to isolated in silico cells at 10 mM glucose. The β-cells at high glucose exhibit a pronounced spike instead of a burst (Fig. (Fig.11 left panel). Second, β-cells exhibit a remaining activity at glucose levels, where isolated in silico cells do not spike at all (Fig. (Fig.11 left panel, blue and green curve). Third, the activity of the β-cells is correlated throughout the islet, i.e., all cells fire in phase (Fig. (Fig.11).
A role of noise was postulated to impact on β-cells in islets (de Vries and Sherman, 2000; Jo et al., 2005). We asked whether it is possible to get desychronized firing of coupled β-cells and the observed transition in a glucose gradient with a deterministic model. The used model for the β-cell electrophysiology (Meyer-Hermann, 2007) does not resolve the intracellular space. However, the position of calcium entry points in the plasma membrane and the gap-junctions are not necessarily the same. The locally increased calcium level might need to diffuse to the gap-junction. The distance between both channels will induce a delay of the gap-junction current with respect to the time of calcium entry. Also, the quick clearance of calcium will reduce the calcium concentration that reaches the gap-junction.
A proper model for these effects would be based on a model that resolves the intracellular space and includes intracellular organelles. We leave this for future investigations and estimated here, whether these effects might change the behavior of connected β-cells. To this end, a phenomeneological time delay of the current adaptation at the gap-junction to the space-averaged ion concentrations of the cell is used [see Eq. 9]. The delay was estimated by the typical size of a β-cell to be given by Eq. 13. In contrast to the in silico islet without delay of gap-junctions, now bursting of all β-cells in the islet is induced (Fig. (Fig.22 upper left panel). The burst are correlated throughout the islet (see Fig. Fig.22 upper left panel). The overall level of reached calcium is reduced (see Fig. Fig.22 upper right panel). However, the transition between an active and an inactive part of the islet with a glucose gradient is not altered, in particular, it did not become sharper (Fig. (Fig.22 upper right panel).
When increasing the delay to 500 ms, the β-cell activity is only weakly correlated (Fig. (Fig.22 lower panels). The nonlinearity of the model becomes a visible property of the model. The transition between an active and an inactive part of the islet is now more pronounced even though not reaching the desired sharp transition. [See also Supplementary Material Movie; the movie shows 3×9 β-cells in a two-dimensional plane connected by gap-junctions. Each cell-cell connection contains 30 gap-junctions with a conductance of 5.1 pS each. The gap-junction currents adapt to the cytosolic state of the cell with a delay of 500 ms. The color codes the membrane potential from low (blue) to high (red).] Note that the activity of β-cells is interrupted in the whole islet with irregular intervals (Fig. (Fig.22 lower left panel)
The membrane potentials of individual cells in an islet with glucose gradient are shown in Fig. Fig.33 for a delay of 50 ms (upper panels) and 500 ms (lower panels). It can be seen that the longer delay induces less activity in the membrane potential for the same glucose level. This confirms that the transition between an active and inactive part of the islet is sharper for longer gap-junction delays.
How is the islet behaving when the gap-junction current is increased? It is expected that strong gap-junction currents will induce a stronger coupling of all β-cells in the islet. Indeed, a 25-fold gap-junction current induces a stronger correlation of β-cell spiking and almost identical calcium waves among all connected β-cells (see Fig. Fig.4).4). No bursts but single spikes are observed. Interestingly, also, the β-cells that see a glucose level of only 2 mM now spike and exhibit calcium waves. Thus, gap-junctions can induce synchronization of islet β-cells down into the part of the islet that senses glucose levels below the β-cell activity threshold.
In a gap-junction knock-out experiment, the gap-junction density is set to zero for all β-cells in the islet. This electrically isolates all β-cells and allows us to consider their behavior as isolated cells for different glucose levels. The corresponding spiking or bursting behavior is depicted in Fig. Fig.5.5. This might be compared to an islet sensing the same glucose gradient but with β-cells connected by gap-junctions in Fig. Fig.33 (using a 50 ms delayed gap-junction current). It can be seen that connected β-cells exhibit activity at glucose levels for which the isolated β-cells in this knock-out experiment do not show any activity. This, demonstrates that gap-junctions have the capability to transport the activity of β-cells to regions of the islet at which subcritical glucose levels are encountered. This confirms the relevance of gap-junctions for the activity of β-cells embedded in islets.
The gap-junction density, so far, was assumed to be identical for all β-cell connections in the islet. We asked about the impact of individual gap-junction densities for different connections. The random choice of the individual gap-junction expression level was done in a simple way.
For all neighbor pairs ij, where r[0,2a] is a random number and a determines the range of variability, which is set to a=0.5 here. In average the expression level is identical to the value ρgap that was previously used for all β-cells.
Using this range of expression levels between 50% and 150% noise is introduced to the system but the general behavior of the system does not change (data not shown). Thus, a stochastic variation in gap-junction currents does not induce major effects onto β-cell activity in islets.
In order to get a sharper transition between the active and the inactive part of the islet, we investigated a stochastic expression of K-ATP-channels. The random expression level of the K-ATP-density is chosen according to the law in Eq. 15 with a=0.5. At spatially constant glucose levels of 10 mM this variation exhibited division of the islet into two distinguished groups of cells. The two groups differed in the degree of calcium activity (see Fig. Fig.66 upper panels).
This observation depends on the random number generator seed value and the grouping can be more or less pronounced compared to the example shown in Fig. Fig.6.6. In particular, the depicted example exhibits high calcium activity in the part of the islet that corresponds to the low glucose part in the case of application of a glucose gradient. Thus, a division of the islet into an active and inactive part could not be expected or might even be suppressed for this random number seed. Surprisingly, this is not the case. Even though the grouping is not very pronounced it is still stronger than without K-ATP stochasticity (data not shown).
Using a different random number generator seed value that does not induce a strong grouping at constant glucose, no bias is set by the random K-ATP-expression profile. Then, the activity division of the islet into two parts became even more pronounced (see Fig. Fig.66 lower panels). As in experiment (Rocheleau et al., 2004), the transition from the active to the inactive part happens within a width of two cells at glucose levels between 5.75 mM and 7.0 mM. The classification is more pronounced for a longer time delay of gap-junction adaptation (compare Fig. Fig.6,6, right panels).
Within the active part of the islet the β-cells now exhibit the behavior expected from experiment (see Fig. Fig.66 lower panels).
The inactive β-cells show a remaining low activity, which is not reflected in the membrane potential by bursts but by single spikes. This is shown for the two β-cells representing the transition from the active to the inactive part in Fig. Fig.66 (lower right panel) at 5.75 mM and 7.0 mM glucose (see Fig. Fig.77).
We have investigated β-cell coupling in islets based on experimental data (Rocheleau et al., 2004) by combining a quantitative mathematical model of β-cell electrophysiology (Meyer-Hermann, 2007) with a recent mathematical model for gap-junctions (Benninger et al., 2008). The interesting experimental observation of the division of islets into an active and an inactive part after application of a glucose gradient was analyzed. Experiments show a balance in the coupling mechanisms sufficient to synchronize stimulated region but not enough to allow propagations into unstimulated region. This result is reproduced by the newly developed model for β-cell electrophysiology including gap-junctions. It is found in silico that coupled β-cells fire in a synchronous way throughout the islet if a spatially constant glucose level is applied. A glucose gradient induces a graded behavior of coupled β-cells. A division of the islet into an active and an inactive part of the islet is not observed and only slightly promoted by a delay of gap-junction adaptation to changed intracellular membrane potential and ion concentrations.
Previous theoretical work for different cells (de Vries, 2001) and for single β-cells (Jo et al., 2005) suggest a role of stochastic expression levels on β-cell coupling. This was further supported by the observation that the number of gap-junctions per cell-cell connection varies by at least 50%. In the present model stochastic variation in the gap-junction expression level did not exhibit a strong effect on the transition between active and inactive parts of the islet in a glucose gradient.
Thus, the mathematical model can only qualitatively describe the behavior found in experiment. The model of coupled β-cells confirms a role of gap-junctions in coordinating electrical activity. As the agreement between theory and experiment was not fully quantitative we used the model to find hints for an explanation of the sharp transition between the active and inactive parts of islets in a glucose gradient. One might consider to add missing currents, or that the way we handle the coupling architecture or mechanism of action is not fully correct. A ~20% decrease in gap-junction conductance was observed at 2 mM compared to 11 mM (unpublished data). Also, a weak voltage gating of the gap-junction (Moreno et al., 2004) may impact on the results.
The most important hint coming from the present model is the reproduction of the experimentally observed sharp transition from active to inactive β-cells in response to stochastic variation in the expression levels of K-ATP channels. A 50% variation leads to a transition with characteristics that are in agreement with experiment. Thus, again noise is found to strongly impact on β-cell activity in islets even though not the noise in gap-junction expression.
The situation is even more improved when assuming a delay of the gap-junctions currents. The time delay of the gap-junction Eq. 9 is formulated as a delay of the adaptation of the membrane potential. This effectively reflects a delay caused by limited ion diffusion because the membrane potential is solely built up by ions. The membrane potential itself will certainly equilibrate faster in the cell than the electrochemical gradient is set-up at the location of the gap-junction. Thus, one should be aware that this is just a phenomenological description of a spatial effect. A concise description of this phenomenon requires a spatial modeling approach including calcium buffering in organelles and explicit diffusion of ions, which is left for the future. The qualitative finding that the time delay of the adaptation of the gap-junction current to changed cell states improves the in silico data with respect to experiment, suggests that spatial effects may play an important role for β-cell bursting in islets of Langerhans.
Experiments were performed as described in detail (Rocheleau et al., 2004). Briefly, islets were starved 1 h at 2 mM glucose, then loaded into a two channel microfluidic device. The media were delivered down two channels: one channel at 2 mM and another channel at 12 mM. NAD(P)H was imaged using two-photon microscopy with 710 nm excitation and 380–550 nm detection. Calcium was imaged through staining the islet for 1–2 h with 4 μM Fluo-4 at room temperature during starvation. Fluo-4 was imaged with 488 nm excitation and 505–550 nm detection.
Solving the set of differential equations: The set of ordinary differential equations Eqs. 1, 3 is solved using a self-developed fourth-order Runge–Kutta algorithm. In contrast to the model without gap-junctions (Meyer-Hermann, 2007), Eq. 9 attributes to each gap-junction its own dynamics. Thus, the delay equation has to be included into the implicit solution method of the set of differential equations. Technically, this is realized by extending the cell specific properties of the β-cell-class in the C++-code by a conductance state of each of its gap-junctions.
The solver is equipped with an adaptive time step control. If a time step reduction by a factor of 2 induces an above threshold deviation of the calculation in any of the time steps, the time step is correspondingly reduced. All simulations were done in a regime of time steps in which this criterion was never met (Δt=0.00001 s). All simulations were performed on a standard PC.
Robustness of simulations: Simulations were done in a two-dimensional in silico islet of 3×9 β-cells fully connected by gap-junctions. Simulations are evaluated by monitoring the electrical activity and the calcium dynamics in the central layer of the islet. The robustness of the results was tested by varying the size and the dimension of the system. On one hand a larger part of the islet with 5×9 beta-cells was used. In addition, the results were counter checked with a three-dimensional islet with 5×5×9 β-cells. It turned out that the choice of the system size and dimension did not significantly impact on the results and left the conclusions untouched. The observed robustness also excludes major effects stemming from the boundary of the modeled islet. For this reason, we have chosen to present the results using the smallest system, which is easier to be presented in the graphs. In a more quantitative modeling approach small deviations stemming from the system’s size and dimension would have to be considered.
The results derived for stochastic variation in transmembrane protein expression levels, where derived using a nonweighted distribution as in Eq. 15. Using a weighted distribution (such as Gaussian) mainly changes the variability range needed to reproduce the same result. In view of the missing experimental data about the variability of protein expression levels in pancreatic β-cells, we have used the simplest nonweighted distribution. The range a was varied. Smaller ranges still exhibit the grouping of β-cells that leads to the discussed sharp transition as observed in experiment. However, the calcium curves of different β-cells in the islet become less and less different with reduced variation range a.
Thus, the grouping of β-cells (in terms of activity) in the pancreatic islet turns out to be a robust feature of islets with variable protein expression levels. In particular, this interesting main result also persists when all ion-conducting proteins are varied with a range of a=0.5 (data not shown).
If the results exhibit an unstable behavior the simulation was repeated with slightly perturbed initial conditions with a different random number generator seed value, and with a smaller time step size. In all presented cases the characteristics of the simulation results remained unchanged while the details of the time course were sensitive to the modifications.
M.M.H. and R.K.P.B. developed the mathematical model. M.M.H. performed the simulations. M.M.H. and R.K.P.B. wrote the article. M.M.H. and R.K.P.B. thank Jonathan V. Rocheleau and David W. Piston who developed the experimental techniques and performed the experiments (Rocheleau et al., 2004) that formed the motivation and basis of the present analysis. M.M.H. thanks Michele Solimena and Jaaber Dehghani for regular and fruitful discussion. FIAS is supported by the ALTANA AG. M.M.H. is supported by the EC-NEST-project MAMOCELL. R.K.P.B. is supported by funding from the NIH under Contract Nos. K99DK85145 (to R.K.P.B.) and R01DK53434 (to D.W.P.).