Branching and oscillations in the epigenetic landscape of cell-fate determination

. The well-known Waddington’s epigenetic landscape of cell-fate determination is not static but varies because of the dynamic gene regulation during development. However, existing mathematical models with few state variables and fixed parameters are inadequate in characterizing the temporal transformation of the landscape. Here we simulate a decision-switch model of gene regulation with more than two state variables and with time-varying repression among regulatory factors. We are able to demonstrate multi-lineage differentiation at different timescales that portrays the branching canals in Waddington’s illustration. We also present a repressilator-type system that activates suppressed genes via sustained oscillations in a flattened landscape, hence providing an alternative scheme for cellular reprogramming. The time-dependent parameters governed by gradient-based dynamics regulate cell differentiation, dedifferentiation and transdifferentiation. Our prediction integrates the theories of branching and structural oscillations in cell-fate determination, which reveals key temporal patterns of cell differentiation and associated diseases, such as cancer. paper, we have shown theoretical predictions of novel developments in


Box 1: Mathematical model
We consider a minimal gene regulatory network (GRN) of the form shown in Figure 1b, where the number of nodes (n) is arbitrary. This GRN is a representation of the network that characterizes decision switches in cell-fate determination [10, 37,43,50,75]. For simplicity, we draw the network in such a way that each node represents a gene regulatory factor (GRF) involved in expressing a specific cell type/phenotype [10,13]. A node in this GRN represents either a specific GRF, or a coarse-grained subnetwork of multiple factors that can be treated as single GRF. An example of a GRN in the form shown in Figure 1b is the coarse-grained mesenchymal transcription network with PPAR-γ, RUNX2 and SOX9 as nodes (Fig. 1c) [18].
One of the simple models that describe the qualitative dynamics of the GRN with multiple GRFs is the following system of differential equations: apply a gradient-based optimization method to the time-varying repression strength γij to drive cell-fate induction in the direction of the nearby steep canal (see Methods). This method is governed by a timescale factor that declines through time [5,12].
The equilibrium points and sustained oscillations of this model lie on the hyperspace 1  , and (2) Investigating the effect of dynamic model parameters provides insight on the upregulation and down-regulation of a gene. For example, it is possible to reverse the fate of an initial condition that leads to the down-regulation of a gene (convergence to low-valued stable equilibrium point; see growth rate g (Fig. 2c). Enhancement of the maximal growth rate β+g, such as by introducing external stimulus, is a direct technique in activating genes [20][21][22][23][24]. Another technique to . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted January 30, 2015. ; https://doi.org/10.1101/007831 doi: bioRxiv preprint drive gene activation is the reduction in the degradation rate ρ (Fig. 2d) [25,26]. Extrinsic and intrinsic stochastic noise also plays significant role in the up-regulation and down-regulation of gene expression (e.g., Fig. 2e), as discoursed in various theoretical and experimental studies [18,[27][28][29][30][31][32][33].
Nevertheless, this basic model of gene expression is often not enough to explain the complexity of cell-fate determination. To advance the theory of gene regulation, it is necessary to include many GRFs because there are key dynamics that are not conceivable in a model with only one variable. The canalization and oscillations in Waddington's epigenetic landscape are observable in an ordinary differential equation model with at least two or with at least three state variables, respectively. The oscillation-induced transdifferentiation and dedifferentiation discussed in this study are only possible with GRNs involving three or more GRFs.

Results.
We observe two significant dynamics in our simulations. The first one is multilineage differentiation via sequentially branching developmental pathways. This sequential branching portrays the canalization in Waddington's landscape at different timescales. The pathways trailed by the differentiating cells depend not only on the structure of the GRN and parameter values but also on the initial condition (see Supplementary Fig. 1). The second one is flattening of the epigenetic landscape which eliminates the deep valleys, resulting in sustained oscillations. The GRN that generates the oscillations (oscillator) can boost the strength of a suppressed GRF in partially or terminally-differentiated cells. This oscillatory behavior is structural in nature (i.e., based on the structure of the GRN) and not due to noise nor delay. Therefore, we are able to integrate into a single theory the theories of branching and structural oscillations in cell-fate determination.
. CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted January 30, 2015. ; https://doi.org/10.1101/007831 doi: bioRxiv preprint Mathematically, the sequential branching in the epigenetic landscape towards different cell types represents convergence to one of the equilibrium points ( Supplementary   Fig. 1). The variations in the interaction coefficients drive changes in the topography of the landscape ( Supplementary Fig. 2) and eventually stabilize at different timescales, hence sequential branching become possible (refer to the timescale factor in Methods). The pathway bifurcates from one branch (a primed state) to multiple branches (lineages), which may further have sub-branches that culminate in branch endpoints (cell-fates); see There are cases where a stable equilibrium point vanishes and a stable limit cycle emerges, especially when the repressive reciprocal interaction among GRFs is asymmetric (Fig. 4a). This attracting limit cycle is generated by an oscillator that perpetually up-regulates suppressed GRFs and down-regulates dominant GRFs, resulting in fluctuating kinetics (Fig.   4b). One example of an oscillator is a repressilator-type network (Fig. 4c), which is similar . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted January 30, 2015. ; https://doi.org/10.1101/007831 doi: bioRxiv preprint but not exactly identical to the repressilator proposed by Elowitz and Leibler [34]. In this repressilator-type network, the strength of repression in one loop is stronger than the strength of repression in the reverse loop (asymmetric interaction). The sustained oscillations generated by this repressilator-type network can arise in a GRN with three or more nodes (odd or even number of GRFs; e.g., Figs. 5a and 5b). The dynamics of this oscillator can be illustrated by an epigenetic landscape with flattened topography, that is, there are no deep valleys in the route of the differentiating cells, and the cells are continually sliding in zigzag canals without endpoint.
The oscillations drive the strengths of the GRFs to have alternating positive and negative rates, whereas the gradient-based optimization method forces the dynamics towards positive rates only. Thus, oscillatory behavior is not optimal in the sense of differentiation towards cell types located at deep valleys in the landscape. We then expect that gradientbased dynamics that are persistent for some period result in damped oscillations (Figs. 5c and 5d). There are cases where the damped oscillations illustrate multipotency (or pluripotency depending on the GRN) which is represented by the equal probabilities of differentiating towards all the considered cell types (Fig. 5c). In some cases, the interaction coefficients vary with different timescales, resulting in partial differentiation and sometimes in the reversal of the status of the initially dominant GRF (Fig. 5d). Note that if the dampening of the oscillations is fast, the initial oscillations can be unnoticeable yet can still activate suppressed genes ( Supplementary Fig. 19).

Discussion.
Various studies have attempted to model the cell differentiation process, but there are still more to uncover in epigenetics. Further theoretical prediction and experimental validation are needed to fully explain cell-fate determination and reprogramming. Varying the efficiency of GRF in expressing a gene (β), the degradation rate (ρ), or the constitutive . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted January 30, 2015. ; https://doi.org/10.1101/007831 doi: bioRxiv preprint growth rate (g) is a straightforward technique in stimulating the activation or deactivation of a GRF and its corresponding gene [6,7,10,18,35]. However, regulating the repression strength of GRFs (γij) has not been explored, and we have shown numerical illustrations where variations in this GRF-GRF repression affect the qualitative behavior of the cell differentiation system (Figs. 3 and 5). We are able to replicate Waddington's model using a single set of equations with many GRFs involved. A GRN with only two nodes generally cannot describe the sequential bifurcation of canals and the oscillations in cell-fate determination. The different timescales involved in gene interaction influence the outcome of cellular regulation [5,[36][37][38][39].
One of the aims of this study is to spur more discussions on non-equilibrium dynamics and oscillations arising from asymmetric systems with multiple GRFs, which can broaden our understanding about the mechanisms of gene regulation. Reversal of the route from differentiated state to pluripotency is previously thought to be impossible but now many dedifferentiation techniques have been proposed and tested by experiments [9,40,41]. We propose another alternative technique for cellular reprogramming, which is by rewiring the GRN to have a repressilator-type network, possibly with the aid of external input and stochastic noise [18,42]. Exogenous inhibitors can be introduced to weaken the repression in one loop (Fig. 4c). Our numerical predictions can help design cellular engineering strategies for generating induced multipotent stem cells (or pluripotent stem cells depending on the GRN [18,43,44]) using the oscillations that can activate silenced genes. Several studies have discussed that oscillating GRF expression is indeed an attribute of progenitor cells [45][46][47], thus supporting our claim. However, note that in reprogramming back to pluripotency, we also need to assure activation of defined factors in the core pluripotency circuitry, such as Oct4, Sox2 and Nanog [18,40].
. CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted January 30, 2015. ; https://doi.org/10.1101/007831 doi: bioRxiv preprint The oscillator motifs (e.g., repressilators) which are part of a larger GRN contribute to the fluctuations observed in gene regulation dynamics. In fact, there are many types of oscillators [48][49][50] was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted January 30, 2015. ; https://doi.org/10.1101/007831 doi: bioRxiv preprint the theory of gene regulatory networks. One experimental approach to demonstrate our numerical predictions is by employing synthetic biology techniques [33,34]. In reality, the temporal transformation in the epigenetic landscape is due to multiple intrinsic and extrinsic factors. Here we only consider changes in the GRF-GRF interaction coefficient γij but we should not disregard that gene regulation consists of the interplay among many factors and processes. For example, GRF-GRF interaction can be regulated not only through γij but also through the modifications in the maximal growth rate (β+g) or through the degradation rate (ρ) [6,7,10,18,35]. Increasing the maximal growth rate or decreasing the degradation rate of a certain GRF enhances the steady state strength of the GRF, which in turn intensifies the repression of the other GRFs. Furthermore, the dynamics observed from empirical data combine the various effect of many parameters. For example, a decline in the strength of a GRF suggests various possible reasons, such as due to an increased degradation rate or due to an increased repression by an antagonist GRF. Hence, we need to interpret data by considering all possible factors. There are other mechanisms not explicitly discussed in this paper, such as spatial pattern of cell differentiation, chromatin remodeling and DNA methylation [62][63][64].
The time-dependent parameter governed by a timescale factor corresponds to the delay activities in gene regulation. Another way of investigating the effect of time delay in cell differentiation (e.g., delay-induced oscillations) is by employing a delay differential equation model (for example, see [65]). Differential equation models with time delay are generally considered as infinite-dimensional systems.
In summary, our simulations predict the following outcomes: First is the sequential branching of lineages in cell-fate determination that portrays differentiation from pluripotent state to transient states (lineage progenitors) towards specialized cell types. Second is the cellular reprogramming driven by oscillations generated by a repressilator-type network, . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted January 30, 2015. Methods. In our simulations, we use the following differential equation model (see Box 1): We assume the following parameter values: β=1, ρ=0.05 and . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted January 30, 2015. ; https://doi.org/10.1101/007831 doi: bioRxiv preprint (5)  The timescale factor is represented by exp(-εit) with decline rate εi, as described in various studies [5,12]. As time progresses (e.g., as cell matures), the timescale factor declines and the value of ui leads to equilibrium. In addition, the dynamic parameter ui is initialized with value ui(0)=0.001 for all i. Let σu=0 and σu=0.01 for deterministic and stochastic simulations, respectively. For simplicity, we approximate the partial derivative in Equation Quasi-potential, cell type probability and network entropy. We suppose that the quasipotential (Φ) of the epigenetic landscape is computed as (6) . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted January 30, 2015. ; https://doi.org/10.1101/007831 doi: bioRxiv preprint where the dt dX i is deterministic [6,17]. This quasi-potential is assumed to characterize the elevation of the landscape [6]. The basal constitutive growth rate g, which can also represent effect of external stimulus, is excluded in our computation of the quasi-potential.
The probability of differentiating to cell type i or the expected proportion of cells . To visualize the canalization in Waddington's illustration, we use three coordinate axes: time, quasi-potential, and cell type probability. We also compute for the network entropy [19,74]    . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted January 30, 2015. ; https://doi.org/10.1101/007831 doi: bioRxiv preprint

Competing financial interests
The authors declare no competing financial interests. Figure 1. A sketch of the epigenetic landscape, and a gene regulatory network (GRN) for cell-fate determination. . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

Figures
The copyright holder for this preprint (which this version posted January 30, 2015. ; https://doi.org/10.1101/007831 doi: bioRxiv preprint . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted January 30, 2015. ; https://doi.org/10.1101/007831 doi: bioRxiv preprint . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted January 30, 2015. ; https://doi.org/10.1101/007831 doi: bioRxiv preprint loop that generates sustained oscillations. This repressilator-type network has auto-activation unlike the repressilator by Elowitz and Leibler [34]. was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted January 30, 2015. ; https://doi.org/10.1101/007831 doi: bioRxiv preprint