Cooperative regulation of cellular identity in systems with intercellular communication defects

The cooperative dynamics of cellular populations emerging from the underlying interactions determines cellular functions, and thereby their identity in tissues. Global deviations of this dynamics on the other hand reflects pathological conditions. However, how these conditions are stabilized from dis-regulation on the level of the single entities is still unclear. Using generic Hodgkin-Huxley type of models, we investigate here how channel defects in single beta-cells affect the intercellular communication, and whether the collective dynamics of the population can be thereby influenced. We identify a stark contrast between the dynamical behavior of homogeneous and heterogenous populations of cells. Whereas in the first case, when all of the cells acquired channel defects, a coupling interval can be found for which the pathological dynamics is stabilized and will be exhibited with ≈ 100% probability, in a heterogenous population, at least N/2 + 1 cells must have the channel defects in order for the pathological dynamics to be stabilizied. The numerical simulations demonstrated however that the 1 ar X iv :1 90 9. 02 94 4v 1 [ qbi o. C B ] 2 6 A ug 2 01 9 probability to exhibit the pathological dynamics in this case is maintained at < 10%, even when the population contains only a single cell that displays physiological dynamics. This results therefore suggest that the physiological bursting dynamics of a population of beta-cells is cooperatively regulated, even when single-cells acquire disfunctional channels that induce intercellular communication defects.


INTRODUCTION
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] to be executed. On the other hand, the collective behavior of homogenous 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-synchronised electrical activity (spiking and bursting) of primary beta-cells in the pancreatic islet is what precedes insulin secretion [4], [5], [6], 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-10s, whereas the duration of the bursting period varies from 0.2 to 5min. In contrast, isolated cells show electrical activity with very different time scales, do not burst and thereby do not secrete insulin effectively [7], [8]. In this case, the opening probability of the potassium channels is too small for the individual cells to present regular spiking signal that conversely only occurs in population of synchronized cells, thereby demonstrating that the intercellular communication in the organized islet structure is a necessary prerequisite for effective functioning of this coupled-cell system. 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 dis-regulation on the level of the interacting elements is however unclear. This question becomes even more prominent considering that it has been experimentally demonstrated that generically, the collective behavior of cells on the level of tissues is generally stable against mutations that occur in single cells, allowing them to retain their identity [9].
Here we use numerical simulations to tackle how a change in the global dynamics of a coupled system can be established by dis-regulation 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 beta-cells [10], [11], [12]. Additionally, we also use a modified version of this model where a potassium-like ion channel with decreased opening probability is introduced [13] 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 disfunction in a form of blocking or inactivation. We use these models to investigate how deviation from the physiological dynamics of the population of coupled beta-cells is triggered depending on the number of cells acquiring disregulation in channel activity. Three cases are therefore analyzed: (1) dynamics of two coupled beta-cells under physiological conditions (bursting activity), (2) dynamical behavior of two coupled cells that have acquired the disfunction (stable steady state), but one or both of them do not manifest it in isolated cells and (3) dependence of the global population dynamics on number of cells with channel defects.
We found that the dynamics of a population of two coupled cells, each characterized with bistability between silent and bursting state, is either equivalent to the same attractor in which both cells are poised when in isolation, or dictated by the physiological bursting dynamics when the uncoupled cells are initially poised in different attractors. In contrast, the bursting dynamics of a heterogeneous 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 N/2 + 1 cells have the channel disfunction. 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 heterogeneous population of beta-like cells reflected through their bursting dynamics is cooperatively regulated, even when a large sub-fraction of this population acquired dis-functional channels that induce intercellular communication defects.

