Tipping the Balance: A Criticality Perspective

Cell populations are often characterised by phenotypic heterogeneity in the form of two distinct subpopulations. We consider a model of tumour cells consisting of two subpopulations: non-cancer promoting (NCP) and cancer-promoting (CP). Under steady state conditions, the model has similarities with a well-known model of population genetics which exhibits a purely noise-induced transition from unimodality to bimodality at a critical value of the noise intensity σ2. The noise is associated with the parameter λ representing the system-environment coupling. In the case of the tumour model, λ has a natural interpretation in terms of the tissue microenvironment which has considerable influence on the phenotypic composition of the tumour. Oncogenic transformations give rise to considerable fluctuations in the parameter. We compute the λ−σ2 phase diagram in a stochastic setting, drawing analogies between bifurcations and phase transitions. In the region of bimodality, a transition from a state of balance to a state of dominance, in terms of the competing subpopulations, occurs at λ = 0. Away from this point, the NCP (CP) subpopulation becomes dominant as λ changes towards positive (negative) values. The variance of the steady state probability density function as well as two entropic measures provide characteristic signatures at the transition point.


Introduction
The generation of heterogeneity in a population is often a consequence of the underlying stochastic nonlinear dynamics.Cell biology offers a large number of examples of population heterogeneity characterized by the coexistence of multiple subpopulations with distinct phenotypic traits [1,2,3].The frequently observed case is that of bimodality with two distinct subpopulations defining the heterogeneity.A well-understood physical basis of such bimodality is as follows.A single positive feedback loop or multiple loops governing the underlying deterministic dynamics create the potential for bistabiity in a specific parameter regime [4].In the state space defined by the concentrations of the key dynamical variables, the two stable steady states are separated by an unstable steady state.The corresponding landscape picture is that of two valleys separated by a hill.The minima of the valleys define the stable steady states, the attractors of the dynamics, and the top of the hill the unstable steady state.The stochastic component of the dynamics generates fluctuations (noise) in the magnitudes of the dynamical variables opening up the possibility of fluctuation-driven transitions from one attractor in the state space to the other.In the landscape analogy, a ball, left to itself, resides at the bottom of one of the valleys.Noisy dynamics, in the form of, say, random kicks imparted to the ball, give rise to a finite probability that the ball crosses the top of the hill and reaches the other valley.The stochastic dynamics smear out the steady states in the form of a steady state probability distribution with two distinct peaks, the case of bimodality.In this case, bistable deterministic dynamics are essential for the observation of bimodality in the stochastic case.In an alternative scenario, the bimodality is purely noise-induced with the deterministic dynamics yielding single stable steady states in the full parameter regime, even in the presence of positive feedback loops.In cell biology, there is now experimental evidence of noise-induced bimodality in the key cellular phenomenon of gene expression [5].One of the earliest examples of purely noise-induced transitions, resulting in bimodality, relates to a model of population genetics, referred to as the genetic model [6,7].In this model, the deterministic dynamics are governed by the rate equation   = 0.5 −  + (1 − ) (1) in which  represents a state variable, with   [0,1], and  is a parameter describing the coupling of the system to the environment with   [−∞, +∞].The model has a realistic interpretation in a population genetics context as well as in the case of a scheme of chemical reactions [7].From (1), the steady state is given by A globally stable unique steady state exists for each value of the external parameter .In the presence of a fluctuating environment, the coupling parameter  acquires the form,  →  +  , where ξ represents a white noise and  2 is the noise intensity.The dynamics now acquire a stochastic character requiring analysis in the framework of stochastic formalisms.The noise intensity serves as another parameter in the system, besides , and one can show that a noiseinduced transition to a bimodal steady state probability distribution occurs for  2 >   2 , a critical noise intensity.One can identify a critical point at   = 0,   2 = 2,   = 0.5.For  2 <   2 , the probability distribution has a single peak centred on   = 0.5.The single peak splits into a pair of peaks above the critical point.As in the case of equilibrium critical behavior, the noiseinduced nonequiibrium transition at the critical point is characterized by a set of critical exponents.The critical behavior has been shown to belong to the mean-field Ising universality class [7].For  ≠ 0, the critical value of the noise intensity, above which bimodality is observed, becomes a function of .
A recent study has shown that biochemical models with positive feedback also exhibit meanfield Ising-like critical behavior [8].The transition to a bimodal distribution is, however, not purely noise-induced if the underlying deterministic dynamics give rise to bistability in a specific parameter regime.In general dynamical models, the dynamics undergo regime changes at the bifurcation points which are analogous to phase transition points in equilibrium models [9,10,11].The bifurcations bringing about discontinuous (continuous) changes at the bifurcation points are akin to first-order (continuous/critical-point) phase transitions.In dynamical systems, the transition at a saddle-node (SN) bifurcation is discontinuous whereas that occurring at a supercritical pitchfork (SP) bifurcation is continuous.The latter bifurcation point has the character of a critical point.The mapping between the mean-field Ising and the nonequilibrium models rests on the equivalent forms of the equation of state and the steady state equation in the vicinity of the respective critical points.The effective thermodynamic quantities (temperature, magnetic field and magnetization) of the equilibrium model are expressed in terms of the effective biochemical parameters of the nonequilibrium model [8,10].The utility of such mappings has been demonstrated in single cell experiments on T cell signaling [8].
A prominent example of population heterogeneity is that of intra-tumour heterogeneity (ITH) in terms of both genotypic and phenotypic variability [12,13].Genotypic changes in the form of mutations are responsible for the formation of tumours, an abnormal conglomeration of cells and tissues.A tumour further harbours subpopulations with heterogeneity at both the genotypic and phenotypic levels.The first type of heterogeneity has been extensively studied whereas the interest in analyzing the second type of heterogeneity is of more recent origin.The study of phenotypic heterogeneity opens up the possibiity for developing targeted therapies against the subpopulations with a dominant role in the growth and spread of cancer.
Phenotypic diversity can have both genotypic and non-genetic origins.There is now considerable experimental evidence that the cells in a tumour are broadly of two types: cancer stem cells (CSCs) and non-cancer stem cells (NCSCs) [12].The CSCs arise when normal stem cells undergo mutations.The two subpopulations have considerable phenotypic differences, e.g., the proliferative potential of the CSC is much higher than that of the NCSC.The protein HER2 ( human epithelial growth factor receptor 2) promotes the quick growth of cancer cells in a type of breast cancer [14,15].The tumour population has two phenotypically distinct subpopulations designated as HER2+ and HER2-.The former (latter) subpopuation has a higher (lower) HER2 protein level in a cell than the level prevailing in a normal cell.The HER2+ subpopulation has a higher growth rate and tends to spread (metastatize) at a faster rate.The dominance of a cancer-promoting subpopulation in a tumour tilts the balance towards the eventual fate, cancer.
The onset of cancer involves a transition from the healthy to the pre-disease state followed by a sudden deterioration to the disease state with its characteristic clinical symptoms [16,17].In theoretical models, the sudden regime shift occurs via a saddle-node bifurcation so that the transition has the character of a first-order phase transition.The variation of an appropriate parameter generates a line of first order phase transitions terminating at a critical point in the phase diagram.A number of quantitative measures have been proposed to provide the early warning signals of an imminent transition to the pre-disease state [18,19,20].Before this state is reached, one can reverse the progressive deterioration by adopting suitable measures like drug treatment.The sudden regime shift to the disease state brings in irreversibility due to the existence of hysteresis in the response curve.In this paper, we explore the critical-like behavior in a model of cancer in which the ITH is defined in terms of two distinct subpopulations.In Section 2, we introduce the model and show that it can be mapped onto the population genetics model with the dynamics as described in (1).The stochastic dynamics of the model are described using the formalism of the Fokker-Planck equation (FPE).The existence of a purely noise induced transition to bimodality at a critical point is demonstrated and the critical exponents are identified as those of the mean-field Ising type.In Section 3, the phase diagram in the  −  2 plane is computed with the phase boundaries separating regions of unimodality from that of bimodality.The computation is based on stochastic considerations arising from the formalism of the FPE.The phase diagram is interpreted in terms of the sudden regime shifts to the disease state.The corresponding hysteresis curve is also exhibited.In Section 4, the transition from the regime in which both the subpopulations have equal probability for dominance to a regime in which one of the subpopulations is dominant, is investigated.The signatures of the transition, which tips the balance towards a specific subpopulation, are obtained via some statistical measures.Section 5 contains some concluding remarks.

