Gain of function mutation in K(ATP) channels and resulting upregulation of coupling conductance are partners in crime in the impairment of Ca 2 + oscillations in pancreatic ß -cells

Gain of function mutations in the pore forming Kir6 subunits of the ATP sensitive K + channels (K(ATP) channels) of pancreatic β -cells are the major cause of neonatal diabetes in humans. In this study, we show that in insulin secreting mouse β -cell lines, gain of function mutations in Kir6.1 result in a significant connexin36 (Cx36) overexpression, which form gap junctional connections and mediate electrical coupling between β -cells within pancreatic islets. Using computational modeling, we show that upregulation in Cx36 might play a functional role in the impairment of glucose stimulated Ca 2 + oscillations in a cluster of β -cells with Kir6.1 gain of function mutations in their K(ATP) channels (GoF-K(ATP) channels). Our results show that without an increase in Cx36 expression, a gain of function mutation in Kir6.1 might not be sufficient to diminish glucose stimulated Ca 2 + oscillations in a β -cell cluster. We also show that a reduced Cx36 expression, which leads to loss of coordination in a wild-type β -cell cluster, restores coordinated Ca 2 + oscillations in a β -cell cluster with GoF-K(ATP) channels. Our results indicate that in a heterogenous β -cell cluster with GoF-K(ATP) channels, there is an inverted u-shaped nonmonotonic relation between the cluster activity and Cx36 expression. These results show that in a neonatal diabetic β -cell model, gain of function mutations in the Kir6.1 cause Cx36 overexpression, which aggravates the impairment of glucose stimulated Ca 2 + oscillations.


Introduction
In rodents, dogs, and humans, plasma glucose levels are maintained within a narrow range by the feedback mechanism between glucose and insulin.When plasma glucose levels are elevated, pancreatic β-cells secrete insulin, which triggers glucose uptake in liver, muscle and adipose tissue cells, and reduces plasma glucose levels [1].Within the pancreas, β-cells are located in endocrine cell clusters called islets of Langerhans.In humans, each islet contains nearly one thousand β-cells [2].In β-cells, insulin secretion is triggered by increased plasma glucose levels, where glucose is transported into the β-cells through GLUT-2 transporters and metabolized to produce ATP.Increased cytosolic ATP blocks ATP sensitive K + channels (K(ATP) channels), which depolarizes the cell membrane.Membrane depolarization increases Ca 2+ influx through voltage gated Ca 2+ channels, and the resulting increase in cytosolic Ca 2+ leads to insulin secretion [3].In mice and healthy humans, insulin secretion is oscillatory, with periods ranging from 5 to 10 min [4,5].Studies show that oscillations increase the effectiveness of the insulin signal on target tissues and play an important role in reducing plasma glucose levels [6,7].Furthermore, in pre-diabetic patients, insulin oscillations are lost long before the emergence of overt diabetes [8], which further emphasizes the importance of oscillatory insulin release from β-cells.
β-cells within the same pancreatic islet are heterogeneous in their molecular, functional, and metabolic properties [9][10][11].In glucose stimulated dispersed β-cells, heterogeneity results in different glucose responsiveness [12] and bursting patterns [13][14][15].In contrast, intact pancreatic islets exhibit robust and coordinated glucose responsiveness and Ca 2+ oscillations that result in oscillatory insulin release [13].At stimulatory levels of glucose, oscillatory insulin secretion results from bursting electrical activity of single β-cells and intra-and inter-islet synchronization [16].The mechanism underlying inter-islet synchronization is not yet determined, but intra-islet synchronization results from the electrical coupling between β-cells mediated by gap junctional coupling through connexin36 proteins (Cx36) [17].Gap junctions are protein structures that form channels between adjacent β-cells and provide direct cell-to-cell communication through transfer of ions and signaling molecules [18].Therefore, in pancreatic islets, gap junctions play a vital role in coordinating the oscillatory release of insulin in response to elevated glucose levels.
In β-cells, K(ATP) channels act as molecular sensors of ATP and play a critical role in glucose stimulated insulin secretion [19,20].Each K (ATP) channel has four Kir6 subunits and four sulfonylurea receptor subunits (Sur1) [21,22].Studies show that Kir6 subunits constitute the primary domain at which ATP reacts with and inhibits the channel activity [23].Gain of function mutations in Kir6.2 subunits of pancreatic β-cells are associated with neonatal diabetes with increased K(ATP) channel activity and reduced glucose stimulated insulin secretion in humans [22,24].There is also evidence suggesting that polymorphisms in the Kir6.2 gene are associated with type 2 diabetes [25].The gain of function mutations in Kir6.2 subunits of K(ATP) channels have also been used to develop diabetic islet models to investigate the role that gap junctional coupling between β-cells plays in insulin secretion [26,27].Previous studies show that Cx36 knockout impairs insulin secretion in both basal and glucose stimulated states in wild-type mouse pancreatic islets [28,29] and insulin secreting cell lines [30].However, in transgenic mouse islets with Ki6.2 gain of function mutation, Cx36 knockout increases glucose stimulated insulin release [31].
In this study, we show that in the NIT-1 β-cell line, gain of function mutations in Kir6.1 result in increased Cx36 expression.To interpret the ramification of this increased coupling strength that results from the Kir6.1 gain of function mutation, we use computational modeling to explore the combined effects of increased K(ATP) conductance and increased coupling conductance between β-cells.We investigate the association between the gap junctional coupling strength and coordinated Ca 2+ oscillations in a heterogeneous β-cell cluster that carries the Kir6.1 gain of function mutations.We analyze the system dynamics at the single-cell and network levels, where we analyze the way in which the dynamical properties of interacting single cells give rise to the emergence of a non-monotonic response to the coupling strength at the cluster level.

Site directed mutagenesis
Kcnj8 (Kir6.1)Mouse Tagged ORF Clone (Origene, Cat No: MR206756) and pCMV6-Entry Mammalian Expression Control Vector (Origene, Cat No: PS100001) were purchased from OriGene.Kcnj8 [G343D] or Kcnj8 [G343D, Q53R] gain of function mutations were generated using QuikChange II XL Site-Directed Mutagenesis Kit (Agilent) according to the manufacturer's protocol.Site-directed mutagenesis in plasmids was verified using Sanger sequencing to confirm the presence of mutations.Mutagenesis primer sequences listed in Supplementary Table S1 were designed using the QuikChange Primer Design Tool.

