From Understanding the Development Landscape of the Canonical Fate-Switch Pair to Constructing a Dynamic Landscape for Two-Step Neural Differentiation

Recent progress in stem cell biology, notably cell fate conversion, calls for novel theoretical understanding for cell differentiation. The existing qualitative concept of Waddington's "epigenetic landscape" has attracted particular attention because it captures subsequent fate decision points, thus manifesting the hierarchical ("tree-like") nature of cell fate diversification. Here, we generalized a recent work and explored such a developmental landscape for a two-gene fate decision circuit by integrating the underlying probability landscapes with different parameters (corresponding to distinct developmental stages). The change of entropy production rate along the parameter changes indicates which parameter changes can represent a normal developmental process while other parameters' change can not. The transdifferentiation paths over the landscape under certain conditions reveal the possibility of a direct and reversible phenotypic conversion. As the intensity of noise increases, we found that the landscape becomes flatter and the dominant paths more straight, implying the importance of biological noise processing mechanism in development and reprogramming. We further extended the landscape of the one-step fate decision to that for two-step decisions in central nervous system (CNS) differentiation. A minimal network and dynamic model for CNS differentiation was firstly constructed where two three-gene motifs are coupled. We then implemented the SDEs (Stochastic Differentiation Equations) simulation for the validity of the network and model. By integrating the two landscapes for the two switch gene pairs, we constructed the two-step development landscape for CNS differentiation. Our work provides new insights into cellular differentiation and important clues for better reprogramming strategies.

section_content: The canonical view of differentiation as an irreversible process has been largely reshaped since the emergence of induced Pluripotent Stem Cells (iPSCs) and other lineage conversions techniques in stem cell biology [1] [2] [3] [4] [5] [6] [7] [8]. The success in inducing a conversion between cellular fates raises several questions [9] : why is a stable mature cell type retrodifferentiable or convertible? Is there a universal principle that can explain cellular development, and is there a fundamental commonality shared by the processes of normal differentiaton, retrodifferentiation and transdifferentiation? What are then the differences among the three processes?

In fact, a first effort to find a general principle traces back to Waddington's pioneering work in embryogenesis which gave rise to his ''epigenetic landscape metaphor'' (Fig. S1 ) for development [10]. Here the landscape metaphor describes differentiation as a ''down-hill'' process, which is about a cell ''rolling'' down from the pluripotent hilltop (the embryonic stem cells) to the lower valleys (the terminal differentiated cells), with multiple bifurcations at the watersheds on the landscape [10]. This metaphor, apparently lacking physical basis in Waddington's time, has long been ignored by experimental biologists until seen a renaissance among them in the recent years [11, 12]. Theorists had revisited this problem at various times. About twenty years after the first revelation of Waddington landscape, Thom proposed the catastrophe theory to explain the branching process in biological system [13]. However, he failed to find a potential function to construct the landscape. Kauffman, in a perspective different from that of Thom, starting from the idea of complex gene regulatory networks (GRNs) proposed that cell types are attractors in GRNs [14, 15]. His work used an efficient mathematical tool -Random Boolen Network (RBN). In parallel, more detailed modeling approaches (like Ordinary Differential Equations or ODEs) were also increasingly applied in modeling gene regulatory circuits. However, detailed studies of differentiation remain scanty and theories from dynamical systems have been applied only recently to the analysis of gene regulatory networks in development, starting with a single binary cell fate branching process [16]. Subsequently, the proposal of ''sequential branching'' model for hierarchical determination of cell fates, implemented both as ODEs or Stochastic Differentiation Equations (SDEs), has led to insights of how gene network dynamics govern pancreas development [17]. The re-discovery of intrinsic stochasticity [18] of gene expression in mammalian cells as well as developments in the theory of stochastic process has led to a first formalization attempt of the ''arrow of time'' (timedirectionality) for cellular differentiation [19]. In parallel, several studies on constructing the potential landscape for different biological systems now finally begin to address Thom's problem [20, 21]. 

In nonequilibrium systems, an explicit potential function generally does not exist, and an intuitive solution is to relate some form of potential to the steady state probability distribution in stochastic system and decompose the force driving the system into the gradient part and the curl part. Using such an approach, Wang constructed the landscape for the cell cycle dynamics and proposed that the curl force is responsible for the dynamics of a limit cycle [22]. Of crucial importance, these theoretical developments led to the first formalization for the one-step of binary branching in Waddington developmental landscape, generating myriad implications in explaining development and reprogramming [23]. However, in these models dimension of the developmental process on the landscape is represented by a hypothetical change of a specific model parameter whose physical validity is not confirmed. Therefore, it is necessary to study the meaning of this operation. In addition to Wang's method, there exists other constructive and ad hoc methods for landscape construction [24]. Even before Wang's construction, Ao proposed a transformation of the SDEs to obtain a potential function for constructing the landscape, similar to finding an effective flux in the Helmholtz-Hodge decomposition [25] [26] [27] [28] [29]. Based on the intuition from Lyapunov theory, Bhattacharya developed a numerical framework for mapping the quasi-potential landscape of a two-gene system [24]. However, none of them has been used to construct the landscape based on a real gene circuits governing the development across two steps of cellular fate branching. 