Single beta-cell models of bursting and silent dynamics
The electrical activity of pancreatic beta-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 [4]. These channels not only regulate membrane potential, ion homeostasis and electrical signaling [14], but a dis-regulation in their activity has been recently also implicated in tumorogenesis and tumor progression [15], [16], [17], [18].
To depict the bursting dynamics of pancreatic beta-cells as observed under physiological conditions (Fig. 1a), as well as to simultaneously model silent dynamics caused upon ion channel dis-regulation (Fig. 1b), we use a model based on the Hodgkin-Huxley formalism [11], [12] with modifications proposed by Stankevich et al. [13]: where V (i) represents the membrane potential of the i-th cell, the functions I Ca (V (i) ), I K (V (i) , n (i) ), I S (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 I K2 (V (i) ) on the other hand defines an additional voltagedependent potassium current. It has been previously demonstrated that a single oscillator, in absence of I K2 (when k (i) = 0, Fig. 1a), is characterized with a bursting attractor that is born via a Hopf bifurcation for V S = −44.7mV . At this parameter value, the equilibrium point looses its stability [19], [20]. Calculating the eigenvalues in this case shows that the equilibrium point is an unstable node (λ 1 , λ 2 , λ 3 ) = (23.138659, −41.832952, 0.089497). Although the bursting attractor is born in the vicinity of the equilibrium point, this point moves away from the bursting attractor for increase in V S . Even more, the two dimensional projection of the phase portrait together with the fast and slow manifolds when V S = −35mV ( Fig. 1c) demonstrates that the periodic trajectories do not intersect the neighborhood of the steady state and the bursting attractor terminates in a homoclininc bifurcation as the trajectory hits the slow manifold [19], [20].
In the presence of Fig. 1b) 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. 1d, corresponding eigenvalues (λ 1 , λ 2 , λ 3 ) = (−37.58, −13.6, −0.24)), estimated for g K2 = 0.2 at the equilibrium point (V 0 ) = (-49.084, 0.0027105, 0.19648) [13]. I K2 (V (i) ) is given in the form: where represents the opening probability of this channel. In contrast to the other potassium channel (Eqs. 3, 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 disfunction, as for instance blocking of the potassium channel or its inactivation [21], and thereby a stable silent state emerges. 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 [13]. Generally, the modified model (k (i) = 0) depicts bistability between physiological, bursting beta-cell dynamics and a pathological, silent state dynamics (Fig. 1d).

Dynamics of two coupled beta-cells under physiological conditions
To investigate the global dynamics of a population of beta-cells, we introduced electrical coupling by adding the following gap-junctional coupling term to the equation for V (i) that describes bi-directional transport of ions between the cells [6]: g c,V is the coupling conductance (coupling strength). We consider only electrical coupling of the cells (Fig. 2a) by adding the following gap-junctional coupling terms to the equations for V in (Eqs. 1). The sum is taken over the whole population, assuming global coupling.
As already noted, if k i = 0 and in absence of coupling, the isolated oscillators display a single stable state -bursting dynamics (Fig. 1c). To infer the effect of coupling on the system's dynamics, we investigated next the bifurcation structure of a minimal model of two coupled cells as a function of the coupling strength g c,V . Bifurcation trees using the membrane voltage variable V (1) were therefore constructed using a Poincare section at the hypersurface n (1) = 0.003 by scanning g c,V ∈ [0, 0.6] from left to right and vice versa (Fig. 2b violet and grey, respectively). Changing the scanning direction of the bifurcation parameter unravelled that for g c,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 chaotic (Fig. 2c,d respectively). This is also reflected in the complex fractal structure of the respective basins of attraction, as exemplified for g c,V = 0.15 in Fig. 2e. Given that the synchronized bursting state lies only on the diagonal of the basin (grey line in Fig. 2e), it is most probable to reach one of the other dynamical states in this parameter region. For g c,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 beta-cells that communicate via gap junctions.

Dynamics of two coupled beta-cells that contain dis-functional potassium channels
For k i = 1 however, the presence of the potassium channel which has a decreased opening probability (Fig. 1b, Eqs. 6,7) ensures that the equilibrium point is also stabilized, such that isolated oscillators display bistability between a bursting and a silent state in absence of coupling (Fig. 1d). 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 beta-cells that are individually characterized with bistability between the two states ( Fig. 3a), and investigated the two possible scenarios: i) both oscillators in the absence of coupling are poised either in the bursting attractor or in the stable node attractor, and ii) one of the oscillators is poised in in the stable node, whereas the other oscillator is poised in the bursting attractor. When 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 scenario ii), when isolated cells 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. 2b has a complex structure for varying coupling strength g c,V (Fig. 3b, top). To identify the existing dynamical solutions in this case, we complemented the bifurcation diagram with direct calculations from semirandom initial conditions (Fig. 3b, bottom). In particular, for each g c,V we calculated 500 time series for the system of two coupled beta-cells such that the initial conditions for one of the oscillators were chosen exactly at the stable node (V 0 ∈ (0, 0.12) and S (2) 0 ∈ (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 6dimensional 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 [22]. Therefore, the combination of both methods sheds light on the dynamics of the coupled beta-cells model under the introduced pathological channel defects.
The bifurcation tree and the direct calculations plot (Fig. 3b) demonstrate that under weak coupling (g c,V < 0.017), the system of two coupled betacells where initial conditions differentially poise the oscillators in the two stable attractors, predominantly displays either bursting or complex, chaotic dynamic. For the g c,V sub-intervals for which bursting or chaotic dynamics was observed, the probability for the coupled system to exhibit steady state solution was less than 10%. That the bursting or the chaotic dynamics are the dominant regimes under weak coupling is also reflected in the corresponding basins of attraction estimated for the respective g c,V values (Fig. 3c,d). In contrast, for a large part of intermediate coupling strength interval (0.017 < g c,V < 0.08), the probability that the coupled system exhibits steady state dynamics was ≈ 100%, with less than 10% variability in small 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, with only small lines of the basin of the bursting attractor (Fig. 3e).
A further increase in g c,V continuously decreases the probability for the full system to exhibit steady state behavior. At g c,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 state.
For high enough coupling strength g c,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. 3f). These results therefore demonstrate that a system of two coupled beta-cells that have disfunctional potassium channels and thereby possibility to manifest a stable steady state, predominantly exhibits 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.

Dynamical behavior of a minimal mixed population model
To model 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 dis-regulation in a form of inactivation or blocking of a potassium channel (Fig. 4a). Bifurcation analysis showed that the mixed population does not display steady state dynamics over the whole g c,V interval (Fig.4b). This is contrary to the previous case (Fig. 3b) where silent state was dominant for intermediate g c,V values when bistability characterizes both cells, although the initial conditions of only one of them were poised in the stable node. 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 [23]. 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. 4c). For intermediate coupling, different patterns of alternating between the small-and the large-amplitude bursting occur (Fig. 4d-f), to finally converge to almost synchronized bursting of both oscillators under strong coupling (Fig. 4g). The analysis therefore demonstrates that in a mixed population of two coupled cells, where only one of them has dis-regulation of a potassium channel, pathological silent dynamics is not observed on the level of the population. This implies that deviations from the physiological bursting dynamics might be possible for a larger sub-population size that has channels dis-regulation.

