Cancer presents a serious threat to human health. The understanding of the cell fate determination during development and tumor-genesis remains challenging in current cancer biology. It was suggested that cancer stem cell (CSC) may arise from normal stem cells or be transformed from normal differentiated cells. This gives hints on the connection between cancer and development. However, the molecular mechanisms of these cell-type transitions and the CSC formation remain elusive. We quantified landscape, dominant paths, and switching rates between cell types from a core gene regulatory network for cancer and development. Stem cell, CSC, cancer, and normal cell types emerge as basins of attraction on associated landscape. The dominant paths quantify the transition processes among CSC, stem cell, normal cell, and cancer cell attractors. Transition actions of the dominant paths are shown to be closely related to switching rates between cell types, but not always to the barriers in between, because of the presence of the curl flux. During the process of P53 gene activation, landscape topography changes gradually from a CSC attractor to a normal cell attractor. This confirms the roles of P53 of preventing the formation of CSC through suppressing self-renewal and inducing differentiation. By global sensitivity analysis according to landscape topography and action, we identified key regulations determining cell-type switchings and suggested testable predictions. From landscape view, the emergence of the CSCs and the associated switching to other cell types are the results of underlying interactions among cancer and developmental marker genes. This indicates that the cancer and development are intimately connected. This landscape and flux theoretical framework provides a quantitative way to understand the underlying mechanisms of CSC formation and interplay between cancer and development. Cancer Res; 75(13); 2607–18. ©2015 AACR.

Major Findings

We developed a landscape and path theoretical framework to investigate the global natures and dynamics for a core CSC gene network. Landscape exhibits four basins of attraction, representing CSC, stem cell, cancer, and normal cell states. We quantified the kinetic rate and paths between different attractor states. We uncovered certain key genes and regulations responsible for determining the switching between different states.

With the Hill function describing the activation or repression regulations, the ordinary differential equations (ODE) have the form as:

formula

F represents the driving force for the change of protein concentration with time. Here, in Eq. A, i = 1, 2, …, 6, so we have a total of six equations. S represents the threshold of the explicit sigmoidal functions, i.e., the strength of the regulatory interaction, and n is the Hill coefficient that determines the steepness of the sigmoidal function (1) or cooperativity of the gene regulatory binding (dimer, trimer, tetramer, etc.). Here, parameters for Hill function are specified as: S = 0.5, n = 4 (see Supporting Information for the determination of parameter values). In addition, k is the self-degradation constant, b is the repression constant, and a is the activation constant. In the above equation, the first term represents the self-degradation, the second term represents the activation from the node j to the node i (m1 represents the number of activations to node i, and this term represents self-activation when i = j), and the last term denotes repression from the node j to the node i (m2 represents the number of repressions to node i, and this term represents self-repression when i = j). Here, we designated the values of parameters uniformly (Eq. A), i.e., all activation strengths a are same and also for the repression strengths b, because so far we have no corresponding information about the regulation strength—or the magnitude of activation and repression parameters—between different genes in the underlying gene network of cancer and development.

When the diffusion coefficient D is small, the moment equations can be approximated to (2, 3)

formula

Here, x, |$\bisigma (t)$|⁠, and A(t) are vectors and tensors, and |$A^{\rm T} (t)$| is the transpose of |${\bi A}({\bi t})$|⁠. The matrix elements of A is Aij = ∂Fi[x(t)]/∂xj(t). From these equations, we can solve |$\bar {\bi x}(t)$|and |$\bisigma (t)$|⁠. Here, we only consider diagonal elements of |$\bisigma (t)$| from mean field splitting approximation. Therefore, the evolution of probabilistic distribution for each variable can be acquired through the mean and variance from Gaussian approximation:

formula

|$\bar x(t)$| and |$\sigma (t)$| are solutions from Eqs. B and C.

  • The DNA binding/unbinding speed is much faster than the protein degradation or synthesis speed, so that the activation or repression regulations are described by the Hill functions.

  • The parameter values for different regulation strengths are uniform.

  • In cells, there are external noise from environments and intrinsic noise from small number of molecules. We assume that the external noise is the major source of noise in our cancer stem cell model. This is due to the fact that the protein production in mammalian cells is abundant.

Cancer presents a serious threat to human health. Conventionally, the cancer has been seen as a disease caused by mutations. It has become more and more clear that cancer is not merely a disease only originating from some single gene mutation, but a disease of state. Cancer should be a particular natural state of cell regulation normally inaccessible, which is hidden under the layers of complex molecular networks (4–9). These molecular networks form a variety of cell types by evolution. The driver gene mutations may play a permissive role by allowing cells to access these hidden but immanent cellular states, but not a constructive role in a stepwise causation of cancer cell phenotype. Cancer would be viewed as an intrinsic state and a hidden default program, only released by a series of mutations.

Cancer states should be naturally emerging functional entities, the result of collective action of all the genes and interactions in the gene regulatory network. They become robust due to the flow of the unstable states in their state space neighborhood (basin of attraction). A single network can give rise to multiple attractors, leading to multistability. Thus, to learn the mechanism of cancer, we need to investigate the underlying gene regulatory networks and associated dynamics, which govern the evolution and behavior of normal and cancer states. There have been increasing numbers of studies on the nature of the network (10) in terms of either deterministic or stochastic chemical kinetics. It often probes only the local properties of the network (4, 11). However, global nature of the system cannot be easily revealed from such analysis directly.

