Cooperative maintenance of cellular identity in systems with intercellular communication defects

The emergent properties of complex biological systems are largely determined by the interactions of the elements that constitute these systems. Synchronous firing of ensembles of interacting neurons enables, for example, 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 collective behavior of homogeneous or heterogeneous cell types including their arrangement and types of contacts determines specific functions of tissues and organs. For instance, it has been experimentally observed that the well-synchronized electrical activity (spiking and bursting) of primary $β$-cells in the pancreatic islet is what drives insulin secretion,^5–7 where the secretion increases with the fraction of time that the cells spend in the spiking state. The duration of the silent phase between two bursts is, in this case, regulated by the rate at which calcium is removed from the cell interior. The spiking oscillations typically display a period of 1–10 s, whereas the duration of the bursting period varies from 0.2 to 5 min. In contrast, although disassociated cells show increased electrical activity at elevated glucose concentrations, this activity is variable from cell to cell, with some cells being electrically “silent” and others continuously “spiking” at any given glucose concentration.^8–12 Thus, synchronization of electrical activity via gap junctions has long been argued as essential for normal glucose-dependent insulin secretion in the intact islet.^13,14 Deviation from the characteristic dynamical behavior of the system, on the other hand, reflects pathological conditions. How such global deviations in the dynamics are stabilized from dysregulation on the level of the interacting elements is, however, unclear. This question becomes even more prominent considering that it has been experimentally demonstrated that in transgenic mice, the collective behavior of $β$-cell populations where approximately $70%$ of the cells express zero functional potassium channels can maintain strong glucose dependence of metabolism, Ca, membrane potential, and insulin secretion from intact islets, although the glucose dependence is lost in isolated transgenic cells.¹⁵ 

Here, we use numerical simulations to tackle how a change in the global dynamics of a coupled system can be established by dysregulation on a subpopulation level. To address this, we employed a generic Hodgkin–Huxley formalism, which on a single cell level describes physiological bursting activity of pancreatic $β$-cells.^16–18 Additionally, we also use a modified version of this model where a potassiumlike ion channel with decreased opening probability is introduced¹⁹ to stabilize a silent state in addition to the bursting activity, thereby enabling bistability on a single-cell level to be established. This silent state mimics ion channel dysfunction in a form of modified opening probability due to nonmonotonic voltage dependence. We use these models to investigate how deviation from the physiological dynamics of the population of coupled $β$-cells is triggered depending on the number of cells acquiring dysregulation in the channel activity. Three cases are, therefore, analyzed: (1) dynamics of two coupled $β$-cells under physiological conditions (bursting activity), (2) dynamical behavior of two coupled cells that have acquired the dysfunction (stable steady state), but one or both of them do not manifest it in isolated cells (heterogeneous population), and (3) dependence of the global population dynamics on number of cells with channel defects (mixed population). We found that the dynamics of a population of two coupled cells, each characterized with bistability between the silent and bursting states but being differentially poised in the stable attractors (heterogeneous case), can be dominated by the silent pathological state for an intermediate coupling range, despite the small basin of attraction of the silent state on a single cell level. In contrast, the bursting dynamics of a mixed population, in which only a portion of cells can exhibit bistability between silent and bursting dynamics and are subsequently poised in the silent state, can be only disrupted if at least $N/2+1$ cells have the channel dysfunction. The probability for the full population to exhibit silent dynamics in this case is, however, less than $10%$⁠. These results thereby demonstrate that the identity of a mixed population of $β$-like cells reflected through their bursting dynamics is cooperatively maintained, even when a large subfraction of this population acquired dysfunctional channels that induce intercellular communication defects.

The electrical activity of pancreatic $β$-cells relies on multiple types of voltage- and ligand-gated ion channels that are permeable to inorganic ions such as sodium, potassium, chloride, and calcium^5,20,21 and regulate membrane potential, ion homeostasis, and electrical signaling.²² 

To depict the bursting dynamics of pancreatic $β$-cells as observed under physiological conditions, as well as to simultaneously model silent dynamics caused upon ion channel dysregulation, we use a model based on the Hodgkin–Huxley formalism^17,18 with modifications proposed by Stankevich and Mosekilde,¹⁹ 

where $V(i)$ represents the membrane potential of the ith cell, the functions $ICa(V(i))$⁠, $IK(V(i),n(i))$⁠, and $IS(V(i),S(i))$ define the three intrinsic currents, fast calcium, potassium, and slow potassium, respectively, such that