Generation of NIT-1 cells with stable overexpression of wild type or gain of function mutant Kcnj8 (Kir6.1)
To generate an in vitro diabetes model, NIT-1 cells were transfected with either pCMV6-Entry Mammalian Expression Control (Origene, Cat No: PS100001) vector or one of Kcnj8 (Kir6.1)Mouse Tagged ORF Clone (Origene, Cat No: MR206756; [WT]), Kcnj8 [G343D] or Kcnj8 [G343D, Q53R] vectors using Lipofectamine 3000 (Thermo Fisher Scientific, Invitrogen) reagent according to the manufacturer's guidelines to develop stable cell lines.NIT-1 cells were then selected in complete culture medium supplemented with Geneticin (G418) at a concentration of 250 μg/mL for at least 2-3 weeks to obtain G418-resistant cells.Then, stable overexpression of Kcnj8 in NIT-1 cells was confirmed at mRNA and protein levels.

Generation of NIT-1 cells with stable inhibition of Gjd2 (Cx36)
Gjd2 is the mouse orthologue of Cx36 in humans.To investigate the effects of Gjd2 knockdown on the control and Kcnj8 overexpressing cells, NIT-1 cells were transfected with either shRNAs targeting Gjd2 or control shRNA vector purchased from OriGene (Cat No: TR500814).Briefly, NIT-1 cells were transfected with vectors using Lipofectamine 3000 (Thermo Fisher Scientific, Invitrogen) reagent according to the manufacturer's protocol.Then, cells were selected in complete media containing Puromycin (2 μg/mL) in addition to Geneticin, which is used to maintain stable overexpression of Kcnj8, to enrich the cells with stable inhibition of Gjd2.Suppression of Gjd2 was verified at mRNA and protein levels.The most effective shRNA construct was used to carry out the further experiments.

RNA extraction, complementary DNA (cDNA) synthesis and quantitative Real-Time PCR (qRT-PCR)
Total RNAs were isolated from NIT-1 cells with varying levels of Kcnj8 or Gjd2 expressions following the EcoPURE Total RNA Isolation Kit (Ecotech Biotechnology, Erzurum, Turkey) standard protocol, and the concentration and purity of the RNA samples were determined using Epoch Microplate Spectrophotometer (BioTek Instruments, Winooski, VT, USA).Complementary DNA (cDNA) synthesis from total RNA samples was carried out with "High-Capacity cDNA Reverse Transcription Kit" (Thermo Fisher Scientific) according to the manufacturer's recommendations.qRT-PCR reactions were prepared using 5x HOT FIREPol EvaGreen qPCR Mix Plus (Solis Bio-Dyne) following the manufacturer's protocols.Primer sequences used in qRT-PCR experiments are listed in Supplementary Table S2.β-actin was used for standardization of gene expressions.qRT-PCR was carried out in a Qiagen Rotor Gene Q real-time thermal cycler (QIAGEN) using standard parameters.All qRT-PCR reactions were performed in triplicates and relative expression levels were calculated using 2 − ΔΔCT method [32].

Fluorometric imaging intracellular Ca 2+
Intracellular free calcium ion concentration ([Ca 2+ ] i ) was measured using Fluo-4 Direct™ Calcium Assay Kit (Invitrogen, Thermo Fisher Scientific) following the manufacturer's instructions.Initially, µ-Dish 35 mm, high Glass Bottom ( ˙Ibidi GmbH, Gräfelfing, Germany) were precoated with 0.1 μg/mL poly-L-lysine as described before [33].NIT-1 cells (3×10 5 ) were seeded in a 35 mm glass-bottom dish.After 72 h of incubation, NIT-1 cells were washed once with phosphate-buffered saline and incubated in DMEM low glucose medium (Ecotech Biotechnology, Erzurum, Turkey) for 1 h.Subsequently, cells were incubated in Fluo-4 Direct calcium reagent solution with a probenecid concentration of 5 mM (Invitrogen) for 30 min 5 % CO 2 at 37 • C and 30 min at room temperature in the dark.Then, cells were stimulated with 20 mM glucose [34] and imaging initiated using LSM 710 Confocal Laser Scanning Microscope (Zeiss, Oberkochen, Germany).Excitation of Fluo-4 Direct calcium reagent solution was provided by the 488-nm line of an argon laser.Images were acquired at 5 sec intervals.Data were analyzed by Leica LAS AF Lite software.

Computational model of a single β-cell
The electrical activity and Ca 2+ dynamics of an individual β-cell are mathematically formulated using a modified version of the model proposed by [35].This 4-variable model is simpler than the more highly developed Integrated Oscillator Model [36] making it more amenable to network simulations and fast/slow analysis.Values of the model parameters and their explanations are given in Supplementary Table S3.The model describes the rate of change of the membrane potential of an individual β-cell as a function of ionic currents as follows: Here, C m represents the electrical capacitance of the cell membrane, I Ca is the voltage-sensitive Ca 2+ current, I K(ATP) is the ATP sensitive K + current, I K is the delayed-rectifying K + current, and I S is a slowlyactivating K + current that was formulated to be analogous to a Ca 2+ activated K + current [37].The ionic currents are as follows: Here, g K(ATP) , g Ca , g K and g s are parameters for the maximum conductances of the associated ion channels, and the values are proportional to the expression density of the associated ion channel in the cell membrane.O K(ATP) is the fraction of open K(ATP) channels.m ∞ is the activation function of voltage sensitive Ca 2+ channels and defined as follows: where the parameter V m sets the half-activation voltage and θ m is a shape parameter.The dynamic variables n and s are the activation variables for the delayed rectifier K + current and slow K + current, respectively.Finally, V K and V Ca are Nernst potentials for K + and Ca 2+ ions, respectively.The dynamic activation variables in the model are defined by the following differential equations: Here τ n and τ s are time constants.n ∞ and s ∞ are steady-state activation functions and they are expressed with the following sigmoid functions of the membrane potential: M. An et al.
Finally, the equation that governs the rate of change of intracellular free Ca 2+ concentration ([Ca 2+ ] i ) is as follows: Here, [Ca 2+ ] i is the intracellular free Ca 2+ concentration, f is the ratio of free-to-bound cytosolic Ca 2+ , and k Ca is the pumping rate of Ca 2+ ions out of the cell per unit time.α is a conversion parameter from Ca 2+ current to Ca 2+ flux.
Model β-cells produce electrical bursting with periods of ~3.5 min (Fig. 1A, black), driven by the sawtooth oscillations in the slow variable s (Fig. 1A, red).Each burst of spikes is accompanied by an increase in the cytosolic Ca 2+ concentration (Fig. 1B) that reflects Ca 2+ influx through voltage gated Ca 2+ channels.