Effect of intercellular communication defects on the dynamics of heterogeneous beta-cells population
To identify whether stabilizing pathological dynamics is a system-size effect, we investigated the dynamics of heterogeneous populations of increasing sizes, starting with a minimal population of N = 3 cells, where (N/2 + 1) of them can exhibit bistability between silent and bursting dynamics. 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. 5a). However, the probability to observe this pathological, silent dynamics in the heterogeneous cell population is significantly smaller than the one in the homogenous population. Direct calculations from semirandom initial conditions as in Fig. 3b showed that for the minimal system of N = 3 cells, this probability is smaller than 1% (Fig. 5b). 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. 5a, 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 population of N = 12 cells, where 11 of them have dis-functional channel increased to approximately 10%, almost over the full range of coupling strengths (Fig. 5c). Despite this increase in the probability to observe the pathological state, the population of heterogeneous cells can still cooperatively regulate the physiological bursting dynamics, even when the majority of the population has acquired dis-regulation in the channel's activity. This is contrary to a homogeneous population of cells with dis-regulation, where a coupling interval exists for which the probability to observe the pathological dynamics is ≈ 100% (Fig. 3b). When the number of cells that lack dis-functional 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 dis-functional channels for the silent state to be stabilized, the numerical simulations demonstrated that when 2 or 3 cells display only bursting dynamics, the minimal heterogeneous population size for which the silent state can be stabilized is N = 5 (Fig. 6a,c) or 7 (Fig. 6b,c) cells, respectively. The stabilization threshold scales with the number of cells displaying only bursting dynamics, being highest when the heterogeneous population has the largest number of cells exhibiting solely bursting dynamics. Similarly to Fig. 5a, a step-wise increase of the population by a cell that displays bistability between the silent and the bursting state also decreases the coupling threshold for stabilization (Fig. 6c).

Discussion
Intercellular communication dictates coordinated dynamics of cellular populations, and thereby determines cellular identity in terms of their function. For example, insulin secretion of pancreatic beta-cells is preceded 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 contrast, isolated cells do not burst and thereby do not secrete insulin effectively. Such dynamical behavior of single cells has been related to the decreased opening probability of the potassium channel, resulting in silent dynamics. We considered here the case when cells can exhibit silent dynamics, even when coupled in a beta-cells population. In particular, we addressed the question 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 a 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, the minimal heterogeneous population was only characterized with bursting dynamics over the full coupling range. This suggested that the pathological state can only be stabilized as a system-size effect in a heterogeneous population of cells, when a significant proportion of cells in the population have channels with decreased probability of opening. The numerical simulations of larger heterogeneous populations on the other hand demonstrated that although the silent state could be indeed stabilized when the fraction of cells with channel dis-functions is increased, the bursting dynamics still remained a dominant dynamical behavior of the population, with 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 beta-cells can be reliably maintained, even in the presence of dis-functional channels on the level of single cells.