The gating variables for $n(i)$ and $S(i)$ are the opening probabilities of the fast and slow potassium currents,

The function $IK2(V(i))$⁠, on the other hand, defines an additional voltage-dependent potassium current. It has been previously demonstrated that a single oscillator, in the absence of $IK2$ [when $k(i)=0$⁠, Fig. 1(a)], is characterized with a bursting attractor that is born via a Hopf bifurcation for $VS=−44.7mV$⁠. At this parameter value, the equilibrium point looses its stability.^23,24 Although the bursting attractor is born in the vicinity of the equilibrium point, this point moves away from the bursting attractor for increase in $VS$⁠. Even more, the two dimensional projection of the phase portrait together with the fast and slow manifolds when $VS=−35mV$ [Fig. 1(c)] demonstrates that the periodic trajectories do not intersect the neighborhood of the steady state and the bursting attractor terminates in a homoclinic bifurcation as the trajectory hits the slow manifold.^23,24 Calculating the eigenvalues in this case shows that the equilibrium point is an unstable node $(λ1,λ2,λ3)=(23.14,−41.83,0.09)$⁠.

In the presence of $IK2(V(i))$ [$k(i)≠0$⁠, Fig. 1(b)] that varies strongly with the membrane potential in the vicinity of the equilibrium point; however, the unstable node is stabilized without affecting the global flow of the model [Fig. 1(d), corresponding eigenvalues $(λ1,λ2,λ3)=(−37.58,−13.6,−0.24)$⁠, estimated for $gK2=0.2$ at the equilibrium point $(V0(1),n0(1),S0(1))=(−49.084,0.0027105,0.19648)$].¹⁹ $IK2(V(i))$ is given in the form,

where

represents the opening probability of this channel. In contrast to the other potassium channel [Eqs. (3) and (5)], which opens with probability $n∞=1$ when the membrane voltage reaches a threshold value, the opening probability of the modified channel will be equal only to 0.5.¹⁹ From physiological point of view, this can be interpreted as an ion channel dysfunction²⁵ due to nonmonotonic voltage dependence, manifested through a stable silent state. Between the stable node and the bursting attractor, there is a rejecting current which enables the system to remain in the stable steady state when starting from initial conditions in its vicinity.¹⁹ Generally, the modified model (⁠$k(i)≠0$⁠) depicts bistability between physiological, bursting $β$-cell dynamics and a pathological, silent state dynamics [Fig. 1(d)].

To investigate the global dynamics of a population of $β$-cells, we introduced electrical coupling by adding the following gap-junctional coupling term to the equation for $V(i)$ in Eq. (1) that describes bidirectional transport of ions between the cells,^7,26

with $gc,V$ being the coupling conductance (coupling strength). In this way, only electrical coupling of the cells [Fig. 2(a)] is considered. The sum is taken over the whole population, assuming global coupling.

As already noted, if $ki=0$ and in the absence of coupling, the isolated oscillators display a single stable state-bursting dynamics [Fig. 1(c)]. To infer the effect of coupling on the system’s dynamics, we next investigated the bifurcation structure of a minimal model of two coupled cells as a function of the coupling strength $gc,V$⁠. Bifurcation trees using the membrane voltage variable $V(1)$ were, therefore, constructed using a Poincaré section at the hypersurface $n(1)=0.003$ by scanning $gc,V∈[0,0.6]$ from left to right and vice versa [Fig. 2(b) violet and gray, respectively]. Changing the scanning direction of the bifurcation parameter unraveled that for $gc,V∈[0,0.19]$⁠, there is a coexistence between a synchronized bursting state and an additional dynamical regime. Numerical simulations for distinct coupling values from this interval showed that this dynamical regime could be periodic or aperiodic [Figs. 2(c) and 2(d), respectively]. This is also reflected in the complex fractal structure of the respective basins of attraction, as exemplified for $gc,V=0.15$ in Fig. 2(e). Given that the synchronized bursting state lies only on the diagonal of the basin [gray line in Fig. 2(e)], it is most probable to reach one of the other dynamical states in this parameter region. For $gc,V>0.19$⁠, however, only synchronized bursting between both oscillators was observed. This analysis, therefore, demonstrates that the physiological, synchronized bursting behavior is the dominant dynamical behavior over a wide coupling range in a population of identical $β$-cells that communicate via gap junctions.