Modeling a heterogeneous β-cell cluster with electrical coupling
The computational model describes a cluster of coupled β-cells with a 5×5×5 cubical arrangement (Fig. 2).In the past, it has been suggested that in pancreatic islets, a small subset of highly connected β-cells (hubs) [38] or highly active β-cells (leaders) [39] acted as pacemakers and they were necessary for coordinated oscillations Ca 2+ .However, computational modeling revealed that such pacemaker subsets were not necessary for coordinated oscillations, where pancreatic islets behaved like a nonlinear average of its constituent cells [40].Therefore, we use a cubical structure for its simplicity in describing electrical coupling between neighboring cells, with a regular nearest neighbor distribution.Cells in the interior of the cube are coupled to 6 neighbors.Cells on a face of the cube are coupled to 5 neighbors, except cells on an edge which are coupled to 4 neighbors, or cells on a corner that are coupled to 3 neighbors.In the model, the electrical coupling between two adjacent cells is defined as the product of the potential difference between two cells and the coupling strength [41].The coupling strength between the model cells would reflect the expression level of the connexin proteins.
The coupling current between cell i and a neighbor cell j is: Where, V i and V j are the membrane potentials of the i th and j th cells, respectively, and g c is the coupling strength between the cells.The total coupling current between cell i and all its neighboring cells is then: where n i is the set of all the cells adjacent to cell i.Finally, this current is added to the eq.( 1) for the rate of change of the membrane potential of the i th cell: In the model, heterogeneity in electrical excitability of the β-cells is introduced by randomly distributing K(ATP) channel densities.For each cell, g K(ATP) is randomly chosen from a normal distribution: where μ is the mean and σ is the standard deviation.

Model integration and computer simulations
The model parameter values and their physiological explanations are provided in the Supplemental Table S3.The system of differential equations was numerically solved in the Python programming environment using the LSODA algorithm, which is freely available within the integration module of the scientific computation package for Python (scipy.integrate1.7.3)[42].The numerical continuation analysis was carried out on the XPP/AUTO dynamical systems analysis and integration toolbox [43], and the resulting data files were exported to Python to generate figures.

Statistical analysis
The data is represented as mean ± standard deviation.Statistical significance is tested using students t-test and significance levels are reported as p-values.

Gain of function mutation in Kir6.1 results in Cx36 upregulation
To develop NIT-1 cells overexpressing the Kir6.1 subunit (Kcnj8) with gain of function mutations, we utilized site directed mutagenesis to produce a G343D mutation alone and in combination with the mutation Q53R on a wild type Kcnj8 clone as described in Methods.After confirmation of directional mutagenesis using Sanger sequencing (data not shown), we transfected NIT-1 cells with either vectors carrying wild type Kcnj8 (Kcnj8 [WT]) or vectors with Kcnj8 gain of function mutations (Kcnj8 [G343D] or Kcnj8 [G343D, Q53R]).Overexpression of Kcnj8 was confirmed at both mRNA and protein levels (Fig. 3A and B).To understand the effect of Kcnj8 overexpression on the connexin36 (Cx36) expression, we measured Gjd2 levels in cells overexpressing Kcnj8 with and without gain of function mutations.Our results showed that cells with WT or mutant Kcnj8 overexpression exhibited a significant Gjd2 upregulation (Fig. 3C and D), where the most prominent Gjd2 upregulation was found in cells carrying both gain of function mutations of Kcnj8 (Fig. 3C and D).To further verify Gjd2 upregulation in cells overexpressing Kcnj8 with gain of function mutations, we used confocal microscopy to capture the Kcnj8 (Fig. 3E, green) and Gjd2 levels (Fig. 3E, red) in mutant and control cells.Our results showed a  S3, with g K(ATP) =150 pS.significantly higher Gjd2 colocalization with mutant Kcnj8 (Fig. 3E, Kcnj8 [G343D, Q53R]) compared to the controls (Fig. 3E, control), further verifying the positive association between the expression of Kcnj8 [G343D, Q53R] and Gjd2.