We have developed a conceptual yet quantitative and physical nonequilibrium landscape theory for cancer to meet the above-mentioned challenges (9). The concept that emerges from the study is as follows: (i) The normal states and cancer states can be quantified as basins of attractions on the cancer landscape, the depth of which represent the associated probabilities; (ii) the global stability, behavior, and function of cancer can be quantified through the landscape topology (barriers between the normal and cancer states and basin depths); (iii) quantification of the paths and transition speed from normal (cancer) state to cancer (normal) state representing the underlying tumor-genesis and reverse processes, leading the hope for prevention and treatment; (iv) the results of global sensitivity analysis based on landscape topography can identify multiple targets, including genes and gene regulations influenced by the environments for cancer. In this approach, the paradigm for cancer study is shifted from exploring local properties (cancer as a disease of individual driver mutation) to quantitatively studying global natures (cancer as a disease of states) of cancer through the landscape of regulatory networks.

However, there is another very important aspect we have not mentioned in the above cancer studies: the development and differentiation. Waddington (12) proposed a picture for development. The differentiation process is viewed as a marble rolling down from the top of the mountain hill (stem cell) to the bottom of the valleys (differentiation). Although Waddington provides a picture for development, the description is rather qualitative and is only at a metaphor level. Recently, quantification of Waddington landscape for development has been achieved (13–16). The stem cell and differentiated states can all be described by the basins of attractors. The developmental and differentiation process can then be quantitatively described as the transition between these attractors. The depths of individual attractor and barrier heights between the attractors become quantitative measures of the degrees of difficulty for differentiation and development. The dominant paths between the attractor states provide the quantification of the process of differentiation and development. From this one can identify the differentiation and reprogramming path critical for the practice of regeneration medicine. The global sensitivity analysis from the landscape topography helps to identify the key genes and regulations critical for the differentiation and reprogramming for targeting.

The hint of the link between development and cancer starts to surface out because of appearance of cancer stem cell (CSC). The CSC studies suggest that tumors grow and develop from a minority population of cells called CSCs (17–19). CSCs are highly tumorigenic cell types and are suggested to play important roles in oncogenesis because of their abilities of initiating tumors and driving metastasis (20, 21). Like normal stem cells, CSCs have the ability of self-renewal and differentiation to multiple cell types. Some researchers believe the frequent comebacks of cancer after chemotherapy or radiotherapy is because of the CSCs not being removed in the above-mentioned procedures (17). Therefore, uncovering the relationship between the cancer and differentiation may provide hope for developing a strategy of combatting cancer. Furthermore, it was postulated that CSCs may arise from normal stem cells, or be transformed from normal differentiated cells because of mutations (18, 22). Meanwhile, some experimental studies have shown the possibility of transformation between cancer cells and cancer stem-like cells (23). However, the molecular mechanisms of these cell-type transitions as well as the CSC formation remain elusive.

In this work, we aim to uncover the relationship and interplay between cancer and development. We start with the underlying gene regulatory network with both cancer and development marker genes as well as the interactions between them. We uncovered potential landscape for cancer and development. Four attractors emerge on the landscape, and characterize individually the CSC, stem cell, cancer, and normal cell types. The cancer, normal, and stem cells emerge as a result of gene regulatory interactions of cancer genes and of developmental genes, respectively. The CSC and its associated transformation to other cell types (cancer, normal, stem) emerge from the interactions between cancer and development marker genes. We also uncovered the transition processes between these cell types by identifying the dominant paths among these four attractors. We found that the transition actions S for the dominant paths are correlated with the kinetic transition rates among these cell types quantified by the first passage time [mean first passage time (MFPT)]. We also identified the key genes and regulations responsible for the cell-type transitions by performing the global sensitivity analysis on the kinetic transitions between the basins of attraction.

Model of gene regulatory network for cancer and development

By integrating some cancer or stem cell–related core circuits (8, 24–29), we constructed a simplified gene network of cancer and development (Fig. 1), which includes 6 gene nodes (Supplementary Table S1) and 16 edges (7 activation regulations and 9 repression regulations). MDM2 and P53 are cancer marker genes, whereas the ZEB and OCT4 are development (stem cell) marker genes. ZEB is also a major transcription factor promoting the epithelial–mesenchymal transition (EMT), which gives rise to the formation of CSCs (30). The miRNA provides regulations to both cancer and development (25, 26, 31). They mediate the interactions between the cancer and developmental genes. Circle nodes represent genes or miRNAs repressing tumors whereas hexagonal nodes represent stem cell marker genes. From the topology of the network, the major tumor repressor gene P53 represses the stem cell marker gene ZEB and OCT4 through two miRNAs, miR-200 and miR-145. There are some pairs of mutual repression links, such as the ones between miR-200 and ZEB, between miR-145 and ZEB, and between miR-145 and OCT4. In addition, oncogene MDM2 represses the expression of P53. There are also self-activations for certain genes, such as for P53, ZEB, and OCT4. Two classes of major maker genes (stem cell maker genes and tumor repressor genes) determine the attractor states (phenotypes) of the system. For example, high ZEB and OCT4 expression denotes stem cell–related states, and low P53 expression denotes cancer-related states. In reality, it might be much more complicated. Lower P53 expression might not necessarily lead to cancer-related states, implied by the heterogeneity in tumor cells. More relevant experimental data are needed for measuring corresponding gene expression level in different cell types, to construct more precise model for CSC.

Starting from the topology of the network, we can write down the ODEs that govern the time evolution of variables (gene expression level), shown in Eq. A.