For $ki=1$⁠, however, the presence of the potassium channel which has a decreased opening probability due to nonmonotonic voltage dependence [Fig. 1(b), Eqs. (6) and (7)] ensures that the equilibrium point is also stabilized such that isolated oscillators display bistability between a bursting and a silent state in the absence of coupling [Fig. 1(d)]. Thus, for initial conditions in the vicinity of the stable node, the physiological bursting behavior is not observed. We, therefore, refer to this silent regime as pathological. To unravel how such channel defects affect the global dynamics of the population, we used again a minimal model of two coupled $β$-cells that are individually characterized with bistability between the two states [Fig. 3(a)] and investigated the two possible scenarios: (i) the initial conditions of both oscillators are poised either in the bursting attractor or in the stable node attractor and (ii) a heterogeneous case where one of the oscillators is poised in the stable node, whereas the other oscillator is poised in the bursting attractor.

When the initial conditions of both oscillators are poised in the same attractor, the solution of the coupled system is trivial—the population behavior is synchronized and identical to the dynamics of the cells in isolation (results not shown). In contrast, for the heterogeneous scenario (ii), where both cells have the same bifurcation structure but populate distinct attractors such that the initial conditions for one cell are poised in the bursting, whereas for the other in the stable node, the bifurcation tree constructed analogously to the one in Fig. 2(b) has a complex structure for varying coupling strength $gc,V$ [Fig. 3(b), top]. To identify the existing dynamical solutions in this case, we complemented the bifurcation diagram with direct calculations from semirandom initial conditions [Fig. 3(b), bottom]. In particular, for each $gc,V$⁠, we calculated 500 time series for the system of two coupled $β$-cells such that the initial conditions for one of the oscillators were chosen exactly at the stable node $(V0(1),n0(1),S0(1))=(−49.084,0.0027105,0.19648)$⁠, whereas the initial conditions for the second oscillator were randomly chosen from a phase space cube in the interval $V0(2)∈(−60,−20)$⁠, $n0(2)∈(0,0.12)$⁠, and $S0(2)∈(0.17,0.2)$⁠. Although the size of this cube includes both the bursting as well as the stable node attractor in the uncoupled case, the basin of the bursting one is dominant. The full set of initial conditions covers the 6-dimensional phase space of the system with 3 degrees of freedom per oscillator densely enough such that one can detect all stable coexisting attractors with a significant basin of attraction. These direct calculations represent a different approach from the bifurcation analysis and are valid for large system sizes as well.²⁷ Therefore, the combination of both methods sheds light on the dynamics of the coupled heterogeneous $β$-cells model under the introduced pathological channel defects.

The bifurcation tree and the direct calculations plot [Fig. 3(b)] demonstrate that under weak coupling (⁠$gc,V<0.017$⁠), the heterogeneous system of two coupled $β$-cells where initial conditions differentially poise the oscillators in the two stable attractors, predominantly displays either bursting or complex, aperiodic dynamic. For the $gc,V$ subintervals for which bursting or aperiodic dynamics was observed, the probability for the coupled system to exhibit steady state solution was less than $10%$⁠. That the bursting or the aperiodic dynamics are the dominant regimes under weak coupling is also reflected in the corresponding basins of attraction estimated for the respective $gc,V$ values [Figs. 3(c) and 3(d)]. In contrast, for a large part of intermediate coupling strength interval (⁠$0.017<gc,V<0.08$⁠), the probability that the coupled system exhibits steady state dynamics was $∼100%$⁠, with less than $10%$ variability in small $gc,V$ ranges. The exemplary basin of attraction plot estimated for $gc,V=0.02$ also depicts the dominance of the stable equilibrium point, with only small lines of the basin of the bursting attractor [Fig. 3(e)]. This is surprising, given the small basin of attraction of the steady state for isolated cells [Fig. 1(d)]. The mechanism leading to this stabilization will be discussed elsewhere (manuscript in preparation).

A further increase in $gc,V$ continuously decreases the probability for the full system to exhibit steady state behavior. At $gc,V≈0.12$⁠, for example, the system visits the bursting attractor and the equilibrium point with equal probability. This implies that in the presence of noise, the full system can switch between bursting physiological and silent pathological states.