Cx36 overexpression aggravates the impact of Kir6.1 gain of function mutations
Our experimental data show that gain of function mutations in K (ATP) channels result in Cx36 overexpression in β-cells.Since Cx36 proteins constitute the gap junctional connections between β-cells, we used computational modeling to understand the role of increased gap junctional conductance on the impairment of Ca 2+ oscillations in a heterogenous β-cell cluster with gain of function mutations in their K (ATP) channels (GoF β-cells).In the model, heterogeneity within a β-cell cluster is introduced by randomly sampling the maximal K(ATP) channel conductance (g K(ATP) ) from a random normal distribution with N (μ,σ) for each cell (Fig. 4A).To simulate a wild-type-like healthy β-cell cluster, we used a lower mean and variance for the g K(ATP) distribution (Fig. 4A, green), whereas we used a rightward shift in the g K(ATP) distribution with a higher variance to represent a GoF β-cell cluster (Fig. 4A, red).Depending on the value of g K(ATP) , the single model β-cells exhibited three different types of behavior.β-cells with low g K(ATP) are labelled as overactive cells (Fig. 4A, yellow area), which exhibit continuous spiking and have high cytosolic Ca 2+ concentration (Fig. 4B).β-cells with moderate g K(ATP) are labelled as intrinsicallybursting cells (Fig. 4A, green area), since in isolation, they exhibit bursting with periods ranging between 2-5 min at stimulatory levels of glucose (Fig. 4C).Finally, β-cells with high g K(ATP) are labelled as intrinsically-silent cells (Fig. 4A, red area), since they remain silent at stimulatory levels of glucose in isolation (Fig. 4D).
To test whether Cx36 expression and the resulting increase in coupling conductance play a functional role in the impairment of [Ca 2+ ] i oscillations in β-cells with gain of function mutations in their K (ATP) channels, we performed model simulations for heterogenous β-cell clusters with varying gap junctional coupling strength.In the model, gap junctional conductance is determined by parameter g c .Therefore, we used higher values of g c to account for Cx36 overexpression.The cluster of model wild type β-cells (Fig. 4A, green distribution) produced coordinated Ca 2+ oscillations with moderate g c (45 pS).This is shown using the population average of [Ca 2+ ] i in Fig. 5A.In the GoF β-cell cluster (Fig. 4A, red distribution), combined with an increase of coupling conductance to 75 pS, the cells were completely silent (Fig. 5B).When g c was reduced to the baseline level (45 pS), [Ca 2+ ] i oscillations were restored in the GoF β-cells cluster (Fig. 5C).This result shows that the observed increase in Cx36 expression as a result of the Kir6.1 gain of function mutation might have a functional role in the impairment of glucose stimulated Ca 2+ oscillations.
To understand the way excitability is related to the gain of function mutations in K(ATP) channels and Cx36 overexpression, we investigated the bursting threshold in a heterogenous β-cell cluster as a function of the mean g K(ATP) value in the cell cluster and the coupling strength between cells.Fig. 5D shows the bursting region for a heterogenous β-cell cluster on the mean g K(ATP) -g c plane.In the blue region, all or a subset of the cells exhibit bursting electrical activity and [Ca 2+ ] i oscillations, whereas in the red area, the entire cluster remains silent (Fig. 5D), which indicates no response to stimulatory glucose.The curve that separates the blue and red regions serves as the mean g K(ATP) -g c bursting threshold.In Fig. 5D, the angle between green and red arrows (θ) shows the relative increase in g c in response to increased mean g K(ATP) .A heterogenous β-cell cluster (Fig. 5D, green circle) can remain within the bursting region despite increased mean g K(ATP) as long as θ is small.However, for large θ, an increase in mean g K(ATP) is accompanied by a profound increase in g c and the cluster becomes silent (Fig. 5D, red circle).This result shows that in a heterogenous β-cell cluster, increased g K(ATP) can be well tolerated to a certain level unless there is a significant increase in g c .Hence, Cx36 overexpression seen in GoF β-cells might play an important role in the impairment of glucose stimulated [Ca 2+ ] i oscillations.
To experimentally test the effect of reduced Cx36 expression on the activity of β-cells with gain of function mutations in Kir6.1, we used shRNAs targeting Gjd2 in Kcnj8 [G343D, Q53R] overexpressing NIT-1 cells (Kcnj8 [G343D, Q53R]/shRNA Gjd2) and measured cytosolic Ca 2+ concentrations using Fluorometric imaging as described in Methods.shRNA targeting reduced Gjd2 expression to the control levels in NIT-1 cells overexpressing Kcnj8 [G343D, Q53R] (Fig. 6A and B, shRNA Gjd2).Fluorometric Ca 2+ measurements were recorded from Kcnj8 [G343D, Q53R] overexpressed control (Fig. 6C and Supplemental video 1) and shRNA Gjd2 targeted NIT-1 cells (Fig. 6D and Supplemental video 2).For each group, between t = 0 − 50 sec, 6 images from Fig. 2. Heterogeneous β-cell population is modeled as a 5×5×5 cubical cell cluster.Electrical excitability of the cells is indicated by the color gradient from red to green.Cells that cannot intrinsically produce bursting are called 'silent' cells, and they are illustrated with shades of red.Cells that can intrinsically produce bursting are called 'bursting' cells, and they are shown with shades of green.Silent cells have a higher K(ATP) channel conductance (g K(ATP) ), which implies a higher K(ATP) channel expression on their plasma membrane.The gap junctional conductance between two adjacent cells is given by g c .zoomed-in regions are illustrated (Fig. 6C and D).In each zoomed-in region, normalized [Ca 2+ ] i signal intensities from indicated single cells (Fig. 6C and D, yellow and magenta markers) are shown for 6 min (Fig. 6E).The cells with Kcnj8 [G343D, Q53R] overexpression exhibited no detectable Ca 2+ oscillations (Fig. 6E, yellow; Supplemental video 1).In contrast, shRNA targeted Gjd2 knockdown restored glucose stimulated Ca 2+ oscillations in Kcnj8 [G343D, Q53R] cells (Fig. 6E, magenta; Supplemental video 2).

Reduced gap junctional coupling impairs electrical activity in cluster of bursting cells
To understand the effect of changes in Cx36 expression on the electrical activity of a β-cell cluster, we investigated the average [Ca 2+ ] i patterns of each cell as a function of coupling strength, g c .Fig. 7 shows the model simulations for the average intracellular Ca 2+ concentration (Fig. 7A) of a wild type β-cell cluster with varying g c levels (75-0 pS).For these simulations, the g KATP levels of the individual β-cells are selected from the green distribution shown in Fig. 4A, which is a balanced distribution among overactive (Fig. 4B), intrinsically-bursting (Fig. 4C) and intrinsically-silent (Fig. 4D) cells.The simulations show that for high g c levels (75-60 pS), the cell cluster generates coordinated [Ca 2+ ] i oscillations with periods around 5 min (Fig. 7A, 0-60 min).For moderate g c levels (45-30 pS) the cell cluster maintains bursting [Ca 2+ ] i oscillations with a reduced amplitude and period (Fig. 7A, 60-120 min).Reducing g c below 30 pS impairs the coupling between the cells and diminishes synchronous oscillations.The box plot in Fig. 7B shows the distribution and 30-min average of the intracellular [Ca 2+ ] i concentration of the  7B, g c =75-30 pS).Reducing g c to 15 pS, the average activity remains high with increased variance (Fig. 7B, purple).Setting g c to 0 completely terminates the coupling between the cells and significantly reduces the activity at both cluster (Fig. 7A, 150-180 min) and individual level (Fig. 7B, brown).These results show that in a wild-type-like model cell cluster that has a balanced g KATP distribution, reducing g c not only impairs the coordination between the cells, it also reduces the average activity.

