Exploring the Mechanisms of Differentiation, Dedifferentiation, Reprogramming and Transdifferentiation

We explored the underlying mechanisms of differentiation, dedifferentiation, reprogramming and transdifferentiation (cell type switchings) from landscape and flux perspectives. Lineage reprogramming is a new regenerative method to convert a matured cell into another cell including direct transdifferentiation without undergoing a pluripotent cell state and indirect transdifferentiation with an initial dedifferentiation-reversion (reprogramming) to a pluripotent cell state. Each cell type is quantified by a distinct valley on the potential landscape with higher probability. We investigated three driving forces for cell fate decision making: stochastic fluctuations, gene regulation and induction, which can lead to cell type switchings. We showed that under the driving forces the direct transdifferentiation process proceeds from a differentiated cell valley to another differentiated cell valley through either a distinct stable intermediate state or a certain series of unstable indeterminate states. The dedifferentiation process proceeds through a pluripotent cell state. Barrier height and the corresponding escape time from the valley on the landscape can be used to quantify the stability and efficiency of cell type switchings. We also uncovered the mechanisms of the underlying processes by quantifying the dominant biological paths of cell type switchings on the potential landscape. The dynamics of cell type switchings are determined by both landscape gradient and flux. The flux can lead to the deviations of the dominant biological paths for cell type switchings from the naively expected landscape gradient path. As a result, the corresponding dominant paths of cell type switchings are irreversible. We also classified the mechanisms of cell fate development from our landscape theory: super-critical pitchfork bifurcation, sub-critical pitchfork bifurcation, sub-critical pitchfork with two saddle-node bifurcation, and saddle-node bifurcation. Our model showed good agreements with the experiments. It provides a general framework to explore the mechanisms of differentiation, dedifferentiation, reprogramming and transdifferentiation.


Introduction
A pluripotent undifferentiated cell can differentiate into types of differentiated cells. Each cell type has a specific regulated gene expression. Cellular differentiation is determined by the underlying gene regulatory network during the process of development, which leads the primary cell into its ultimate fate-a particular phenotype. Induced pluripotent stem (iPS) cells provide the opportunity to obtain pluripotent stem cells which potentially have therapeutic uses [1,2]. Recently many studies have been reported that one type of cells can be converted to another type of functional cells directly [3][4][5][6][7]. This is a big step forward in the cell biology since there is no need to create iPS cells first for cell type switching, skipping many intermediate steps. This direct reprogramming technology is called the lineage reprogramming. Thus an adult cell can be reprogrammed directly to new cells as lineage switching. The lineage switching through direct transdifferentiation without going through the iPS state might be applied to regenerative medicine with less risk of cancer. However, it is still challenging to quantify the mechanisms of the differentiation, dedifferentiation, reprogramming and transdifferentiation [3][4][5][6][7][8][9][10][11].
The concept of ''epigenetic landscape'' was first introduced by Waddington in 1940s [12] The quantifications of the Waddington potential landscape for the process of cell differentiation have been explored recently [13][14][15][16][17]. Different valleys represent different cell phenotypes (cell fates) on the cell development potential landscape [13][14][15][16][17]. Waddington visualized the undifferentiated state as the local maximum and differentiated states as the local minimum on the landscape [12]. In our landscape picture, the undifferentiated state and differentiated state are both local minima in certain regions of the landscape. Undifferentiated state has relatively low expressions of differentiation mark genes while differentiated state has at least one high expressions of differentiation mark genes. In addition, Waddington believed the differentiation is a downhill process driven by the funneled landscape gradient. In our picture, the differentiation can occur with several different mechanisms, through funneled landscape, through stochastic fluctuations and the probability fluxes even when the landscape is not funneled towards the differentiated states, and through induction.
For development and differentiation system, we represent a cell as a chemical system having given genomic makeup, with each and every possible phenotype as a potential ''state'' [18,19]. This is very much analogous to the notion of a polypeptide, as a chemical molecule, can have many different possible ''conformational states'', although each individual protein molecule has only a particular state at a given moment in time. This chemical definition of ''the system'' is important. Imagine that proteins are defined only through biological functions; then different conformations of a polypeptide will be considered as ''different molecules.'' Then the notion of spontaneous conformational change would not make sense. Indeed, there are still cell biologists who think different cells from the same person as different cells; rather than as a ''same chemical system in different states'' [18,19]. The process of the cell development can be viewed as the system moving from one valley (primary or stem cell phenotype) through bifurcation to another valley (differentiated cell phenotype) on the potential landscape. And the transdifferentiation process can be viewed as the system escaping from one stable differentiated valley to another differentiated valley through certain paths on the potential landscape shown in Figure 1(A). The differentiated cells (S A ) can switch to another lineage cell type (S B ) through an explicit pluripotent stable state (S C ). Indirect transdifferentiation mechanism which requires an initial dedifferentiation step S A ?S C ?S B shown in Figure 1(A). It illustrates a differentiated cell (S A ) reprogrammed back to a pluripotent state (S C ) with less differentiated, and then can be re-differentiated to another type of differentiated cell (S B ) [3,5,6]. This is a possible strategy of pluripotent lineage reprogramming while the enhancement of efficiency is required. The underlying process is a transdifferentiation involving a stepwise dedifferentiation. In addition to indirect transdifferentiation, there is another lineage reprogramming approach: the direct transdifferentiation mechanism as S A ?S O ?S B shown in Figure 1(A). Direct transdifferentiation is a mechanism of converting one type of differentiated cells to another type of differentiated cells without undergoing through a pluripotent state or progenitor cell type. The differentiated cells (S A ) down regulate their own cell-specific genes (X 1 ) and activate the target cell-specific genes (X 2 ), thus they can switch to another lineage cell type (S B ) through an explicit intermediate stable state (S O ) or a series of indeterminate states [3][4][5]8,9]. In our study, the intermediate state is defined as an intermediate stable state with low or medium pluripotency and having very low expressions of the differentiation mark genes, while a series of indeterminate states are defined as a series of unstable states with low or medium pluripotency and very low expressions of differentiation mark genes in the course of lineage switching. Sridharan et al [20] showed that partially reprogrammed cells as an intermediate stage of the reprogramming process can switch to the completely reprogrammed iPS state. Thus the states of partially reprogrammed cells may exist along the paths from a differentiated state S A or S B to iPS state S C . The research by Mikkelsen [21] showed that partially reprogrammed cells can be trapped at a common intermediate state. Thus the states of partially reprogrammed cells may exist along the paths from a differentiated state S A to another differentiated state S B through an intermediate S O or indeterminate states. These intermediate state and indeterminate states may have certain expressions of stem cell marker genes and thus can be viewed as partially reprogrammed cells. This is supported by the observation that fibroblast cells specific genes are efficiently silenced and the embryonic reprogramming is not fully induced in partially reprogrammed cells [20]. We believe that different experimental and environmental conditions can lead to quite different results and change the topological structure of the potential landscape [20,21]. The partially reprogrammed cells may be trapped in certain regions in the gene expression space.
In this study, we term direct transdifferentiation as transdifferentiation and indirect transdifferentiation requiring an initial dedifferentiation or reprogramming step as dedifferentiation. The goal of regenerative medicine can potentially be realized through the processes of differentiation, dedifferentiation, reprogramming and transdifferentiation [4]. Here we use cell type switchings short for the terms ''differentiation, dedifferentiation, reprogramming, and transdifferentiation''. Recent advances have shown that there are three possible driving forces for cell type switchings: (1) Stochastic Fluctuations. Cells choose their pathways of differentiation stochastically in the process of development without apparent regards to environment or history [22]. Some studies in cell development reveal that intrinsic stochasticity is an important mechanism for development [22]. The extrinsic fluctuations are also expected to play a role in cell development. Thus the fluctuations can be a driving force for the processes of cell type switchings. (2) Gene Regulation. cell type switchings can be achieved by the change of regulation strengths of their lineage specific genes in many studies [6,8,9,14,15]. (3) Induction. Lineage specific cells can be reprogrammed to a pluripotent state through over-expressions of some defined transcription factors [23,24]. Transfection of certain cell specific genes into the primary cells, and over-expressions of the target lineage specific genes as well as certain stem cell-associated genes can induce the processes of cell type switchings.
Given the three driving forces for cell fate decision making, it is still challenging on how to quantify the processes of cell type switchings on the landscape, and how to connect them to experiments. These processes of cell type switchings are controlled by their underlying gene regulatory network. The lineage-specific transcription factors play a critical role in the processes of cell type switchings. In this study, we explored a simple cell differentiation network module with autoregulation and mutual antagonism between transcription factors (lineage-specific genes) [15,17], which exists in many cell differentiation processes, shown in Figure 1(B). The lineage-specific genes can strongly instruct the cellular lineage choice. The circuit is composed of a pair of self activating autoregulation and mutual inhibiting cross-antagonism cell-specific genes X 1 and X 2 [15,17]. In iPSC or ESC (embryonic stem cell), pluripotent genes are often highly expressed, and most lineage related genes are off. However, there are examples of gene regulatory circuits with the same architecture in our study which control binary decisions at branch points of cell differentiation in multi-potent cells. Such mutual antagonism gene circuit modules (where the self activation can also be indirect) in binary branch points of cell lineage commitment can often be found. A lot of studies have explored the primed multipotent common myeloid progenitor (CMP) can differentiate to either myeloid cell or erythroid cell in blood cell formation by mutual antagonism interaction of transcription factor gene GATA1 and PU:1 shown in Figure 2(A) [25,26]. GATA1 and PU:1 are both self-activated. In the genetic regulation of the inner cell mass/trophectoderm lineage decision, Oct4 represses expression of Cdx2, and Cdx2 represses expression of Oct4 to allow the segregation of inner cell mass and trophectoderm lineages [27,28]. Oct4 and Cdx2 are mutual inhibited and self-activated [27,28] shown in Figure 2(A). In the genetic regulation of the epiblast/primitive endoderm lineage decision, antagonism between Nanog and Gata6 results in segregation of primitive endoderm and epiblast within the inner cell mass [27,29,30] shown in Figure 2(A). Nanog and Gata6 are also both self-activated [30]. These three circuits all can be viewed as X 1 and X 2 in our network.
We will study this key network module to uncover the underlying functional mechanisms of cell type switchings. The phase diagram in Figure 1(C) suggests that the system can have five different phase regions, each of which has different underlying landscapes with different distribution of valleys. Furthermore, we show how stochastic fluctuation, gene regulation and induction induce the cell type switchings. The potential landscape and flux both direct the processes of cell type switchings. Probability flux provide a curling force breaking the detailed balance and lead the biological paths of cell type switchings to be deviated from the paths obtained by steepest descent gradient of the landscape. The forward and backward paths of cell type switchings are irreversible, without passing through the saddle point. Furthermore, the flux can become the main driving force for cell type switching when the landscape is not biased towards the specific processes [16,31]. Barrier height and dynamic transition speed are used to quantify the global stability of the landscape topography. The stability here represents the ability for a cell to stay at a certain cell type state against certain fluctuations. In practice, the fluctuations in some cases maybe small but never zero. We uncover and classify four mechanisms of cell type switchings: super-critical pitchfork bifurcation, sub-critical pitchfork bifurcation, sub-critical pitchfork with two saddle-node bifurcation, and saddle-node bifurcation.

