Entropy, Ergodicity and Stem Cell Multipotency

Populations of mammalian stem cells commonly exhibit considerable cell-cell variability. However, the functional role of this diversity is unclear. Here, we analyze expression fluctuations of the stem cell surface marker Sca1 in mouse hematopoietic progenitor cells using a simple stochastic model and find that the observed dynamics naturally lie close to a critical state, thereby producing a diverse population that is able to respond rapidly to environmental changes. We propose an information-theoretic interpretation of these results that views cellular multipotency as an instance of maximum entropy statistical inference.

Populations of mammalian progenitor cells commonly exhibit considerable cell-cell variability. However, the functional role of this diversity is unclear. Here, we analyze the dynamics of the stem cell surface marker Sca1 in mouse hematopoietic progenitor cells using a simple stochastic model and propose that multipotency is an example of maximum entropy statistical inference: in the absence of defined instructions, individual cells remain maximally non-committal with respect to their molecular identity, thereby generating a maximally diverse population that is able to respond optimally to a wide range of unforeseen future environmental changes. These results indicate that cell-cell variability is central to collective decision-making by mammalian progenitor cells and suggest that regulation of cellular dynamics is exerted at the population, rather than individual cell, level. Clonal populations of unicellular organisms often exhibit phenotypic diversity, which confers selective advantage under adverse environmental conditions. Wellknown examples include antibiotic bacterial persistence, the lysis-lysogeny switch of λ-phage, competence development and sporulation of B. subtilis, and lactose uptake by E. coli [1]. The ubiquity of this phenomenon indicates that it is a generic, evolvable, mechanism that facilitates collective cellular decision-making by enabling robust, rapid responses to diverse environmental changes. Recently, stochastic fluctuations in expression of important marker proteins have been seen to generate functional diversity within multipotent mammalian stem cell populations, suggesting a similar role for cell-cell variability in cellular decision-making in higher organisms [2]. These observations have motivated speculation that functional multipotency (the ability to differentiate along a number of distinct cellular lineages) is a collective property of stem and progenitor cell populations, reflective of fitness constraints imposed at the population, rather than individual cell, level [3]. This perspective is appealing since such regulated cell-cell variability in principle allows a cellular population to remain primed to respond quickly to a range of different differentiation cues while remaining robust to cell loss. However, convincing demonstrations of the potency of individual stem cells appear to argue strongly against such a collective view (for example, single long-term repopulating hematopoietic stem cells are able to fully reconstitute the blood system of lethally irradiated adult mice and small numbers of pluripotent stem cells are able to rescue development of genetically compromised embryos [4]). Thus, it is still unclear how population-level and cell-intrinsic regulatory programs interact to control mammalian stem and progenitor cell dynamics.
Here we propose a theoretical framework that reconciles these disparate observations, which views cellular multipotency as an instance of maximum entropy statis-tical inference. In this view, individual cells satisfy any minimal regulatory constraints imposed upon them (such as basic metabolic requirements, etc.) yet, in the absence of defined instructions, are maximally noncommittal with respect to their remaining molecular identity, thereby generating a diverse population that is able to respond optimally to a range of unforeseen future environmental changes. Thus, rather than viewing the stem/progenitor state as an attractor of the underlying molecular regulatory dynamics (i.e. associating cellular identities with well-defined, stable, patterns of gene expression -a common modeling assumption, that has received some experimental validation for differentiated cell types [5]) we suggest that individual stem/progenitor cells are characterized by fundamental uncertainty in their molecular state and their populations exhibit variability in accordance with this intrinsic uncertainty. However, since this model exchanges the attractor hypothesis at the single cell level for an ergodicity assumption for the underlying stochastic processes, each individual cell has the latent potential to assume every identity within the population, and thereby retains the regenerative capacity of the entire population. As this view is fundamentally stochastic, its corollary is that regulation of cellular dynamics is exerted at the level of probabilities (i.e. at the population level), rather than at the individual cell level.
In order to illustrate this perspective we consider here the non-equilibrium expression dynamics of the stem cell surface marker Sca1 (stem cell antigen 1) in populations of multipotent eml mouse hematopoietic progenitor cells. It has previously been shown that Sca1 levels fluctuate stochastically in eml cells in culture, with extrinsic 'transcriptome-wide' noise driving transitions between Sca1 high and Sca1 low states, which transiently prime individual cells for erythroid and myeloid differentiation respectively and generate a characteristically bimodal Sca1 expression distribution within the population at equilibrium (see Fig. 1, bottom panel and Ref. [6]). However, the underlying mechanisms by which these stochastic fluctuations are regulated are not known. In the absence of this knowledge we assume here that the intracellular dynamics of the Sca1 expression level z(t) are described by a generic stochastic differential equation: where ξ(t) is a standard one-dimensional white noise process [ ξ(t) = 0 and ξ(t)ξ(s) = δ(t − s)] and d(z) accounts for fluctuations in Sca1 levels due to both intrinsic sources (i.e. noise in the molecular processes involved in Sca1 production/decay, such as transcription, translation, translocation and degradation etc.) and extrinsic sources (i.e. fluctuations in upstream regulators and uncontrolled environmental noise). Rather than model Sca1 levels directly it is convenient to introduce a reaction coordinate x(z) such that the Fokker-Planck equation for the probability density ρ(x, t) has the form with scalar potential ψ(x) and diffusion coefficient σ. Such a transformation, which maps the original dynamics to those of a Brownian particle in a one-dimensional potential field, may always be achieved by application of Itō's lemma. The stationary solution of Eq. (1) is the Boltzmann-Gibbs distribution This solution exists so long as ψ(x) grows sufficiently rapidly as |x| → ∞ that the partition function Z remains finite. In this case, the dynamics are ergodic and the free energy where E(ρ) and S(ρ) are the energy and entropy functionals respectively, is a Lyapunov functional for the dynamics. Thus, in order to model Sca1 dynamics phenomenologically we need only chose an appropriate reaction coordinate x and form for the potential ψ(x).
Since protein fluctuations typically scale linearly with expression level, a natural choice for the reaction coordinate is x = log z, as has been taken elsewhere [7]. A parsimonious choice for the potential, which allows for observed bimodality without introducing large numbers of free parameters, is: where h is a positive even integer [8]. Intuitively, this is a simple model of a positive-feedback based bistable switch of the kind that commonly regulate cell fate decisions [9]. The stationary distribution p ∞ (x) is then characterized by four nonnegative dimensionless parameters: h, α = α 0 /α 1 , γ = βK/α 1 and θ = σβ 2 /α 2 1 . Along with the diffusion coefficient σ, which sets the timescale for the dynamics, these parameters fully describe the behavior of Eq. (1), which constitutes our model of Sca1 dynamics.
Parameter estimates were obtained by model fitting using maximum likelihood estimation to evolving Sca1 expression distributions obtained experimentally using flow cytometry starting from pre-selected populations of Sca1 low, mid, and high expressing cells as they converge to equilibrium in culture over a period of 18 days (obtained in Ref. [6]). Despite the simplicity of this model, an excellent agreement with the experimental data was observed for the temporal dynamics of Sca1 reconstitution from all three pre-selected subpopulations, using the It has previously been argued, based upon analysis of changing proportions of cells in the Sca1 high and low states, that the observed dynamics are characterized by slow 'sigmoidal' relaxation towards the stationary state [6]. Since a constant probability flux across a barrier naturally leads to exponential relaxation, it was suggested that these dynamics indicate deviation from expected first-order kinetics, possibly due to regulation of Sca1 fluctuations by cell-cell communication or autocrine signaling. However, it is apparent that such recourse is not needed since the experimental system is, by construction, initially far from equilibrium, and therefore far from the regime in which first-order kinetics apply. Rather, in accordance with standard reaction-rate theory, the dynamics are characterized by an initial transient period during which local equilibrium is first established within each potential well, before transitions between wells occur [10]. Examination of the free energy (which is a natural way to assess convergence to equilibrium [11]) shows that this separation of timescales naturally generates the observed convergence dynamics without the need to include additional regulatory mechanisms in the model (see Fig. 3). These results indicate that the observed Sca1 expression dynamics are well described by a simple ergodic process in which individual cells behave independently to one another with respect to Sca1 fluctuations.
This ergodic property is useful since it allows inference of the behavior of individual cells from the population dynamics. While stochastic excursions into Sca1 high and low states have previously been seen to transiently confer different lineage biases to individual progenitor cells in culture, the timescales upon which these switches occur at the single cell level are not known. Thus, the distribution of first passage times (FPTs) out of the Sca1 high and low states are of particular interest. Defining the ranges of Sca1 high and low expression as H = (x 0 , ∞) and L = (−∞, x 0 ) respectively, where x 0 is the intermediate maxima in ψ(x), the FTP T (x) out of S for a cell initially at x ∈ S (where S ∈ {L, H}) may be obtained from the backward Fokker-Planck equation associated with Eq. (1). Denoting G(x, t) = P (T (x) ≥ t) we solve: with initial conditions G(x, 0) = 1 for x ∈ S and boundary conditions G(x 0 , t) = G(±∞, t) = 0, from which the FPT distributions F S (x, t) = P (T (x) < t) = −∂G/∂t for S ∈ {L, H} may be obtained. Conventionally, the FPT distribution F S (x, t) is evaluated from the local minima x S of ψ(x) in S, since this is the state of highest probability. Alternatively, we can weight each initial position within S according to the probability that the cell is at this position at equilibrium. We thus define the expected FPT distribution with respect to the Gibbs measure, 1] is the weight of the population in S. Numerical approximations to F S (x S , t) and F S (t) are shown in Fig. 2. These distributions yield mean FPTs of 62/63 hours for the low state and 1560/1496 hours for the high state using F S (x S , t) and F S (t) respectively. It is significant that these timescales are substantially longer than the eml cell cycle time (approx. 10 -14 hours [12]), since they suggest that Sca1 switching is not simply a consequence of the cell-cycle. Rather, by setting the expected length of time that a pair of cells initially at the same position (e.g. daughter cells from the same cell division) will forget their common origin -and therefore the expected length of time that their identities will be coupled -Sca1 switching appears to encode an elementary form of epigenetic memory that endows individual cells with a transient functional identity.
Since the rate of switching is considerably slower than the rate of cell division it allows the formation of robust communities of cells that maintain the same characteristics though many divisions, and are therefore able to adopt a temporarily stable functional phenotype. Yet, by allowing slow mixing between the communities Sca1 fluctuations also safeguard long-term cell-cell variability and ensure that a robustly heterogeneous population, able to respond to a range of environmental challenges and resilient to removal of cellular sub-types, is maintained.
These results indicate that regulated fluctuations in Sca1 levels may be an intrinsic feature of eml cells in culture since they provide a mechanism by which the population hedges its bets against unforeseen future environmental challenges and thereby retains the capacity to differentiate along both erythroid and/or myeloid lineages as required. If this is the case, then it is natural to ask if  Figure 3. Convergence to equilibrium with respect to the free energy (left). Exponential convergence was observed from all three initial conditions for large time, in accordance with Eq. (1). Entropy of the stationary distribution over the αγplane (right). The empirically observed distribution is marked with a magenta cross and the maximum entropy distribution is marked with a green circle. The empirical distribution lies in the small region of the αγ-plane that is both close to critical and in which the resulting distribution is of high entropy. The colorbar shows percentiles (color online).
the experimentally observed stationary Sca1 distribution is optimal for this purpose; that is, if it is maximally variable in some appropriately defined way. To investigate this, we first note that since the Boltzmann-Gibbs density p ∞ (x) is the minimizer of the free energy F (ρ), it is the most non-committal way to assign probabilities subject to the constrains imposed upon the dynamics by ψ(x) (i.e. it may be seen as an instance of maximum entropy statistical inference) [13]. As each choice of model parameters defines a different potential, which places different constraints upon the dynamics, we may assess the extent to which Sca1 fluctuations optimize population diversity by determining how close the observed equilibrium Sca1 distribution is to the maximum entropy distribution in the ensemble P = {p ∞ (x) | h ∈ 2Z + , α, γ ∈ R + } [14].
Since the coefficient h is, informally, a measure of the sensitivity of the underlying switch to input stimulus, it primarily affects the curvature of the potential around the local minima x 0 (where present) and does not have a strong effect on the entropy. However, by governing a codimension 2 bifurcation that determines whether the underlying switch is in a monostable or bistable state, α and γ can affect the entropy of ρ ∞ (x) considerably. Fig. 3 shows how the entropy of p ∞ (x) varies over the biologically relevant bistable region of the αγ-plane. It can be seen that the point estimate for the entropy of the experimentally observed Sca1 distribution is remarkably close to the maximum entropy distribution in P. Furthermore, it is notable that the observed values of α and γ are also remarkably close to critical. It has long been suggested that such criticality may emerge naturally in biological systems via self-organizing evolutionary processes without the need for fine-tuning (i.e. as an attractor of the evolutionary dynamics), since critical states provide the dual benefits of stability and adaptability [15]. Here, close proximity to criticality is beneficial since it ensures that a heterogeneous population is produced, yet, by minimizing the height of the potential barrier between the Sca1 high and low states, mixing between these populations occurs on a biologically feasible time-scale. These results suggest that Sca1 levels are regulated by a critically self-organizing evolutionary process in which the fitness of the population depends upon a trade-off between robustness of cell-cell variability and the ability to respond rapidly to environmental changes. In summary, our results suggest that variability within mammalian progenitor cell populations has a central role in cellular decision-making and may be seen as an instance of maximum entropy statistical inference. The full mechanisms of cellular collective decision-making in vivo are undoubtedly more complex, and certainly involve complex interactions between different types of stem, progenitor and differentiated cells [16]. The collective decision-making ability of such heterogeneous populations, both in culture and in vivo, remain to be determined.