Transcriptional Dynamics of the Embryonic Stem Cell Switch

Recent ChIP experiments of human and mouse embryonic stem cells have elucidated the architecture of the transcriptional regulatory circuitry responsible for cell determination, which involves the transcription factors OCT4, SOX2, and NANOG. In addition to regulating each other through feedback loops, these genes also regulate downstream target genes involved in the maintenance and differentiation of embryonic stem cells. A search for the OCT4–SOX2–NANOG network motif in other species reveals that it is unique to mammals. With a kinetic modeling approach, we ascribe function to the observed OCT4–SOX2–NANOG network by making plausible assumptions about the interactions between the transcription factors at the gene promoter binding sites and RNA polymerase (RNAP), at each of the three genes as well as at the target genes. We identify a bistable switch in the network, which arises due to several positive feedback loops, and is switched on/off by input environmental signals. The switch stabilizes the expression levels of the three genes, and through their regulatory roles on the downstream target genes, leads to a binary decision: when OCT4, SOX2, and NANOG are expressed and the switch is on, the self-renewal genes are on and the differentiation genes are off. The opposite holds when the switch is off. The model is extremely robust to parameter changes. In addition to providing a self-consistent picture of the transcriptional circuit, the model generates several predictions. Increasing the binding strength of NANOG to OCT4 and SOX2, or increasing its basal transcriptional rate, leads to an irreversible bistable switch: the switch remains on even when the activating signal is removed. Hence, the stem cell can be manipulated to be self-renewing without the requirement of input signals. We also suggest tests that could discriminate between a variety of feedforward regulation architectures of the target genes by OCT4, SOX2, and NANOG.

model but can also be used to compare directly to experimental results, since the model uses measurable equilibrium reaction rates. We illustrate this by first evaluating the transcription rate using the statistical mechanics approach and then comparing to the equilibrium approach for a simple system. Consider an example with a single gene that is autoregulatory, i.e. the transcription rate is enhanced when it is bound by its own protein product, X. There are then four possible states of the system. With the logic shown in Table S1, the transcription rate is given by [X] (S1) In Eq. (S1), δ i are related to the binding free energies, and [R] is the RNAP concentration.
Denote by G 0 , G X , G R and G XR the four states of the gene; unbound, bound by the transcription factor X, bound by the RNAP and bound by both X and RNAP. These states are subject to a normalization The following reaction scheme defines the network At thermodynamic equilibrium, one has where K i = k if /k ib are the equilibrium constants. From Eqs. (S3,S4) and Eq. (S2), we compute the fractional probability P of the gene being bound by RNAP, which is functionally similar in form to Eq. (S1). Both methods are equivalent as the only physical requirement is thermal equilibrium for these reactions. Although the statistical mechanics approach is more intuitive, the equilibrium approach allows us to define a reaction scheme, which can be described in terms of measurable kinetic constants. It is also clear that such a scheme is required to carry out a stochastic simulation.