Activity of a heterogenous cluster of intrinsically silent cells is a nonmonotonic function of the coupling strength
To understand the effect of coupling strength on a cell cluster that exhibits gain of function mutation in K(ATP) channels, we investigated the average [Ca 2+ ] i patterns of a heterogenous cell cluster with a rightshifted g KATP distribution (Fig. 4A, red).Fig. 8 shows the model simulations for the average intracellular Ca 2+ concentration of the β-cell cluster with varying g c (75-0 pS) (Fig. 8A).For high g c (75-60 pS), the cell cluster remains silent with low [Ca 2+ ] i levels (Fig. 8A, 0-60 min), except  S3 were used.S3 were used.
for a single burst when the coupling conductance is reduced to 60 pS.Reducing g c to 45-30 pS, the cell cluster generates coordinated bursting (Fig. 8A, 60-120 min) with high overall activity (Fig. 8B green and red box plots).Nevertheless, in this range, some cells remain intrinsically silent and only produce sub-threshold oscillations due to the intercellular coupling.Reducing g c to 15 pS, the cell cluster continues to oscillate, but with reduced average amplitude (Fig. 8A, 120-150 min).Finally setting g c to 0 completely terminates the coupling between the cells and most oscillations are lost (Fig. 8A, 150-180 min).The box plot in Fig. 8B shows that the mean [Ca 2+ ] i across the population first rises when the coupling conductance is decreased, and then falls when the coupling is terminated altogether.Also, the variance in average [Ca 2+ ] i across the population increases with decreasing coupling strength.Thus, unlike the case of the model wild-type cells, in the mutant cells, the average [Ca 2+ ] i across the population has an inverted u-shaped relation with the coupling strength.

Recruitment and silencing are responsible for the inverted u-shape average Ca 2+ distribution in a cluster of bursting cells
To understand the inverted u-shaped relation between g c and the activity of the model mutant cell cluster, we analyzed the evolution of the system dynamics for different coupling conductances in a simple two-cell model, which comprised an intrinsically-bursting (β b ) and an intrinsically-silent model β-cell (β s ).For the burster, g K(ATP) is set to 150 pS, whereas for the silent cell g K(ATP) is set to 200 pS.The remaining parameter values are identical for each cell.In the two-cell model, the voltage equations for β b and β s are: ) The remaining variables are as described in Methods, for each cell.
In the two-cell model, setting g c to 0 isolates the cells and β b generates bursting oscillations (Fig. 9A and B, green), whereas β s remains silent (Fig. 9A and B, orange).When g c is increased to 25 pS, β s begins to burst in phase with β b , but with a shorter active phase (Fig. 9C and D).Because of this recruitment, the average [Ca 2+ ] i level of the two-cell population increases by nearly 25 % (Fig. 9B and D, magenta; 0.08 vs 0.1 μM).When the coupling strength is doubled, to 50 pS, the effect is to stop both cells from bursting and remain silent instead (Fig. 9E and F).This results in the lowest average [Ca 2+ ] i (Fig. 9F, magenta; 0.06 μM).This behavior, of recruitment of silent cells when the coupling strength is increased from 0, and then silencing of the system when the coupling becomes too strong, is what occurs in the 125-cell network and gives the inverted u-