I. The model of cell fate network
We start with gene circuit module for typical differentiation. The gene regulatory circuit for cell fate decision has two mutual repression and self-activation lineage-specific transcription factors: X 1 and X 2 shown in Figure 1(B). It is more complete to consider three or more gene system. But the challenge is that a network with more genes requires more parameters to describe and therefore much bigger search space to explore exhaustively for uncovering the underlying mechanisms. Furthermore, with more genes, it is more difficult to visualize the results. The two gene system we considered is the simplest to exhaustively and effectively explore the underlying mechanism in parameter space [15][16][17]25]. We would like to use this model to explore the basic underlying mechanisms. The dynamics of this circuit is described by a set of two-variable ordinary differential equations below, with the rate of expression change for these two genes: where x 1 and x 2 are the time-dependent expressions of the two cell-specific transcription factors X 1 and X 2 [15,17,25]. Parameter a 1 and a 2 are the self activation strength of the transcription factors X 1 and X 2 respectively. b 1 and b 2 are the strength of the mutual repression for transcription factors X 1 and X 2 respectively. k 1 and k 2 are the first-order degradation rate for X 1 and X 2 respectively [15,17,25]. S represents the threshold (inflection point) of the sigmoidal functions, i.e., the minimum concentration needed for appreciable changes, and n is the Hill coefficient which represents the cooperativity of the regulatory binding and determines the steepness of the sigmoidal function. For simplicity, we do not include studies of all the different parameters of S and n in the main text. We included the studies in the supporting information. We show the phase diagrams for varying these parameters in Figure S1 in File S1. We can see varying these parameters can also lead to bi-stable states or tri-stable states and also the phase transitions. In the main text, the parameters for Hill function and degradation rate for X 1 and X 2 are specified as: S~0:5,n~4, and k~k 1~k2~1 :0 [15][16][17]25]. In this section, we assume the symmetric situation a~a 1~a2 and b~b 1~b2 .
Although the values of parameters can be different in organisms under different circumstance, the mathematical model here describes a simple yet representative motif gene circuit, and these values (S~0:5,n~4,k~k 1~k2~1 :0) are used in many previous studies [15][16][17]25].
1. The phase of cell fate network. To explore the dynamics under different conditions mimicking by different choice of parameters, we showed the phase diagram in Figure 1(C). If we can keep the mutual repression strength b fixed and the self activation a at various levels mimicking the actual developmental process where expression levels of transcription factor change [15](e.g. The expression level of transcription factor Klf 4 can be viewed as the effective self activation a at various levels mimicking the actual developmental process [32]. Because Klf 4 is not required for the maintenance of undifferentiated state of ES cells [32]. Furthermore, the expression level of Klf 4 decreases gradually after induced differentiation [32].), the cells are attracted to different differentiated and undifferentiated states. There are five regions in the parameter phase space in Figure 1(C). Region I with lower self activation a and mutual repression b has only one stable state S O with lower equal levels of the expressions of two lineage specific genes X 1 and X 2 shown in Figure 1(A). This is an intermediate state phase with lower lineage specific genes in the process of transdifferentiation [4]. Region II with higher mutual repression b and lower self activation a has two stable states shown in Figure 1(A): S A which represents the differentiated state with higher expression of X 1 and lower expression of X 2 , S B which represents another differentiated state with lower expression of X 1 and higher expression of X 2 . Region III with lower mutual repression b and relative higher self activation a has three states: S A ,S B and S O . Region IV with higher mutual repression b and self activation a has three states: S A ,S B and S C which represents a pluripotent state with medium equal expressions of X 1 and X 2 in the process of dedifferentiation which can also be viewed as the process of reprogramming. Region V with lower mutual repression b and higher self activation a has all the four stable states: S A ,S B ,S C and S O . By changing the parameters of self activation a and mutual repression b, we can induce the initial differentiated cell to another differentiated cell in region II through the region III or region I by transdifferentiation (the yellow solid line), or through the region IV by dedifferentiation (the red dash line). In regions II, III, IV and V, there also exist tansdifferentiation within each. We will explore the dynamics of gene regulatory network for cell fate decision making process resulted from three driving force of stochastic fluctuations, gene regulation and induction through the instructive changes in details via the corresponding landscape topography for cell development.

