Network Features and Dynamical Landscape of Naive and Primed Pluripotency

Abstract

Although the broad and unique differentiation potential of pluripotent stem cells relies on a complex transcriptional network centered around Oct4, Sox2, and Nanog, two well-distinct pluripotent states, called "naive" and "primed", have been described in vitro and markedly differ in their developmental potential, their expression profiles, their signaling requirements, and their reciprocal conversion. Aiming to determine the key features that segregate and coordinate these two states, data-driven optimization of network models is performed to identify relevant parameter regimes and reduce network complexity to its core structure. Decision dynamics of optimized networks is characterized by signal-dependent multistability and strongly asymmetric transitions among naive, primed, and nonpluripotent states. Further model perturbation and reduction approaches reveal that such a dynamical landscape of pluripotency involves a functional partitioning of the regulatory network. Specifically, two overlapping positive feedback modules, Klf4/Esrrb/Nanog and Oct4/Nanog, stabilize the naive or the primed state, respectively. In turn, their incoherent feedforward and negative feedback coupling mediated by the Erk/Gsk3 module is critical for robust segregation and sequential progression between naive and primed states before irreversible exit from pluripotency.

Figure 1

The pluripotency model: network architecture and cell fate data. (A) Pluripotency displays specific in vitro and in vivo transition patterns among a naive/mESC pluripotent state, N; a primed/mEpiSC pluripotent state, P; and a nonpluripotent differentiated state, D. (B) The regulatory network model is based on a selected set of signaling and transcription factors interacting with each other through various regulatory mechanisms (arrow legend box). (C) A set of perturbation experiments is used as target experimental dataset for the model optimization procedure. These experiments describe the fate outcome j′ = N,P,D of cells that are initially in a state j = N,P,D and then subjected to specific culture conditions (LIF, Activin, FGF) or to upregulation/downregulation (±) of specific factors. (D) Typical relative mRNA levels of pluripotency factors (Oct4-Sox2, Nanog, Esrrb, and Klf4) are associated with naive/mESC, primed/mEpiSC, and nonpluripotent states, which are determined from experiments (Fig. S4) and used for model optimization. To see this figure in color, go online.

Figure 2

Model optimization through computational evolution. (A) The statistical distribution of the global error score Φ is computed over 10³ randomized-parameter models. (B) An example of a computational evolutionary trial showing the decreasing score Φ from 0.8 to 0.2 and the parameter set of the fittest model individual of the family of the ith generation. The generation 0 is composed of N_I = 20 model individuals k of parameter $P_i^k={10}^{{\eta }_k+{\eta }_i}$, where η_i,k values are uniformly distributed between [−1:1]. (C) The statistical distribution of the global error score Φ is computed for the fittest individuals of 10³ evolutionary trials. (D) Given here are condition-specific scores ϕ_i^k of a set of optimal models satisfying Φ < 0.21: model-averaged value (lower panel) and fraction of models with ϕ_i^k > 0.5 (upper panel). (E) Three main classes of optimal models (green, class I; blue, class II; red, class III) are evidenced by the distance matrix of their condition-specific scores (Eq. 5). (F) The three model classes differ in their worst-condition-specific scores. To see this figure in color, go online.

Figure 3

Parameter distribution and core topologies of optimized models. The systematic analysis of the $i=\mathrm{1,49}$ optimized parameters is made for the k = 1,29 optimal models segregated into the three classes identified in Fig. 2E (green, class I; blue, class II; red, class III). (A) The two left panels show parameter values P_i^k (k points in the ith line) and their class-averaged values 〈P_i〉_k⊆class. The two right panels show parameter sensitivities γ_i^k (Eq. 4) and their class-averaged values 〈γ_i^k〉_k⊆class, where the dashed line shows the threshold γ₀ = 0.1 between relevant and irrelevant regulations (Fig. S6). (B) The three classes of optimal model differ in their topology, as evidenced by the distance matrix of model parameter sensitivities D_kl(log γ_i) (Eq. 5). (C) Core pluripotent circuitries for each model class are drawn by including only regulations whose removal results in a significant increase of error score Φ (〈γ_i〉_k⊆class > γ₀ in (A)). Model classes share 15 common regulations (thick network edges). To see this figure in color, go online.