Two-cell dynamics explained through fast/slow analysis
To understand the dynamical mechanism underlying the recruitment and silencing effects in the two-cell network we employ a fast/slow analysis [44,45].The method separates the system into fast and slow subsystems, analyzing each separately, and stitching them together to form an understanding of the full dynamics.On the fast timescale, the slow variables can be treated as parameters of the fast subsystem, and the dynamics of the fast subsystem can be explored as those parameters are changed.In our β-cell model, the fast variables are voltage (V), the activation variable for the voltage-dependent delayed rectifying K + current (n) and the cytosolic Ca 2+ concentration ([Ca 2+ ] i ), which act on the order of milliseconds.In the model, s is the only slow variable, which acts on the order of minutes.Since [Ca 2+ ] i does not affect other variables in the system, we simply ignore this variable in the following fast/slow analysis.
In the two-cell model, setting g c to 0 dissociates the cells.Therefore, the dynamics of each model cell can be investigated separately.In Fig. 10, we perform a fast/slow analysis for each cell, but describe the steps of the analysis only for β b .We start the analysis by generating a bifurcation diagram of the fast subsystem of β b , using s as a bifurcation parameter (Fig. 10A).Although negative s values are physiologically meaningless, they are valid in the differential equations, and allow for completion of the fast-subsystem bifurcation diagram.For small values of s, the system is at a depolarized stable steady state (Fig. 10A, upper red branch).Increasing s, the system goes through a Hopf Bifurcation (HB), where the branch of stable steady states loses stability, and periodic solutions emerge (blue).These periodic solutions represent  S3 were used for a cluster of 125 cells.] i for all cells in the β-cell cluster for different g c .For parameters other than g K(ATP)) and g c , standard values reported in Table S3 were used for a cluster of 125 cells.continuous spiking in the fast subsystem.The branch of unstable steady states (black) turns at a saddle node bifurcation (SNB2).This branch turns at another saddle node bifurcation (SNB1) and regains stability (lower red branch).The upper and lower blue branches that emerge at the Hopf bifurcation show the maximum and minimum values that V takes on during the continuous train of action potentials.The periodic solutions are terminated at a homoclinic orbit (HC), where the lower periodic branch intersects the branch of unstable saddle points.The bifurcation diagram, which we refer to as a `z-curve', summarizes the dynamic behavior of the fast subsystem as a function of the slow  S3 were used for both cells.In each panel, the red and black curves represent the stable and unstable stationary solutions of the fast subsystem, respectively.The blue curves represent the minimum and maximum voltage values of periodic solutions for the fast subsystem.The magenta curve is the s-nullcline defined by s ∞ (Eq.10).Insets show the voltage time course of the corresponding cell, which are also superimposed onto the bifurcation diagrams shown on the s-V axes.HB: Hopf bifurcation; SNB: saddle node bifurcation; HC: homoclinic orbit.For parameters other than g K(ATP)) and g c , standard values reported in Table S3 were used for both cells.variable, s.This diagram shows that the fast subsystem is at rest either at a depolarized or at a hyperpolarized stable steady state, at the high or low ends of the s spectrum, respectively.However, for moderate values of s, periodic solutions emerge, and the fast subsystem generates periodic action potentials.Furthermore, the fast subsystem is bistable for a range of s values between 0.35 and 0.5, where stable periodic solutions and lower voltage stable steady state solutions coexist.Now we incorporate the slow dynamics into the picture to understand how the cell switches between active and silent phases during a burst.For this purpose, we superimpose the projection of the s-nullcline (magenta), which is given by the s ∞ function (Eq.10), onto the z-curve.The z-curve is then considered as a generalized nullcline of the fast subsystem.Therefore, the intersection of the s ∞ curve and the z-curve determine the steady state of the full system.If the intersection takes place on the branch of lower stable steady states, the equilibrium is stable and the system remains at a hyperpolarized silent state.On the other hand, if the intersection takes place deep in the periodic branch, the equilibrium is unstable and the system generates continuous spiking.Finally, if the intersection takes place on the middle branch of unstable steady states, the equilibrium is unstable and the system can switch periodically between spiking active phases and hyperpolarized silent phases following the timescale of the slow variable.That is, the system can burst.The inset in Fig. 10A shows the bursting oscillations generated by the active cell when the s-nullcline intersects the z-curve between the left saddle node bifurcation (SNB1) and the homoclinic bifurcation (HC).This bursting trajectory (green) is also superimposed onto the bifurcation diagram.At the silent phase of a burst, the phase point is on the branch of stable steady states.On the right side of the s-nullcline, the phase point moves to the left because ds/dt is negative in this region (Eq.8).Therefore, during the silent phase of a burst, the phase point follows the branch of stable steady states until SNB1, where it jumps to the periodic branch.When on the periodic branch, spikes are produced and this is the active phase of a burst.Notice that during the active phase, the phase point remains on the left side of the s-nullcline, where ds/dt is positive (Eq.8).Therefore, during the active phase of a burst, the phase point moves to the right until it reaches the homoclinic bifurcation that ends the periodic branch, where it falls back to the branch of lower stable steady states and completes a burst cycle.In Fig. 10B, we perform the same analysis for the silent cell (β s ).For β s , g K(ATP) is set to 200 pS, which shifts the bifurcation diagram to the left (compared to the bifurcation diagram of β b ).Consequently, the s-nullcline intersects the zcurve on the branch of lower stable steady states, so β s remains at the hyperpolarized silent state (inset).
Increasing g c from 0 to 25 pS, the two cells interact with one another according to the relation defined by Eqs.16 and 17, and β s is recruited to generate bursts of action potentials (Fig. 9C).To understand the way bursting in β s depends on β b , we generate a bifurcation diagram of β s using β b 's voltage variable (V b ) as the bifurcation parameter, and setting s s to its average value over the burst period.This allows us to observe the way V s switches between the active and the silent phases during a burst with the help of β b .The bifurcation diagram has a reflected z-curve shape, or an `s-curve', with the same sequence of bifurcations as in Fig. 10 but in the opposite order (Fig. 11A).
To view how bursting in β b recruits β s to burst, we simplify the coupled system by removing the spikes from β b by setting the delayed rectifier activation variable n b to its voltage-dependent equilibrium n b,∞ (V b ).This converts β b into a relaxation oscillator whose fast variable, V b , has the shape of a square wave (Fig. 11B).With this change, β s is still converted to a burster when coupled to β b (Fig. 11C).The burst trajectory is then superimposed onto the s-curve, with different segments color coded to reflect different portions of the β b relaxation oscillation (Fig. 11A).When V b is in a low phase (Fig. 11A, cyan), the β s phase point is on the lower branch of the s-curve, so the cell is in a silent phase.When V b transitions from a low to a high state (Fig. 11B, orange), the β s phase point quickly moves to the right side of SNB2 (Fig. 11A, orange), and enters the basin of attraction of the stable periodic spiking solutions.During the plateau of the relaxation oscillation V b slowly declines (Fig. 11B, green), moving the phase point of β s leftward along the spiking branch (Fig. 11A, green), which is the active phase of a burst V s .As V b drops to a low state of the relaxation oscillation (Fig. 11B, magenta), the β s phase point accelerates its leftward motion (Fig. 11A, magenta).It moves past the homoclinic bifurcation and is attracted to the lower branch of fast-subsystem steady states.The system thus enters a silent phase and completes a burst cycle.
An interesting feature of the coupled two-cell system is that at an intermediate coupling strength, the intrinsically silent cell was recruited to burst, while at larger coupling strengths the intrinsically bursting cell  S3) were used for both cells.
was silenced (Fig. 9).How did this silencing happen?We explain this through a fast/slow analysis in Fig. 12, with strong gap junctional coupling (g c = 50 pS).To generate the fast subsystem bifurcation diagram for β b we set V s to its hyperpolarized equilibrium value.The zcurve is left shifted compared to that of Fig. 10A, so that now the s bnullcline intersects the bottom branch of stable steady states, creating a stable full-system equilibrium in β b .We use this equilibrium value to generate the z-curve for β s and once again there is an intersection of the slow nullcline with the bottom stable branch of the z-curve, creating a stable full-system equilibrium for β s .Therefore, with this strong coupling the asymptotic state for both cells is a hyperpolarized equilibrium.