For high enough coupling strength $gc,V$⁠, however, the basin of attraction of the steady state is significantly decreased, leaving the bursting attractor as the most common solution of the system [Fig. 3(f)]. These results, therefore, demonstrate that a heterogeneous system of two coupled $β$-cells that have dysfunctional potassium channels and thereby possibility to manifest a stable steady state, predominantly exhibiting physiological bursting dynamics. This is true for weak and strong coupling strengths, even when the initial conditions of one of the oscillators correspond to its positioning in the stable node. Only for intermediate coupling values, the silent pathological state is the dominant dynamical regime, even for small steady state basin of attraction for single cells.

To model a physiologically more realistic scenario, we investigate next the dynamics of a minimal mixed system of two coupled cells, where only one of them displays bistability between bursting and a silent state, and thus a channel function dysregulation [Fig. 4(a)]. Thus, in contrast to the previous scenario, the coupled cells have a different bifurcation structure.

Bifurcation analysis showed that the mixed population does not display steady state dynamics for any $gc,V$ value in the given interval [Fig. 4(b)]. This is contrary to the previous case [Fig. 3(b)] where silent state was dominant for intermediate $gc,V$ values when bistability characterizes both cells, although the initial conditions of only one of them were poised in the stable node (coupled heterogeneous $β$-cells). A typical dynamics of the coupled mixed population model is mixed-mode oscillations, defined as complex oscillatory patterns where small-amplitude high-frequency oscillations alternate with large-amplitude long-period cycles.²⁸ Increasing the coupling strength between the two mixed oscillators induces transitions between different manifestations of the mixed-mode oscillations. For example, under low coupling strength, one oscillator displays small-amplitude, whereas the other large-amplitude bursting dynamics [Fig. 4(c)]. For intermediate coupling, different patterns of alternating between the small- and the large-amplitude bursting occur [Figs. 4(d)–4(f)], to finally converge to almost synchronized bursting of both oscillators under strong coupling [Fig. 4(g)]. The analysis, therefore, demonstrates that in a mixed population of two coupled cells, where only one of them has dysregulation of a potassium channel, pathological silent dynamics is not stabilized on the level of the population. This implies that deviations from the physiological bursting dynamics might be possible for a larger subpopulation size that has channels dysregulation.

To identify whether stabilizing pathological dynamics is a system-size effect, we investigated the dynamics of mixed populations of increasing sizes, starting with a minimal population of $N=3$ cells, where at least $(N/2+1)$ of them can exhibit bistability between silent and bursting dynamics. For simplicity, we initially assume global coupling between all cells. The numerical simulations demonstrated that, in this case, there is a threshold of coupling strength for which the silent state can be stabilized on the population level [Fig. 5(a)]. However, the probability to observe this pathological, silent dynamics in the mixed cell population is significantly smaller than the one in the heterogeneous population [Fig. 3(b)]. Direct calculations from semirandom initial conditions as in Fig. 3(b) showed that for the minimal system of $N=3$ cells, this probability is smaller than $1%$ [Fig. 5(b)]. Subsequent increase of the population size by a cell that can display bistability between a silent and bursting dynamics in turn lowered the coupling threshold for which the pathological, silent state was observed [Fig. 5(a), bottom]. The simulations also demonstrated that the probability to observe the pathological state increased with the increase of fraction of cells that have the channel defects. For example, for a mixed population of $N=12$ cells, where $11$ of them have dysfunctional channel, the probability of observing the silent dynamics increased to approximately $10%$⁠, almost over the full range of coupling strengths [Fig. 5(c)].

Physiological studies of $β$-cells activity in the islets of Langerhans have, however, shown that injection of a dye through a microelectrode in a single cell usually labels a small number of neighboring cells,^29,30 presumably by passing through the gap junctions. It has been, therefore, suggested that direct local coupling gives rise to the emergent dynamical behavior of $β$-cells.¹⁷ We have, therefore, repeated the above simulations when the nearest-neighbor coupling was considered, where each cell is coupled with 6 neighboring cells. Restricting the coupling to local interactions did not result in qualitatively different behavior than the one observed for global coupling. In this case, for example, for a mixed population of $N=12$ cells where 11 have a dysfunctional channel, the probability to observe the silent dynamics is slightly smaller than $10%$ [Fig. 5(d)], which was observed for the globally coupled case. Thus, despite the increase in the probability to observe the pathological state with increased system size, the population of mixed cells, being globally or locally coupled, can still cooperatively regulate the physiological bursting dynamics, even when the majority of the population has acquired dysregulation in the channel’s activity. This is contrary to the dynamics of the population where each cell has dysregulation, in which case a coupling interval exists for which the probability to observe the pathological dynamics is $∼100%$ [Fig. 3(b)].