Figure 4

Multistability domains and asymmetric transitions between naive and primed states. Color code represents the three classes of optimal models (green, class I; blue, II; red, class III). (A) Normalized expression/activity levels X^j_∞,j of the i protein species are associated with the stable states j = N,P,D (top, naive state with S_LIF = 1; middle, primed state with S_Act/_FGF = 1; bottom, differentiated state with S_LIF = S_Act/FGF = 0). (B) Distinct multistable subdomains between these three states are depicted by their average expression levels of pluripotency factors 〈X_i^k〉_{i = Oct, Nan, Esr, Klf} as a function of S_LIF (and S_Act/_FGF = 0) and S_Act/FGF (and S_LIF = 0). (Lower panel) An example of multistable property with two bistability domains (shading) and two saddle-node bifurcation points (circles) is illustrated for one optimal model. (C) Dynamical trajectories of optimal models in a 2D section of the concentration space are shown for distinct transition scenarios between naive and primed pluripotent states. (Middle panel, insets) A set of trajectories from a heterogeneous naive state $\left({\overrightarrow{X}}_0={\overrightarrow{X}}_{\infty }^N+\overrightarrow{\eta }\right)$ seems attracted to a slow manifold through the primed state. (D) Primed-to-naive reversion efficiencies ϵ_i^k = 1/θ_i^k are given for the strategy i (overexpression of a given pluripotency factor or 2i, with or without LIF) and optimal model k, where θ_i is the critical level of upregulation or downregulation α_i or of S_LIF that is required for successful reversion. To see this figure in color, go online.

Figure 5

Identification of a modular network through systematic perturbation approach. (A) The approach consists in quantifying the effect of diverse network motif perturbations on the signaling thresholds and expression levels of the naive and primed states. (B) Model-average alterations 〈|ΔY^k|〉_k of expression levels ($Y={\langle X_{\infty ,i}\rangle }_{i=\text{Oct},\text{Nan},\text{Klf},\text{Esr}}\left(S_j=1\right)$ and signaling thresholds (Y = S^C_j) of the naive state (left, j = LIF) and the primed state (right, j = Act/FGF) are computed upon the deletion of a single or multiple (^∗) regulatory links (a ↔ b: a → b ∪ b → a; a → b,c : a → b ∪ a → c). (C) Functional modular circuits are inferred from (B) regulating the naive and primed states. (Black arrows) Activatory regulation; (red circles) inhibitory regulation. (Solid lines) Effective regulation; (dashed line) ineffective regulation. To see this figure in color, go online.

Figure 6

A low-dimensional modular model reproduces the pluripotent dynamics of detailed models. (A) Minimal circuit corresponding to Eq. 6 consists of regulatory modules x_{i = N,I,P}, extrinsic signals s_{i = L,A}, and regulatory parameters a_ii′. (Black arrows) Activatory regulation (red circles) inhibitory regulation. (B) A signal-dependent multistability among three discrete states similar to Fig. 4B is found for 381 parameter sets from a sample of 10⁶ random parameter sets for which log₁₀(a_ii′) values are uniformly distributed between [−1:1] except a_NL = 1 and a_PA = 1. (See Fig. S7A for the quantitative multistability criteria used for parameter set selection.) (C) Three examples of canalized trajectories where heterogeneous initial conditions around the naive-like state (dark shaded circle) are quickly attracted to a slow 1D manifold that passes near the primed-like state (light shaded circle) toward the differentiated-like state (blank circle). (D and E) The 381-parameter set satisfying the multistability criteria in (B) shows specific distribution profiles for each parameter with peaks for high values (D), and significant positive correlations between the parameters a_NN and a_IN (lower panel of E) and between the parameters a_NN and Δλ (upper panel of E). Δλ is a measure of stability of the naive state (Fig. S7, B and C). To see this figure in color, go online.