Self-consistent mean field approximation

The time evolutions of dynamical systems are governed by the diffusion equations. Given the system state P (X1, X2, …, Xn, t) (X1, X2, …, Xn represent the concentration or number of molecules of species), we expect to have n-coupled differential equations, which are hard to solve analytically. Applying a self-consistent mean field approach (14, 32–35), we split the probability into the products of individual ones: |${P}\left({X_1, X_2, \ldots, X_n, t} \right)\sim\prod\nolimits_i^n {{P}\left({X_i, t} \right)}$| and solve the probability self-consistently. This effectively reduces the dimensionality of the system from MN to M × N, and hence makes the problem computationally tractable.

However, for the multidimensional systems, it is still difficult to solve diffusion equations directly. We start from moment equations and assume specific probability distribution based on physical arguments, i.e., we assume some specific connections between moments. In principle, we can construct the probability distribution, if we know all moments. In this work, we assume Gaussian distribution as an approximation, which means we need two moments, mean and variance.

When the diffusion coefficient D is small, the moment equations can be approximated to Eqs. (B) and (C). Further, the steady-state probability distribution can be obtained by Gaussian approximation from Eq. D.

The probability obtained above corresponds to one fixed point or basin of attraction. If the system has multistability, then there are several probabilistic distributions localized at every basin of attraction with different variances. Therefore, the total probability is the weighted sum of all these probability distributions. We can extend the current results to the multidimensional system by considering the total probability as the product of the individual probability for each variable from the mean field splitting. Finally, once we have the total probability, we can construct the potential landscape by: U(x) = −lnPss(x).

Potential landscape and kinetic paths of the cancer and development network

Based on the self-consistent approximation, we obtained the steady-state probabilistic distribution of the gene regulatory network of cancer and development. Furthermore, we quantified the potential landscape by U = −ln(Pss) (9, 13–16, 33, 34). Here, Pss represents the probability distribution of the steady state, and U is the dimensionless potential. We projected the landscape in six dimensions (ZEB, Oct4, MDM2, P53, miR-145, and miR-200) to the two dimensions, one representing the degree of cancerization (P53) and the other representing the degree of development or stemness (ZEB) for visualization. By adjusting self-activation regulation sa and repression regulation b, we can obtain a steady-state solution with four stable states or basins of attractions, as shown in Fig. 2. The blue color region represents lower potential or higher probability, and the red color region represents higher potential or lower probability. Here, the four basins of attraction, respectively, represent four cellular phenotypes. Some recent work (23) has characterized three different cell types including cancer stem cell, luminal and basal cells in breast cancer cells, and the interconversion between these different cell types, from both experimental and theoretical perspectives. It was estimated that the CSC could occupy about 5% of the whole breast cancer cell population in some case. Therefore, we believe that CSC should be a stable cell type quantified by a basin of attraction on landscape. From the landscape point of view, the transition of cellular types can be understood as a ball rolling from one basin to another by surmounting certain barriers (saddle points).

The 6 nodes gene regulatory network of cancer and development we used here is just a simplification of a larger and more realistic network. So, we made some assumptions to determine the cell types in terms of the expression level of marker genes. We assume that high P53/miR-200/miR-145 expression represents normal cell states, whereas low P53/miR-200/miR-145 expression represents cancer-related cell states. Some experimental work showed that miR-200 and miR-145 are effective marker genes for cancer phenotypes (36, 37). In the same way, we treat ZEB and OCT4 as the markers for stemness, which means high ZEB/OCT4 expression represents stem cell related state and low ZEB/OCT4 expression represents non-stem–related state. Here, we assume that ZEB and OCT4 are the markers to distinguish CSC state and cancer state. A more precise model, adding more CSC marker genes (such as CD44, CD24) to the network, might help to elucidate more details. As a result, the gene expression pattern in terms of P53 and ZEB can form four different states or cell types, which are high P53/high ZEB as stem cell state, low P53/high ZEB as CSC state, low P53/low ZEB as cancer state, and high P53/low ZEB as normal cell state.

Some previous work showed that P53 plays an important role in suppressing cellular reprogramming, and therefore provides a barrier for the transition from normal stem cells to CSCs (38–40). Also, P53 prevents dedifferentiation by repressing the expression of CD44, one of the better known CSC marker genes (41). OCT4 has been shown to be highly expressed in both normal stem cells and CSCs (42, 43). ZEB is a major transcription factor promoting the EMT, which gives rise to the formation of CSCs (30). As such, we believe that our assumptions of exploiting P53 and ZEB/OCT4 as the marker genes to determine the cell types are reasonable.

We also calculated the dominant kinetic transition paths between different attractors (path integral methods), which are shown on the landscape. The yellow (from CSC to normal) and magenta (from normal to CSC) lines represent the dominant kinetic paths with arrows denoting the directions of the transitions. We can see that the transition paths are irreversible, reflected by the fact that forward and reverse paths are not identical.

Figure 2 provides a quantitative landscape representing the development of stem cells and CSCs, as well as the transition between stem cells and CSCs. Here, the downhill on the landscape in ZEB characterizes the developmental direction. We can see that the four attractors individually representing CSCs, stem cells, cancer cells, and normal cells are located at different layers on the developmental landscape. For example, stem cells and CSCs are on the top layer of the landscape in ZEB direction. The downhill process from the stem cell attractor to the normal cell attractor in ZEB direction reflects the differentiation process of stem cells. Cancer is located at the intermediate layer below the stem cell and above the normal cell layers in the developmental direction (ZEB). As we see, the stem cells are at the top layer and the normal cells are at the bottom layer in the developmental direction (ZEB). So the cancer cells can be seen as side track attractors or off pathway traps in the normal developmental direction analogous to misfolding scenario in protein folding process (44). In other words, cancer process can be viewed as a side or parallel track (pathway) relative to the normal developmental process (pathway).