When the number of cells that lack dysfunctional channels in the heterogeneous population is, however, increased, we find that the coupling threshold for which stabilization of the silent state is reached is subsequently increased.

Indeed, in accordance with the previous finding that at least $(N/2+1)$ cells must have dysfunctional channels for the silent state to be stabilized, the numerical simulations for globally coupled systems demonstrated that when two or three cells display only bursting dynamics, the minimal mixed population size for which the silent state can be stabilized is $N=5$ [Figs. 6(a) and 6(c)] or $7$ [Figs. 6(b) and 6(d)] cells, respectively. The stabilization threshold scales with the number of cells displaying only bursting dynamics, being the highest when the mixed population has the largest number of cells exhibiting solely bursting dynamics. Similar to Fig. 5(a), a stepwise increase of the population by a cell that displays bistability between the silent and the bursting state also decreases the coupling threshold for stabilization [Figs. 6(c) and 6(d)]. The local coupling did not affect the qualitative behavior of the mixed population in this case as well [Figs. 6(c) and 6(d)].

Intercellular communication dictates coordinated dynamics of cellular populations and thereby determines cellular identity in terms of their function. For example, regulated insulin secretion of pancreatic $β$-cells is driven by well-synchronized bursting activity of the population, where the secretion increases with the fraction of time the cells spend in the spiking state. In turn, disassociated $β$-cells also show increased electrical activity at elevated glucose concentrations; however, the activity varies from cell to cell such that some cells are electrically “silent,” whereas others are continuously “spiking” at any given glucose concentration.^10–12 Using transgenic mice where the potassium channels are functionally “knocked out” in $70%$ of the $β$-cells, it has been, however, demonstrated that the intact islets still show glucose-dependent insulin secretion,¹⁵ suggesting that the functional identity of $β$-cells can be dynamically maintained on a population level due to cell–cell coupling.

Here, we investigated the questions: how channel defects on the level of single cells that alter the intercellular communication can affect the global dynamics of the population and under which conditions the dynamics will deviate from the physiological one. An altered dynamics would, in turn, affect the cellular function and thereby their identity. The bifurcation analysis of minimal homogeneous or heterogeneous populations of $N=2$ coupled cells, where, respectively, either both cells have channel 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 pathological silent state (⁠$∼100%$⁠) for a given coupling range despite the small basin of attraction on a single cell level, the minimal heterogeneous population was characterized only with bursting dynamics over the full coupling range. These findings resemble previously observed cooperative effects in heterogeneous populations of coupled $β$-cells. In particular, it has been demonstrated that bursting on a population level can be restored due to cell–cell coupling, although single cells can only spike continuously or be silent.^9,31 Similar effects have been also observed for coupled van der Pol–FitzHugh–Nagumo units.³² Although in these cases the nature of the bistability between the limit cycle and the steady state was different than in the model used here, what is equivalent is that all of the cells had an equivalent bifurcation structure, where the heterogeneity was introduced by the parametric positioning of the system in one of the stable attractors.

Our findings, therefore, suggested that the pathological state can only be stabilized as a system-size effect in a mixed population of cells. In particular, we coupled cells with different bifurcation structures, cells that only exhibit bursting and cells that have bistability between bursting and silent state, and investigated the effect of the size of the subpopulation which has channel defects on the cooperative dynamics of the full system. For a mixed population of two cells, a coupling range was identified where the silent mode of operation was dominant. The numerical simulations of larger mixed populations, on the other hand, demonstrated that although the silent state could be indeed stabilized when the fraction of cells with channel dysfunctions is increased, the bursting dynamics still remained a dominant dynamical behavior of the population, with the prevalence of $∼80$%–$90%$⁠. Thus, the physiological population dynamics and the thereby associated cellular function, such as insulin production in the case of bursting dynamics of $β$-cells can be reliably maintained, even in the presence of dysfunctional channels on the level of single cells.

A.K. conceptualized the study. N.S. performed the numerical simulations and bifurcation analysis with the help of A.K. N.S. and A.K. interpreted the results and wrote the manuscript.

N.S. acknowledges the DAAD (Project No. 57440915) for the financial support of the scientific visit to MPI-MP in Dortmund, 2019.

The authors declare that they have no conflict of interest.