Here we first quantified the dynamic landscapes for the canonical two-gene fate switch motif under different parameter changes and demonstrated that some change can correspond to a development process with its monotonic increase of entropy production rate (EPR). We then extended previous work to quantify the transdifferentiation paths in addition to the previously quantified differentiation and retrodifferentiation paths on the landscape by using the theory of least action path. We further applied our theory to cell differentiation in the central nervous system (CNS) as an example to construct the two-step developmental landscape. The CNS was chosen because of the extensive researches about neuronal cell differentiation, which have accumulated abundant experimental data on gene regulatory circuit suited for modeling [30]. Intensive theoretical efforts in developmental neurobiology also have persisted for the past thirty years [31]. A considerable proportion of those work have concentrated on signaling cascade and morphogenesis in the neural development. However, GRN has dominant influence when it comes to the differentiation from nervous precursor to three CNS cell types (neuron, astrocyte and oligodendrocyte). Thus here we mainly model gene expression dynamics underlying CNS differentiation, first through qualitative assessment of a core regulatory network, followed by comparing the stochastic simulation (which was also used to explore the optimal cocktails for reprogramming from astrocyte or oligodendrocyte to neuron) with microarray data to confirm the model and network's validity. By gluing the landscape for two fate-switch motifs, we finally constructed the two-step landscape which has never been done before to our knowledge. Our work provides new insights to general principles about the mechanisms of retrodifferentiation and transdifferentiation.
section_name: Introduction

section_name: Materials and Methods

section_content: Recently a large range of cellular fate-switch regulator pairs, characterized by their cross-inhibition and self-activation, are confirmed [32]. This common gene pair is increasingly regarded as a general network design for the cellular fate bifurcation [33] (Fig. 1a ). In order to introduce the developmental landscape idea and demonstrate how a landscape can formulate the global developmental (also the reprogramming) processes mathematically and physically, we studied the canonical two-gene module (Fig. 1a ) as our first step. To characterize dynamics of this two-gene network, we assumed that autocatalysis and mutual-inhibition of genes are independent (their effects are thus additive), and thus used the following standard kinetic model [16, 23] :

in each of the above two nonlinear equations the first term represents the contribution to rate of change in gene expression from autocatalysis where the transcription and translation processes are lumped into one step; the second term represents that of inhibition from another gene and the last term first-order degradation. The first and second terms are implemented in the form of Hill equation with n = 4, but do not necessarily imply the cooperation between regulators [34]. Basically, we set a 1 ~a2 ~a~1,b 1 ~b2 ~b~1,k 1 ~k2 ~k~1, and S = 0. 5 to ensure simplicity and symmetry of the model. The rate of change in gene expression can also be regarded as a force driving the gene network. In vector form, we obtain:

In order to obtain a global and intuitive understanding of the development process described by the above kinetic equations of two fate-switch genes, a developmental landscape under the concept of Waddington landscape is needed [23, 35]. Until recently, there are three different approaches available for landscape construction [23, 24, 26]. In the following we derived the formula for landscape construction from Wang's framework [23], which relates the quasipotential to the steady state probability of a biological stochastic system. The mathematical basis of two other methods [24, 26] and the corresponding two-gene landscapes are included in the supporting information (Fig. S2 ). 

In nonequilibrium system with Nw1 dimensions, it is not easy to obtain a potential function U which satisfies +U~F(x): Since the intrinsic and extrinsic noise exists universally in biological system, an intuitive solution is to relate some form of potential (quasi-potential) to the steady state probability distribution. The above equations under noise influence can be related to a continuous stochastic model, whose probability distribution evolution is determined by the diffusion equation. 

According to the probability conservation law (The probability (P(x 1,x 2,t)) change at point (x 1,x 2 ) equates to the sum of the flux J(x 1,x 2,t) flows in/out of that point), we obtain the following equation:

where probability flux is defined as J~FP{D+P based on the following Fokker-Planck (FP) equations:

where D is the diffusion matrix which is set as a constant (which is basically 0. 1) in our simulation for simplification. In steady state, the LPss Lt ~{+ : J ss ~0 but probability flux J ss itself doesn't need to be zero, therefore, F~J ss =P ss zD : +P ss =P ss ~Jss =P ss z({D : +({ ln P ss (x,t??))):

Following the Boltzmann law in equilibrium system, the steady state probability can be converted into the dimensionless quasipotential according to:

Then,

The force driving the nonequilibrium system is thus decomposed into two parts: the gradient term F g (relates to the steady state probability) and the curl term F c (relates to the divergence-free flux). 

Even though the potential in nonequilibrium or equilibrium system relates to steady state probability distribution in a similar way, the dynamics of the nonequilibrium system is also determined by the curl force F c in addition to the gradient force F g. The obtained quasi-potential is then used to map stage-specific and whole dynamic landscape (Fig. 1b, c & Fig. S3A, B ). Following Wang's method, the change of parameters is used to represent the dynamic process. By fixing parameters a=b=k=S for F 1,F 2, we calculated U, the quasi-potential for each point in the gene expression state space of x 1,x 2 from the above equation, and then constructed the stage-specific landscape, which is defined because the parameters in the model are constants thus we can relate the landscape to a dynamic stage. The stage-specific landscape was further manipulated to construct the whole dynamic landscape. We rotated anticlockwise the x 1,x 2 coordinates and the landscape 45u, followed by extracting the lowest U for each value at the new coordinate X (X ~x1 cos p 4 {x 2 sin p 4 ) which was then used to represent the potential for a specific stage. We then systematically changed one of four parameters to represent the whole dynamic process and quantified the corresponding stage-specific landscape. After combining a serial of that set of potential, we mapped out the whole dynamic landscape. The EPR (Entropy Production Rate) is used to explore the relationship between change of a=b=k=S and the development process (Fig. 2d, e & Fig. S3C, D ). An important characteristic of a self-organization system (including the living system) is that it subjects to the principle of maximum entropy production rate (MEPR), that is the EPR will increase as the system evolutes into more ordered states (corresponding to the process that cells develop into specific mature cell types) [36]. Therefore, we assumed that the monotonically increase of EPR for parameter change in the landscape construction can be related to a developmental process. The formula of EPR is defined as [37] :

where k B is the Boltzman constant. In the dimensionless condition, we can use k B T~1. For the dynamic landscape where the EPR increases monotonically along the change of parameters on the landscape construction, they are regarded as the ''verified development landscapes''. Different cell types can be associated with different valleys (attractors) on the landscape [14]. Under the principle of least action, the dominant paths between the stem cell type valley and terminal cell type valleys can be quantified [23]. These dominant paths are regarded as the biological development (paths from stem cell valley to two terminal cell valleys), retrodifferentiation (paths from terminal cell valleys to stem cell valley) and transdifferentiation (paths between two terminal cell valleys) processes (Fig. 2 ). The dominant paths were obtained by minimizing the following action of the paths, for which we used the simulated annealing algorithm [23, 38]. 

The theoretical details about how to derive the above equation can be referred to Wang's papers [23, 38]. To find the transdifferentiation paths, two searching methods were used. The first one was operated under the symmetrical change of parameter a 1 and a 2, while the other was operated under asymmetrical change of parameter a 1 and a 2, which means that we fixed a 1 and changed a 2 if we wanted to transdifferentiate the cell from x 1 high expressed cell type to x 2 high expressed cell type, or vice versa. 

We further studied the effects of noise over the landscape topography and the dominant paths by varying the noise ratio from 0. 1 to 1 (Fig. 3 ).
section_name: Two-gene Landscape

section_content: The commitment of neuron lineage depends on the lateral inhibition, in which the prefated neural progenitors inhibit the neurogenesis of the neighboring cells through the Delta{Notch pathway. Proneural genes, like Mash 1, in the downstream of this pathway, are crucial for neural induction [39]. Their inhibitors, such as Hes 5, on the other hand, are glia inducers [39]. However, the key regulators which are in charge of the astrocyte-oligodendrocyte fate switch have not been fully identified, excepting limited reports like decreasing of Scl leads to Olig 2 high expression, and results in the formation of oligodendrocyte. For the purpose of modeling, we need to construct a minimal network for CNS differentiation by integration of scattered evidences available. We focus on the master genes, and do not consider the signaling molecules and influences from cellular environment in our model, as stated in the section of introduction. By doing so, we collected genes responsible for CNS fate decision and lineage markers genes. In particular, we started the network construction from finding Mash 1, Brn 2, Zic 1, Hes 5 and Olig 2's upstream regulators or targets because they are the five factors used in the lineage conversion from fibroblast to functional neurons, implying their importance in the CNS differentiation [5]. 

The curated minimal network responsible for the differentiation of CNS is shown in Fig. 4a which is characterized by the coupling of two-step fate decision. References supporting this network are summarized in Tab. S1.
section_name: Minimal Network Responsible for CNS Differentiation

section_content: We model the network in Fig. 4a as previous work [17]. The formula of gene expression change rate in the SDEs (Eq. 11) for the neural differentiation network contains three terms. The first term describes the gene's interaction with other genes (or itself), which is implemented in the form of a customary Hill equation. To reduce the computational complexity for high dimension system while focus on the qualitative behavior, we minimize the number of rate constants and combine the activation and repression effects into a single term. The expression of positive regulator will appear in both of the denominator and numerator of the regulated gene's equation; while for a negative regulator, it will appears only in the denominator. We multiply their expression to model the cooperative regulations between different genes. a,a e are used to represent different levels of production rate. The second term of the equations reflects first order degradation. To capture the gene expression stochasticity, we also added a Gaussian white noise as the third term. The noise is defined as elsewhere,

where the first equation means that the ensemble average of the noise is zero and the second equation describes the independence of the noises between different time points, implying a Wiener process. The set of SDEs for the twelve genes in Fig. 4a is defined as:

Feedback from terminal cell to CNS progenitor gene Pax 6 :

The parameters of the above equations in our simulation are set as: c~0:02,n~4,g~0:1,a~4,a s ~2:2,a e ~3,k~1. Several other sets of parameters were implemented and similar qualitative results were obtained. To represent the negative feedback of the mature terminal cells back to the precursor gene Pax 6, for which there is functional but no molecular evidence, we enforced an exponential decrease of Pax 6 in our model (the first equation for the change of Pax 6). However, we obtained qualitatively similar conclusion when the differentiation process was modeled by negative feedback from terminal cells to the x 1 (the second equation for the change of Pax 6, Fig. S5A, B ). The equations for the two three-gene motifs (x 1 {(x 2,x 6 ) and x 6 {(x 7,x 8 )in Fig. 4a ) are the core of the CNS development system, which was further used to constructed the CNS developmental landscape. Verification of the model by SDEs simulation. After searching the GEO [40], three datasets were retrieved as the reference for the stochastic simulation (GDS3442, GSE10796, GDS2379 [41] [42] [43], which study neurogenesis, astrogenesis and oligodendrogenesis, respectively). For each dataset, we first normalized the expression levels of specifically chosen genes by their maxima in the time series. By doing so, we can quantitatively demonstrate the expression pattern during neurogenesis which can be used to compare the simulated gene expression dynamics (Fig. 5a & Fig. S6-7A ). We also compared gene expression level of marker genes during three different CNS cell type genesis (Fig. 5b & Fig. S6-7B ). We made single cell simulation (one run of simulation) with Eq. 10 and compared the simulation dynamics with that of the corresponding genes from microarray datasets (Fig. 5c, Fig. S5 & Fig. S6-7C ).
section_name: Model the GRN Dynamics of CNS Differentiation

section_content: It is impossible to visualize a landscape of minimal CNS differentiation network which has twelve genes. Therefore we made a reasonable simplification: we mapped the landscape based on gluing two binary fate-decision modules' landscapes since they are the master regulatory pairs responsible for CNS cell types' fate commitment and the downstream genes don't have impact on the upstream genes due to the fact that we ignored the feedback from the downstream genes in our simplified model. Each fate-switch landscape was constructed for the corresponding three-gene module in a similar way as the previous two-gene developmental landscape construction (the remaining nine dimensions were used as the parameters of the module's equations). 

In constructing the whole landscape, the expression level of Pax 6 (x 1 ) was used to represent the development stage (Fig. 6a ), consequently it can be regarded as the bifurcation parameter of binary fate commitment. We first visualized the bifurcation diagram for the motifs of Mash 1{Hes 5 and Scl{Olig 2 and identified the ranges of bifurcation parameter (x 1 ) for the first and the second bifurcations (Fig. 6a ). The range of expression level of Pax 6 for the first bifurcation is about from 4. 5 to 2. 0 while the second bifurcation is about from 2. 0 to 0. 9. The cut-off values of 2. 0 and 0. 9 are the values where the first and the second bifurcations reach the most remarkable region (Line I and II on Fig. 6a ), respectively. Therefore, we assumed that at Pax 6~2:0 the neuron matures and glia forms from which the second-step differentiation starts while at Pax 6~0:9 the astrocyte and oligodendrocyte mature and the second-step differentiation stop. In order to construct the two-step CNS developmental landscape, we then mapped the first and the second-step switch genes' landscape according to the parameter ranges of bifurcation. Specifically, we used the decrease of expression level, from 4. 5 to 2. 0, of Pax 6 or x 1 (similar to the parameter a in two-gene landscape) to represent the first step differentiation, then calculated the quasi-potential of Mash 1 and Hes 5 state space at different x 1 by solving corresponding FP equation and finally mapped the firststep developmental landscape (Fig. 6b ). Similarly, we used the range from 2. 0 to 0. 9 of x 1 to represent the second-step differentiation. We then found the attractor states of Hes 5 high expression at each x 1 value and substituted the expression value of Hes 5 at each attractor into the equations of the second-step fate determinants (Scl and Olig 2) and finally mapped the second-step developmental landscape after solving the FP equation (Fig. 6c ). 

We used an ad hoc coordinate system to integrate the two landscapes into one. When the ordinate coordinate is between 4. 5 and 2. 0 (expression level of Pax 6), the X coordinate corresponds the transformed (similar to the two-gene landscape transformation) gene expression of Mash 1 and Hes 5 and the corresponding region (Region A in Fig. 6d ) represents the first-step landscape. While the ordinate coordinate is between 2. 0 and 0. 9, the left part of X coordinate corresponds to the extending value of the previous coordinate while x 1 ~2:0 (Region B in Fig. 6d ). The right part corresponds to the compressed Scl and Olig 2 coordinate system (Region C in Fig. 6d, see explanation below). 

In order to glue the first-step development landscape (region A in Fig. 6d ), the extending neuronal landscape (region C) and second-step development landscapes (region B), three specific processes are involved. Firstly, the attractors on the joining of the two landscapes (while a~2:0) are aligned. In practice, the average height of the second-step landscape (including that of the joining attractor) can be adjusted through the size of mesh used in calculating the underlying probability with smaller size a larger average probability (P) per unit area and thus lower height (u~{ln P) for the landscape (The probability used in calculating the landscape amounts to 1). By adjusting the size of mesh, we obtained a similar potential for the joining attractors. Secondly, while fixing the joining attractor, the second-step landscape can be stretched or compressed horizontally to maximally matching the boundaries near the two joining attractors since the X-coordinate is ad hoc and that of first and second-step doesn't have direct relationship. Lastly, for the ''bumps and cliff'' remains between the joining of regions A, B and C, we first cut off parts of the boundary of the second-step landscape with high altitude and then used interpolation to create a smooth transition between different regions. This is warranted for that the boundaries with very high altitude have very low probability (which in reality is possibly of nonexistence) thus can be disregarded and kept as spaces for the interpolation. On the other hand, the barrier between different valleys obtained after interpolation is enough to describe the stability of cell types and the difficulty of transdifferentiation between neuron and astrocyte (or oligodendrocyte).
section_name: Central Nervous System Developmental Landscape

section_name: Results

section_content: We generalized Wang's framework to further quantify the landscape under other parameter change and to explore whether or not the change can relate to a development process. 

We found that, while fixing other parameters, varying of parameters: b (a~0) (the measure of inhibition), k (the measure of degradation strength), S (the reaction threshold constant), in addition to parameter a (the measure of self-activation strength) in Wang's case [23], we can map similar landscapes (Fig. 1a, b and Fig. S3A, B ). The b-change landscape is different from others in that the landscape (Fig. S3A ) undergoes a bifurcation from one attractor to two attractors analogous to type I bifurcation in ref. [16]. Previously Wang's work [23] demonstrated that the change of parameters leads to the bifurcation. Is this reasonable to represent the developmental process on landscape construction? Do they have certain physical relevance with the development process? We were thus motivated to explore the relationship between the change of parameters and the entropy production rate (EPR, Eq. 8). In Fig. 1d, e, we show that the EPR increases monotonically along the corresponding change of a=k in constructing the landscape. Therefore, we demonstrate that the change of a=k in the landscape construction may be applicable to characterize the developmental process. When parameter b increases from 0 to 1. 5, however, the EPR does not always increase monotonically (Fig. S3C ). The starting non-monotonic region may offer a physical explanation on the idea that type-I bifurcation in ref. [16] is possibly not suitable for cell differentiation. Interestingly, the EPR decreases monotonically along the change of parameter S from 0 to 1. 5 (Fig. S3D ). We speculated that this decrease of S may correspond to the carcinogenesis because that, according to the MEPR, the system will evolve from terminal expression states to a promiscuous expression state which is a well-known phenomenon in the carcinogenesis (Fig. S3B ).
section_name: Two-gene Developmental Landscape

section_content: The center valley (A in Fig. 2a ) on the verified development landscape corresponds to stem cells while the side valleys (B and C in Fig. 2a ) correspond to two different terminal cell types. We quantified the dominant paths between valleys of different cell types. Similar to previous results [23], we first found that the development (paths from center valley A to marginal valleys B or C in Fig. 2a ) and retrodifferentiation paths (paths from marginal valleys B or C to center valley A) over the landscape are distinct and irreversible. We then extended previous work [23] to quantify the dominant paths between two terminal cell types on the landscape. We assume them as biological transdifferentiation paths (paths between two marginal valleys B and C). By implementing the two searching methods introduced in the section of Material and Method, we reached a consistent dominant transdiferentiation path between two mature cell types (a~0), which indicates that the optimal approach to transdifferentiate needs to maintain the parameter staying at a~0 during the whole process. The transdifferentiation path is direct and just transverses from the starting valley on the landscape and directly moves towards the target valley without firstly moving towards the stem cell valley (Fig. 2 ). Moreover, the corresponding dominant paths between valleys B, and C over the developmental stage-specific landscape (parameter a is fixed to 0) are also nearly identical (Fig. 2b ). These results drive us to speculate that the trandifferentiation may be reversible. However, at certain condition, we found that the transdifferentiation paths over developmental landscape will gradually move forward to the stem cells at first, then downward to the target cell type (data not shown). Applying these discoveries into the biological situation, it would imply that the optimal transdifferentiation strategy may depend on the source and target cell types as well as the the microenvironment where they locates in: some conversions need to firstly converse into the precursors then to the target cell types, while in extreme situation (here a~0 or no autocatalysis of gene expression) a direct phenotype conversion would be optimal. 

We then varied the intensity of noise for whole dynamic landscape and found that the shape of landscape and the direction of dominant paths is apparently influenced by the intensity of gene expression noise (Fig. 3 ). When the noise of gene expression increases, the landscape becomes flatter and the barrier between different valleys becomes lower, thus the cell types become less stable. Accordingly, we found that the reprogramming and normal development paths become increasingly straight when the noise of the system increases from 0. 1 to 1. Biologically, it indicates that the development, reprogramming processes could be accelerated when the gene expression stochasticity is magnified. In contrast to Waddington's perspective about a static development landscape, this result could provide a dynamic view of the landscape [10, 12, 35]. 

Similar to previous studies, we further found that the mean first passage time (MFPT) for escaping from one valley decreases as noise increases [22] (Fig. S4 ), in line with the idea about the influence of the magnitude of developmental fluctuations due to the stochasticity of gene expression in ref. [18, 44, 45] (The equation of MFPT and the corresponding graph is included in the supporting information).
section_name: Quantify the Dominant Developmental, Retrodifferentiation and Transdifferentional Paths on Two-gene Landscape

section_content: Previous construction of development landscape only demonstrated on a general two-gene module which doesn't consider multiple bifurcation of a specific biological system. In order to demonstrate the feasibility of the quantification of a landscape for a complex development process, we attempted to construct the developmental landscape for central nervous system development, starting from establishing a minimal gene regulatory network. In order to prove the validity of regulatory network as well as the model, we first implemented SDEs simulation and compared the simulated gene expression dynamics with microarray datasets. We also used the set of SDEs to simulate the transdifferentiation from oligodendrocyte/astrocyte to neurons and found that the simulation result matches conclusions drawn from the transdifferentiation paths of the two-gene landscape model. The two-step landscape was then constructed based on the confirmed SDEs of the coupled two fate switch circuits. 

The minimal CNS differentiation network is shown in Fig. 4a. We notice that this network consists of the same two motifs (Fig. 4b ), where there is a common upstream regulator positively regulating two mutual-inhibiting genes, a typical motif found in extensive network analysis in morphogenesis [46]. This motif is also used to explain the mechanism for stemness maintenance and reprogramming.
section_name: Core CNS Developmental Network

section_content: We found that, when assigning Pax 6 a high expression level and other genes low expression level at the beginning of the simulation, genes expressed sequentially in the network hierarchy in our stochastic simulation model, and captured the ordered developmental dynamics. Our result shows that the gene expression levels of Mash 1 and Hes 5 will first increase to a certain threshold, then follow a random switch between the expression of Mash 1 and Hes 5 in the opposite direction which respectively commit to neuron or glia cellular fate, indicating a pitchfork-shape bifurcation in development (Fig. S5A ). A similar bifurcation also appears for Scl and Olig 2 gene expression dynamics when the second step development occurs and cells commit to astrocyte or oligodendrocyte (Fig. S5A ). 

We found that our results qualitatively match the public microarray data. Taking single cell simulation of neuron differentiation as an example, the simulation agrees with public microarray data (GDS3442) [41] in which the progression from undifferentiated neuron precursor at E9:5 to neurogenesis at E11:5 and to maturation of neurons at E13:5 was measured (Fig. 5 ). In the microarray data, the expression of Pax 6 decreases constantly, while that of Mash 1 and Hes 5 first increases from E9:5 to E11:5. Then, Hes 5 decreases dramatically from E11:5 to E13:5 while Mash 1 decreases only slightly in this time period. Our simulation shows the similar expression dynamics with a clear bifurcation between Mash 1 and Hes 5 (Fig. S5A ). Our simulation is further confirmed by the microarray data from E9:5 to E13:5, which shows that expression of markers of other cell types is much lower comparing to the expression of Tuj 1. The stochastic simulation of the gene expression dynamics for astrogenesis and oligodendrogenesis also qualitatively agrees with the corresponding datasets of microarray experiments, as shown in the supporting information (Fig. S6-7 ). 

We next explored the optimal ''cocktails''of genes to be ecotypically expressed and their order under manipulation for converting astrocyte/oligodendrocyte to neuron [17]. By testing a serial of simulation, we found that the conversion could occur with higher efficiency when we activated sequentially Mash 1, Pax 6, followed by inhibiting Scl or Olig 2 and Hes 5. We further found that oligodendrocyte/astrocyte will firstly dedifferentiate into CNS precursor then redifferentiate into neurons during our transdifferentiation simulation (Fig. S5C ). With strong enforced activation of Mash 1, a direct transdifferentiations into neuron occurs in the simulation (data not shown). Our model also suggests that the cell types can automatically inter-converse between three cell types under large noise influence (data not shown). All the above three results agree with conclusion from the two-gene landscape model.
section_name: Confirmation of the Minimal Network and Model for CNS Differentiation by Stochastic Simulation

section_content: Based on the above stochastic simulation confirmed model and network, we herein constructed the CNS developmental landscape. Of our knowledge, this is the first physical realization for constructing a two-step developmental landscape for a biological system. 

Fig. 6b shows the first-step developmental landscape for Mash 1 and Hes 5 where the center valley corresponds to the CNS precursors, the left valley corresponds to the neuron cell type and the right valley corresponds to the glia cell type. Fig. 6c shows the second-step developmental landscape for Scl and Olig 2 while the center valley corresponds to the glia, the left valley corresponds to the astrocyte and the right valley corresponds to the oligodendrocyte. Using an ad hoc coordinate system and gluing technique (see Material and method), we obtained the two-step development landscape by integrating the first and the second-step developmental landscapes (Fig. 6d ). Our two-step CNS developmental landscape is similar to the original Waddington landscape (Fig. S1 ) with a two-step sequence of bifurcations that generates multiple valleys. 

Under the same framework for the quantification of dominant paths over the two-gene landscape, we further explored the dominant paths for the development from the CNS progenitors to neuron/glia and from glia to astrocyte/oligodendrocyte, and the corresponding retrodifferentiation paths. Similar to the results of those from the two-gene landscape, the normal development paths on the CNS developmental lanscape are distinct from the retrodifferentiation paths. In addition, by using the searching method applied in the two-gene network (the Pax 6 corresponds to the parameter a in two-gene landscape), we found the dominant transdifferentiation paths on the landscape need firstly move back to the CNS progenitor then downward the target terminal cell types (data not shown). These discoveries are similar to the results from the stochastic simulation where the reprogramming cocktails of genes will firstly reprogram the source cell types into CNS progenitors then into the target cell types.
section_name: Central Nervous System Developmental Landscape

section_content: The studies on the two-gene landscape provide a few interesting insights into the mechanism of development, retrodifferentiation as well as transdifferentiation. 

We demonstrated that the development and retrodifferentiation paths are irreversible because the dominant paths don't follow gradient since there is an additional curl force determining the dynamics of the system [23]. Such an irreversibility of forward and backward paths may be a common phenomenon in the nonequilibrium systems, like the living organisms. This conclusion can be extended to embrace recent findings that the genomic gene expression dynamics during developmental fibroblast formation is distinct from that in the iPSCs induction [47]. 

Intriguingly, a direct transdifferentiation path without passing through the medial precursor states appears on the landscape at certain development stage (a~0). Meanwhile, in certain range of parameter a, the transdifferentiation path will firstly move toward the stem cell valley and then to the target cell valley. These results may well associate with previous ideas (as also shown in our transdifferentiation simulation) that the transdifferentiation can happen through firstly dedifferentiating into precursors, then redifferentiating into the target cells. At extreme situation, a direct phenotypic conversion is applicable. It can be related to recent findings that, under additional signal induction, the direct conversion between fibroblast and neural progenitor has a higher efficiency than that by first reprogramming to the stem cells, followed by instructing to target cell types [48]. 

Our work demonstrates the importance of gene-expression stochasticity in development. With higher noise term in the SDEs system, the cell types can automatically inter-converse (data not shown); similarly, a larger diffusion matrix D will flatten the landscape, and lead to straight dominant developmental/retrodifferentiational paths and convenient transition between different valleys (cell types) over developmental landscape. These results reflect previous considerations that the reprogramming efficiency may be improved when biological noise processing mechanism (Proteasome or Wnt signal pathway, et al. ) is enhanced [49]. 

In practise of iPSCs, the somatic cells are induced by adding ecto-expressing factors (OKSM, Oct 4, Klf 4, Sox 8, Myc). Previous studies [23] for explaining the reprogramming by two-gene model (Fig. 1a ) with a parallel increase of parameter a 1 =a 2 in Eq. 1 didn't explicitly consider the ecto-expression of pluripotent genes. In that explanation all the two fate genes will finally return to low expression level. The enforced expression of OSKM in reprogramming cells may be safely explained by an improved perspective. The mutually-inhibited genes in the two-gene model are possibly a part of the complex differentiation network responsible for different somatic cell types fate commitment. They are further activated by the pluripotency network genes including the OKSM. A simplified model for this idea is shown in Fig. 4b as a three-gene model. In this model, we found that, when we enforce the expression of upstream regulator, the side terminal cells can be indeed reprogrammed into states of center stem cells (Fig. S8 ). 

Our work is different from that Wang's in considerable ways. We extended Wang's study to construct the b=k=S-change landscape and provided physical explanation why a=k change can be used to map a developmental process while b=S change cannot. We also quantified the transdifferentiation paths over the two-gene landscape and found they are identical and direct (parameter a doesn't change). In contrast to Wang's one-step landscape, we also constructed the first two-step Waddington development landscape based on a manual-curated network of CNS development. 

There are several relevance and differences between our CNS developmental landscape and the original Waddington landscape. Similar to the original landscape, our CNS developmental landscape has two-step bifurcation with multiple valleys. Each valley corresponds to a specific CNS cell type (Fig. 6d ). Like what shown on the Waddington landscape, the development on CNS landscape starts from the center valley and ends with the formation of three CNS cell type after two-step bifurcation. However, developmental landscape in Waddingtonnian sense is static [10] while our landscape is dynamic due to the noise influence. The dominant transdifferentiation paths over CNS developmental landscape do not correspond to a direct phenotypic switch as in the two-gene landscape. This may result from the fact that the bifurcation on the CNS developmental landscape increases from one stable steady states to two stable states (Fig. 6a ) while that of the two-gene developmental landscape decreases from three stables states to two stable steady states (Fig. 1b, c ). 

Previous studies model the irreversibility of development by certain feedback loops [49, 50]. Here we assumed an exponential decrease of Pax 6 to represent the feedback from all the mature cell types. But we also considered supposed direct feedback from marker genes (Tuj 1, Aldh1L 1 and Sox 8) to Pax 6 and obtained similar results. In hematopoiesis, the gene expression dynamics of multiple-step two-gene motifs is investigated by stochastic simulation and bifurcation analysis in ref. [49, 51]. However, we herein present a physical model for CNS development, thus provide a much fundamental and global understanding of cellular development which cannot be captured by stochastic simulation and bifurcation analysis.
section_name: Discussion

section_content: Figure S1 The original developmental/epigenetic landscape proposed by Cornad Waddington [10]. The original epigenetic landscape. Here the landscape concept refers to the recurring metaphor describing differentiation as a ''down-hill'' process of a cell ''rolling'' down from the high pluripotent hilltop (the embryonic stem cells) to the lower valleys (the terminal differentiated cells). (EPS) Table S1 Crucial genes and regulatory relationship for CNS development by manually curation. In this table, the Source regulator column shows the upstream genes and the Target gene shows their corresponding targets, both of which are crucial genes underling CNS development. The relationship between regulators and targets is listed in the third column, either positive (upregulation) or negative (down-regulation). If the regulation is proposed, we denote a ''FALSE'' item in the fourth column. For the confirmed regulation, we list the reference (Pubmid) for readers' check, while for proposed regulation, the listed reference evidenced the importance of the target genes in the formation of CNS. The data in this table is used to draw the Figure 4a in main text. (PDF)

Supporting information S1 This combined fourteen pages supporting file includes four sections about the mathematical basis for Ao's and Sbhattacharya's approaches to landscape construction, calculation of MFPT under different noise influence, the confirmation of astrogenesis and oligodendrogenesis simulations by comparing with microarray data and the proposal for explaining the stemness maintenance by the three-gene model, respectively. (PDF)
section_name: Supporting Information

section_name: Acknowledgments

section_content: We are indebted to helps from Sui Huang, Ping Ao, Qiang Zhang, Sudin Sbhattacharya, Joseph X. Zhou, Jin Wang, Yuan Wang, Jiemin Wong during various stages of the manuscript preparation. We also thank for the insightful comments from two anonymous reviewers which significantly improved the quality of this manuscript.
section_name: Acknowledgments

section_content: This work was supported by the National 973 Key Basic Research Program ( 2010CB945401, 2012CB910400 ); the National Natural Science Foundation of China ( 31171264, 31071162, 31000590 ); and the Science and Technology Commission of Shanghai Municipality ( 11DZ2260300 ). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

section_name: Author Contributions

section_content: Conceived and designed the experiments: XJQ SSD TLS. Performed the experiments: XJQ SSD. Analyzed the data: XJQ SSD. Contributed reagents/materials/analysis tools: XJQ SSD. Wrote the paper: XJQ SSD TLS.
section_name: Author Contributions