In addition, the landscape for the gene network of cancer and development indicates not only one direction for development (the route from stem cells to normal cells or from CSCs to cancer cells) in ZEB, but also another direction for cancerization (the route from stem cells to CSCs or from normal cells to cancer cells) in P53. In fact, both developmental and cancerization directions are important in characterizing the whole process of cancer and development. The developmental process and cancerization process are intimately coupled. The presence of CSC cell type provides a hint of such connection.

The landscape picture in Fig. 2 only represents a special case for a set of fixed parameter values, where four stable states all exist. In the realistic situations, the parameters (e.g., regulations strengths and production or degradation rates) might not be constant and change with time as well as environments. As such, we explored the change of landscape topography when self-activation regulation sa and repression regulation b change (Fig. 3). With self-activation sa decreased (first column), the landscape changes from four stable states coexisting to a bistable state (CSC state and cancer state), and finally to a monostable cancer state. This process of landscape change caused by the decrease of self-activation provides a possible mechanism for cancer formation. That is, starting from a stem cell state (Fig. 3A), the cell transforms to a CSC state (Fig. 3B) with the decrease of self-activation (could be caused by certain mutations) in the cancerization direction, which finally differentiates to the cancer state (Fig. 3C, a monostable cancer state) in the developmental direction. This also suggests that the cancerization process and cellular differentiation process are connected. As reflected by the change of the landscape, the tumor-genesis can start from pluripotent stem cells, which turn into CSCs because of the accumulation of mutations. Next, the CSCs can differentiate to cancer cells owing to the unlimited differentiation ability of CSCs. For the horizontal direction (repression constant b increases), the major change of landscape is that the cancer state gradually disappears or loses its stability and the normal state becomes more stable. Therefore, increasing the repression strength b provides a possible way to make cancer cells come back to normal cells.

Bifurcation diagrams for repression/self-activation strength and corresponding landscape

To investigate the influence of regulation strengths on system dynamics, we show the bifurcation diagram for repression constant b (Fig. 4A) and corresponding landscape (Fig. 4B). To obtain the bifurcation diagram, we first calculated the nullclines separately for (dxi/dt = 0, i ≠ P53) and (dxi/dt = 0, i ≠ ZEB) on the P53 and ZEB plane at fixed parameter values. By calculating the cross-points of the two nullclines, we can find the stable fixed points and unstable fixed points of the system. Repeating this process through changing parameters, we can acquire a bifurcation diagram for corresponding parameters.

We can see that with the repression strength b increased, the cell experiences a transition from a monostable state to a bistable state, which can also be seen by the landscape change from one basin of attraction representing cancer cell state (Fig. 4B, b = 0.1) to two basins of attraction representing CSC state and normal cell state coexist (Fig. 4B, b = 0.8). This indicates that by increasing repression strengths of regulations, cells are inclined to make a transition from a cancer state to a normal state.

We also explored the bifurcation diagram for the self-activation regulation strength sa1 (self-activation for P53) and sa2 (self-activation for ZEB and OCT4) and the corresponding landscapes, which are shown in Fig. 5. As self-activation sa1 increases (Fig. 5A), the cell transforms from a monostable state to a bi-stable state and finally to another monostable state. The corresponding landscape changes from one attractor (CSC state) to two coexisting attractors (CSC and normal state coexist), and finally to a monostable normal state attractor (Fig. 5B–D).

It was shown that P53 has a role of suppressing self-renewal and inducing differentiation after DNA damage in embryonic stem cells (ESC; refs. 38 and 45). Here, regarding the CSC state as an abnormal state (DNA damage), our landscape view provides an intuitive and quantitative explanation for the afore-mentioned role of P53 as a tumor repressor. With P53 activated (sa1 increases), the landscape changes gradually from a dominant CSC attractor to a dominant normal state attractor, which demonstrates the role of P53 repressing cancer and inducing differentiation of ESCs with DNA damage (preventing the formation of CSCs). In contrast, when P53 is inactivated (sa1 decreases), the landscape changes gradually from a normal cell state to a CSC state, which provides an explanation for the fact that a differentiated cell can gain mutations that drive the formation of cancer (38).

Similarly, as the self-activation strength sa2 (self-activation of ZEB and OCT4) increases, the cells gradually transform from a normal cell state to a CSC and normal coexisting state and finally to a dominant CSC state (Fig. 5E–H shows the bifurcation diagram and corresponding landscapes). This indicates that by enhancing activation of ZEB and OCT4, the normal cell types can be transformed to CSC cell types. This is reasonable because the activation of ZEB and OCT4 switches on the stemness of the system, and therefore makes tumor-genesis possible.

Transition action, potential barrier, and MFPT

