The potential landscape of genetic circuits imposes the arrow of time in stem cell differentiation

Abstract

Differentiation from a multipotent stem or progenitor state to a mature cell is an essentially irreversible process. The associated changes in gene expression patterns exhibit time-directionality. This "arrow of time" in the collective change of gene expression across multiple stable gene expression patterns (attractors) is not explained by the regulated activation, the suppression of individual genes which are bidirectional molecular processes, or by the standard dynamical models of the underlying gene circuit which only account for local stability of attractors. To capture the global dynamics of this nonequilibrium system and gain insight in the time-asymmetry of state transitions, we computed the quasipotential landscape of the stochastic dynamics of a canonical gene circuit that governs branching cell fate commitment. The potential landscape reveals the global dynamics and permits the calculation of potential barriers between cell phenotypes imposed by the circuit architecture. The generic asymmetry of barrier heights indicates that the transition from the uncommitted multipotent state to differentiated states is inherently unidirectional. The model agrees with observations and predicts the extreme conditions for reprogramming cells back to the undifferentiated state.

Figure 1

Dynamics of the canonical gene regulatory circuit of two mutually opposing transcription factors that positively self-regulate themselves. (A) Circuit architecture for the two genes X₁ and X₂. (B) Bifurcation diagrams indicating the stable position of S(x₁, x₂) where x₁ = x₂ for the symmetric case (vertical axis), during the symmetric change of a = a₁ = a₂ over the indicated range of values (horizontal axis), for the other parameter values b₁ = b₂ = 1, k₁ = k₂ = 1, and S = 0.5, n = 4. (C) Force field in the X₁ – X₂ state space for two parameter values for parameter a on both sides of the respective critical point in the bifurcation diagram. (D and F) Steady-state probability distribution P_ss(S) calculated from the Fokker-Planck equation (Eq. 1) as function of the parameter a in panel D or the noise parameter D in panel F. Colors indicate the probability P as shown in the color bar. (E and G) The corresponding quasipotential landscape where the elevation of the landscape (quasipotential) represents –ln(P(S)).

Figure 2

Dynamical behavior of the probability P(S, t) and flux vectors during fate decision of a stem cell in the multipotent state S^∗_C. P(S, t) is evaluated during bifurcation from tristable to the bistable regime as the parameters a₁ and a₂ are decreased according to a₁ ∼ exp(–λ₁^∗t) and a₂ ∼ exp(–λ₂^∗t) with λ₁ = 0.01 and λ₂ = 0.015. In panel A, the initial state is near the central attractor S^∗_C, P(S = (0.3, 0.3), t = 0) = 1, whereas in panel B, the initial state is near the attractor S^∗_A, i.e., P(S(1.0, 0.0), t = 0) = 1. Other parameters are the same as Fig. 1.

Figure 3

Relative barrier height as a function of parameter a accounts for directionality of state transitions around the bifurcation. (A) The bifurcation diagram for same parameters as in Fig. 1B with large arrows representing the transitions across the respective barriers U_SC and U_SA that separate the central stem cell attractors S^∗_C and the differentiated cell attractors S^∗_A and S^∗_B. (B) Computed heights of the barriers U_SC and U_SA as a function of a = a₁ = a₂ (for noise level D = 0.05). Here a_crit = bifurcation point, a′= value of a at which the relative barrier heights reverse. (C) Sections through the potential landscape illustrating the barriers at the three indicated values of a (dashed arrows). (D) The transitions mean first passage times τ for $S_C^{\ast }\to S_A^{\ast }$ (τ_CA) and S^∗_A → S^∗_C (τ_AC). (E) Direct comparison of U and τ as function of a: log-scale relative barrier heights or transition times for the transitions S^∗_A → S^∗_C and S^∗_C → S^∗_A, respectively. Other parameters are the same as Fig. 1.

Figure 4

Dependence from the noise level D of barrier heights U and state transition times τ (A and B) and their equivalence (C and D), regime a = 1. The values τ_AC, τ_AB, and τ_CA denote the transition times (given a noise level D) for noise-driven transitions between the respective attractors, and U_SA and U_SC represent the barrier heights separating the respective attractors. Here a₁ = a₂ = 1 (other parameters the same as Fig. 1).