Dynamical systems theory of cellular reprogramming

In cellular reprogramming, almost all epigenetic memories of differentiated cells are erased by the over-expression of a few genes, resulting in regaining pluripotency, the potential for differentiation. Considering the interplay between oscillatory gene expression and slower epigenetic modiﬁcations, such reprogramming is perceived as an unintuitive, global attraction to the unstable manifold of a saddle, which represents pluripotency. The universality of this scheme is conﬁrmed by the repressilator model and by gene regulatory networks that are randomly generated and extracted from embryonic stem cells.

In the development of multicellular organisms, cells with identical genomes differentiate into distinct cell types.This cellular differentiation process has often been explained as balls falling down the epigenetic landscape, as originally proposed by Waddington [1]: balls start from the top of the landscape, and as development progresses, they fall into distinct valleys, which correspond to differentiated cell types.In modern biology, such landscapes are believed to be formed by epigenetic regulation, including DNA and chromatin modifications [2][3][4][5].For pluripotent cells, these modifications are small, whereas each differentiated cell type has a different epigenetic modification pattern [6][7][8][9].Cells with pluripotency, such as embryonic stem (ES) cells, are located in the vicinity of the first branching point into the valleys, because they can easily differentiate into different types of cells with just slight stimuli [10].
In 2006, a seminal study by Takahashi and Yamanaka reported that differentiated cells can regain pluripotency only by overexpressing few genes (so-called the four Yamanaka factors) without direct manipulations of epigenetic modifications.This was termed as reprogramming of induced pluripotent stem (iPS) cells [11].The reprogramming is often described as "climbing" the epigenetic landscape [12][13][14].This hypothesis, however, has two problems that still need to be addressed: (1) cells have many degrees of freedom, with expression and epigenetic modifications of many genes, whereas reprogramming manipulation involves only few degrees of freedom.How is it possible?;moreover, (2) if the initial pluripotent state is represented by the top of the landscape, it is not a stable point.Thus, how can reprogramming robustly make the cells head toward such an "unstable" state?
Theoretically, these issues should be resolved based on dynamical systems theory.The interplay between fast gene regulation and slow epigenetic dynamics shapes the epigenetic landscape, and differentiated cells are represented by different attractors [15][16][17].Therefore, upon the reprogramming operation, cellular states starting from different attractors first converge into a unique pluripotent state, which is not stable, from which states * kaneko@complex.c.u-tokyo.ac.jp move toward various attractors afterwards.At first glance, these requirements seem to be incompatible; an unstable state (e.g., repeller) is not attracted from different initial conditions.Hence, to satisfy these requirements, the pluripotent cell is expected to be represented at least by a saddle that is attracted from many directions and departs only along unstable directions (manifolds), which represent the cell differentiation process, leading to attractors of different destinations.To regain pluripotency by reprogramming, cellular states must be placed on the stable manifold of the saddle by common manipulations from different attractors.Such manipulation, however, would require fine-tuned control.In contrast, reprogramming is mediated by the overexpression of a few common genes across a variety of differentiated cell types.Therefore, some dynamical systems concept beyond just a saddle is needed.
A recent experiment provides some clues on this subject.Temporal oscillations in DNA methylation and corresponding gene expression levels are observed during cellular differentiation [18,19].In fact, gene expression oscillations have also been reported during somitogenesis and in embryonic stem cells [20][21][22], whereas its relevance to cell differentiation has been theoretically investigated for decades [23][24][25][26][27][28][29].Recalling possible significance of oscillatory dynamics, it is reasonable to consider that if there is an oscillation of fast gene expression around the saddle point of the slow epigenetic dynamics, global attraction to it from broad initial conditions may be attained beyond its stable manifold.As the oscillation dynamics are extended beyond the stable manifold of a saddle, global attraction to the vicinity of the saddle may be possible by taking advantage of the interplay between fast expression and slow modification dynamics.
Herein, we verified this possibility by using a dynamical system model with a gene regulatory network (GRN) and epigenetic modification.We consider a cell model in which the cellular state was represented by the expression x i and epigenetic modification level θ i for each gene i, with i = 1, 2, . . ., N .Gene expression dynamics, with faster time scales, are governed by GRN with mutual activation or inhibition by transcription factors [30][31][32][33][34], whereas slower epigenetic dynamics change the feasibility of gene expression, which follows the gene expression patterns.We assumed the epigenetic feedback reinforement, meaning that as more a gene is expressed (silences), the more feasible (harder) to express.This hypothesis was based on the experimental observations on the Trithorax (TrxG) and Polycomb (PcG) group proteins, two of the essential epigenetic factors for cellular differentiation [35].Specifically, we adopted: In Eq. (1a), gene expression shows an on-off response to the input by adopting the function F (z) = tanh(βz), whereas β = 40 [36].If J ij is positive (negative), gene j activates (inhibits) gene i, whereas J ij is set to 0 if no regulation exists.External input I i (t) is applied only during reprogramming manipulation to flip the expression of the gene i.For simplicity, I i (t) takes a constant non-zero value when gene i is overexpressed for reprogramming manipulation and zero otherwise.
In Eq. (1a), −θ i works as the threshold of the expression of the gene i, which represents the epigenetic modification status (when there is no epigenetic modification, it takes zero).Eq. (1b) represents the epigenetic feedback regulation.Following the experimental observation of positive epigenetic feedback [35,[37][38][39][40][41][42], we adopted this simple form as its specific form has not yet been confirmed [29,[43][44][45].Here, τ denotes the characteristic timescale for epigenetic modifications, which is assumed to be sufficiently larger than 1; the change in epigenetic modification is much slower than that of gene regulatory dynamics [46][47][48].
Recalling the relevance of oscillatory dynamics, we chose a GRN in which oscillatory dynamics were generated for appropriate θ i values (specifically at θ i ∼ 0).First, we adopted a repressilator model as a minimal model (see Fig. S4a [49]), consisting of three genes that repress the expression of the next gene in a cyclic manner [50].Specifically, we chose J 21 = J 32 = J 13 = −g = −0.4 in Eq. (1a).
The expression of x i in this model showed a limit-cycle oscillation when θ i was close to zero.Thus, for the epigenetic modification to change θ i following Eq.(1b), the states were differentiated into three fixed-point attractors {θ 1 , θ 2 , θ 3 } = {−1, 1, 1}, {1, −1, 1}, {1, 1, −1}, after first approaching a straight line θ 1 = θ 2 = θ 3 , as shown in Fig. 1a [51] (see also Fig. S5a[49]).In these fixed points, dθ i /dt = 0 were satisfied; that is, the differentiation of expression x i was embedded into the epigenetic modification θ i .Now, we considered "reprogramming."Starting from one of the differentiated fixed points, we added external input I i (t) to invoke the transient oscillation again (black dotted line in Fig. 1b).Later, I i (t) was set at zero.After reprogramming manipulation, they approached a line with θ 1 = θ 2 = θ 3 around the origin, and then deviated from the line to one of the three fixed points (Fig. 1b), in the same manner as the differentiation process.During this reprogramming process, memory of the differentiated states was erased.Once the oscillation in x was recovered, the approach to the straight line and deviation from it always followed (Fig. 1c).
As shown in Fig. 1cd, the straight line to which all trajectories converged agreed with the unstable manifold v u (Fig. 2a).Of course, attraction to the v u axis was natural if the initial conditions were restricted to the stable manifold for {θ 1 , θ 2 , θ 3 } = {0, 0, 0}.However, we observed an attraction toward the unstable axis over a wide range of initial conditions for {θ i }, which supports the oscillation of x i .Furthermore, the magnitudes of the eigenvalues for the stable and unstable eigenvectors were of the same order (Re(λ u ) = 0.31, Re(λ s1 ) = Re(λ s2 ) = −0.66,see Fig. S5c [49]).Thus, the reprogramming dynamics shown in Fig. 1b could not be explained just by the linear stability.
To elucidate whether the nonlinear effect suppresses the instability along the v u axis, we computed d θu /dt.As shown in Fig. 2c, d θu /dt was drastically reduced from the linear case.We also computed (d θs1 /dt, d θs2 /dt) for a certain θu value (that is, the flow structure in the ( θs1 , θs2 ) plane, sliced along the θu axis), which showed that θs1 = 0 changed from stable to unstable at θu = θth u (∼ 0.4) (see Fig. 2d and Fig. S6 [49]).Up to θu < θth u (∼ 0.4), θ i in ( θs1 , θs2 ) plane was attracted to the θu axis.By further increasing θu beyond θth u , θ i departured from the θu axis rotating in the ( θs1 , θs2 ) plane, leading to differentiation into three distinct fixed points.
Next, we considered how attraction to the v u from the ( θs1 , θs2 ) plane was lost at the θu = θth u .By considering θu as a parameter, the direction of flow in the ( θs1 , θs2 ) plane toward the v u was determined by the sign of ∂Θ s1 /∂ θs1 = ∂ xs1 /∂ θs1 − 1 ≡ νΘ s1 (we defined Θ s1 as a projection on v s1 ).As shown in Fig 2b, with the increase in θu , the bifurcation point from the limit cycle to the fixed point approached the θu = 0 line.Hence, by slightly changing θs1 , x reached fixed points.Accordingly, ∂ xs1 /∂ θs1 increase beyond one, so that ∂Θ s1 /∂ θs1 became positive at θu , approaching θth u ∼ 0.4, as shown in Figs.2d and 3b.
Thus, we have unveiled how attraction to the unstable manifold is achieved by slow epigenetic fixation of the oscillation of fast gene expression in the repressilator model.Following this picture, reprogramming is possible by forcing the cells to return to the oscillatory state.Then, the cell is attracted to a pluripotent state with low epigenetic modification θ i ∼ 0, from which differentiation to distinct cell types with specific θ values follows.
To verify the generality of this reprogramming scheme, we examined several GRN models with more degrees of freedom.As discussed in [52], differentiation from oscillatory states is often observed in GRNs (e.g., 20% of randomly generated GRNs show oscillatory dynamics for N = 10).An example is shown in Fig. S8a [49].From a differentiated state, we overexpressed three genes to regain oscillatory expression (black line in Fig. S8a [49]).Later, global attraction to unstable manifold also occurred as discussed above.Then, the cell states branch to distinct fixed point states again (blue line in Fig.
S8a [49]).In these cases, the original pluripotent state with θ = 0 was an unstable fixed point, with one positive eigenvalue for the Jacobi matrix of θ dynamics (Fig. S8d [49]), as in the repressilator model.Even though the degrees of freedom increased, the unstable manifold is one-dimensional, and the attraction to the manifold occurred from a higher-dimensional state space.This implies that reprogramming manipulation requires only partial degrees of freedom compared with the total number of genes.In fact, overexpression of three genes is sufficient for reprogramming in GRN models with N = 10, as far as we have investigated.
The present mechanism also works in a model extracted from GRN for an ES cell [53], as a core network with five genes (Nanog, Oct4, Gata6, Gata4, and Klf4 ) [29] (see Fig. S4b).Oct4, Sox2, and Klf4 are known as factors to induce reprogramming.The model involves a negative feedback loop, as in the repressilator, in addition to positive feedback regulation.In this five-gene model, x i and θ i oscillate in the region near the origin, and then differentiation to three fixed points progresses as in the repressilator case (three lines in Fig. 4a), whereas θ i = 0 for all i represents a saddle point with one unstable manifold and four stable manifolds, as shown in Fig. S9b.After overexpression of Oct4, Nanog, and Klf4 in one of the differentiated cell types for a certain time span (black dotted line in Fig. 4b) [54], the epigenetic state θ i approached the unstable manifold for the unstable fixed point θ i = 0, leading to recovery of pluripotency (blue line in Fig. 4b).
In this letter, we have shown that oscillatory gene expression dynamics with slow epigenetic modifications lead to cellular reprogramming by overexpression of only few genes.The global attraction to the unstable manifold of the saddle point explains the reprogramming process.Now, the return to the top of the landscape by reprogramming, which is seemingly unstable, is explained by the strong attraction toward the unstable manifold of the saddle, and suppressed instability along with the unstable manifold, owing to the approach of limit-cycle of bifurcation to fixed points.The memory of the cellular state before reprogramming manipulation was erased through this reprogramming process.
Moreover, regain of oscillation was found to be the main requirement for reprogramming, whereas elaborate manipulations to induce a cellular state into specific states is not necessary.This explains the role of oscillations in gene expression in pluripotent cells [21] and epigenetic modification through the differentiation process [18], as well as it explains how reprogramming is possible by overexpressing just few genes among thousand of that [11,14].The timescale separation between fast expression dynamics and slow epigenetic modification feedback required is also consistent with previous observations [46,47].In future studies, experimental support is necessary, as well as theoretical analysis of slow-fast dynamical systems [55,56].
and fixed-point analysis around {θ 1 , θ 2 , θ 3 } = {0, 0, 0}.The repressilator model obviously has symmetry against a transformation in 1 Eq. (S1) can be written as the following circulant matrix: Here, we considered the following matrix C: ∂Θ/∂θ can be expanded as: Matrix C can be diagonalized by a following 3 × 3 discrete Fourier transform (DFT) matrix F : where ω j = exp(2πij/3) with j = 0, 1, 2 (i, used here, is an imaginary unit).Thus, the circulant matrix ∂Θ/∂θ can be diagonalized using the DFT matrix F .The jth eigenvector v j is now written as: Each eigenvalue λ j corresponding to the eigenvector v j is Accordingly, we obtained Now, we estimated the value of a 0 , a +1 , a −1 .a 0 = ∂ xi /∂θ i is given by the change in gene expression against the change in epigenetic modification of the same gene (see Fig. S4a).In our model, we adopted a positive feedback relationship between gene expression and epigenetic modification of the gene i; a 0 is positive.In contrast, a +1 is given by the change in gene expression against the change in epigenetic modification of the repressing gene; a +1 is negative.Similarly, a −1 is positive.In addition, the change in gene expression was amplified throughout the negative repressilator loop; thus, a 0 , a +1 , a −1 should satisfy Indeed, from numerical calculations, the values of a 0 , a +1 , a −1 followed these estimations (a 0 = 0.66, a +1 = −1.08, a −1 = 1.72 for g = −0.4).From these estimations, we obtained the relation of eigenvalues as: In the present case (a 0 = 0.66, a +1 = −1.08, a −1 = 1.72 for g = −0.4),λ 0 is positive, and λ 1 , λ 2 is negative: Thus, θ = 0-saddle point consists of an eigenvector v 0 with λ 0 > 0 and eigenvectors v 0 , v 1 with Re(λ 1 ) = Re(λ 2 ) < 0. To analyze and plot in real θ space, we adopted with a normalization factor.The symmetry of the repressilator against the transformation in 1 → 2 → 3 → 1 is not broken as long as θ 1 = θ 2 = θ 3 is satisfied.Therefore, the direction of the eigenvector v u always coincides with that of the unstable manifold.
B. Linear stability analysis for x dynamics with fixed θ We discussed x dynamics of repressilator model with fixed θ as parameter.Full model with fixed θ can be written as Thus, Jacobi Matrix ∂X i /∂x j can also be written as First, we studied the stability of the fixed point for {θ 1 , θ 2 , θ 3 } = {0, 0, 0} to elucidate the condition for oscillation in x dynamics with no epigenetic modification.In this case, {x 1 , x 2 , x 3 } = {0, 0, 0} satisfies the fixed-point condition.
By inserting x i = 0, θ i = 0 for all i in Eq. (S17), we obtained: From Eq. (S18), eigenvalue λ is given by In the repressilator model(g > 0), the first value of (S19) is always negative.The fixed point is unstable for where the limit-cycle attractor appeared as a result of Hopf bifurcation, to specific g > 0.05 for β = 40.
Next, we considered the bifurcation from the limit cycle to a fixed point by changing the parameter θ.From Eq. (S17), the eigenvalues of the Jacobi matrix {λ} are given by the solution of the following equation: Here, we plotted −βg/ cosh 2 (−βg(x − z)) as a function of x in Fig. S1.The 1/ cosh 2 function has a tall spike at x = z for a sufficiently large β.Considering that we adopted the tanh function, sigmoidal function with range{−1, 1}, the value of x i is restricted to {−1, 1}.This suggests that if even one of θ i /g(i = 1, 2, 3) in Eq. (S21) exceeds one, and the second term in Eq. (S21) is close to zero, and the solution to Eq. (S21) is only determined by (−1 − λ) 3 , then λ = −1, corresponding to a stable fixed point.Considering that Eq. (S19) gives the solution corresponding to an unstable fixed point, the bifurcation point from the limit cycle to the fixed point was determined as follows: (a)   The eigenvalues λ k of the ∂Θi/∂θj matrix given by the fixed-point analysis of Eq. ( 2) around θi = 0 for all i.In Fig. 4, we adopted eigenvectors vs 1 , vs 3 with different Re(λ k ).

FIG. 2 .
FIG. 2. (a) Stream plot of (d θu/dt, d θs 1 /dt) in ( θu, θs 1 ) space according to Eq. (2).Red (blue) line represents the direction of the eigenvector vu(vs 1 ).(b) Attractors in the x-space for each fixed θ value.For the green (blue)-colored region, the attractor in the x space was the limit cycle (fixed point).Each trajectory shows the dynamics in the x1 − x2 space for ( θu, θs 1 ) (both x1 and x2 axes for all figures are ranged within {−1.05, 1.05}).(c) d θu/dt plotted as a function of θu.For comparison, we plotted d θu = λu θu(red dotted line).(d) Degree of attraction to vu. ν, in the text, is plotted as a function of θu.If it is negative (positive), {θ} was attracted to (departed from) vu.See Fig. S6[49] for more detail.

FIG. 4 .
FIG. 4. (a) Cell differentiation and (b) reprogramming in the five-gene model (a) Three orbits starting from the vicinity of the saddle point θi = 0 for all i (black dotted point), reached three distinct cell types.(b) From differentiated cell types (red point), we added external input Ii(t) to Nanog, Oct4, and Klf4 for a certain time span (black dotted line).After such reprogramming manipulation, we set Ii(t) = 0.The cell state then spontaneously approached the saddle point, and the reinitiated the differentiation progression again (blue line).τ = 10 3 .See Fig. S9[49] for more detail.

1 FIG
FIG. S1. −βg/ cosh 2 (βg(x − z)) as a function of x (β = 40, g = 0.4).This function has a tall spike at x = z.Considering the range of the tanh function, this term could be regarded as zero if |z| > 1 for a sufficiently large β.
FIG. S8.Random gene regulatory network model with N = 10.We randomly generated the gene regulatory matrix Jij and selected one sample that showed oscillation dynamics for xi(t) around θi = 0. Jij was generated according to the normal distribution, with an average of zero and a standard deviation of 1/ √ N .(a) Cell differentiation and reprogramming.For plotting, we adopted a principal component analysis (PCA) of {θi}(i = 1, 2, . . ., 10) over 1,000 trajectories starting from θi ∼ 0. The trajectory {θi(t)} was plotted by adopting the first, second, and third PCA modes.(Red lines) Trajectories starting from θi ∼ 0, five different fixed points are reached after transient oscillation dynamics.(Black dotted line) From one of the differentiated cell types, we added external input Ii(t) to three genes as reprogramming manipulation.(Blue line) After reprogramming manipulation, the cell state was attracted to the unstable manifold and differentiated again.(b) Time development of cell differentiation and reprogramming (Fig.S8a.The upper bar shows the differentiation and reprogramming process, similar to Fig.1ab.(c) Time development of cell differentiation (red region in Fig.S8b) for a short time window.(d) The eigenvalues λ k of ∂Θi/∂θj matrix given by fixed point analysis of Eq. (2) around θi = 0 for all i.
FIG. S9. (a)Time development for each trajectory of cell differentiation (yellow, green, and red line in Fig.4a) and reprogramming (black dotted and blue lines in Fig.4b).(Yellow, green, and red) Three trajectories starting from θi ∼ 0 are plotted, which reach three different attractors.(Black) From the differentiated state, we added an external input Ii(t) to Nanog(x1), Oct4 (x2), and Klf4 (x5) as reprogramming manipulation.(Blue) After reprogramming manipulation, the cellular state was first attracted to the unstable manifold, and then differentiated again.(b) The eigenvalues λ k of the ∂Θi/∂θj matrix given by the fixed-point analysis of Eq. (2) around θi = 0 for all i.In Fig.4, we adopted eigenvectors vs 1 , vs 3 with different Re(λ k ).