Discussion
Gain of function mutations in the pore forming Kir6 subunits of ATP sensitive K + channels of pancreatic β-cells is the major cause of neonatal diabetes mellitus in humans [26,46,47].In this study, we showed that a gain of function mutation in Kir6.1 causes overexpression of connexin36 (Cx36) protein in insulin secreting mouse β-cell lines (Fig. 3).Using computational modeling, we showed that such an increase in Cx36 expression plays an important role in the impairment of glucose stimulated Ca 2+ oscillations in β-cell clusters with gain of function mutations in their K(ATP) channels (GoF β-cells).Our results showed that without an increase in Cx36 expression, a moderate gain of function mutation in K(ATP) channels might not be sufficient to diminish glucose stimulated [Ca 2+ ] i oscillations as seen in neonatal diabetic β-cell lines (Fig. 5).Our fluorometric [Ca 2+ ] i imaging data verified that reducing Cx36 expression in NIT-1 mouse β-cells with Kir6.1 gain of function mutations restored glucose stimulated cytosolic [Ca 2+ ] i oscillations in some cells (Fig. 6).Our modeling results showed that clusters of wild type β-cells responded differently from clusters of GoF β-cells to changes in gap junctional coupling strength.In a cluster of wild type β-cells, overall cluster activity exhibited a monotonically decreasing trend with reduced coupling strength between the cells (Fig. 7), whereas in a GoF β-cell cluster, there was an inverted u-shaped non-monotonic relation between cluster activity and the coupling strength (Figs.8,9).A fast/slow analysis of a two-cell system demonstrated the origin of the u-shaped relationship (Figs.[10][11][12]. Gain of function mutation Kir6.2 subunits of the K(ATP) channels (coded by the Kcnj11 gene) have been used for developing diabetic human β-cell models before [48].Recently, Remedi et al. showed that Kir6.1 subunits (coded by the Kcnj8 gene) are also expressed in mouse pancreatic β-cells [49].More importantly, both single Kcnj8 [G343D] and double Kcnj8 [G343D, Q53R] gain of function mutations were demonstrated to cause either glucose intolerance or diabetes, respectively, similar to the effect of Kir6.2 GoF in mouse β-cells.Remedi et al. also showed that only 2 out of 10 Kcnj8 [G343D] mice showed diabetes with dramatically elevated fasting blood glucose levels and a marked impairment in glucose tolerance [49].However, all Kcnj8 [G343D, Q53R] mice developed severe diabetes with profoundly impaired glucose tolerance, as a result of diminished plasma insulin levels in response to glucose challenge, pointing to the better potential of Kcnj8 [G343D, Q53R] double gain of function mutation as a model to study diabetes [49].Previous studies showed that gain of function mutations in Kir6 subunits of K(ATP) channels diminish glucose stimulated [Ca 2+ ] i oscillations and insulin release [31,50].These studies also showed that Cx36 knockout results in increased islet activity and insulin release in these cells.Our results confirm these findings and show that in a cell cluster that exhibits a gain of function mutation in K(ATP) channels, Cx36 knockdown results in increased activity with loss of coordinated oscillations (Fig. 6).These results show that when coupling is strong, in a heterogenous cell cluster, intrinsically silent cells suppress the intrinsically bursting cells and diminish [Ca 2+ ] i oscillations.When coupling between the cells is removed, intrinsically bursting cells produce [Ca 2+ ] i oscillations, so the overall level of insulin secretion would be higher than with the strong coupling.However, the current study takes this one step further and shows that the overall activity of a cell cluster is a non-monotonic inverted u-shaped function of the coupling strength between the cells.Indeed, with coupling of intermediate strength, intrinsically bursting β-cells can recruit the silent cells to become bursters, thereby increasing the level of insulin secretion from the population.This result indicates the existence of an optimum level of coupling between the cells that maximizes insulin secretion from the GoF cells.
Our experimental results show that β-cells with gain of function mutation in Kir6.1 overexpress Cx36, and they do not generate [Ca 2+ ] i oscillations at stimulatory glucose levels (Fig. 6C and E; Supplemental video 1).However, targeted knockdown of Cx36 using shRNAs restores oscillations in a fraction of the cells (Fig. 6D and E; Supplemental video 2).This finding is consistent with previous studies [31,50] showing an increased insulin release and activity in β-cells with Kir6.2 gain of function mutation after Cx36 knockout.Although the cells produce out of phase oscillations, we note that the in vitro model we use does not allow formation of 3-dimensional large β-cell clusters.Therefore, it is not possible to conclude whether dyscoordination between β-cells result from very low Cx36 expression or the 2-dimensional morphology of the in vitro cell clusters.Indeed, islet architecture plays an important role in coordination of [Ca 2+ ] i oscillations in pancreatic β-cells [51].The islet architecture is one of the major determinants of the average number of Fig. 12. Fast/slow analysis for the case of strong coupling, g c = 50 pS.(A) The intrinsically bursting β-cell β b has a stable rest state at a hyperpolarized voltage when V s is at its equilibrium value, due to intersection of the s-nullcline (magenta) with the stable lower branch (red) of the z-curve.(B) The intrinsically silent cell β s is at a stable hyperpolarized equilibrium when V b is at its equilibrium value.For parameters other than g K(ATP)) and g c , standard values (Table S3) were used for both cells.neighbors a cell interacts with [52].Another limitation of the current study is that the mathematical model only accounts for the ionic current mediated through gap junctional connections and ignores the transfer of other metabolites and signaling molecules, which may profoundly affect [Ca 2+ ] i oscillations [53,54].This effect could be explored by employing a more comprehensive β-cell model that accounts for cellular energy metabolism [55] as well as paracrine and autocrine pathways [56].However, this would increase the computational cost significantly and make dynamical analysis less tractable.Previously, such a modeling paradigm was used to model a small cell cluster, and it was shown that increased Ca 2+ permeability between cells improved coordination within the islets even when electrical coupling was weak [57].However, due to the significantly high computational cost the study was limited to a small 3×3×3 β-cell network.
Our results show that reduced Cx36 expression improves Ca 2+ oscillations in Kir6.1 GoF NIT-1 cells.These findings suggest Cx36 as a potential target for neonatal diabetes.Previous studies show a relation between Cx36 and epilepsy in humans and animals [58].Quinine, an anti-malaria drug, is suggested as a potential treatment in these patients to relieve epileptic seizures by targeting Cx36 [59].This family of drugs might be a potential compound that targets Cx36.However, mefloquine, another anti-malaria drug that also blocks Cx36 channels, was shown to inhibit depolarization-induced insulin release in mouse pancreatic islets, but not glucose-induced insulin secretion [60].Since mefloquine also blocks both K(ATP) and L-type Ca 2+ channels, the mechanism of its net effect on insulin release is not clear [60].Nevertheless, these studies suggest the potential of this family of drugs as a means for blocking Cx36.Our results show a non-monotonic inverted U-shaped relation between gap junctional coupling strength and cluster activity.This implies that any therapy that targets Cx36 must be well-optimized to induce an improving effect.
In conclusion, our results show that in an in vitro mouse β-cell model, gain of function mutation in Kir6.1 results in a significant increase in Cx36 expression.Furthermore, computational simulations and fast/slow analysis suggest that in neonatal diabetes mellitus, gain of function mutations in K(ATP) channels and the resulting increase in Cx36 might be partners in crime in diminishing glucose stimulated [Ca 2+ ] i oscillations.The experimental data provided by the current study and earlier work also suggest that Cx36 might be a potential target for the treatment of neonatal diabetes.

