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.


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: The state variable Xi represents the strength/concentration of the GRF involved in expressing the gene towards the i-th cell type. The parameters β>0, ρ>0, g≥0 and γij≥0 are the efficiency of GRF in expressing the corresponding gene, degradation rate, basal constitutive growth rate, and time-varying interaction coefficient associated with the inhibition of Xi by Xj, respectively. The kinetics of GRF auto-activation is a sigmoidal increasing function that is negatively influenced by the strength of the other GRFs. This model is originally proposed by Cinquin and Demongeot [10,13] but their model has symmetric parameters unlike Box Eq. 1 which admits asymmetric repression coefficient γij.
We assume that the goal of gene regulation is to maximize the strength of GRFs so that the outcome is moving towards the steepest canal and possibly towards the deepest valley in the epigenetic landscape, which is expected in cell-fate determination. Here we 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 Fig. 2a) by increasing the value of gene expression efficiency β ( Fig. 2b) or by adding external stimulation that augments the effect of the basal constitutive 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 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.
Mathematically, the sequential branching in the epigenetic landscape towards different cell types represents convergence to one of the equilibrium points (Supplementary 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 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 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].
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]. The oscillations generated by various oscillator motifs when combined may produce fluctuations similar to stochastic noise. However, note that the combined large and small oscillations are sometimes not entirely random noise, especially when the detected fluctuations are part of the GRN structure and not coming from stochastic sources [29,51].
There are theoretical perspectives that interpret cell-fate determination as a chaos-driven process [52]. Another perspective is that the oscillator motif of a larger GRN can be used for artificial transdifferentiation by generating oscillations that can prompt cell-lineage switching, similar to what stochastic fluctuations can do [29,53]. In fact, direct transdifferentiation between related cell types branching from one lineage can be more straightforward compared to dedifferentiation to pluripotency [54].
From our simulations, we formulate some conjectures. (i) The dedifferentiation caused by abnormal oscillators (e.g., aberrant repressilator-type network) play a role in the existence of cancer stem cells and mutator phenotype [30,[55][56][57][58]. Abnormal changes in the structure of GRF-GRF interaction, such as abnormal timescale factor and abnormal weakening of repression links, can lead to disease. Indeed, partially reprogrammed cells and excessive plasticity can cause cancer [30,54,59,60]. Moreover, it is also possible that these oscillators play a part in epigenomic reprogramming and influence transgenerational epigenetic inheritance [61]. Abnormalities in the GRF-GRF interaction could be passed-on to offspring. (ii) We can reprogram cells back to pluripotency by regulating the wiring of the GRN. This implies that there are no unique reprogramming factors, and we can reprogram cells using any regulatory factor as long as it can lower the "gravity" of the epigenetic landscape [54]. In this paper, we have shown theoretical predictions of novel developments in 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]. can be helpful in drug discovery [71].
Methods. In our simulations, we use the following differential equation model (see Box 1): We assume the following parameter values: β=1, ρ=0.05 and  (5) is used for finding relative optimum and is similar to the trait dynamics frequently used in evolutionary biology [72,73]. The trait dynamics in fitness landscape (a concept in evolutionary biology) is similar to the cell-fate determination dynamics in Waddington's epigenetics landscape [58]. 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 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] defined by This network entropy measures the potency of population of cells. Population of multipotent and pluripotent cells have high network entropy because they can proliferate to heterogeneous cell types.
Numerical method. The ordinary and stochastic differential equations are solved using Runge-Kutta 4 and Euler-Maruyama with 0.01 as step size, respectively. For supplementary mathematical discussions about the differential equation model (Box Eq. 1), refer to related literatures [10,18,35].

Notes on abstraction.
One of the advantages of the model in Box 1 is that it is straightforward, thus any peculiar dynamics can be clearly interpreted. However, the model is an abstraction of the cell differentiation process. Hence, we focus on the qualitative dynamics rather than on exact quantitative values. Likewise, the qualitative dynamics arising from the model are possible to arise in a more complex system that contains the minimal GRN (Fig.   1b) as a subnetwork.
The authors declare no competing financial interests.   loop that generates sustained oscillations. This repressilator-type network has auto-activation unlike the repressilator by Elowitz and Leibler [34].