Parameter dependence of the stem cell box switch behavior
Next we investigate how the model parameters determine the dynamics of the stem cell box.
More precisely, we probe the parameter range for which one obtains a switch-like behavior.
In other words, is the switch-like behavior robust? The curves of the steady state values of the concentrations of OCT4, SOX2, NANOG and the OCT4-SOX2 complex in Figs. 3 and 4 in the main text as functions of the environmental factor A + and signal B − exhibit bistable switch-like behavior with two saddle-node bifurcations (SN). These curves were obtained using the parameter sets given in the Table S2. We study the parameter dependence by asking as to how the switch-like behavior changes as the parameters around the nominal parameter set displayed in Table S2 are varied. Thus, we follow the SN of the bifurcations in Fig. 3 as functions of the various parameters. In effect, these are two-parameter bifurcation plots. We therefore plot the locus of the thresholds of the environmental factor A 0 + , at which the switch is turned on, as a function of the parameters. In each of the plots, that we will subsequently display (Figs. S1-S4) the y-axis represents the particular parameter value being varied, and the x-axis the thresholds A 0 + at which the SN occurs. We obtain either of two cases: For each parameter value, there are two points corresponding to a bistable switch with two turning points or a single threshold corresponding to an irreversible bistable switch. If the threshold value is zero, which happens for extremely low parameter values, it corresponds to the disappearance of the switch-like behavior altogether.
Binding strengths of the OCT4-SOX2 and OCT4-SOX2-NANOG complexes. In studying the parameter dependence of the binding strengths of the protein complexes on the three genes, we use the fact that the nature of the regulation on OCT4 and SOX2 is similar.
Hence, it is sufficient to study the parameter dependence of the binding at the NANOG gene and either of the OCT4 and SOX2 genes. We therefore consider binding strengths a 3 and a 2 which determine how strongly OCT4-SOX2-NANOG and OCT4-SOX2 bind to OCT4, respectively, and e 2 and e 1 , the binding strengths of the OCT4-SOX2-NANOG and OCT4-SOX2 complexes on the NANOG gene respectively (see Eq. (1) in the main text.
From Fig. S1, one observes the threshold at which the switch-like behavior occurs reduces as the free energy of binding decreases, or as a 3 , a 2 , e 2 and e 1 increase. This is intuitively obvious, since the stronger that these complexes bind to OCT4, SOX2 and NANOG, the higher the transcription rate, and hence a smaller threshold of the environmental factor is required to turn the switch on. Although qualitatively similar results are obtained for the binding strengths of the OCT4-SOX2-NANOG and OCT4-SOX2 complexes on OCT4 and NANOG, the threshold A 0 + of the environmental factor, which is required to turn the switch on, is much smaller for the NANOG gene, as compared to the binding at OCT4 and SOX2.
This suggests that the OCT4-SOX2-NANOG complex binding to NANOG, can be much weaker as compared to its binding at OCT4 and SOX2, and yet have significant effect on the dynamics. This is because the regulation at NANOG is autoregulatory which turns out to be crucial in kick-starting the stem cell box switch. As can be seen From fig. S1, the values for a 3 , a 2 , e 2 and e 1 range over almost three orders of magnitude, i.e. from .01 to 10, for which the switch-like behavior is exhibited 1 . Hence, the switch is extremely robust.
There is an additional interesting feature in the behavior of a 3 , e 2 , the binding strengths of the OCT4-SOX2-NANOG complex on the OCT4 (SOX2) and NANOG genes respectively. Formation and degradation of the OCT4-SOX2 complex. In Fig. S2 we show the formation, dissociation and degradation rates, k 1c , k 2c and k 3c respectively, of the OCT4-SOX2 complex as functions the threshold A 0 + , at which the switch is turned on. As k 2c and k 3c increase, A 0 + also increases since the SOX2 and OCT4 transcription rates must be increased to account for loss of the OCT4-SOX2 complex. The opposite effect occurs for k 1c , since as it increases, the rate of complex formation is increased and hence a lower threshold is required to turn the switch on. For k 1c > 0.5, the system becomes a irreversible switch.
Degradation rates of OCT4, SOX2 and NANOG proteins. As can be seen from Fig.   S3, increasing the degradation rates of OCT4, SOX2 and NANOG effectively increase the threshold of the switch, since a higher activation from the environmental factor is required to account for loss of proteins through degradation. Degradation of the NANOG protein is particularly sensitive in that the threshold dramatically increases on increase of its degradation rate. The switch-like behavior is lost when the degradation rate of NANOG reaches g 3 > 1.5 and degradation rate for OCT4/SOX2 reaches g 1 > 5. The degradation rate of NANOG is much more sensitive to any of the other parameters due to its strong impact upon the dynamics. Being autoregulatory, the degradation rate of NANOG is crucial in determining its steady state value.
Basal transcription rates of the OCT4, SOX2 and NANOG genes. The basal transcription rates have the effect of reducing the threshold of activation of the switch, since transcription occurs even for a smaller environmental factor A + . Hence, as can be seen from Parameter values for Figures 10, 11 in the main text. The same parameters as in Table S2 are used to generate Figs. 10 and 11 in the main text except for the changes shown in Table S3. To obtain Fig. 11 we first assume that the stem cell box is on (steady state values of all the species for A + = 100). A is then set to zero (removal of the environmental factor) and the steady states of the system are evaluated as function of signal B − .

Equations and Parameters for Target Gene Regulation
To study the dynamics of the target genes, due to regulation by OCT4-SOX2 and NANOG, we initially ignore the feedback from NANOG to OCT4 and SOX2. Hence, we assume the input to the system, OCT4-SOX2, to be constant. The equations determining the regulation of the target gene product are then given by where [T G] refers to the concentration of the target gene product. The parameters used for Fig. 6 are given in Table S4 (with f 3 = 0). The curves in Fig. 6, which refer to the simple feedforward network where NANOG does not have autoregulation are obtained by simulating the above equations for e 2 = 0, f 2 = 0, and with the same values as in Table S4, for the remaining parameters. In Fig. 6C the expression for the self-renewal gene is plotted for the parameters in Table S4, with the change, h 2 = 0.001.
The equations, which are used to describe the regulation of the self-renewal target genes (TG) with OCT4-SOX2 and NANOG as activators and differentiation TG's with OCT4-SOX2 and NANOG as repressors, are given below. These follow from the truth table given in Table 3B and Table 3C respectively.
The above equations are used in conjunction with the equations for the stem cell box, to generate the curves in Fig. 8 in the main text.

Parameters for Integrated Model
In the integrated model we combine Eqs. (1,2), from which we generate Fig. 7 in the main text. NANOG now feeds back to OCT4 and SOX2. For the combined model the parameters for the OCT4, SOX2 and NANOG genes are the same as in Table S2. The target gene parameters are given in Table S5. For Fig. 7B in the main text, which shows the output of self-renewal genes, the parameter h 2 = 0.001 is used, and the remaining parameters are the same as in Table S5.
For Fig. 8, in which we study the dependence of the self-renewal (Fig. 8A), and differentiation ( Fig. 8B) TG products as a function of the input environmental factor A + , we use the parameters given in Table S6 (these correspond to Eq. S7.). The target gene parameters are given in Table S5.
The parameters used to generate Fig. 9 are the following: The target gene parameters for all sets of curves are given in Table S5. For Fig. 9A, the stem cell box parameters are the same as in Table S1, except for the changes given in Table S7. For Fig. 9B, which displays the TG expression for strong binding of NANOG to the OCT4 and SOX2 genes, we use the parameters given in Table S8. The only parameter which is changed is h 2 = 0.01.
For Fig. 12, which shows the irreversible bistability for NANOG over-expression, we use the same parameters as in Table S2, except for the following changes, η 5 = 35, η 6 = 0.035. The target gene parameters are given in Table S5.
In all of the figures generated which are referred to in this section the parameter f 3 = 0.
Dependence of Target Gene expression assuming an incoherent feedforward architecture.
Next we comment on the dependence of the target gene expression assuming an incoherent feedforward architecture, on network parameters. The motivation is to study the response of the amplitude filter curve as a function of the binding strengths of the transcription factors at the target gene. We use the parameters given in S4. From Eq. (2) in the main text, for the rate of change of the target gene product, the two key parameters are g 1 , h 2 , the strength of binding of the OCT4-SOX2 complex, and the OCT4-SOX2-NANOG complex to the target genes, respectively. We noted earlier that as the free energy of binding of the OCT4-SOX2-NANOG to the target gene increases, or as h2 decreases, the target gene gets increasingly expressed. This is seen in the third panel in Fig. S5, which shows that as h 2 decreases, the target gene repression increases and we loose the sharp cutoff of the steady state value of the protein level with respect to OCT4-SOX2 protein levels. In the middle panel, g 1 , is varied, which determines the rate of transcription of the target gene.
Increasing the amount of the input OCT4-SOX2 complex has the effect of simply increasing the total amount of target gene protein level. Finally, in the uppermost panel, we vary e 2 , the binding affinity of the OCT4-SOX2-NANOG complex to the NANOG gene, which essentially tunes the amount of autoregulation. Decreasing the amount of autoregulation leads to less expression of NANOG, and hence less repression of the target genes.

Tables
Table S1: Logics and rates for a single autoregulatory gene.              Figure S1: Two parameter bifurcation plots for binding strengths of protein complexes to the genes. The y-axis represents the value of the parameter, in this case a 3 , a 2 , e 2 and e 1 , and the x-axis the threshold A 0 + of the environmental factor, at which the switch turns on.  Figure S2: Two parameter bifurcation plots for formation and degradation rates, k 1c , k 2c , k 3c of the OCT4-SOX2 complex. Increasing k 1c beyond 0.5, leads to a irreversible switch.    Figure S5: Target gene expression curves when varying parameters. Top panel: e 2 represents the strength of binding of OCT4-SOX2 and NANOG onto the NANOG gene, and therefore it represents the amount of autoregulation due to positive feedback. As we decrease the feedback, NANOG expression reduces, and hence the expression of the differentiation target genes increases. Middle panel: g 1 determines how strongly OCT4-SOX2 binds to the target gene. Hence for lower values of g 1 , the target gene expression decreases. Lower Panel: h 2 represents the repression due to the OCT4-SOX2, NANOG complex binding to the target gene gene. Reducing the affinity of binding of the complex leads to a release of the expression of the differentiation target genes.