Model of Tumour Heterogeneity 1.1 Similarity with Population Genetics Model
We consider a tumour population the heterogeneity of which is defined in terms of two distinct subpopulations, say, HER2+ and HER2-subpopulations or cancer-promoting (CP) and noncancer-promoting (NCP) subpopulations.Let  1 () and  2 () be the sizes (size is given by number of cells) of the NCP and CP subpopulations respectively at time  with the total population size () =  1 () +  2 ().Each cell in a subpopulation can undergo a symmetric cell division to yield a pair of daughter cells belonging to the same subpopulation resulting in the growth of the subpopulation.Let  1 and  2 be the growth rate constants associated with the NCP and CP subpopulations respectively.These rate constants define effective growth rates when the cell death due to apoptosis is included in the rates.A cell in a subpopulation can also divide asymmetrically with one daughter cell belonging to the parent subpopulation and its companion to the other subpopulation.Through this process of interconversion, the other subpopulation increases in size [14] whereas the size of the parent subpopulation remains the same.The rate constants of these processes are  1 (asymmetric division of a cell in the NCP subpopulation) and  2 (asymmetric division of a cell in the CP subpopulation).The dynamical equations capturing the above described processes are [15]: The symmetric division event giving rise to two daughter cells belonging to the other subpopulation are not considered explicitly as such events are of rare occurrence [15].The dynamical equations in (3) and ( 4) are similar in form to those appearing in a number of studies dealing with the problem of phenotypic heterogeneity [15,21,22,23].We now define  1 () and  2 () to be the fractions of the total population belonging to the NCP and CP subpopulations respectively at time  with  1 () =  1 () and  2 () = 1 −  1 ().Expressing  2 () in terms of  1 (), (3) and ( 4) can be combined into a single equation for  1 , namely where   =  1 −  2 and the rate constants  1 and  2 are taken to be equal [15],  1 =  2 =   .The theoretical results, for the evolution of the subpopulation fractions  1 and  2 as a function of time and with different initial conditions, fit the experimental data quite well [15] conferring validity on the tumour model.Through a simple rescaling of the time,   = 2  , and introducing the dimensionless parameter,  =   2  , ( 5) is rewritten in the form We point out that the dynamical equation ( 6) has the same form as that of the population genetics model [6,7].One can draw parallels between the two models by noting that the population genetics model considers a constant population of haploid individuals in which two alleles, with frequencies  and 1 −  in the population, compete for a specific gene locus.The frequencies undergo changes due to two mechanisms: mutation between the two alleles and natural selection.The parameter β in ( 1) is a selection coefficient with dependence on the state of the environment.A positive (negative) value of β favours the allele with frequency −  in (1) represents the change in frequency brought about by random mutations.In the cancer model, described by ( 3) and ( 4), the two mechanisms for the change in frequency are the two-way transitions between the phenotypes, analogous to mutations, and the unequal growth rates of the subpopulations, favouring one subpopulation over the other, as in natural selection.The crucial difference between the two models is that in the genetic model, the total population size is conseved whereas in the case of the tumour model the total size () is a function of time.The dynamical variables,  1 and  2 , representing the fractions of the total population that belong to the NCP and CP subpopulations respectively, can, however, attain finite steady state values due to a cancellation of the time-dependent factors in the numerator and the denominator of and , in the limit of large times (attainment of the steady state).This is possible due to the process of interconversion between the phenotypes via asymmetric cell division.In Appendix A, we show this explicitly for the tumour model (with dynamics as in ( 3) and ( 4)) under steady state conditions.The steady state value of  1 is given by the expression in (2) with λ replacing β, confirming the validity of the reduced tumour model in (6).From the general to the reduced tumour model, the number of parameters decreases from four to one.