Fig. 3 .
Fig. 3. (A) Relative expression of Kcnj8 in cells with stable overexpression of wild type or mutant Kcnj8 [G343D, Q53R] at the mRNA level (B) Relative expression of Kcnj8 in cells with stable overexpression of wild type or mutant Kcnj8 [G343D, Q53R] at the protein level (C) Relative expression of Gjd2 in cells with stable overexpression of wild type or mutant Kcnj8 [G343D, Q53R] at the mRNA level (D) Relative expression of Gjd2 in cells with stable overexpression of wild type or mutant Kcnj8 [G343D, Q53R] at the protein level (E) Confocal microscopy images of Kcnj8 and Gjd2 in cells with stable overexpression of control plasmid or mutant Kcnj8 [G343D, Q53R].The DAPI staining is used to determine the location of the cell nuclei.Control refers to the cells stably transfected with pCMV6-Entry Mammalian Expression Control vector.*p < 0.05, **p < 0.01, ***p < 0.001

Fig. 5 .
Fig. 5. (A) Population average [Ca 2+ ] i level in a heterogenous β-cell cluster with g K(ATP) chosen from the distribution N(150 pS, 45 pS) and g c =45 pS.(B) Average [Ca 2+ ] i level in a cell cluster with g K(ATP) chosen from N(200 pS, 60 pS) and g c =75 pS.(C) Average [Ca 2+ ] i level in a cell cluster with g K(ATP) chosen from N(200 pS, 60 pS) and g c =45 pS.(D) Bursting region in the plane of the mean of g K(ATP) and coupling strength g c .The green arrow indicates an increase in mean g K(ATP) , whereas the red arrow indicates a combined increase in both mean g K(ATP) and g c .For parameters other than g K(ATP)) and g c , standard values reported in TableS3were used.

Fig. 7 .
Fig. 7. Intracellular Ca 2+ concentration of a wild-type-like model β-cell cluster for different g c levels.The K(ATP) channel conductance for cells in the cluster is sampled from the green distribution in Fig. 4A.(A) The average [Ca 2+ ] i time course of the β-cell cluster for different g c levels (75-0 pS).The g c level is reduced by 15 pS steps after 30 min.(B) The distribution and mean of the 30 min average of [Ca 2+ ] i for all cells in the β-cell cluster for different g c .For parameters other than g K (ATP)) and g c , standard values reported in TableS3were used for a cluster of 125 cells.

Fig. 8 .
Fig. 8. Intracellular Ca 2+ concentration of a β-cell cluster with gain of function mutation K(ATP) channels for different g c levels.The g K(ATP) values for cells in the cluster are sampled from the red distribution in Fig. 4A.(A) The average [Ca 2+ ] i time course of the β-cell cluster for different g c levels (75-0 pS).The g c level is reduced by 15 pS steps each 30 min.(B) The distribution and mean of the 30 min average [Ca 2+ ] i for all cells in the β-cell cluster for different g c .For parameters other than g K(ATP)) and g c , standard values reported in TableS3were used for a cluster of 125 cells.

Fig. 9 .
Fig. 9. (A, C, E) Voltage and (B, D, F) [Ca 2+ ] i time courses of an intrinsically-bursting β-cell (β b , green) and an intrinsically-silent β-cell (β s , orange) for varying coupling strengths (g c =0, 25 and 50 pS) in a two-cell model.In panels B, D and F, the magenta line shows the average value of the two cells.For parameters other than g K(ATP)) and g c , standard values reported in TableS3were used for both cells.

Fig. 10 .
Fig. 10.Fast/slow analysis of the dynamics of the intrinsically bursting (β b ) (A) and the intrinsically silent (β s ) (B) β-cell.V b and V s show the membrane potential of the β b and β s , respectively.In each panel, the red and black curves represent the stable and unstable stationary solutions of the fast subsystem, respectively.The blue curves represent the minimum and maximum voltage values of periodic solutions for the fast subsystem.The magenta curve is the s-nullcline defined by s ∞ (Eq.10).Insets show the voltage time course of the corresponding cell, which are also superimposed onto the bifurcation diagrams shown on the s-V axes.HB: Hopf bifurcation; SNB: saddle node bifurcation; HC: homoclinic orbit.For parameters other than g K(ATP)) and g c , standard values reported in TableS3were used for both cells.

Fig. 11 .
Fig. 11.The dynamics of the coupled two-cell system for g c = 25 pS.(A) The bifurcation diagram of the fast subsystem of the intrinsically silent β-cell (β s ) using the bursting β-cell's (β b ) membrane potential (V b ) as the bifurcation parameter.To generate the bifurcation diagram, the slow variable of β s (s s ) is set to its average value over a burst period (0.3).Black arrows represent direction of the flow.(B) The voltage time course of β b , where the activation variable of the delayed rectifier K + current (n b ) is set to its steady state (n b,∞ ), so that β b becomes a relaxation oscillator.(C) The bursting voltage time course of β s for s s =0.3.The time courses and trajectories are color-coded according to segments of the β b relaxation oscillation.SNB: Saddle node bifurcation; HC: Homoclinic orbit; HB: Hopf Bifurcation.For parameters other than g K(ATP)) and g c , standard values (TableS3) were used for both cells.