When the parameters are specified as a = 0.5, b = 0.8, and sa1 = sa2 = 0.5, the system exhibits a bistable behavior (CSC state and normal state coexist). To quantify the landscape topography, we define barrier height as the potential difference between each attractor minimum and the corresponding saddle point on landscape. We define USCSC as the potential difference between the CSC state and the saddle point, UsaddleUCSC, and USN as the potential difference between normal state and the saddle point, UsaddleUN. Here, UCSC and UN denote individually potential at the minimum for the CSC attractor and the normal state attractor, and Usaddle denotes the potential at the saddle point between these two basins. In the bistable regime, by changing sa1, we explored the correlation among the transition action for the dominant path from one attractor to the other, the barrier, and the MFPT, shown in Fig. 6.

We see that as the regulation strength sa1 (self-activation for P53) decreases, the transition action S for the path from normal cell attractor to CSC attractor increases (Fig. 6B). This is reasonable because the activation of P53 inclines the bistability toward the normal state, making the transition from the normal state to CSC state more difficult, and thus leads to the increase of the corresponding transition action S and decrease of the chance of transition. Similarly, the transition action S for the dominant path from CSC attractor to normal attractor decreases as sa1 increases (Fig. 6B). This leads to easier transition from CSC attractor to normal attractor. In addition, Fig. 6A shows that the barrier for CSC attractor and the barrier for normal attractor both increase with sa1 increased. This indicates that two stable states both become more stable as sa1 increases.

Figure 6C shows the actions for the dominant paths between normal state attractor and CSC attractor versus the corresponding barrier. We can see that the barrier and the transition action are not always positively correlated (Fig. 6C), but the relative change of transition actions and barriers are correlated and consistent [Fig. 6D shows ΔBarrier (USNUSCSC) versus ΔAction (SN→CSCSCSC→N)]. This behavior arises because the dynamics of the system is determined by two driving forces, one from the landscape and the other from the curl flux originating from the non-equilibriumness of the system (33, 46). The barrier only reflects the effects of landscape, whereas the transition action reflects how hard it is to make transitions between different attractors, which involves the effects of both the landscape and the curl flux.

In addition, we calculated the MFPT between CSC attractor and normal attractor (see Supplementary Fig. S1 for the results between CSC attractor and cancer attractor). The MFPT was calculated based on the transition action obtained from path integrals (47). The MFPTs here are determined by two factors: exponential part and prefactor part. Here, we concentrate on the exponential part for the MFPT through the action. Because the MFPT is exponentially related to the action, this is the main contribution to the MFPT. On the other hand, the prefactors are less sensitive than the exponential to the MFPT and they are set to be equal and constant (set it equal to 1). The quantification of prefactors requires detailed information on the exact degradation and production rates for different genes, which are currently not available. Therefore, the MFPT results we obtained here are only relative and only represent the main trend with respect to the changes. The correlation between MFPT and transition action as well as between MFPT and barriers are separately shown in Fig. 6C and E. We found that the MFPT is correlated with the transition action S (inset in Fig. 6C). This is reasonable because the transition action and the MFPT both reflect how difficult it is to make transitions between different attractors. Therefore, they are consistent. On the other hand, the MFPT is not always correlated to the barrier height (Fig. 6E), because of the similar reasons for the relationship between barrier and action. We noticed that the relative change of MFPT and barrier is correlated as shown in Fig. 6F, where ΔMFPT (MFPTN→CSC − MFPTCSC→N) versus ΔBarrier (USNUSCSC) is plotted.

Clinically, we often see that after the chemotherapy or radiotherapy, the cancer returns. In our landscape picture, although the cancer state populations can be destroyed (the population represented by the probability for cancer attractor becomes much smaller), the CSC state may still exist. After a finite period of times, the evolution of the population will lead to the recurrence of cancer as shown in Fig. 7A and B. This picture corresponds to after the treatment (chemotherapy or radiotherapy), cancer cell populations are destroyed. But after some time, the cancer returns. We see that from landscape perspective (U = −lnPss, Pss is the steady-state probability distribution), destroying the cancer population may not completely kill the cancer as if the landscape is not be significantly altered (Fig. 7C and D), the cancer state is still stable and subsequent evolution dynamics will eventually fill in cancer state population even initially it is destroyed. To destroy a cancer state, one has to alter the underlying landscape such that the cancer state and the CSC state are no longer a stable and preferred state. To realize this, the more effective way seems to be either modifying the genes or the regulation interactions among genes in the underlying cancer gene regulatory network. The former might be achieved by the genetic manipulations and the latter might be achieved by genetic, epigenetic, and environmental changes.

Global sensitivity analysis