Critical-point Transition to Bimodality
Taking cues from the population genetics model [6,7], the results for a critical point transition to bimodality in the cancer model are as follows.The stochastic dynamics are governed by the differential rate equation: where  and 1 −  represent the fractions of the total population belonging to the NCP and CP subpopulations respectively and () = (1 − ).For convenience, the dimensionless time   is represented by  itself.During the evolution of a tumour, significant changes occur in the tissue microenvironments, to which the tumour cells are exposed, so that the parameter λ, representing the coupling of the system to the environment, changes progressively as the tumour evolves and is also subjected to fluctuations.The average value of  changes from positive to negative as the growth rate constant  2 , associated with the CP population, exceeds  1 , the growth rate constant of the NCP subpopulation.A positive (negative) value of λ favours the NCP (CP) subpopulation with fractional value  ((1 − )).The fluctuations in the dynamics are taken into account by replacing λ by  + .Let (, ) denote the probability of finding  in the interval (,  + ) at time .In the Ito formalism, the stochastic differential rate equation ( 7) yields the Fokker-Planck equation (FPE): where () = 1 2 −  + ( − 1).Note that the noise has the character of a multiplicative noise as it depends on the state of the system.In the case of additive noise, there is no state dependence and the FPE has the form as in (8) with () = 1.From ( 8), the steady state probability density function (PDF) is given by where  represents the normalization constant.The extrema,   , of the PDF ( The significance of the extrema lies in the fact that they are the suitable indicators of a transition [7].In a unimodal to a bimodal distribution, the number of maxima of   () change from one to two.The maxima have a natural interpretation in terms of the 'phases' of the system since they denote the most probable values preferentially observed in experiments.Knowing   (), one can define a stochastic potential, () as From ( 9), the explicit expression for () is The maxima of   () correspond to the valleys of the potential whereas a minimum becomes a hilltop.In the case deterministic dynamics, the deterministic potential, (), is defined as . The minima of () are the stable steady states of the dynamics.In the case of additive noise, the deterministic and stochastic potentials are identical modulo an inessential constant so that the minima of the stochastic potential are the macroscopic steady states.This is also true in the case of the multiplicative noise if the noise intensity  2 is low.If the noise intensity is sufficiently large, the stochastic potential may develop new features like the appearance of new extrema.This is clear from (10) in which the first two terms contribute to the deterministic steady state equation and the last term represents the effect of the external noise.In the case of additive noise, the last term does not contribute as () = 1 and the extrema of the steady state PDF coincide with the deterministic steady states (in the case of the tumour model there is one unique steady state, shown in (2), with β replaced by λ).In the ball-in-a landscape picture, the effect of the additive noise is to jiggle the ball around in the potential landscape with the potental retaining its shape.The valleys and the hilltops in the deterministic case retain their identity in the presence of aditive noise.In the case of multipicative noise and low noise intensity, the number and the nature of the steady states do not change from those in the deterministic case.For large noise intensity, the possibility occurs that the solutions of (10) do not match in number and position with those in the deterministic case.Now the ball not only jiggles around in the landscape but the shape of the potential also changes from, say, a single-valley structure to one with two valleys.
Returning to (10), the roots of the equation, defining the extrema of   () are for  = 0: From ( 13), it is clear that a purely noise-induced transition occurs at a critical value   2 = 2 of the noise-intensity from a unimodal (a single maximum  0 ) to a bimodal (two maxima  + ,  − and a minimum  0 ) steady state PDF   ().Figure 1 shows the plots of   () versus  for  2 <   2 ,  2 =   2 and  2 >   2 .The PDF is flat-topped at the critical noiseintensity, characteristic of a critical-point transition.

Critical Exponents
The critical exponents for the purely-noise-induced transition in the population genetics model are calculated using some standard procedure [7].For biochemical models with positive feedback and bistable deterministic dynamics, Bose and Ghosh [10] have shown that the deterministic steady state equation serves as a starting point for the calculation of the critical exponents if the noise is additive in character.The method is generalised to the case of multiplicative noise by starting with (10), the solutions of which are the extrema of the steady state PDF.The proposed method is applicable irrespective of the character of the noise, additive, multiplicative or with both the components present.We rewrite (10) as In the absence of noise, the equation reduces to the deterministic steady state equation.We briefly sketch the method of derivation of the critical exponents [10].The idea is to reduce the equation (  ) = 0 to the form of the equation of state of the mean-field Ising model close to its critical point expressed as [24] ℎ −  − 1 3  3 = 0 (15) In (15),  represents the average magnetization per spin,  = is the reduced temperature with   being the critical temperature and ℎ is the reduced magnetic field.The order parameter of the transition is the spontaneous magnetization m (ℎ = 0) which has a zero value for  >   and a non-zero value for  <   .The critical point is given by  = 0, ℎ = 0.The equation of state further corresponds to the normal form of the steady state equation of an imperfect SP bifurcation [25]: with the bifurcation point,  = 0,  = 0, separating a regime of monostability from that of bistability.
To proceed with the derivation, the expression (  ) = 0 is Taylor expanded around a point   to the third order with the proviso that  ´ˊ(  ) = 0, so that a second-order term is absent in the Taylor expansion as in the case of the magnetic equation of state (15).The expansion has the form Comparing with (15), one can write down the following relations between the thermodynamic and dynamic quantities: For the model under consideration, At the critical point, the order parameter  = 0, i.e.,   =   .Also,  = 0 and  2 = 2.
Thus, (10) has a triple root   =  0 =  ± =   = 1 2 . In terms of the dynamical quantities, one gets (  ) = 0,  ˊ(  ) = 0 at the critical point, i.e., the SP bifurcation point.One can check that  ˊ(  ) changes sign at the critical noise intensity   2 = 2 and  ˊˊˊ(   ) is always negative.These features are consistent with the change in sign of the magnetic quantity, in (18), at the critical temperature.The noise intensity plays the role of temperature in the nonequilibrium model.The critical behaviour is identical to that of the mean-field Ising model with the parameter λ playing the role of the magnetic field.The order parameter  =  + − 0.5 exhibits power-law singularity close to the critical point, ~ ( 2 −   2 ) 1/2 with the critical exponent having the value

Nonequilibrium Model with 𝜆 ≠ 0
The critical point transition, belonging to the mean-field Ising universality class, requires the parameter λ to be zero.When λ has a finite value, the transition from unimodality to bimodality occurs when the noise-intensity exceeds a critical value dependent on the value of λ.The bifurcation curves are obtained in a parametric form since the parameter λ cannot be written explicitly as a function of  2 .The region of bistability is enclosed within a pair of bifurcation boundaries with each boundary representing a line of SN bifurcations, equivalent to first-order phase transitions [25].The computation of the boundaries is carried out in the following manner.We start with (14) and note that for a SN bifurcation to occur, the conditions to be satisfied are (  ) = 0 and  ˊ(  ) = 0.The two conditions yield the expressions From ( 20) and ( 21), one obtains the parametric relations The bifurcation curves are defined by ( 22) and ( 23).For 0 <   < 1, the points ( 2 (  ), (  )) are plotted in the ( 2 − ) plane (Figure 2).Note that the variable   represents an extremum (maximum) of the steady state PDF so that the phase diagram has a stochastic character.One identifies a region of bimodality separating two regions of unimodality.In the bimodal (unimodal) region, the steady state PDF,   (), has two peaks (one peak).The critical point is depicted by the point (2, 0).The termination of the bifurcation curves at this point is reminiscent of a line of first-order transitions terminating at a critical point in the cases of equilibrium thermodynamic phases transitions like the liquid-gas and the paramagneticferromagnetic phase transitions [24].In the deterministic case, the region of bimodality disappears and one has only monostability in the full parameter regime.In the stochastic case, the bimodality is accompanied by hysteresis [6,7] the plots of which are computed from ( 14) for both   and (1 −   ), the population fractions, corresponding to the two subpopulations.Figure 3 shows the hysteresis curves ( 2 = 4) and one finds a reflection symmetry, between the curves, across the vertical line at  = 0.

Cancer as a Phase Transition
The heterogeneity in the tumour population is described in terms of the two subpopulations: NCP and CP.Examining the expression for the steady state PDF,   (), in where  and 1 −  represent the fractions of the NCP and CP subpopulations, respectively, of the total population.We first consider the case  = 0.This is the situation when, on an average, neither of the subpopulations is favoured by the environment.When the noiseintensity  2 is less than the critical value   2 = 2, the steady state PDF,   (), has a unique maximum at  0 = 1 2 which coincides with the deterministic stable steady state.
When the noise intensity exceeds the critical value,   () becomes bimodal with the peak positions at  ± (see (13)).One notes the relation  + = 1 −  − with  − <  + .For each subpopulation, the fractional value  has an equal probability (Figure 1) to peak at a low or a high value.Since   () =   (1 − ), if  for the NCP subpopulation is low (high), that of the CP subpopulation is high (low).One designates the  = 0 state as a state of balance between the subpopulations since each of them has equal probability to have a larger fractional value.The peak positions  − and  + move respectively to 0 and 1 respectively as the noise intensity  2 → ∞.
A simple analogy with a population of colour-coded chemical reactants provides a clear physical picture of the population heterogeneity in the two distinct regimes:  2 <   2 and  2 >   2 for  = 0 [7].The dynamics of the system of reactants are as given in (1).There are two subpopulations of reactants of yellow and blue colour with the fractional values  and 1 −  respectively.For  2 <   2 , neither subpopulation dominates and the reacting system exhibits a flickering green colour.Above the critical noise-intensity, the reacting system will either be predominantly yellow or predominantly blue with equal times, on an average, spent in both the states.In terms of the cancer model, the population is mostly NCP or CP, each possibility occurring with equal probability.When  has a finite value, the balance tilts in favour of the NCP (CP) subpopulation for  > 0 ( < 0). Figure 4 shows the cumulative distribution function (CDF) of the random variable  versus  with the CDF defined as ( * ) gives the probability that the fractional value  of the NCP subpopulation is less than or equal to  * .The curves labeled "" and "" correspond to   (0.5) and 1 −   (0.5) respectively.From Figure 4, it is clear that the NCP subpopulation becomes the dominant one as λ becomes positive (Curve  lies above Curve ), the dominance becoming more prominent as λ increases towards more positive values.For a sufficiently large value of λ, the CDF reaches its maximal value of 1.For the CP subpopulation, the plots  and  represent the CDFs 1 −   (0.5) and   (0.5) respectively.The subpopulation becomes more and more dominant as  acquires more negative values.Figure 5 shows the evolution of   () as λ changes from positive to negative values and with  2 = 4.0.Invoking the relation in (24), the same set of plots represent the evolution of   (1 − ) but with the sign of λ changed.In the region of bimodality and for  = 0, the NCP and CP subpopulations are in a balanced state.The probability that the NCP (CP) subpopulation becomes more dominant is enhanced as λ increases towards positive (negative) values.In the regions of unimodality (  = 1.5,  = −1.5), the steady state PDF has a single peak with the NCP (CP) subpopulation dominating over the other subpopulation for positive (negative) λ value.Within the region of bimodality, the opportunity exists to reverse the dominance of the CP population through appropriate strategies including drug therapy but once the bifurcation boundary (Figure 2) is crossed into the region of unimodality, irreversibility in the form of hysteresis (Figure 3) sets in.The sudden deterioration from the pre-disease to the disease state at the bifurcation point effectively constitutes a 'point of no return'.Figure 6 shows the evolution of the stochastic potential () as a function of λ in a few representative cases with  2 = 4.0.From Figure 6, one finds that in the case of  = −1.5, a very small fraction of the total population belongs to the NCP subpopulation so that the CP subpopulation constitutes the major component of the tumour population.The relative stability between the NCP and CP subpopulations shifts with changes in the parameter λ.
The key point emerging from our study is that the progression of a tumour population from the normal to the cancer state can be understood in terms of the critical-point and first-order phase transitions [26,27,28] which occur in the nonequlibrium (dynamical phase transitions).In the normal state, the NCP and the CP subpopulations are kept in balance with both the subpopulations equally favourable.This is specially so for  = 0 and the noise-intensity  2 <   2 .A critical-point transition to a region of bimodality occurs at the critical noise-intensity (Figure 1) with the steady state PDF developing two peaks around the states,  − and  + = 1 −  − .The fractional values,  − and  + occur with equal probability so that the subpopulations continue to remain in balance.The critical-point transition into the region of bimodality opens up the possibility of one subpopulation outcompeting the other.This asymmetric situation or tilting of the balance is possible when the parameter  has a non-zero value with the NCP (CP) subpopulation becoming dominant for positive (negative) values of λ.
When the growth rate constant of the CP subpopulation exceeds that of the NCP subpopulation, the parameter λ becomes negative and the pre-disease state is reached at the bifurcation boundary in Figure 2. Up to the pre-disease state, the deterioration towards the disease state is reversible and the opportunity for damage control exists.In the region of unimodality and for negative values of λ, the disease state is reached with the CP subpopulation achieving irreversible dominance.This is clear from the hysteresis curves ( Figure 3) which show that to return to the original branch, once a bifurcation has taken place, an "overshoot" of the control parameter value is requred.In the region of bimodality, the upward and downward jumps occur at the same parameter values, i.e., the transitions are reversible.

Quantitative Signatures of the Onset of Dominance
Cell-fate transitions, as in the cases of cell differentiation and the sudden deterioration in the progression of complex diseases like cancer, share similarities with phase transitions [17,26,29,30].Rapid advances in single-cell techniques have made it possible to obtain network-wide gene expression data at single-cell resolution for a large ensemble of cells.The large-scale data provide signatures of dynamical bifurcations like the pitchfork and SN bifurcations, coinciding with experimentally observed cell-fate transitions.The pre-disease state in the route from the normal to the disease state marks the onset of irreversibility as in the case of the SN bifurcation.In literature the transition to irreversibility is referred to as the tipping point and sometimes as the "critical-point" transition [20] though the latter nomenclature is reserved in the physics literature for continuous transitions.The method of dynamical network biomarker (DNB) has been developed to detect the early warning signal (EWS) of the approach to the pre-disease state, in the proximity of a bifurcation point [17,18].The DNB consists of a dominant group of genes in the gene regulatory network which provides the EWS in terms of the following features of the experimental data: (i) the coefficient of variation (CV) of the expression of a DNB gene increases significantly as the bifurcation point is approached, (ii) an increased covariation between the members of the DNB and (iii) the rapid decrease in the correlation between a DNB and a non-DNB element of the network.Appropriately constructed criticality indices capture the increased covariation in the vicinity of the tipping point [17,18].The DNB method has been applied successfully to high-throughput experimental data towards the early detection of the pre-disease state in a number of complex diseases including cancer.
In recent studies, some specific entropic measures have been proposed to detect the predisease state from the data on reference samples and a single case sample.The concept of entropy is fundamental in statistical mechanics and information theory [31,32].The entropy of a random variable provides a quantitative measure of the uncertainty associated with the random variable.On an average, it is a measure of the amount of information required for a description of the random variable.Various entropic measures have been proposed so far to analyse single-cell datasets obtained in high-throughput experiments.One of these measures is that of relative entropy, a measure of the distance between two probability distributions.In the case of discrete random variables, the relative entropy, in terms of the Kullback-Leibler (KL) distance [31,32], is defined as where () and () represent the probability mass functions.We take the base of the logarthm to be 2.In the case of continuous random variables, the summation is replaced by an integral with (), () representing the PDFs.The KL distance provides a measure of the loss in information when an approximation () is used for the true mass function ().The KL distance is always positive and has the value zero only if () = ().It is, however, not a bonafide distance measure as it is not symmetric and does not satisfy the triangle inequality.The Jensen-Shannon divergence is defined to be [33,34]   (||) = 1 2   (||) + 1 2   (||),  = 1 2 ( + ) (27) One can check that   (||) =   (||) (symmetry) and 0 ≤   (||) ≤ 1.The zero (one) value indicates perfect (zero) overlap between the distributions.Also, the quantity √  (||) defines a metric distance satisfying the triangular inequality.The definitions ( 26) and ( 27) hold true when  and  represent the cumulative probability.Based on the last definition, incorporated in a single-sample-based Jensen-Shannon Divergence (sJSD) method, a criticality index ICI (inconsistency index) has been computed from the real datasets of a few complex diseases including some specific types of cancer [34].The criticality index exhibits a sharp peak as a signature of a tipping point transition from a pre-disease to a disease state.
In our two-subpopulation model of tumour evolution, the key dynamical variables are the fractions  and (1 − ) of the NCP and CP subpopulations, respectively, in the total population.In this case, as already discussed, the NCP (CP) subpopulation becomes progressively more dominant as the value of λ increases towards more positive (negative) values.In the following, taking cues from the EWS of the approach to the pre-disease state [17,18], we provide three quantitative signatures of the transition from the state of balanced subpopulations,  = 0, to the emergence of a dominant subpopulation.We confine our attention to the parameter regime,  2 >   2 .Figure 7 shows the plot of the variance of the steady state probability distribution (( 9)) versus the parameter  for  2 = 4.The approach to the state of balance from both the sides is indicated by a rise in the magnitude of the variance with the peak value occurring in the state of balance itself.If initially the NCP subpopulation is dominant, an increase in the variance would indicate the approach to the state of balance.The state of balance effectively serves as a tippng point from the dominance of the NCP subpopulation to that of the CP subpopulation.The deterioration can be reversed as long as the system is in the region of bimodality with irreversibility setting in as the bifurcation boundary is crossed.), all the cells in the tumour population belong to either the NCP or the CP subpopulation.The situation is characterized by the term "contradictory information (CI)" [35] in contrast to the case of "minimum information", a feature of the uniform distribution.In the latter case, the random variable  has an equal probability to assume any value in a range of values whereas in the former case the variable has a choice between two extreme values with equal probability.Another interesting feature of the bipolar limit is that the transitions between the peaks at  = 0 and 1 are extremely rare so that the stabilization of one phenotype occurs with respect to the other phenotype [7].When  ≠ 0, the growth rates of the subpopulations become unequal but the so-called "unfit" phenotype may still be stabilized.This happens because even if the system is not at the bipolar limit, the transitions between the probability peaks are still infrequent.