Super-critical and sub-critical pitchfork bifurcation
versus saddle-node bifurcation in cell fate network. We explored the bifurcation for cell fate decision network for different conditions. When mutual repression regulation parameters b~b 1~b2 increase with small self activation regulation a~a 1~a2~0 :2, the phase diagram has a super-critical pitchfork bifurcation which is a second order phase transition [33,34] shown in Figure 3(A). The solid lines represent stable fixed points while the dash lines represent unstable fixed points. We can see a stable state S O becomes an unstable state and splits into a pair of new stable states S A and S B at the critical point [33,35]. As the self activation regulation strength a increases, the phase diagram changes to a new form of sub-critical pitchfork with two saddlenode bifurcation which is a first order phase transition shown in Figure 3(B) as a~a 1~a2~0 :5. The initial state S O is mono stable at lower mutual repression b, then a pair of new stable states S A and S B (two saddle-node bifurcations) emerge at somewhere far away from the initial state S O as mutual repression b increases. After the critical point of sub-critical pitchfork, the center initial stable state S O at the center becomes unstable, only the two new stable states S A and S B are left in the phase space. Super-critical pitchfork bifurcation represents a type of ''second-order transition'' in physics [36]. The difference between super-critical pitchfork bifurcation and sub-critical pitchfork bifurcation is that: super-critical pitchfork bifurcation represents one stable equilibrium splits into two stable equilibrium and a unstable equilibrium while sub-critical pitchfork bifurcation represents two unstable equilibrium and a stable equilibrium merge into an unstable equilibrium. Thus super-critical pitchfork bifurcation differs from the sub-critical one in that two new stable equilibrium S A and S B , when they appear, already have a significant distance away from the middle stable equilibrium S O . But the two stable fixed points and the two unstable fixed points in sub-critical pitchfork with two saddle-node bifurcations are both symmetric in x 1 {x 2 {b three dimensional space, while they are not symmetric in x 1 {b two dimensional space shown in Figure 3(B). These two bifurcations shown in Figure 3(A) and (B) are similar to the picture described in Waddington's epigenetic landscape [12].
The phase diagrams shown in Figure 4(A), Figure 5(A) and Figure 6(A) are saddle-node bifurcations. A saddle-node bifurcation denotes a collision and disappearance of two equilibria rather than a pitchfork bifurcation [33,35]. The saddle-node bifurcation is a first order phase transition [33,34]. We can see that the initial valley S A does not split into new valleys as the description of Waddingtons epigenetic landscape (a pitchfork bifurcation) [35]. New valleys S B and S C or S O are born at somewhere far from the existing valley S A in the state space. It is anther way of creating or eliminating the valleys from the potential landscape besides a pitchfork bifurcation [35]. The cell moves to the new valley S C or S O and S B in sequence under fluctuations since its own valley disappears in another saddle-node bifurcation. We have already explored another form of bifurcation for cell fate network as self activation a decreasing with b~1:0 in our previous study [15,17]. The phase diagram was drawn on the intrinsic potential landscape as the black lines in Figure 1(D) which is a sub-critical pitchfork [33,34] at the phase transition point (a~0:774).
We would like to explore these mentioned non-equilibrium phase transition under fluctuations and gene regulation. We might monitor the expressions of the differentiation marker genes in time and obtain the correlation functions. The singularity of the selfcorrelation function indicates the first order phase transition (saddle-node bifurcation) and the continuity of that shows the second order phase transition [33,34]. Thus we might distinguish these mechanisms of cell type switchings. We will explore these mechanisms of four bifurcations through our potential landscape theory in details below.
3. Intrinsic potential landscape. We obtained the intrinsic potential landscape w 0 (see the section of Methods) with Lyapunov properties to quantify the global stability by solving the zero fluctuation limit Hamilton-Jacobi equation and the associated intrinsic flux velocity in the zero noise limit [37]. The population potential landscape of cell development can be obtained through the exploration of the underlying probability dynamics, by solving the Fokker-Planck diffusion equation (see the section of Methods) [15]. The population potential landscape U is related to steady state probability distribution P ss through {lnP ss under fluctuations. The intrinsic potential landscape is quantified at the zero noise limit while the population potential landscape is quantified under finite fluctuations. Both show the global properties of the cell developmental process. Although intrinsic potential landscape gives less information (only at zero noise limit) about the network than population potential landscape, it can be used to quantify the global stability due to its nature of being a Lyapunov function [37]. We can illustrate two-dimensional potential landscape (the coordinates x 1 and x 2 ) to one dimension. One dimensional cross section coordinate x links S A side minimum through S C middle minimum to S B side minimum. x represents the gene expression levels, xv0 shows gene X 1 is dominant while xw0 shows X 2 is dominant. If the self activation strength a decreases relatively slowly, relative to gene regulation in development, the potential landscape can be viewed as a succession of one dimensional potential slice. Figure 1(D) shows the intrinsic potential landscape for normal cell differentiation development process from pluripotent state (S C ) to differentiated states (S A and S B ) and the pluripotent reprogramming process from differentiated states (S A and S B ) to pluripotent state (S C ). We can see the intrinsic potential landscape w 0 can be used to quantify the Waddington's picture and has almost the same shape with the population potential landscape [15].
The red dash lines and the yellow solid line shown in Figure 1(D) schematically described the lineage reprogramming process: dedifferentiation and transdifferentiation, respectively. The dedifferentiation process shows that differentiated state S A follows a step backward to a pluripotent state S C and then is induced to re-differentiate to another differentiated state S B . While the transdifferentiation process shows that differentiated state S A converts directly to another differentiated state S B through certain intermediate stable state or not. Much work has been done on lineage reprogramming and progress has been made in manipulating the key regulator gene to convert cell lineages [3][4][5][6]8,9]. The understanding of the underlying mechanism is still challenging. We will discuss the possible mechanisms of these lineage reprogramming process in detail using this simple gene regulatory circuit.
We can see that when self activation a~a 1~a2 is strong with higher mutual repression b~b 1~b2~1 :0, the valley of the central pluripotent state S C is much deeper and the system is attracted to this valley shown in Figure 1(D). As the strength of self activation a decreases, the valleys of side differentiated attractors S A and S B become deeper while the central pluripotent state S C becomes weaker. When the strength of self activation a approaching to zero, the central state S C becomes a ridge and therefore it is not stable while the side states S A and S B become stable. This result of intrinsic potential landscape with global Lyapunov property of global stability shows the similar mechanism with the result obtained from exploring the population potential landscape [15][16][17].
In order to quantify the stability of each state from the potential landscape topography, we can apply barrier height to measure the relative weights between different stable states. We showed barrier height of intrinsic potential landscape versus the strength of self activation a in Figure 7A. We set Dw 0S A~w0S {w 0A and Dw 0S C~w0S {w 0C , where w 0S is the value of the intrinsic potential landscape at the saddle point between state S A and state S C , w 0A represents the minimum value of the intrinsic potential landscape at differentiated state S A while w 0C represents the value of that at pluripotent state S C . Barrier height Dw 0S C decreases as a decreases, and state S C vanished after the phase transition critical point a c1~0 :774, where the system transits from three stable states (S A ,S B ,S C ) to two stable states (S A ,S B ). It implies that the attraction of state S C becomes shallower. Barrier height of Dw 0S A increases first, then decreases. It shows that the attraction of the differentiated state S A (S B ) becomes deeper first, then becomes weaker after another critical point around a c2~0 :4. So differentiated states S A and S B at a c2 are more stable. Figure 7B  A cell is a non-equilibrium open system with exchanges of energy and information from the outside environment. This leads to dissipation which is determined by both potential landscape and flux. The dissipation can give another global physical characterization of the non-equilibrium system. Non-equilibrium system dissipates both energy and entropy in steady state, where the entropy production rate is equal to heat dissipation rate. The heat dissipation rate is formulated as HDR~Ð F : Jdx [13,[37][38][39][40], which increases first then decreases as self activation a decreases as shown in Figure 7(D). This indicates that larger area of the dominant probability flux leads to more heat dissipation because the system needs to consume more energy [37]. The system   II. The mechanisms of cell type switchings 1. Stochastic Fluctuations. The cell type switchings at a given stage of development with different symmetric self activation a at fixed mutual repression b. The stochastic or inductive cell development can often be influenced by the external environment. We showed the paths of state transitions in cell development on the intrinsic potential landscapes for different self activation a with fixed mutual repression b~1:0 due to stochastic fluctuations shown in Figure 8. We can see the green lines represent the reprogramming or dedifferentiation paths from differentiated state S A or S B to pluripotent state S C while the red lines represent the differentiation paths from pluripotent state S C to differentiated state S A or S B shown in Figure 8(A)(B) when self activation a is relative stronger and the system has three stable states. Its worth pointing out that a green path from differentiated state S A to pluripotent state S C connected to a red path from pluripotent state S C to another differentiated state S B can provide a possible mechanism of the process of dedifferentiation first and then redifferentiation shown in Figure 8(A)(B). We also showed that both the green and the red lines represent the transdifferentiation paths from one differentiated state to another differentiated state shown in Figure 8(C)(D) when self activation a is relative weaker and the system has only two stable states, just as a toggle switch. The intrinsic flux velocity ((J ss =P ss ) D?0 ) represented by purple arrows are perpendicular to the negative gradient of intrinsic potential ({+w 0 ) represented by the white arrows in Figure 8 (see the section of Methods).
The cell type switchings processes at a given stage of development with symmetric changing mutual repression b while fixing self activation a. We considered the potential landscape changing under fluctuations with varying mutual repression parameter b~b 1~b2 at a given state with fixed self activation a~0:9. Any given cell may take a completely different route back to their pluripotent state in principle. Certain sequence of stages can emerge in the process of cell type switchings [4]. In experiments, if there are several pathways, one can collect the statistics and find out the relative probabilities of each path, giving the quantification of the path weights. In modeling, path integral weights are  Figure 10. These processes are fluctuation or induction induced transition. The purple lines represent the paths from state S B to state S A while the black lines represent the paths from state S A to state S B [15,37]. We can see there are two dominant paths with the same color for transdifferentiation from a certain differentiated state to another differentiated state in each sub figures, one path is through intermediate state S O while the other path is through pluripotent state S C . We also found the two different colored development paths between each two states follow quite different routes. It is irreversible between the forward dedifferentiation and the backward dedifferentiation paths through the pluripotent state S C , and between the two transdifferentiation paths through intermediate state S O or without an explicit intermediate state. This illustrates the irreversibility of the developmental paths which can be verified from the ongoing and future dynamical experiments.
The path weight represents the probability of each route for cell type switchings. It can be used to predict the probability of different routes for cell type switchings. The path probability can be obtained by the action A(x) for cell development (See methods for details). We labeled A P C as the action of the path through state S C , and A P O as the action of the path through state S O . Figure 9(C) showed the logarithm of dedifferentiation path probability through state S C divided that of transdifferentiation through state S O decreases as mutual repression strength b becomes weaker. This showed that the dedifferentiation path probability through state S C decreases or the transdifferentiation path probability through state S O increases as mutual repression strength becomes weaker.
The purple arrows represent the direction of the probability flux J while the black arrows represent the direction of the negative gradient of population potential landscape {+U shown in Figure 10. We can see the flux is almost perpendicular to the negative gradient of the population potential landscape [13,37]. The dynamics of transdifferentiation and dedifferentiation processes are determined by both gradient landscape and probability flux. Probability flux provides a curling force breaking the detailed balance, and leads the system to stay at the non-equilibrium state. The gradient force attracts the system into stable valleys. The potential landscape and flux both direct the processes of cell type switching. Flux can lead a system to move on even a relatively flat landscape, e.g., the limit cycle attractor, thus ''flux-directed differentiation'' and ''down-hill-directed differentiation (Waddington)'' both can occur in cell development. ''down-hill-directed differentiation (Waddington)'' leads to the exponential waiting of barrier crossing while ''flux-directed differentiation'' gives a much more precise timing. Flux also can lead the biological paths of cell type switchings to be deviated from the paths obtained by steepest descent gradient, and the corresponding paths of cell type switchings are irreversible. We would like to point out additional flux can emerge from epigenetics of slow (non-adiabatic) transcription and translation regulations [41] often encountered in eukaryotic cells. The flux generated by the slow time scales can lead to the new mechanism of differentiation and reprogramming [31,42]. The competition of barrier crossing and slow binding can lead to optimal speed of cell type switching. [31,42,43].
It is worth noting that even though state S O disappears in Figure 10(C), there still exist transdifferentiation paths through a series of indeterminate states near (0,0) position. This provides the possible mechanism of two ways of lineage reprogramming. We labeled the saddle point between state S A and state S C as s 1 while the saddle point between state S A and state S O as s 2 . In    We also can explore MFPT by t*exp(A(x)) [44]. Importantly, MFPT is also useful to characterize stability of the network for changing the regulations represented by the self activation a and mutual repression b under a small but fixed fluctuations (during the regulation changes or induction) mimicking the real environments. Figure 9(D) showed MFPT along dedifferentiation and transdifferentiation paths versus mutual repression strength b. We can see that the transdifferentiation path through state S O becomes more preferred than dedifferentiation path through state S C , and MFPT becomes shorter for transdifferentiation path through state S O as mutual repression strength decreases. In other words, transdifferentiation process is easier (harder) and the dedifferentiation process is harder (easier) when mutual repression is weaker (stronger).

Gene Regulation.
Decreasing self activation a 1 and increasing self activation a 2 induce the transdifferentiation process from state S A to state S B with lower mutual repression strength b. The instructive change of landscape via varying regulation strengths is another important mechanism in action for cell development. Down regulating the lineage specific gene for initial primary differentiated cell and up regulating the lineage specific gene for final target differentiated cell can induce transdifferentiation or dedifferentiation. We explored this mechanism below with changes in decreasing self activation a 1 for gene X 1 and increasing self activation a 2 for gene X 2 .
Self activation strength can be set for describing the time evolution of the self activation regulation parameters as: a 1~a0 exp({l 1 t) [25] which continuously decreases in time (down-regulates cell specific gene X 1 for differentiated state S A ) and another self activation regulation strength a 2~a0 ½1{exp({l 2 t) which continuously increases in time (upregulates cell specific gene X 2 for target differentiated state S B ) in cell developmental process due to the influences of the regulations of other genes. l 1 and l 2 are the rates for the decrease of self activations a 1 and a 2 . We assumed the same value of l~l 1~l2~1 0 {4 for simplicity for the latter calculations. At this value of l, self activation strength a 1 and a 2 decrease relatively slowly compared with regulation dynamics of gene X 1 and X 2 . Thus the dynamics is a slow non-equilibrium relaxation process. a 0 is the scaled value of self activation a 1 and a 2 [25].
We explored the transdifferentiation mechanism below with decreasing self activation a 1 and increasing self activation a 2 with lower mutual repression strength b. Figure 4(A) shows the saddlenode bifurcation phase diagram for decreasing self activation strength a 1 with lower b~0:2 and smaller a 0~1 :0. Figure 4(B) shows barrier height versus decreasing self activation a 1 with D~0:005. We defined the saddle point between state S A and state S O as s 1 , and the saddle point between state S B and state S O as s 2 . Barrier height is defined as: Ba sij~Usi {U j , where U si is the potential value of the i saddle point, and U j is the minimum at valley S j . Barrier height can quantify the degree of global robust and stability at a valley. We can see the cell stays at the monostable differentiated state S A at the beginning of the transdifferentiation. As self activation a 1 decreases, an intermediate state S O emerges. Valley S A is much deeper than valley S O due to barrier height Ba s1A of valley S A being higher than that of Ba s1O . It means the differentiated state S A is more preferred and more attractive than intermediate state S O . The system is preferred to stay at state S A with gene X 1 being dominant. As self activation strength a 1 becomes weaker and self activation a 2 becomes stronger, the valley of state S A becomes shallower while the valley of state S O becomes deeper. Then, the valley of state S O is more attractive than that of state S A since barrier height Ba s1 A is lower than barrier height Ba s1O , and gene X 1 and X 2 are both at lower expressions. After state S A disappears, the cell is driven into intermediate state S O . As self activation strength a 1 decreases further, the other differentiated state S B emerges, and barrier height Ba s2B becomes higher than barrier height Ba s2O . Finally, the cell is forced into state S B . This process interprets the mechanism of transdifferentiation from state S A to state S B through an intermediate state The above results showed the dynamics at certain stage of transdifferentiation. We can also explore the continuous dynamics controlled by the set of equations below: where l~10 {4 and a 2~a0 {a 1 . The continuous time dynamics of down-regulating gene X 1 and up-regulating gene X 2 is shown in Figure 4(C) with a 0~1 :0,D~0:005,b~0:2 using Eq.2. We obtained the transdifferentiation paths on the four dimensional potential landscape. The purple path is from state S A to state S B while the green path is the reverse transition both through intermediate state S O . It implies that the system with small mutual repression strength b (large inhibition) prefers the transdifferentiation path through intermediate state S O . Although transdifferentiation process does not seem to occur naturally, it has been observed in many experiments. For example, the exocrine cells in adult mice can transdifferentiate into b-cells using defined factors for direct reprogramming without passing through a pluripotent state but through an unnatural intermediate state [4,8,9]. Figure 5(A) shows the phase diagram of saddle-node bifurcation under a 0~1 :0 and mutual repression b~0:3. Figure 5(B) shows barrier height versus self activation a 1 with D~0:005. We defined barrier height as Ba si~Us {U i , where U s is the potential value of saddle point s between state S A and state S B , and U i is the potential value at state S i . Figure 5(C) shows the paths and the landscape for continuous dynamics using Eq.2 with a 0~1 :0,D~0:005,b~0:3. We can see the cell stays at differentiated state S A with higher barrier height Ba s A at first, then the landscape valley tilts the cell from state S A to state S B , barrier height Ba s B becomes higher than barrier Ba s A and the valley of state S A eventually disappears. Finally, valley S B becomes deeper. The weights of these two valleys exchange at the end of transdifferentiation process [35]. This process interprets the mechanism of transdifferentiation from state S A to state S B directly without through a specific intermediate state but through a series of indeterminate states. This result can be used to explain the mechanism that the enforced expressions of C=EBP with endogenous PU:1 can reprogram B cell into macrophages [4,5]. B cell specific marker is CD19 while the macrophage specific genes is Mac{1. The gene regulatory circuit is shown in Figure 2(B). B cell commitment factor Pax5 can up-regulate many B cell specific genes (such as CD19). The macrophage commitment factor C=EBP can up-regulate many macrophage cell specific genes (such as Mac{1) and down-regulate B cell specific genes (such as Pax5) [4,5]. Transcription factor PU:1 is needed in the process of transdifferentiation. The gene PU:1 has the property of autoactivation. Mikkola's work indicated that C=EBP and Pax5 act in mutual antagonisms [5,45]. The dashed lines for the autoactivation indicate uncertainty in Figure 2(B). Thus we can reduce the gene regulatory circuit in to two markers of CD19 and Mac{1 similar as our mutual antagonistic and self activation X 1 and X 2 [4,5]. C=EBPs inhibit B cell commitment transcription factor (B cell-specific genes) which down-regulates B cell marker CD19 (X 1 ) in B cell, and co-activate macrophage specific genes which up-regulates its target marker Mac{1 (X 2 ) in macrophages. B cells pass through a series of indeterminate states with lower expressions of B cell-specific genes CD19 (X 1 ) and macrophagespecific genes Mac{1 (X 2 ), which does not seem to undergo an initial dedifferentiation [4,5]. Figure 11(A)(B) show the logarithms of MFPT versus barrier heights using the same parameters in Figure 4(B) and Figure 5(B) respectively. We can see the time spent from one state to another and barrier height have the relationship as: t*exp(Ba). It implies that the harder the system is out from one valley with higher barrier height, the longer the escape time is.
We also explored the behavior for the system when regulation is not symmetrical, not only for the case when self-activation strength a 1 is not equal to self-activation strength a 2 , but also for the case when self-activation strength a 2 is not changing synchronously with self-activation strength a 1 . In Figure S2 in File S1, we showed the potential landscape of continuous dynamics with self-activation strength a 2 set as a constant (a 2~0 :65) while the self-activation strength a 1 continuously decreases. The other parameters are diffusion coefficient D~0:005, mutual inhibition strength b 1~b2~0 :2. We can see the cell may stay at differentiated state S A at first since the basin of differentiated state S A is lower than differentiated state S B when self-activation strength a 1~1 :5, then the landscape basin tilts the cell from the differentiated state S A to the intermediate state S O , and the basin of state S A eventually disappears. Finally, the basin S B becomes deeper, and the system shifts from the intermediate state S O to the differentiated state S B . The green path is from state S A to state S B while the purple path is the reverse transition from state S B to state S A both through the intermediate state S O . In Figure S3 in File S1, we showed the potential landscape of continuous dynamics with self-activation strength a 2 set as a constant (a 2~0 :1) while the self-activation strength a 1 continuously decreases. The other parameters are diffusion coefficient D~0:005, mutual inhibition strength b 1~b2~0 :2. We can see the cell may stay at differentiated state S A at first when self-activation strength a 1~1 :0, then the cell shifts from the differentiated state S A to the intermediate state S O , and eventually the basin of state S A disappears. The green path is from state S A to state S O while the purple path is the reverse transition from state S O to state S A . Here, intermediate state S O may represent the partially reprogrammed cells.
Decreasing self activation a 1 and increasing self activation a 2 induce dedifferentiation process from state S A to state S B with higher mutual repression strength b. We assumed self activation a 1 and a 2 at relatively higher average scaled values with a 0~4 :0 and relative larger mutual repression strength b~0:5 to induce the initial cell undergoing through a balanced pluripotent state [24]. Figure 6(A) showed the saddle-node bifurcation phase diagram for decreasing self activation strength a 1 with a 0~4 :0 at different time. Figure 6(B) showed the barrier height versus the parameter a 1 with D~0:005. We defined the barrier height as Ba i~Ucmax {U i , where U cmax is a constant relative maximum value of population potential landscape and U i is the minimum value of population potential at valley S i . We can see the system begins with a monostable differentiated state S A with higher expression of cell-specific gene X 1 and lower expression of cell-specific gene X 2 . As parameter a 1 decreases, a saddle node bifurcation emerges, giving rise to another differentiated state S B with lower expression of cellspecific gene X 1 and higher expression of cell-specific gene X 2 . Initially, barrier height Ba A of valley S A is much higher than that of valley S B , thus valley S A is much more stable than valley S B . So the system prefers to stay at differentiated state S A . As a 1 becomes weaker and the corresponding a 2 becomes stronger, two self activations a 1 ,a 2 for two cell specific mutually exclusive genes X 1 ,X 2 are over-expressing balanced (relative higher expression), another stable pluripotent state S C with medium expressions of gene X 1 and X 2 emerges, and the potential landscape has three valleys. Valley S A quantified by barrier height Ba A is deeper than valley S C and valley S B quantified by barrier height Ba C and Ba B at the beginning of valley S C emerging. As self activation a 1 decreases and a 2 increases further, valley S C and valley S B become deeper while valley S A becomes shallower. Barrier height Ba C is higher than Ba A and Ba B at a 1~a2 . Therefore, the system with differentiated state S A shifts to under pluripotent state S C as a process of dedifferentiation. A recent experimental studies [24] proposed a model for the coupled pluripotency module (selfactivation of Oct4 and Sox2) and for the differentiation module with mutual antagonism between the MEs (mesendodermal) and ECTs (ectodermal) shown in Figure 2(C). MEs inhibit the activation between Oct4 and Sox2, then Oct4 can only activates gene MEs, and inhibits gene ECTs [24]. This process can be viewed as MEs have the effect of self activation. Thus, this module can be reduced to two mutual antagonism gene MEs and ECTs with indirect self activation as our gene regulatory circuit of X 1 and X 2 shown in Figure 2(C). It implies that higher self activation strength a 1 and a 2 being balanced can lead the differentiated cell back towards the pluripotent cell. As self activation a 1 keeps on decreasing and a 2 keeps on increasing, the valley of the other differentiated cell state S B becomes deeper than that of pluripotent cell state due to barrier height Ba B being higher than Ba C and Ba A . Eventually, the valleys of S A and S C disappear at their saddle-node bifurcation [35]. Thus the cell leaves the pluripotent cell state S C and is forced to enter into the other differentiated cell state S B . The results showed the mechanism of dedifferentiation and redifferentiation. This mechanism can be used to explain many studies of cell dedifferentiation process during tissue regeneration both in vitro and in vivo [6]. For example, Pax5 is essential for initiating B cell commitment and is continuously required to maintain B cell lineage commitment [6,7,45]. Pax5 deletion can convert committed B cells into hematopoietic progenitors with pluripotency [6,7,45]. It is partly similar as the circuit in Figure 2(B) if we substitute Mac{1 into other lineage specific genes. Pax5 deletion means down-regulating the B cell specific genes (such as BCs) as the effect of self activation a 1 . This gene regulation can lead B cells (S A ) to dedifferentiate to hematopoietic progenitors (S C ). Then these cells can re-differentiate to T cell, macrophage or granulocyte (S B ) under appropriate culture conditions, such as the T-cell-deficient circumstance to reconstitute T cell development [6,7]. The appropriate culture conditions can be achieved by up-regulating the target cell genes as the effect of another self activation a 2 .
The population potential landscape at different developmental stage of decreasing self activation parameter a 1 after the relaxation process to a steady state among x 1 and x 2 is shown in Figure 6(C) using Eq.2. The green line represents the dedifferentiated path from differentiated state S A to another differentiated state S B through pluripotent state S C . The purple line represents the backwards dedifferentiated path from differentiated state S B to another differentiated state S A also through pluripotent state S C . We can see the irreversible paths on the four dimensional population potential landscape due to non-zero flux. The dedifferentiated landscape and the paths can be quantitatively described for predictions.
Decreasing mutual repression strength b induces differentiation and dedifferentiation process from state S C to state S A (S B ) with certain self activation a. Figure 3(A) shows the phase diagram of super-critical pitchfork bifurcation under self activation a~a 1~a2~0 :2 while changing mutual repression strength b. We can see the potential landscape of continuous dynamics shown in Figure 3(C) using Eq.2 with D~0:005, self activation a~a 1~a2~0 :2 and decreasing mutual repression b is similar to Waddington's epigenetic landscape [12,35]. A cell valley can form from an undifferentiation state around S O . S O can be viewed as a stem cell state with lower expressions of differentiation gene markers while S C can be viewed as the stem cell with medium expressions of the stem cell markers [17,46]. When decreasing mutual repression strength b, the initial valley splits into two other valleys and the initial valley becomes a ridge [35]. The cell will choose one valley as its fate. Figure 3(B) also shows the phase diagram of another form of sub-critical pitchfork with two saddle-node bifurcation under larger self activation strength a~a 1~a2~0 :5 when increasing mutual repression b. The continuous potential landscape shown in Figure 3(D) using Eq.2 with D~0:005, self activation a~a 1~a2~0 :5 and decreasing mutual repression b is also similar to Waddington's epigenetic landscape [12,35] except the surrounding of the critical point. Around the critical point, there coexist three stable states S O , S A and S B . We also quantified the paths on the potential landscapes. The purple lines represent the differentiation paths from undifferentiation state to differentiation state while the green lines represent the dedifferentiation or reprogramming paths. We can see the paths are irreversible even in the pitchfork bifurcation due to the existence of flux. This mechanism can describe the autonomous cell fate specification [47]. Stem cells must fulfill two tasks of self-renewal and generation of differentiated cells. In symmetric cell division, each stem cell can divide to generate either two daughter stem cells or two differentiated cells symmetrically while in asymmetric cell division, each stem cell splits to one daughter stem cell and one daughter differentiated cell [48]. The pitchfork bifurcation in this study can represent an asymmetry event that a polarized mother stem cell splits into two daughter cell M or N with different expressions of X 1 or X 2 . If daughter cell M with a very low value of X 1 or X 2 , it might fall into differentiated state, while daughter cell N with relative higher expression of X 1 or X 2 still stays at the pluripotent state [35,48]. The asymmetric cell division usually occurs early in embryogenesis [35,47].
3. Induction of over expression. Cell fate is influenced by inductive stimulus from a group of surrounding cells [23,35]. Over-expressions of defined transcription factors can induce one cell type to another cell type which does not depend on gene regulations. This has been achieved in practice using over expression of stem cell marker transcription factors. In our previous gene circuit studies of cell fate decision making for stem cell differentiation and development [15,17], the two genes in the network are both differentiation markers. The idea is that the specific differentiation markers when imbalanced will give differentiation of one cell fate or the other (two side basins S A and S B ). A more balanced differentiation marker setup (between the two) will lead to iPS stem cell state (center basin S C ). Our theoretical work [15,17] has predicted the possibility of the seasaw mechanism (balance or imbalance) of reprogramming. That is over-expressing both the concentrations of differentiation marker genes in a balanced way can induce and force differentiated cells into iPS stem cells or pluripotent cells [15,17].
The two mutually exclusive differentiation markers MEs (S A ) and ECTs (S B ) shown in Figure 2(C) with balanced over-expressions of key transcription linage specific factors can induce the lineage cell into pluripotent state (S C ) instead of the stem cell markers Oct4 and Sox2 for pluripotency of reprogramming as a ''seesaw model'' [24]. Our theoretical work [15,17] has already predicted this possibility of expressing differentiation markers for reprogramming and the seasaw mechanism suggested in their work [24].
Significant efforts have been made towards the experimental converting fibroblasts (S A ) to cardiomyocytes (S B ) by induction of over-expressing key genes. It is reported that direct transdifferentiation can be achieved by over-expressing gene Gata4, Tbx5 and Mef 2c from fibroblasts (S A ) to cardiomyocytes (S B ) [11]. Gata4, Tbx5 and Mef 2c are core transcription factors during early heart development and can co-activate other cardiac gene expression [11]. So Gata4, Tbx5 and Mef 2c can be viewed as cardiomyocyte cell specific gene C M s which have self activation. Over-expressing gene Gata4, Tbx5 and Mef 2c (C M s) can transdifferentiate fibroblasts to cardiomyocytes not through a pluripotent state [11]. Another experiment showed that the indirect transdifferentiation can be achieved with an initial dedifferentiation from fibroblasts (S A ) through pluripotent precursor-Cardiac progenitor (S C ) by over-expressing some stem cell markers Sox2,c{Myc,Oct4 and Klf 4, and then be induced to cardiomyocytes (S B ) [10].

Conclusions
In this study, we applied our potential and flux framework to explore the mechanisms of cell developmental processes of differentiation, dedifferentiation, reprogramming and transdifferentiation. The potential landscape of two gene regulatory circuit shows that the system has four stable valleys at specific regulation regions, two differentiated state S A and S B , one pluripotent state S C , and an intermediate state S O . Our work provides a quantitative basis for explaining the mechanisms of the transition among the four states. Barrier height based on the population potential landscape or the intrinsic potential landscape can quantify the stability of the attractors and the efficiency of switching among the attractors. We can acquire the dynamical transition rate of the system from one valley of attraction to another by MFPT for escape and the dominant paths for dedifferentiation and transdifferentiation via the path integral method. We can see the paths of cell type switchings are irreversible due to non-zero probability flux.
In this study, we have discussed three driving forces: stochastic fluctuations, gene regulation and induction, which can lead to cell type switchings. The cell type switching driven by stochastic fluctuations is a spontaneous transition, gene regulation is much like a non-autonomous varying of time-dependent landscape, and induction is a condition of initial value re-setting process with no apparent paths. The fluctuations maybe small in some cases but never zero. When exploring the stochasticity, we used fixed set of the values of self activation and mutual repression regulation parameters a and b. We not only discussed the possibility of cell type switching through stochastic dynamics but also other two mechanisms including the induction and regulation changes. We also explored the different dynamics with different sets of the parameter a and b. For gene regulation, we varied the parameters a and b for regulating the cell type switchings. For induction, we did not change any parameters. Instead, we just gave the cell an initial set (condition) with over-expression of its lineage specific gene. We quantitatively investigated the mechanism of cell type switching through the induction without the change of the underlying landscape and through the changes in regulations leading to the changes of the underlying landscape topography. Furthermore, these two types of cell type switchings driven by gene regulation and induction are not spontaneous transitions only due to fluctuations, but a controlled process under either the changes in regulations with regulation-dependent potential landscape or the induction with fixed potential landscape.
We found that the topography of the global potential landscapes is strongly correlated to the self activation strength and the mutual repression strength of the transcription factors. Dedifferentiation can be induced by the core regulators of pluripotent genes using in iPS or the synergistic effect of lineage specifiers in specification of differentiated cells. We can adjust two self activation strength a 1 and a 2 to be relative larger to force the differentiated cell to a pluripotent cell with higher inhibition strength b, and then redifferentiate the pluripotent cell to our target differentiated cell type [24]. This process can be viewed as an initial epigenetic activation phase representing the redifferentiation after a temporal overexpression of pluripotent reprogramming factors to a pluripotent state [4,24,49]. Somatic cells can be transdifferentiated by temporal over-expressions of pluripotent reprogramming transcription factors. Transdifferentiation can be induced by down regulating the lineage specific marker gene (X 1 ) of the original differentiated cell (decreasing self activation a 1 ) while activating another lineage specific marker gene (X 2 ) of the final differentiated cell (increasing self activation a 2 ) at relative lower inhibition strength b, through an intermediate state or a series of indeterminate states. This process can be viewed as lineageinstructive transcription representing the induction of lineage specific gene for the target differentiated cells [4,49]. This gives us a new understanding that the topography of underlying potential landscapes in cell development dynamics determines the feasibility and efficiency of cell type switchings.
We also classified the mechanisms of pitchfork bifurcations depicted Waddington's epigenetic development landscape including super-critical pitchfork bifurcation, sub-critical pitchfork bifurcation, sub-critical pitchfork with two saddle-node bifurcation, and saddle-node bifurcation depicted the transdifferentiation landscape [35]. We uncovered a pitchfork bifurcation of Waddington's epigenetic landscape and the irreversible paths (caused by the non-equilibrium flux) between differentiation and reprogramming. We also uncovered the saddle-node bifurcation landscape. Saddle-node bifurcation can give the explanation of possible mechanisms of dedifferentiation and transdifferentiation processes and can further explain the irreversibility of the paths for differentiation, dedifferentiation, reprogramming and transdifferentiation processes as hysteresis loop even without the presence of the non-equilibrium flux. We noticed a special kind of sub-critical pitchfork with two saddle-node bifurcations also shares the certain features with saddle-node bifurcation (hysteresis loop) and certain features of pitchfork bifurcation (Waddington's landscape).
Importantly, we uncovered some novel mechanisms as a starting point to decipher the mysterious code of the cell type switchings. Our theory can be used to guide the designs of the differentiation, dedifferentiation, reprogramming and transdifferentiation processes.

Methods
Quantifying non-equilibrium potential landscape, flux, non-equilibrium thermodynamics and the paths Fluctuations exist widely in biological systems [13,[50][51][52][53][54][55]. The dynamics in noisy fluctuating environments can be formulated as: _ x x~F(x)zf. F(x) is the deterministic force, where x is the vector representing different concentrations in state space. f is Gaussian noise term and its autocorrelation function is vf(x,t)f(x,0)w~2D(x)d(t), where D(x) is diffusion coefficient matrix. We set D(x)~DG(x), where D is the diffusion coefficient representing the level of noise strength while G is the scaled diffusion matrix described the anisotropy phenomenon. We can explore the corresponding Fokker-Planck diffusion equation [56,57] for probability distribution P(x,t): LP=Lt{ + : (FP{DG(x)+P). In this study, we set G(x) as a unit matrix for simplicity. The probability flux J is defined as: J~FP{D+P. In steady state, the force decomposition is shown as: F~{D+UzJ ss =P ss [13,50,51].
We obtained the Lyapunov function w 0 as the intrinsic potential from the zero fluctuation limit Hamilton-Jacobi equation(HJE) [37,58]: F : +w 0 z+w 0 : G : +w 0~0 by a numerical method -level set method using the Mitchell's level-set toolbox [59]. The force decomposition in zero fluctuation limit is shown as: F~{G : +w 0 z(J ss =P ss )D D?0 . From the Hamilton-Jacobian equation above, we can obtain (J ss =P ss )D D?0 : +w 0~0 [37,51,60,61]. We also can obtain the mean first passage time from the following equation [56]: F : +tzD+ 2 t~{1. t i j represents the mean first passage time from state i to state j.
The path integral approach we used is shown as below. The path probability starts from initial state x i at t~0, and end at the final state of x f at time t. The path integral formula is shown as [15,44]: P(x f ,tDx i ,0)~ð Dx exp½{ : (dx=dt{F(x)))~ð Dx exp½{A(x)Ð Dx exp½{ Ð L(x(t))dt, where L(x(t)) is the Lagrangian and A(x) is the action for each path [15,44]. The path integral over Dx represents the sum over all possible paths connecting x i at time t~0 to x f at time t. The exponent factor gives the weight of each specific trajectory path. The probability from initial state to the final state is equal to the sum of all possible paths with different weights. Every dynamical path doesn't contribute to the same weight and each path is exponentially weighted. Therefore, the path integrals can be approximated with a set of dominant paths while the other subleading path weights can be neglected for their relative small values. We can find the dominant paths with the optimal weights through minimization of the action A(x) or Lagrangian L(x(t)). Thus, we can identify the optimal paths which give more contribution to the weight as biological paths or cell type switching pathways in our study.