We also performed a global sensitivity analysis to identify the key regulations governing the relative stability of different attractors. For simplicity, we used the bistable system to carry out the sensitivity analysis by exploring the change in transition action S (Supplementary Figs. S2A and S2B) and barrier heights (Supplementary Figs. S2C and S2D) after changing regulation strengths. We use A1, …, A7 to represent seven activation links (Supplementary Table S2), and R1, …, R9 to represent nine repression links (Supplementary Table S3). We can see that the key activation regulations include A1(P53 → P53), A4(P53 → miR-145), A5(ZEB → ZEB), A6(OCT4 → OCT4), and A7(P53 → MDM2) (Supplementary Fig. S2A), and the key repression links include R1(MDM2 −∣ P53), R6(miR-145 −∣ ZEB), and R9(miR-145 −∣ MDM 2) (Supplementary Fig. S2B). These results confirm the role of the P53-MDM2 negative feedback loop (A1, A7, R1, R9) on cancer formation (48), and the role of miRNAs on cancer and development. It has been suggested that MDM2 dysregulation caused by downregulation of miR-143 and miR-145 contributes to epithelial cancer development (28). For example, when the regulation strength of A7(P53 → MDM 2) is increased, the transition action from the CSC attractor to normal attractor (SCSC→N) increases and the transition action from the normal attractor to CSC attractor (SN→CSC) decreases. This means that the normal attractor becomes less stable relative to the CSC attractor. This is reasonable since that MDM2 is an oncogene, and therefore when MDM2 is activated the cancer is more likely to appear (the CSC attractor becomes more stable relative to the normal attractor). The results from the global sensitivity analysis also confirm the role of ZEB as an oncogene (A5, R6). For example, when the self-activation of ZEB (A5) is strengthened, the transition action from CSC attractor to normal attractor (SCSC→N) increases and the transition action from normal attractor to CSC attractor (SN→CSC) decreases. This indicates that the CSC attractor becomes more stable and the normal attractor becomes less stable relatively, which confirms the role of ZEB as an oncogene by promoting EMT transition (30, 49). We can find that the results from barriers (Supplementary Figs. S2C and S2D) are basically in line with the ones from transition actions (Supplementary Figs. S2A and S2B). But they are not completely consistent, for example the regulation A3 and A6. This is because of the effects of flux, as we have discussed above for Fig. 6.

In this work, we developed a landscape and flux theoretical framework to quantify landscape topography, paths, and rates of switching between the cell types for a cancer and developmental core network under fluctuations. We constructed a core gene regulatory network for the cancer and development, and uncovered its potential landscape. The landscape for the cancer and development exhibits four stable basins of attraction at specific regulation regimes, which individually represent the stem cell, CSC, cancer, and normal cell types. The existence of CSC as a basin of attraction on the landscape and associated transformations to other cell types (cancer, stem, normal) illustrate the correlation between cancer (cancerous) and development (stemness). We also quantified the dominant transition paths among these attractor states to uncover the underlying mechanism and the detailed processes of cell-type transitions. When the self-activation regulation strength sa1 increases (P53 activated), the landscape changes gradually from a CSC attractor to a normal attractor. This provides a quantitative explanation for the role of P53 as a tumor repressor, for preventing the formation of CSC through suppressing self-renewal and inducing differentiation.

We showed that the weight factors (W = exp[−S]) or transition actions S for the dominant paths are closed related to the MFPT. This suggests that the transitions among different attractors are orchestrated by the least action principle and the magnitude of the transition action determines the difficulty for the transition to occur. Least action principle, as a principle in modern physics and mathematics, states that nature often finds the most efficient way through optimal path (with least action leading to optimal weight) from one point to another. Here, our results show that optimal paths for switching cell types in the organisms (among normal, cancer, stem, and CSCs) can be derived from the least action principle.

In addition, the barrier is not always consistent with the transition action or MFPT because of the presence of another driving force of the system, the curl flux (46). The transition action for the dominant path, the MFPT, and the potential barrier provide quantitative measures for the relative stability of different attractors, as well as the transitions among them.

Recent evidences suggest that cells that undergo EMT gain stem cell-like properties, thus giving rise to CSCs (30, 49). This suggests ZEB, an EMT marker gene, also promotes the formation of CSCs. In our landscape picture, when a cell switches from cancer attractor to CSC attractor or from normal cell attractor to CSC attractor, both of these processes are accompanied by an EMT process, marked by the activation of ZEB gene. When the cancer cells or normal cells undergo an EMT, not only this initiates the metastasis, but also it promotes the formation of the CSCs, making this disease fatal. So, our results also suggest that the ZEB and relevant miRNA (miR-200 and miR-145) may serve as potential anticancer targets by blocking EMT process or formation of the CSC. In this work, we used Hill function to represent all activation and repression regulations for simplicity. As such, the model can be further complemented by considering the detailed miRNA-mediated regulation dynamics (8, 50).

Of note, the current model is only a simplified model with a core gene network to describe CSC system. However, the interactions between cancer cells and stem cells are complex and many other relevant genes could play roles. Therefore, it is expected that a more complete and detailed CSC network constructed, with more marker genes and regulation pathways, may reveal more detailed mechanisms for CSC. This will help our understanding of the mechanisms about cancer and development as well as the associated dynamics. It will also help in developing new strategies for cancer prevention and treatment.

No potential conflicts of interest were disclosed.

Conception and design: C. Li, J. Wang

Development of methodology: C. Li, J. Wang

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): C. Li, J. Wang

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): C. Li, J. Wang

Writing, review, and/or revision of the manuscript: C. Li, J. Wang

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): C. Li, J. Wang

Study supervision: J. Wang

C. Li thank Dr. Li Xu and Mrs. Ning Liu for helpful discussions.