Conclusions
Population heterogeneity, in terms of distinct phenotypes, is a widely observed feature in cellular populations.The heterogeneity is mostly driven by noise, of both intrinsic and extrinsic origins, and has a significant role in cellular decision-making processes as well as in the implementation of the bet-hedging strategies for the survival of populations [1,2,3].In the presence of noise, the dynamics have a stochastic component and the attractors of the dynamics in the state space describe the "stable" states of the system.An attractor corresponds to a cloud of points in the state space which reduces to a single point or a cycle in the limit of deterministic dynamics.The cell population dynamics are subjected to two effectively opposing influences, one attracting the dynamics towards the cloud centre or to a cycle of states ( the "drift" term appearing as the first term on the right-hand-side of the FPE ( 8)) and the other giving rise to diffusion in the state space ( the "diffusion" or "noise" term appearing as the second term on the right-hand-side of the FPE).The idea of attractors as different cell types was demonstrated in a pioneering study by Kauffman [36] and forms the central concept in the analysis of experimental observation-based data [30,37,38].In the context of a tumour population, Kauffman first proposed the hypothesis that the cancerous state could be an attractor of the tumour dynamics [39].
It is now well-recognized that "dramatic" changes in the tissue microenvironments are associated with the progression of tumours [12,13] and these microenvironmental changes play an active role in tumorigenesis, specifically, in promoting the dominance of the CP subpopulation.One principal result of our study is to establish the equivalence between the two-subpopulation dynamical model of tumour progression, as described in ( 3) and ( 4), and the well-known population genetics model with dynamics as described in (1).The mapping enables one to reduce the number of parameters from four ( 1 ,  2 ,  1 ,  2 ) to one, namely, λ.The parameter λ takes care of the coupling to the environment and is expected to undergo large changes in its average value as a tumour progresses to the cancerous state.As shown in (6), the steady state subpopulation fractions, (  ,  = 1, 2), are governed by the single parameter λ, which is a ratio of the parameters   and 2  .The parameter is further subjected to considerable fluctuations due to the random variations in the tissue microenvironments.
Because of the equivalence, some well-known results of the stochastic population genetics model are of validity in the case of the two-subpopulation tumour model.The most important of these is the existence of a noise-induced transition from unimodality to bimodality when the noise intensity exceeds a critical value.The maxima of the steady state PDF, which correspond to the most probable states, define the "phases" of the system.We have developed a procedure, based on the knowledge of the extrema of the steady state PDF(( 14)), to explore the phase diagram of the NCP-CP model exhibiting both critical-point and first-order phase transitions in the nonequilibrium.The critical exponents, associated with the critical-point transition  = 0,  2 = 2, belong to the mean-field Ising universality class.The evolution of the steady state PDF and the stochastic potential as the parameter λ changes (Figures 5 and 6) capture the progression of the tumour population from a "healthy" to a "disease" state.In the former (latter) state, the NCP (CP) subpopulation becomes dominant.Some of the concepts and methodologies of the statistical mechanical description of phase transitions, like the state transition from a unimodal to a bimodal frequency distribution through a flattened unimodal profile (Figure 1) and the depiction of the state-transitions via potential landscapes and hysteresis curves (Figures 3 and 6), have made their appearance in recent studies on cancer.In such studies, cancer is interpreted in terms of a state-transition in the time series of transcriptome (gene expression) data with the transition points (bifurcation/phase transition points) identified as the points in state space at which dynamical regime shifts take place [17,28,34,40,41] .
In the NCP-CP model, apart from the usual bifurcation transitions of the SP and SN type, there is another transition from a state of balance to a state of dominance in terms of the NCP and CP subpopulations.This transition can be detected through the changes in the variance of the steady state PDF and the entropic measures like the JS divergence and the CPE as the parameter λ changes.The use of information-theoretic measures to detect the pre-disease state at which the irreversible transition to the disease state occurs is now a standard tool in the hands of data analysts.Some strategies for escaping the cancer attractor have been proposed by Huang and Kauffman [42].A noise-based strategy of limited application till now arises from the experimental evidence of noise modulators, which are chemical compounds and biochemical complexes [43,44].These modulators (noise enhancers) do not change the mean expression level and by regulating the noise intensity can bring about an escape from a specific, say cancer, attractor.The stochastic NCP-CP model demonstrates the key role played by multiplicative noise in the progression of a tumour population to the cancerous state in which the CP subpopulation becomes dominant.The model, the key variables of which are the subpopulation fractions, is simple enough to be analytically tractable and may serve as a starting point for the characterization of the ITH taking into account the random variability of the tissue microenvironment.The assumptions underlyng the tumour model, with the defining equations as in ( 3) and ( 4), have considerable experimental support [14,15].The proposal of a purely noise induced transition brought about by multiplicative noise, as demonstrated in the case of the population genetics model, is realizable in a number of experimental systems [7].
The link between the experimentally motivated two-subpopulation tumour model and the population genetics model lays open the possibility of designing experiments to probe how the tumour microenvironment noise tilts the phenotypic composition of the tumour in favour of the cancer-promoting subpopulation.
solution of the equation[7]

