Cooperative regulation of cellular identity in systems with intercellular communication defects
CCooperative regulation of cellular identity insystems with intercellular communicationdefects
Nataliya Stankevich , and Aneta Koseska September 24, 2019 Department of Applied Cybernetics, St. Petersburg State University,Saint-Petersburg, Russia Department of Radio-Electronics and Telecommunications, Yuri GagarinState Technical University of Saratov, Saratov, Russia Department of Systemic Cell Biology, Max Planck Institute of MolecularPhysiology, Dortmund, Germany
Abstract
The cooperative dynamics of cellular populations emerging fromthe underlying interactions determines cellular functions, and therebytheir identity in tissues. Global deviations from this dynamics onthe other hand reflects pathological conditions. However, how theseconditions are stabilized from dis-regulation on the level of the sin-gle entities is still unclear. Here we tackle this question using genericHodgkin-Huxley type of models that describe physiological burstingdynamics of pancreatic beta-cells, and introduce channel dis-functionto mimic pathological silent dynamics. The probability for patholog-ical behavior in beta-cell populations is ∼ ∼ a r X i v : . [ q - b i o . CB ] S e p hysiological bursting dynamics of a population of beta-cells is co-operatively regulated, even under intercellular communication defectsinduced by dis-functional channels of single-cells.Keywords: Beta-cells dynamics; disfunctional channels; heteroge-neous cell population. The emergent properties of complex biological systems are largely deter-mined by the interactions of the elements that constitute these systems.Synchronous firing of ensembles of interacting neurons enables for exam-ple, complex collective tasks such as Hebbian learning and memory [1], [2]or wake-sleep cycles [3], [4] to be executed. On the other hand, the col-lective behavior of homogeneous or heterogeneous cell types including theirarrangement and types of contacts determines specific functions of tissuesand organs. For instance, it has been experimentally observed that the well-synchronised electrical activity (spiking and bursting) of primary beta-cellsin the pancreatic islet is what precedes insulin secretion [5], [6], [7], wherethe secretion increases with the fraction of time that the cells spend in thespiking state. The duration of the silent phase between two bursts is in thiscase regulated by the rate at which calcium is removed from the cell interior.The spiking oscillations typically display a period of 1-10s, whereas the du-ration of the bursting period varies from 0.2 to 5min. In contrast, isolatedcells show electrical activity with very different time scales, do not burst andthereby do not secrete insulin effectively [8], [9]. In this case, the openingprobability of the potassium channels is too small for the individual cells topresent regular spiking signal that conversely only occurs in population ofsynchronized cells, thereby demonstrating that the intercellular communica-tion in the organized islet structure is a necessary prerequisite for effectivefunctioning of this coupled-cell system. Deviation from the characteristicdynamical behavior of the system on the other hand reflects pathologicalconditions. How such global deviations in the dynamics are stabilized fromdis-regulation on the level of the interacting elements is however unclear.This question becomes even more prominent considering that it has been2xperimentally demonstrated that generically, the collective behavior of cellson the level of tissues is generally stable against mutations that occur insingle cells, allowing them to retain their identity [10].Here we use numerical simulations to tackle how a change in the globaldynamics of a coupled system can be established by dis-regulation on a sub-population level. To address this, we employed a generic Hodgkin-Huxleyformalism which on a single cell level describes physiological bursting activityof pancreatic beta-cells [11], [12], [13]. Additionally, we also use a modifiedversion of this model where a potassium-like ion channel with decreasedopening probability is introduced [14] to stabilize a silent state in additionto the bursting activity, thereby enabling bistability on a single-cell level tobe established. This silent state mimics ion channel disfunction in a form ofblocking or inactivation. We use these models to investigate how deviationfrom the physiological dynamics of the population of coupled beta-cells istriggered depending on the number of cells acquiring disregulation in channelactivity. Three cases are therefore analyzed: (1) dynamics of two coupledbeta-cells under physiological conditions (bursting activity), (2) dynamicalbehavior of two coupled cells that have acquired the dis-function (stablesteady state), but one or both of them do not manifest it in isolated cellsand (3) dependence of the global population dynamics on number of cellswith channel defects. We found that the dynamics of a population of twocoupled cells, each characterized with bistability between silent and burstingstate, is either equivalent to the same attractor in which both cells are poisedwhen in isolation, or dictated by the physiological bursting dynamics whenthe uncoupled cells are initially poised in different attractors. In contrast,the bursting dynamics of a heterogeneous population in which only a portionof cells can exhibit bistability between silent and bursting dynamics and aresubsequently poised in the silent state, can be only disrupted if at least N/ Results
The electrical activity of pancreatic beta-cells relies on multiple types ofvoltage- and ligand-gated ion channels that are permeable to inorganic ionssuch as sodium, potassium, chloride and calcium [5],[15], [16]. These chan-nels not only regulate membrane potential, ion homeostasis and electricalsignaling [17], but a dis-regulation in their activity has been recently alsoimplicated in tumorogenesis and tumor progression [18], [19], [20], [21].To depict the bursting dynamics of pancreatic beta-cells as observed un-der physiological conditions, as well as to simultaneously model silent dy-namics caused upon ion channel dis-regulation, we use a model based on theHodgkin-Huxley formalism [12], [13] with modifications proposed by Stanke-vich et al. [14]: τ dV ( i ) dt = − I C a ( V ( i ) ) − I K ( V ( i ) , n ( i ) ) − k ( i ) I K ( V ( i ) ) − I S ( V ( i ) , S ( i ) ) τ dn ( i ) dt = σ ( n ∞ ( V ( i ) ) − n ( i ) ) τ S dS ( i ) dt = S ∞ ( V ( i ) ) − S ( i ) (1)where V ( i ) represents the membrane potential of the i-th cell, the functions I C a ( V ( i ) ), I K ( V ( i ) , n ( i ) ), I S ( V ( i ) , S ( i ) ) define the three intrinsic currents, fastcalcium, potassium and slow potassium respectively, such that: I C a ( V ( i ) ) = g C a m ∞ ( V ( i ) )( V ( i ) − V C a ) (2) I K ( V ( i ) , n ( i ) ) = g K n ( V ( i ) − V K ) (3) I S ( V ( i ) , n ( i ) ) = g S S ( V ( i ) − V K ) (4)The gating variables for n ( i ) and S ( i ) are the opening probabilities of thefast and slow potassium currents: ω ∞ ( V ( i ) ) = [1 + exp V ω − V ( i ) Θ ω ] − , ω = m, n, S (5)4igure 1: Dynamical structure of the beta-cell models (Eqs. 1-7).
Schematic representation of a beta-cell in a. physiological conditions (blue:Ca-channel, green:K-channel with physiological opening probability) and b. with an additional potassium channel with decreased opening probability(red). Fast and slow manifolds (black solid and dashed lines, respectively),attractors and their basins of attraction (blue: bursting attractor, red: stablenode) are shown for c. k ( i ) = 0, corresponding to Fig. 1a and d. k ( i ) = 1,Fig. 1b. Stable and unstable nodes are marked with red and blue circles,respectively. Parameters: τ = 0 . τ S = 35, σ = 0 . g Ca = 3 . g K = 10 . g K = 0 . g S = 4 . V Ca = 25 . V K = − . V m = − . θ m = 12 . V n = − . θ n = 5 . V S = − θ S = 10 . V P = − . θ p = 1 . I K ( V ( i ) ) on the other hand defines an additional voltage-dependent potassium current. It has been previously demonstrated thata single oscillator, in absence of I K (when k ( i ) = 0, Fig. 1a), is charac-terized with a bursting attractor that is born via a Hopf bifurcation for V S = − . mV . At this parameter value, the equilibrium point looses itsstability [22], [23]. Although the bursting attractor is born in the vicinity ofthe equilibrium point, this point moves away from the bursting attractor forincrease in V S . Even more, the two dimensional projection of the phase por-trait together with the fast and slow manifolds when V S = − mV (Fig. 1c)demonstrates that the periodic trajectories do not intersect the neighborhoodof the steady state and the bursting attractor terminates in a homoclinincbifurcation as the trajectory hits the slow manifold [22], [23]. Calculatingthe eigenvalues in this case shows that the equilibrium point is an unstablenode ( λ , λ , λ ) = (23 . , − . , . I K ( V ( i ) ) ( k ( i ) (cid:54) = 0, Fig. 1b) that varies strongly withthe membrane potential in the vicinity of the equilibrium point however, theunstable node is stabilized without affecting the global flow of the model(Fig. 1d, corresponding eigenvalues ( λ , λ , λ ) = ( − . , − . , − . g K = 0 . V (1)0 , n (1)0 , S (1)0 ) = (-49.084,0.0027105, 0.19648)) [14]. I K ( V ( i ) ) is given in the form: I K ( V ( i ) ) = g K p ∞ ( V ( i ) )( V ( i ) − V K ) (6)where p ∞ ( V ( i ) ) = [ exp V ( i ) − V p Θ p + exp V ( i ) p − V Θ p ] − (7)represents the opening probability of this channel. In contrast to the otherpotassium channel (Eqs. 3, 5) which opens with probability n ∞ = 1 when themembrane voltage reaches a threshold value, the opening probability of themodified channel will be equal only to 0.5. From physiological point of view,this can be interpreted as an ion channel disfunction, as for instance blockingof the potassium channel or its inactivation [24], and thereby a stable silentstate emerges. Between the stable node and the bursting attractor there isa rejecting current which enables the system to remain in the stable steadystate when starting from initial conditions in its vicinity [14]. Generally, themodified model ( k ( i ) (cid:54) = 0) depicts bistability between physiological, burstingbeta-cell dynamics and a pathological, silent state dynamics (Fig. 1d).6 .2 Dynamics of two coupled beta-cells under physio-logical conditions To investigate the global dynamics of a population of beta-cells, we intro-duced electrical coupling by adding the following gap-junctional couplingterm to the equation for V ( i ) that describes bi-directional transport of ionsbetween the cells [7], [25]: I C ( V ( i ) ) = (cid:88) j ∈ Γ i g c,V ( V ( i ) − V ( j ) ) (8)with g c,V being the coupling conductance (coupling strength). We con-sider only electrical coupling of the cells (Fig. 2a) by adding the followinggap-junctional coupling terms to the equations for V in (Eqs. 1). The sumis taken over the whole population, assuming global coupling.As already noted, if k i = 0 and in absence of coupling, the isolated os-cillators display a single stable state - bursting dynamics (Fig. 1c). To inferthe effect of coupling on the system’s dynamics, we investigated next thebifurcation structure of a minimal model of two coupled cells as a functionof the coupling strength g c,V . Bifurcation trees using the membrane voltagevariable V (1) were therefore constructed using a Poincar´ e section at the hy-persurface n (1) = 0 .
003 by scanning g c,V ∈ [0 , .
6] from left to right and viceversa (Fig. 2b violet and grey, respectively). Changing the scanning direc-tion of the bifurcation parameter unraveled that for g c,V ∈ [0 , . g c,V = 0 .
15 inFig. 2e. Given that the synchronized bursting state lies only on the diago-nal of the basin (grey line in Fig. 2e), it is most probable to reach one ofthe other dynamical states in this parameter region. For g c,V > .
19 how-ever, only synchronised bursting between both oscillators was observed. Thisanalysis therefore demonstrates that the physiological, synchronized burstingbehavior is the dominant dynamical behavior over a wide coupling range ina population of identical beta-cells that communicate via gap junctions.7igure 2:
Behavior of two coupled beta-cells with bursting dynamics.a.
Schematic representation of the model. b. Bifurcation trees estimatedfrom Poincar´e sections at the hypersurface n (1) = 0 .
003 when g c,V is scannedin both directions (violet: forward, gray: backward). k ( i ) = 0 , i = 1 , c. periodic or d. chaoticattractors are also shown for g c,V = 0 .
03 and 0 .
15 respectively. e. Estimatedbasins of attraction of the coexisting bursting attractors for g c,V = 0 .
15 whenstarting from the following initial conditions: n (1)0 = 0 . S (1)0 = 0 . n (2)0 = 0 . S (2)0 = 0 . .3 Dynamics of two coupled beta-cells that containdis-functional potassium channels For k i = 1 however, the presence of the potassium channel which has a de-creased opening probability (Fig. 1b, Eqs. 6,7) ensures that the equilibriumpoint is also stabilized, such that isolated oscillators display bistability be-tween a bursting and a silent state in absence of coupling (Fig. 1d). Thus, forinitial conditions in the vicinity of the stable node, the physiological burstingbehavior is not observed. We therefore refer to this silent regime as patholog-ical. To unravel how such channel defects affect the global dynamics of thepopulation, we used again a minimal model of two coupled beta-cells that areindividually characterized with bistability between the two states (Fig. 3a),and investigated the two possible scenarios: i) the initial conditions of bothoscillators are poised either in the bursting attractor or in the stable nodeattractor, and ii) one of the oscillators is poised in the stable node, whereasthe other oscillator is poised in the bursting attractor.When the initial conditions of both oscillators are poised in the same at-tractor, the solution of the coupled system is trivial - the population behavioris synchronized and identical to the dynamics of the cells in isolation (resultsnot shown). In contrast, for scenario ii) , when isolated cells populate distinctattractors such that the initial conditions for one cell are poised in the burst-ing whereas for the other in the stable node, the bifurcation tree constructedanalogously to the one in Fig. 2b has a complex structure for varying couplingstrength g c,V (Fig. 3b, top). To identify the existing dynamical solutions inthis case, we complemented the bifurcation diagram with direct calculationsfrom semi-random initial conditions (Fig. 3b, bottom). In particular, for each g c,V we calculated 500 time series for the system of two coupled beta-cellssuch that the initial conditions for one of the oscillators were chosen exactlyat the stable node ( V (1)0 , n (1)0 , S (1)0 ) = ( − . , . , . V (2)0 ∈ ( − , − n (2)0 ∈ (0 , .
12) and S (2)0 ∈ (0 . , . g c,V < . g c,V sub-intervals for which bursting or chaoticdynamics was observed, the probability for the coupled system to exhibitsteady state solution was less than 10%. That the bursting or the chaoticdynamics are the dominant regimes under weak coupling is also reflected inthe corresponding basins of attraction estimated for the respective g c,V values(Fig. 3c,d). In contrast, for a large part of intermediate coupling strengthinterval (0 . < g c,V < . ∼ g c,V ranges. The exemplary basin of attraction plot estimated for g c,V = 0 .
02 also depicts the dominance of the stable equilibrium point, withonly small lines of the basin of the bursting attractor (Fig. 3e). This is sur-prising, given the small basin of attraction of the steady state for isolatedcells (Fig. 1d). The mechanism leading to this stabilization will be discussedelsewhere (manuscript in preparation).A further increase in g c,V continuously decreases the probability for thefull system to exhibit steady state behavior. At g c,V ≈ .
12 for example,the system visits the bursting attractor and the equilibrium point with equalprobability. This implies that in the presence of noise, the full system canswitch between bursting physiological and silent pathological state. For highenough coupling strength g c,V however, the basin of attraction of the steadystate is significantly decreased, leaving the bursting attractor as the mostcommon solution of the system (Fig. 3f). These results therefore demon-strate that a system of two coupled beta-cells that have disfunctional potas-sium channels and thereby possibility to manifest a stable steady state, pre-dominantly exhibits physiological bursting dynamics. This is true for weakand strong coupling strengths, even when the initial conditions of one of theoscillators correspond to its positioning in the stable node. Only for interme-diate coupling values the silent pathological state is the dominant dynamicalregime, even for small steady state basin of attraction for single cells.10igure 3: Behavior of two coupled beta-cells with disfunctionalpotassium channels. a.
Schematic representation of two the system. b. (Top) Bifurcation tree estimated from Poincar´e sections at the hypersurfaceand n (1) = 0 .
003 constructed for fixed initial conditions: V (1)0 = − . n (1)0 = 0 . S (1)0 = 0 . V (2)0 = − . n (2)0 = 0 . S (2)0 =0 . g c,V . k ( i ) = 1 , i = 1 , c. g c,V = 0 . d. g c,V = 0 . e. g c,V = 0 . f. g c,V = 0 .
12. For c.-f. , the initial conditions of the other variables werefixed in steady state n (1)0 = 0 . V (2)0 = − . n (2)0 = 0 . S (2)0 = 0 . .4 Dynamical behavior of a minimal mixed popula-tion model To model physiologically more realistic scenario, we investigate next the dy-namics of a minimal mixed system of two coupled cells, where only one ofthem displays bistability between bursting and a silent state, and thus adis-regulation in a form of inactivation or blocking of a potassium channel(Fig. 4a).Bifurcation analysis showed that the mixed population does not displaysteady state dynamics over the whole g c,V interval (Fig.4b). This is contraryto the previous case (Fig. 3b) where silent state was dominant for intermedi-ate g c,V values when bistability characterizes both cells, although the initialconditions of only one of them were poised in the stable node. A typical dy-namics of the coupled mixed population model is mixed-mode oscillations, de-fined as complex oscillatory patterns where small-amplitude high-frequencyoscillations alternate with large-amplitude long-period cycles [27]. Increasingthe coupling strength between the two mixed oscillators induces transitionsbetween different manifestations of the mixed-mode oscillations. For ex-ample, under low coupling strength, one oscillator displays small-amplitudewhereas the other large-amplitude bursting dynamics (Fig. 4c). For inter-mediate coupling, different patterns of alternating between the small- andthe large-amplitude bursting occur (Fig. 4d-f), to finally converge to almostsynchronized bursting of both oscillators under strong coupling (Fig. 4g).The analysis therefore demonstrates that in a mixed population of two cou-pled cells, where only one of them has dis-regulation of a potassium channel,pathological silent dynamics is not stabilizied on the level of the population.This implies that deviations from the physiological bursting dynamics mightbe possible for a larger sub-population size that has channels dis-regulation. To identify whether stabilizing pathological dynamics is a system-size ef-fect, we investigated the dynamics of heterogeneous populations of increasingsizes, starting with a minimal population of N = 3 cells, where ( N/ Behavior of a minimal mixed population model. a.
Schematic representation of minimal population of two coupled hetero-geneous cells. b. Bifurcation tree estimated from Poincar´e sections atthe hypersurface n (1) = 0 .
003 when g c,V is scanned in both directions. k (1) = 0 , k (2) = 1 and other parameters as in Fig. 1. Exemplary time se-ries of characteristic dynamical regimes for: c. g c,V = 0 . d. g c,V = 0 . e. g c,V = 0 . f. g c,V = 0 . g. g c,V = 0 . Behavior of heterogeneous beta-cells population. a. (Top)Schematic representation of networks of increasing size, where only a singlecell does not have a channel dis-regulations (red color of node responds to k ( i ) = 1, purple color to k ( i ) = 0) and (bottom) a corresponding plots of cou-pling strength thresholds for steady state stabilization in dependence on thesize of the population ( N ). Probability of appearing of co-existing attractorsin dependence on the coupling strength for the minimal population b. of N = 3 coupled cells, as well as for c. N = 12. Other parameters as in Fig. 1.14igure 6: Behavior of heterogeneous beta-cells population as de-pending on the number of cells without channel dis-regulations.a,b.
Schematic representation of networks of interacting cells (red color for k ( i ) = 1, purple color for k ( i ) = 0) and c. corresponding plots depicting thethresholds of coupling strength for which the steady steady state is stabilizedas a function of the size of the population ( N ). Other parameters as in Fig. 1.15umerical simulations demonstrated that in this case there is a threshold ofcoupling strength for which the silent state can be stabilized on the popu-lation level (Fig. 5a). However, the probability to observe this pathological,silent dynamics in the heterogeneous cell population is significantly smallerthan the one in the homogenous population. Direct calculations from semi-random initial conditions as in Fig. 3b showed that for the minimal systemof N = 3 cells, this probability is smaller than 1% (Fig. 5b). Subsequentincrease of the population size by a cell that can display bistability betweena silent and bursting dynamics in turn lowered the coupling threshold forwhich the pathological, silent state was observed (Fig. 5a, bottom). Thesimulations also demonstrated that the probability to observe the pathologi-cal state increased with the increase of fraction of cells that have the channeldefects. For example, for a population of N = 12 cells, where 11 of themhave dis-functional channel increased to approximately 10%, almost over thefull range of coupling strengths (Fig. 5c). Despite this increase in the proba-bility to observe the pathological state, the population of heterogeneous cellscan still cooperatively regulate the physiological bursting dynamics, evenwhen the majority of the population has acquired dis-regulation in the chan-nel’s activity. This is contrary to a homogeneous population of cells withdis-regulation, where a coupling interval exists for which the probability toobserve the pathological dynamics is ∼ N/ N = 5 (Fig. 6a,c) or 7 (Fig. 6b,c) cells, respectively.The stabilization threshold scales with the number of cells displaying onlybursting dynamics, being highest when the heterogeneous population hasthe largest number of cells exhibiting solely bursting dynamics. Similarly toFig. 5a, a step-wise increase of the population by a cell that displays bista-bility between the silent and the bursting state also decreases the couplingthreshold for stabilization (Fig. 6c). 16 Discussion
Intercellular communication dictates coordinated dynamics of cellular pop-ulations, and thereby determines cellular identity in terms of their func-tion. For example, insulin secretion of pancreatic beta-cells is preceded bywell-synchronized bursting activity of the population, where the secretionincreases with the fraction of time the cells spend in the spiking state. Incontrast, isolated cells do not burst and thereby do not secrete insulin ef-fectively. Such dynamical behavior of single cells has been related to thedecreased opening probability of the potassium channel, resulting in silentdynamics. We considered here the case when cells can exhibit silent dynam-ics, even when coupled in a beta-cells population. In particular, we addressedthe question how channel defects on the level of single cells that alter theintercellular communication can affect the global dynamics of the populationand under which conditions the dynamics will deviate from the physiologi-cal one. An altered dynamics would in turn affect the cellular function andthereby their identity.The bifurcation analysis of a minimal homogeneous or heterogeneous pop-ulations of N = 2 coupled cells, where respectively either both cells havechannel defects and thereby can exhibit silent state or only one of them,demonstrated very different dynamical behavior on the population level.Whereas the homogeneous population exhibited dominance of the patholog-ical silent state ( ∼ ∼ − cknowledgements N.S thanks the DAAD (Project Nr. 57440915) for the financial support ofthe scientific visit to MPI-MP in Dortmund, 2019.
Author contributions
A.K. conceptualized the study. N.S. performed the numerical simulationsand bifurcation analysis with help of A.K. N.S. and A.K. interpreted theresults and wrote the manuscript.
Conflict of interest
The authors declare that they have no conflict of interest.
References [1] E. K. Miller and T. J. Buschman, Cortical circuits for the control ofattention, Current Opinion in Neurobiology 23, 216 (2013).[2] R. Follmann, E. E. Macau, E. Rosa, and J. R. Piqueira, Phase oscillatorynetwork and visual pattern recognition, IEEE Transactions on NeuralNetworks and Learning Systems 26, 1539 (2014).[3] S. Daan, D. Beersma, and A. A. Borbely, Timing of human sleep: re-covery process gated by a circadian pacemaker, American Journal ofPhysiology-Regulatory, Integrative and Comparative Physiology 246,R161 (1984).[4] N. Holmgren Hopkins, P. Sanz-Leon, D. Roy, and S. Postnova, Spikingpatterns and synchronization of thalamic neurons along the sleep-wakecycle, Chaos 28, 106314 (2018).[5] F. M. Ashcroft and P. Rorsman, Electrophysiology of the pancreatic β -cell, Progress in Biophysics and Molecular Biology 54, 87 (1989).186] T. R. Chay, Effects of extracellular calcium on electrical bursting andintracellular and luminal calcium oscillations in insulin secreting pan-creatic beta-cells, Biophysical Journal 73, 1673 (1997).[7] P. Rorsman and F. M. Ashcroft, Pancreatic β -cell electrical activity andinsulin secretion: of mice and men, Physiological Reviews 98, 117 (2017).[8] E. Gylfe, E. Grapengiesser, and B. Hellman, Propagation of cytoplasmicca2+ oscillations in clusters of pancreatic β -cells exposed to glucose, CellCalcium 12, 229 (1991).[9] P. Smolen, J. Rinzel, and A. Sherman, Why pancreatic islets burst butsingle beta cells do not? The heterogeneity hypothesis, Biophysical Jour-nal 64, 1668 (1993).[10] M. J. Bissell and W. C. Hines, Why dont we get more cancer? A pro-posed role of the microenvironment in restraining cancer progression,Nature Medicine 17, 320 (2011).[11] A. L. Hodgkin and A. F. Huxley, A quantitative description of membranecurrent and its application to conduction and excitation in nerve, TheJournal of Physiology 117, 500 (1952).[12] A. Sherman, J. Rinzel, and J. Keizer, Emergence of organized bursting inclusters of pancreatic beta-cells by channel sharing, Biophysical Journal54, 411 (1988).[13] A. Sherman and J. Rinzel, Rhythmogenic effects of weak electrotoniccoupling in neuronal models, Proceedings of the National Academy ofSciences 89, 2471 (1992).[14] N. Stankevich and E. Mosekilde, Coexistence between silent and burst-ing states in a biophysical Hodgkin-Huxley-type of model, Chaos 27,123101 (2017).[15] M. Braun, R. Ramracheya, M. Bengtsson, Q. Zhang, J. Karanauskaite,C. Partridge, P. R. Johnson, and P. Rorsman, Voltage-gated ion channelsin human pancreatic ββ