This work was supported by National Science Foundation.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1.
Huang
S
,
Guo
Y
,
May
G
,
Enver
T
. 
Bifurcation dynamics of cell fate decision lineage-commitment in bipotent progenitor cells
.
Dev Biol
2007
;
305
:
695
713
.
2.
Van Kampen
NG
. 
Stochastic processes in chemistry and physics
. 1st ed.
Amsterdam
:
North Holland
; 
1992
.
3.
Hu
G
. 
Stochastic forces and nonlinear systems
.
Hao
B
,
editor
.
Shanghai
:
Shanghai Scientific and Technological Education Press
; 
1994
.
4.
Kauffman
S
. 
Differentiation of malignant to benign cells
.
J Theor Biol
1971
;
31
:
429
51
.
5.
Huang
S
,
Ernberg
I
,
Kauffman
S
. 
Cancer attractors: a systems view of tumors from a gene network dynamics and developmental perspective
.
Semin Cell Dev Biol
2009
;
20
:
869
76
.
6.
Gatenby
RA
,
Vincent
TL
. 
An evolutionary model of carcinogenesis
.
Cancer Res
2003
;
63
:
6212
20
.
7.
Ao
P
,
Galas
D
,
L
H
,
Zhu
X
. 
Cancer as robust intrinsic state of endogenous molecular-cellular network shaped by evolution
.
Med Hypotheses
2008
;
70
:
678
84
.
8.
Lu
M
,
Jolly
MK
,
Levine
H
,
Onuchic
JN
,
Ben-Jacob
E
. 
MicroRNA-based regulation of epithelial-hybrid-mesenchymal fate determination
.
Proc Natl Acad Sci U S A
2013
;
110
:
18144
9
.
9.
Li
C
,
Wang
J
. 
Quantifying the underlying landscape and paths of cancer
.
J R Soc Interface
2014
;
10
:
20140774
.
10.
Davidson
EH
. 
A genomic regulatory network for development
.
Science
2002
;
295
:
1669
78
.
11.
Chen
KC
,
Csikasz-Nagy
A
,
Gyorffy
B
,
Val
J
,
Novak
B
,
Tyson
JJ
. 
Kinetic analysis of a molecular model of the budding yeast cell cycle
.
Mol Biol Cell
2000
;
11
:
369
91
.
12.
Waddington
CH
. 
The strategy of the genes: a discussion of some aspects of theoretical biology
.
London
:
Allen and Unwin
; 
1957
.
13.
Wang
J
,
Zhang
K
,
Xu
L
,
Wang
EK
. 
Quantifying the Waddington landscape and biological paths for development and differentiation
.
Proc Natl Acad Sci U S A
2011
;
108
:
8257
62
.
14.
Li
C
,
Wang
J
. 
Quantifying cell fate decisions for differentiation and reprogramming of a human stem cell network: landscape and biological paths
.
PLoS Comput Biol
2013
;
9
:
e1003165
.
15.
Wang
J
,
Xu
L
,
Wang
EK
,
Huang
S
. 
The potential landscape of genetic circuits imposes the arrow of time in stem cell differentiation
.
Biophys J
2010
;
99
:
29
39
.
16.
Xu
L
,
Zhang
K
,
Wang
J
. 
Exploring the mechanisms of differentiation, dedifferentiation, reprogramming and transdifferentiation
.
PLoS ONE
2014
;
9
:
e105216
.
17.
Marotta
LLC
,
Polyak
K
. 
Cancer stem cells: a model in the making
.
Curr Opin Genet Dev
2009
;
19
:
44
50
.
18.
Lobo
NA
,
Shimono
Y
,
Qian
D
,
Clarke
MF
. 
The biology of cancer stem cells
.
Annu Rev Cell Dev Biol
2007
;
23
:
675
99
.
19.
Dalerba
P
,
Cho
RW
,
Clarke
MF
. 
Cancer stem cells: models and concepts
.
Annu Rev Med
2007
;
58
:
267
84
.
20.
Sampieri
K
,
Fodde
R
. 
Cancer stem cells and metastasis
.
Semin Cancer Biol
2012
;
22
:
187
93
.
21.
Reya
T
,
Morrison
SJ
,
Clarke
MF
,
Weissman
IL
. 
Stem cells, cancer, and cancer stem cells
.
Nature
2001
;
414
:
105
11
.
22.
Klonisch
T
,
Wiechec
E
,
Hombach-Klonisch
S
,
Ande
SR
,
Wesselborg
S
,
Schulze-Osthoff
K
, et al
Cancer stem cell markers in common cancers–therapeutic implications
.
Trends Mol Med
2008
;
14
:
450
60
.
23.
Gupta
PB
,
Fillmore
CM
,
Jiang
G
,
Shapira
SD
,
Tao
K
,
Kuperwasser
C
, et al
Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells
.
Cell
2011
;
146
:
633
44
.
24.
Polytarchou
C
,
Iliopoulos
D
,
Struhl
K
. 
An integrated transcriptional regulatory circuit that rein- forces the breast cancer stem cell state
.
Proc Natl Acad Sci U S A
2012
;
109
:
14470
5
.
25.
Xu
N
,
Papagiannakopoulos
T
,
Pan
G
,
Thomson
JA
,
Kosik
KS
. 
MicroRNA-145 regulates OCT4, SOX2, and KLF4 and represses pluripotency in human embryonic stem cells
.
Cell
2009
;
137
:
647
58
.
26.
Brabletz
S
,
Bajdak
K
,
Meidhof
S
,
Burk
U
,
Niedermann
G
,
Firat
E
, et al
The ZEB1/miR-200 feedback loop controls Notch signalling in cancer cells
.
EMBO J
2011
;
30
:
770
82
.
27.
Chang
CJ
,
Chao
CH
,
Xia
W
,
Yang
JY
,
Xiong
Y
,
Li
CW
, et al
p53 regulates epithelial-mesenchymal transition and stem cell properties through modulating miRNAs
.
Nat Cell Biol
2011
;
13
:
317
23
.
28.
Zhang
J
,
Sun
Q
,
Zhang
Z
,
Ge
S
,
Han
Z
,
Chen
W
. 
Loss of microRNA-143/145 disturbs cellular growth and apoptosis of human epithelial cancers by impairing the MDM2-p53 feedback loop
.
Oncogene
2012
;
32
:
61
9
.
29.
Jolly
MK
,
Huang
B
,
Lu
M
,
Mani
S
,
Levine
H
,
Ben-Jacob
E
. 
Towards elucidating the connection between epithelial-mesenchymal transitions and stemness
.
J R Soc Interface
2014
;
11
:
20140962
.
30.
Mani
SA
,
Guo
W
,
Liao
MJ
,
Eaton
EN
,
Ayyanan
A
,
Zhou
AY
, et al
The epithelial-mesenchymal transition generates cells with properties of stem cells
.
Cell
2008
;
133
:
704
15
.
31.
Burk
U
,
Schubert
J
,
Wellner
U
,
Schmalhofer
O
,
Vincan
E
,
Spaderna
S
, et al
A reciprocal repression between ZEB1 and members of the miR-200 family promotes EMT and invasion in cancer cells
.
EMBO Rep
2008
;
9
:
582
9
.
32.
Sasai
M
,
Wolynes
P
. 
Stochastic gene expression as a many-body problem
.
Proc Natl Acad Sci U S A
2003
;
100
:
2374
9
.
33.
Wang
J
,
Li
C
,
Wang
EK
. 
Potential and flux landscapes quantify the stability and robustness of budding yeast cell cycle network
.
Proc Natl Acad Sci U S A
2010
;
107
:
8195
200
.
34.
Li
C
,
Wang
J
. 
Landscape and flux reveal a new global view and physical quantification of mammalian cell cycle
.
Proc Natl Acad Sci U S A
2014
;
111
:
14130
5
.
35.
Zhang
B
,
Wolynes
PG
. 
Stem cell differentiation as a many-body problem
.
Proc Natl Acad Sci U S A
2014
;
111
:
10185
90
.
36.
Park
SM
,
Gaur
AB
,
Lengyel
E
,
Peter
ME
. 
The miR-200 family determines the epithelial phenotype of cancer cells by targeting the E-cadherin repressors ZEB1 and ZEB2
.
Genes Dev
2008
;
22
:
894
907
.
37.
Dip
N
,
Reis
ST
,
Srougi
M
,
Dall'Oglio
MF
,
Leite
KR
. 
Expression profile of microRNA-145 in urothelial bladder cancer
.
Int Braz J Urol
2013
;
39
:
95
102
.
38.
Aloni-Grinstein
R
,
Shetzer
Y
,
Kaufman
T
,
Rotter
V
. 
p53: the barrier to cancer stem cell formation
.
FEBS Lett
2014
;
588
:
2580
9
.
39.
Hong
H
,
Takahashi
K
,
Ichisaka
T
,
Aoi
T
,
Kanagawa
O
,
Nakagawa
M
, et al
Suppression of induced pluripotent stem cell generation by the p53–p21 pathway
.
Nature
2009
;
460
:
1132
5
.
40.
Kawamura
T
,
Suzuki
J
,
Wang
YV
,
Menendez
S
,
Morera
LB
,
Raya
A
, et al
Linking the p53 tumour suppressor pathway to somatic cell reprogramming
.
Nature
2009
;
460
:
1140
4
.
41.
Godar
S
,
Ince
TA
,
Bell
GW
,
Feldser
D
,
Donaher
JL
,
Bergh
J
, et al
Growth-inhibitory and tumor- suppressive functions of p53 depend on its repression of CD44 expression
.
Cell
2008
;
134
:
62
73
.
42.
Kumar
SM
,
Liu
S
,
Lu
H
,
Zhang
H
,
Zhang
PJ
,
Gimotty
PA
, et al
Acquired cancer stem cell phenotypes through Oct4-mediated dedifferentiation
.
Oncogene
2012
;
31
:
4898
911
.
43.
Trosko
JE
. 
From adult stem cells to cancer stem cells
.
Ann N Y Acad Sci
2006
;
1089
:
36
58
.
44.
Onuchic
JN
,
Wolynes
PG
. 
Theory of protein folding
.
Curr Opin Struct Biol
2004
;
14
:
70
5
.
45.
Lin
T
,
Chao
C
,
Saito
S
,
Mazur
SJ
,
Murphy
ME
,
Appella
E
, et al
p53 induces differentiation of mouse embryonic stem cells by suppressing Nanog expression
.
Nat Cell Biol
2004
;
7
:
165
71
.
46.
Wang
J
,
Xu
L
,
Wang
EK
. 
Potential landscape and flux framework of non-equilibrium networks: robustness, dissipation and coherence of biochemical oscillations
.
Proc Natl Acad Sci U S A
2008
;
105
:
12271
6
.
47.
Wang
J
,
Zhang
K
,
Wang
EK
. 
Kinetic paths, time scale, and underlying landscapes: A path integral framework to study global natures of nonequilibrium systems and networks
.
J Chem Phys
2010
;
133
:
125103
:
1–13
.
48.
Weinberg
RA
. 
The biology of cancer
. 1st ed.
Beijing
:
Science Press
; 
2009
.
49.
Singh
A
,
Settleman
J
. 
EMT, cancer stem cells and drug resistance: an emerging axis of evil in the war on cancer
.
Oncogene
2010
;
29
:
4741
51
.
50.
Lu
M
,
Jolly
MK
,
Onuchic
J
,
Ben-Jacob
E
. 
Toward decoding the principles of cancer metastasis circuits
.
Cancer Res
2014
;
74
:
4574
87
.