Figure 2 .
Figure 2. Parametric plots of bifurcation curves (22) and (23) in the ( −  2 ) plane.The curves separate a region of bimodality from the regions of unimodality.

Figure 3 .
Figure 3. Hysteresis plots for   and 1 −   versus λ with  2 = 4.The two plots correspond to the two phenotypically distinct subpopulations.
, one finds that for  = 0 the PDF has exchange symmetry.It remains invariant under the transformation  → 1 − , i.e., when the subpopulation fractions are interchanged.When λ has finite values, the exchange symmetry is broken and one finds the relation   (, ) =   (1 − , −)(24)

Figure 5 .
Figure 5. Evolution of the steady state PDF   () versus  as the parameter decreases from positive to negative values with  2 = 4.0.The same set of plots represent the evolution of   (1 − ) versus  but with the sign of λ changed.

Figure 6 .
Figure 6.Stochastic potential () ((12)) versus  with  2 = 4.0.From Figure6, one finds that in the case of  = −1.5, a very small fraction of the total population belongs to the NCP subpopulation so that the CP subpopulation constitutes the

Figure 8
Figure8shows the plot of the Jensen-Shannon divergence,   (||) ((27))[34], versus the parameter λ.The steady state probability distributions  and  are as given in(9) with  2 = 4,  ≠ 0 for the distribution  and  2 = 2,  = 0 for the distribution Q.The latter distribution is the flat-topped critical probability distribution shown in Figure1.Thus the JS divergence provides a measure of how different the distribution  is from the critical distribution as the parameter λ is varied.From the Figure it is clear that the closest similarity is for λ = 0. Away from this value, the difference between the two distributions become more and more dissimilar as λ increases towards more positive/negative values.Figure9shows the plot of   (||) versus the parameter λ but now  2 = 4 for both the distributions  and  but the λ-values have opposite signs.Since the NCP and CP subpopulations become dominant for the positive and negative values of λ respectively, the JS divergence captures the dissimilarity of the two situations, one favouring the NCP and the other favouring the CP subpopulation.

Figure 9 . 2 and 2 .
Figure 9.The plot of   (||) versus λ.For both the distributions  and ,  2 = 4 but the λ values are opposite in sign.