Exemplary Final Term Paper 4

Introduction

Background

When embryonic stem cells divide, they either remain stem cells in a state of self-renewal or differentiate to become any other type of cell . Certain transcription factors in these cells control this process through self-regulation . Since other types of cells can be obtained from stem cells, understanding what causes the cells to differentiate will allow scientists to control the differentiation process . Controlling this process will allow stem cells to be used to create more stem cells or to create lineages of specific types of cells .

Previous experiments have shown that numerous transcription factors affect the rate of transcription of the self-renewal and differentiation genes . Transcription factors are proteins that either increase or inhibit the transcription of genes through a variety of mechanisms . Three of these transcription factors, referred to as OCT4, SOX2, and NANOG, were demonstrated to be able to control the renewal or differentiation of stem cells . Each of these transcription factors regulates the other considered transcription factors in a feedforward network, where transcription factors promote each other as well as encourage the transcription of the target genes . The transcription factors are also found within the target genes, further demonstrating the interdependence of this network . Therefore, all of these transcription factors play a role in preventing or promoting the differentiation of the stem cells . Additionally, OCT4 and SOX2 form the OCT4-SOX2 heterodimer, which can form a complex with NANOG also, both of which also act as separate transcription factors affecting the production and repression of the other transcription factors .

The network of transcription factors and their subsequent effect on the transcription of the target genes are controlled by several signals that either promote or prevent the transcription of the target genes related to differentiation or self-renewal . In "Transcriptional Dynamics of the Embryonic Stem Cell Switch" by Vijay Chickarmane et al., the effects that the signals have on the concentrations of the transcription factors and the transcription of the target genes are modeled . Since there are many different signals that affect transcription, this paper will simplify the signals to four signals, referred to as A⁺, A^-, B⁺, and B^- . A⁺ and A^- regulate OCT4 and SOX2 and B⁺ and B^- regulate NANOG . This model is based solely on A⁺ and B^-, since different combinations can cover all possibilities . There are two possible ways to create this model: a coherent model or an incoherent model . In the coherent model, both the complex and the other transcription factor activate the self-renewal genes and repress the differentiation genes . For the incoherent model, one transcription factor complex activates, or promotes the transcription of, both the self-renewal and differentiation genes while another transcription factor represses, or prevents the transcription of, the self-renewal genes weakly and the differentiation genes strongly . Both models reasonably describe the effect the transcription factors have on the transcription of the target genes.

Hypothesis

By creating a model, we will test whether the concentration of outside signals can switch the stem cell from a state of self-renewal to differentiation or the reverse and that this behavior exhibits bistability. Further, by varying parameter values in the model, the ability to cause the stem cell to undergo irreversible switches using high concentration of signals for a brief time will be made possible by either changing the binding and repressive strength between OCT4-SOX2 and the transcription factors or instead increasing the transcription rate of NANOG . Through this model, the effects of signal concentration and parameter values will become more clear, allowing the process of self-renewal and differentiation to be better understood . The model will also make it possible to understand what must be examined to determine if the system follows the coherent or the incoherent model .

Results and Extension Overview

All figures from the paper by Chickarmane et al. were successfully reproduced in Mathematica with some quantitative differences . Qualitatively, almost all figures are equivalent, with only one difference in one of the supplementary figures showing a two-parameter bifurcation plot. Based on these results, the hypothesis that the system undergoes a bistable switch for changing concentrations of the two signals is correct, as well as the hypothesis that by changing the binding and repressive strength of OCT4-SOX2 or increasing the basal transcription rate of NANOG causes the switch to become irreversible.

The extension of the original model included the creation of the corresponding time plots which were not included in the paper by Chickarmane et al., as well as including a component of stochastic noise in the concentration of the three transcription factors . The time plots are important to demonstrate how the findings of the steady-state plots apply to the change in concentrations over time and are biologically relevant. Stochastic noise was included in the extended equations because of the findings of multiple papers that there is noise present in the transcription of transcription factors, therefore having noise present in their concentration . By comparing the time plots, one can see that, even with the presence of noise in the concentrations, the findings of the steady-state plots still apply to this biological system.

Model Description

We model the transcription of target genes which regulate the self-renewal or differentiation of stem cells, which are affected by outside signals which increase or decrease the concentrations of transcription factors OCT4, SOX2, and NANOG and the complexes they form. Four differential equations are used to model the changing concentrations of the three transcription factors and the OCT4-SOX2 heterodimer, which are affected by each of the other factors as well as signals A⁺ and B^-. The transcription of the target genes is then described by one of two models: an incoherent model involving one differential equation to describe either self-renewal and differentiation, depending on parameter choice, and a coherent model involving two differential equations, one describing the self-renewal target gene product concentration and the other the differentiation target gene product concentration.

Incoherent and Coherent Models

The incoherent and coherent models are both reasonable ways to represent this system.

Coherent and Incoherent Feedforward Motifs
(A) Coherent.
(B) Incoherent.
Figure and caption from Chickarmane et al.

The incoherent model involves a feedforward loop where OCT4-SOX2 is an activator for transcription of the self-renewal and differentiation genes and NANOG is a weak repressor for the self-renewal genes and a stronger repressor for the differentiation genes . This is incoherent because OCT4-SOX2 increases the concentration of NANOG and the self-renewal target genes while NANOG inhibits the concentration of the self-renewal target genes . The structure of this loop only allows the system to respond to persistent signals, since the strengths of activation and repression from OCT4-SOX2 and NANOG are delicately balanced .

The coherent model involves a feedforward loop with OCT4-SOX2 and NANOG acting as activators for the self-renewal genes and repressors for the differentiation genes . Since OCT4-SOX2 increases the concentration of NANOG and both are activators for the self-renewal genes, a short signal can cause this system to strongly react due to the cascading effect .

Assumptions

It is known that NANOG, the OCT4-SOX2 heterodimer, and the OCT4-SOX2-NANOG complex positively regulate OCT4, SOX2, and NANOG . OCT4-SOX2 and NANOG regulate the target gene transcription, in which the OCT4-SOX2 heterodimer must first bind to the DNA and then recruit NANOG in order for transcription to occur efficiently . Therefore, the combined effects of OCT4, SOX2, and NANOG regulate transcription.

Additionally, it is assumed that the signals A⁺ and B^- can accurately model all possible combinations of signals, since A⁺ activates OCT4 and SOX2 and B^- represses NANOG . This model also assumes that the signals and transcription factors do not directly interact, since the binding region of the signals is far from the transcription factor binding sites .

Equations

Transcription Factors

The four equations below govern the concentration of the three transcription factors, OCT4, SOX2, and NANOG, and the heterodimer OCT4-SOX2. The equations for OCT4, SOX2, and NANOG are based on the Shea-Ackers rate equation model, such that part of the equation is based on a ratio between the transcription rate of the factor based on outside factors over the degradation rate. These equations also include terms that model the effect of constant degradation rates and formation and dissociation of the OCT4-SOX2 complex from OCT4 and SOX2. The equation for OCT4-SOX2 is based on the formation, dissociation, and degradation rates of OCT4-SOX2 based on constants multiplied by the concentration of the complex. Each term is explained in a table following the equation.

OCT4

$\frac{d[O]}{dt} = \frac{\eta_1 + a_1 [A ^ + ] + a_2 [OS] +a_3 [OS] [N]}{1 + \eta_2 + b_1 [A ^ + ] + b_2 [OS] + b_3 [OS] [N]} - \gamma_1 [O] - k_{1c} [O] [S] - k_{2c} [OS]$

Table 1

$\eta_1$	$a_1 [A ^ + ]$	$a_2 [OS]$	$a_3 [OS] [N]$	$1 + \eta_2$	$b_1 [A ^ + ]$	$b_2 [OS]$	$b_3 [OS] [N]$	$\gamma_1 [O]$	$k_{1c} [O] [S]$	$k_{2c} [OS]$
Basal Transcription Rate	Transcription Rate from A+	Transcription Rate from OCT4-SOX2	Transcription Rate from OCT4-SOX2-NANOG Complex	Basal Degradation Rate	Degradation Rate from A+	Degradation Rate from OCT4-SOX2	Degradation Rate from OCT4-SOX2-NANOG Complex	Constant Degradation of OCT4	Formation Rate of OCT4-SOX2	Dissociation Rate of OCT4-SOX2

SOX2

$\frac{d[S]}{dt} = \frac{\eta_3 + c_1 [A ^ + ] + c_2 [OS] + c_3 [OS] [N]}{1 + \eta_4 + d_1 [A ^ + ] + d_2 [OS] + d_3 [OS] [N]} - \gamma_2 [S] - k_{1c} [O] [S] - k_{2c} [OS]$

Table 2

$\eta_3$	$c_1 [A ^ + ]$	$c_2 [OS]$	$c_3 [OS] [N]$	$1 + \eta_4$	$d_1 [A ^ + ]$	$d_2 [OS]$	$d_3 [OS] [N]$	$\gamma_2 [S]$	$k_{1c} [O] [S]$	$k_{2c} [OS]$
Basal Transcription Rate	Transcription Rate from A+	Transcription Rate from OCT4-SOX2	Transcription Rate from OCT4-SOX2-NANOG Complex	Basal Degradation Rate	Degradation Rate from A+	Degradation Rate from OCT4-SOX2	Degradation Rate from OCT4-SOX2-NANOG Complex	Constant Degradation of SOX2	Formation Rate of OCT4-SOX2	Dissociation Rate of OCT4-SOX2

NANOG

$\frac{d[N]}{dt} = \frac{\eta_5 + e_1 [OS] + e_2 [OS] [N]}{1 + \eta_6 + f_1 [OS] + f_2 [OS] [N] + f_3 [B ^ - ]} - \gamma_3 [N]$

Table 3

$\eta_5$	$e_1 [OS]$	$e_2 [OS] [N]$	$1 + \eta_6$	$f_1 [OS]$	$f_2 [OS] [N]$	$f_3 [B ^ - ]$	$\gamma_3 [N]$
Basal Transcription Rate	Transcription Rate from OCT4-SOX2	Transcription Rate from OCT4-SOX2-NANOG Complex	Basal Degradation Rate	Degradation Rate from OCT4-SOX2	Degradation Rate from OCT4-SOX2-NANOG Complex	Degradation Rate from B-	Constant Degradation of NANOG

OCT4-SOX2

$\frac{d[OS]}{dt} = k_{1c} [O] [S] - k_{2c} [OS] - k_{3c} [OS]$

Table 4

$k_{1c} [O] [S]$	$k_{2c} [OS]$	$k_{3c} [OS]$
Formation Rate of OCT4-SOX2	Dissociation Rate of OCT4-SOX2	Degradation Rate of OCT4-SOX2

Target Genes

There are two sets of equations that are used to describe the concentration of the target gene product. The first is based on the incoherent model, which integrates the self-renewal and differentiation genes concentrations into one equation. The second is based on the coherent model and uses two separate equations.

The incoherent model assumes that OCT4-SOX2 acts as an activator for both the self-renewal and differentiation genes and that NANOG weakly represses the self-renewal genes and strongly represses the differentiation genes. The equation is based on a rate equation model similar to equations above, and the concentration is affected by OCT4-SOX2 and OCT4-SOX2-NANOG, with a degradation rate proportional to the concentration of the product. This model uses changes in parameters to model either the differentiation or the self-renewal genes.

In the coherent model, OCT4-SOX2 and NANOG are activators for the self-renewal genes and repressors for the differentiation genes. The two equations governing the self-renewal or differentiation genes for the coherent model are very similar to the incoherent model, but, unlike in the incoherent model, OCT4-SOX2-NANOG positively affects the transcription of the self-renewal target genes and OCT4-SOX2 only has a negative effect on the differentiation target genes. Each term is further explained in the table following each equation.

Incoherent Target Genes

$\frac{d[TG]}{dt} = \frac{\eta_7 + g_1 [OS]}{1 + \eta_8 + h_1 [OS] + h_2 [OS] [N]} - \gamma_4 [TG]$

Table 5

$\eta_7$	$g_1 [OS]$	$1 + \eta_8$	$h_1 [OS]$	$h_2 [OS] [N]$	$\gamma_4 [TG]$
Basal Transcription Rate	Transcription Rate from OCT4-SOX2	Basal Degradation Rate	Degradation Rate from OCT4-SOX2	Degradation Rate from OCT4-SOX2-NANOG Complex	Constant Degradation of Target Genes

Coherent Self-Renewal Target Genes

$\frac{d[TG_{self-renewal}]}{dt} = \frac{\eta_9 + m_1 [OS] + m_2 [OS] [N]}{1 + \eta_{10} + n_1 [OS] + n_2 [OS] [N]} - \gamma_5 [TG_{self-renewal}]$

Table 6

$\eta_9$	$m_1 [OS]$	$m_2 [OS] [N]$	$1 + \eta_{10}$	$n_1 [OS]$	$n_2 [OS] [N]$	$\gamma_5 [TG_{self-renewal}]$
Basal Transcription Rate	Transcription Rate from OCT4-SOX2	Transcription Rate from OCT4-SOX2-NANOG Complex	Basal Degradation Rate	Degradation Rate from OCT4-SOX2	Degradation Rate from OCT4-SOX2-NANOG Complex	Constant Degradation of Self-Renewal Target Genes

Coherent Differentiation Target Genes

$\frac{d[TG_{differentiation}]}{dt} = \frac{\eta_{11}}{1 + \eta_{12} + q_1 [OS] + q_2 [OS] [N]} - \gamma_6 [TG_{differentiation}]$

Table 7

$\eta_{11}$	$1 + \eta_{12}$	$q_1 [OS]$	$q_2 [OS] [N]$	$\gamma_6 [TG_{differentiation}]$
Basal Transcription Rate	Basal Degradation Rate	Degradation Rate from OCT4-SOX2	Degradation Rate from OCT4-SOX2-NANOG Complex	Constant Degradation of Differentiation Target Genes

Extension Equations

For the extension, the equations for OCT4-SOX2 and the target genes are equilvalent to those above and the equations for each of the three transcription factors newly include Gaussian white noise. This noise accounts for the inherent randomness in the transcription level of these proteins.

ξ(0, 1) gives a value from the Gaussian distribution centered at 0 with a amplitude of 1, which was chosen based on the range of amplitudes of 0.3 to 12 based on outside literature . This noise is applied to each of the three transcription factors since there is inherently noise in the transcription of each.

OCT4

$\frac{d[O]}{dt} = \frac{\eta_1 + a_1 [A ^ + ] + a_2 [OS] +a_3 [OS] [N]}{1 + \eta_2 + b_1 [A ^ + ] + b_2 [OS] + b_3 [OS] [N]} - \gamma_1 [O] - k_{1c} [O] [S] - k_{2c} [OS] + \xi(0, 1) [O]$

Table 8

$\eta_1$	$a_1 [A ^ + ]$	$a_2 [OS]$	$a_3 [OS] [N]$	$1 + \eta_2$	$b_1 [A ^ + ]$	$b_2 [OS]$	$b_3 [OS] [N]$	$\gamma_1 [O]$	$k_{1c} [O] [S]$	$k_{2c} [OS]$	$\xi(0, 1) [O]$
Basal Transcription Rate	Transcription Rate from A+	Transcription Rate from OCT4-SOX2	Transcription Rate from OCT4-SOX2-NANOG Complex	Basal Degradation Rate	Degradation Rate from A+	Degradation Rate from OCT4-SOX2	Degradation Rate from OCT4-SOX2-NANOG Complex	Constant Degradation of OCT4	Formation Rate of OCT4-SOX2	Dissociation Rate of OCT4-SOX2	Formation Noise from Gaussian White Noise with Mean 0 and Amplitude 1

SOX2

$\frac{d[S]}{dt} = \frac{\eta_3 + c_1 [A ^ + ] + c_2 [OS] + c_3 [OS] [N]}{1 + \eta_4 + d_1 [A ^ + ] + d_2 [OS] + d_3 [OS] [N]} - \gamma_2 [S] - k_{1c} [O] [S] - k_{2c} [OS] + \xi(0, 1) [S]$

Table 9

$\eta_3$	$c_1 [A ^ + ]$	$c_2 [OS]$	$c_3 [OS] [N]$	$1 + \eta_4$	$d_1 [A ^ + ]$	$d_2 [OS]$	$d_3 [OS] [N]$	$\gamma_2 [S]$	$k_{1c} [O] [S]$	$k_{2c} [OS]$	$\xi(0, 1) [S]$
Basal Transcription Rate	Transcription Rate from A+	Transcription Rate from OCT4-SOX2	Transcription Rate from OCT4-SOX2-NANOG Complex	Basal Degradation Rate	Degradation Rate from A+	Degradation Rate from OCT4-SOX2	Degradation Rate from OCT4-SOX2-NANOG Complex	Constant Degradation of SOX2	Formation Rate of OCT4-SOX2	Dissociation Rate of OCT4-SOX2	Formation Noise from Gaussian White Noise with Mean 0 and Amplitude 1

NANOG

$\frac{d[N]}{dt} = \frac{\eta_5 + e_1 [OS] + e_2 [OS] [N]}{1 + \eta_6 + f_1 [OS] + f_2 [OS] [N] + f_3 [B ^ - ]} - \gamma_3 [N] + \xi(0, 1) [N]$

Table 10

$\eta_5$	$e_1 [OS]$	$e_2 [OS] [N]$	$1 + \eta_6$	$f_1 [OS]$	$f_2 [OS] [N]$	$f_3 [B ^ - ]$	$\gamma_3 [N]$	$\xi(0, 1) [N]$
Basal Transcription Rate	Transcription Rate from OCT4-SOX2	Transcription Rate from OCT4-SOX2-NANOG Complex	Basal Degradation Rate	Degradation Rate from OCT4-SOX2	Degradation Rate from OCT4-SOX2-NANOG Complex	Degradation Rate from B-	Constant Degradation of NANOG	Formation Noise from Gaussian White Noise with Mean 0 and Amplitude 1

Method

This model is replicated in Mathematica and is used to create a series of one- and two-parameter bifurcation plots. This is done by setting each of the equations equal to zero and using Mathematica to solve the system for lines of equilibria using Solve. These lines are then plotted by varying inputs, including A⁺, B^-, and OCT4-SOX2, thus creating bifurcation plots. These diagrams show the bistable behavior of the system and how varying parameter values changes whether or not the system is capable of bistability.

When creating the time plots without noise for the extension, the sets of differential equations and initial conditions were solved numerically using NDSolve. These were then plotted with varying time. The time plots with noise were created using ItoProcess, which gives Gaussian noise, to create the random process for each equation and then evaluated using RandomFunction set to give 5 plots. Each of these was plotted with varying time.

State Variables and Parameters

All state variables, parameters, and constants are listed in the two tables below, along with descriptions and typical initial value or range of values. If a parameter has a range of values, this is either because the parameter is set to several values to simulate different conditions or the parameter is varied over that range to create a two-parameter bifurcation plot. The initial values were used for the extension, since initial values are not used to create bifurcation diagrams.

Table 11: State Variables

Variable	Description	Initial Value
[O]	Concentration of OCT4	5 or 10
[S]	Concentration of SOX2	5 or 10
[N]	Concentration of NANOG	5 or 10
[OS]	Concentration of OCT4-SOX2	0
[TG]	Concentration of the Target Gene Product	0

Table 12: Parameters and Constants

Parameter	Description	Range of Values
[A+]	Concentration of Signal A+	0 - 600
[B-]	Concentration of Signal B-	0 - 120
k1c	Formation Rate of OCT4-SOX2 Complex	0 - 10
k2c	Dissociation Rate of OCT4-SOX2 Complex	0 - 10
k3c	Degradation Rate of OCT4-SOX2 Complex	0 - 10
a1, a2, a3	Binding Strength between OCT4 and A+, OCT4-SOX2, and OCT4-SOX2-NANOG, respectively	0 - 10
b1, b2, b3	Repression Strength between OCT4 and A+, OCT4-SOX2, and OCT4-SOX2-NANOG, respectively	0.00055 - 0.011
c1, c2, c3	Binding Strength between SOX2 and A+, OCT4-SOX2, and OCT4-SOX2-NANOG, respectively	0 - 1
d1, d2, d3	Repression Strength between SOX2 and A+, OCT4-SOX2, and OCT4-SOX2-NANOG, respectively	0.00055 - 0.011
e1, e2	Binding Strength between NANOG and OCT4-SOX2 and OCT4-SOX2-NANOG, respectively	0 - 10
f1, f2, f3	Repression Strength between NANOG and OCT4-SOX2, OCT4-SOX2-NANOG, and B-, respectively	0 - 0.05
g1	Binding Strength between the Target Genes and OCT4-SOX2	0.1 - 6
h1, h2	Repression Strength between the Target Genes and OCT4-SOX2 and OCT4-SOX2-NANOG	0 - 0.5
m1, m2	Activation Strength between the Target Genes and OCT4-SOX2 and OCT4-SOX2-NANOG for Self-Renewal	0.1
n1, n2	Repression Strength between the Target Genes and OCT4-SOX2 and OCT4-SOX2-NANOG for Self-Renewal	0.001 - 0.01
q1, q2	Repression Strength between the Target Genes and OCT4-SOX2 and OCT4-SOX2-NANOG for Differentiation	0.001 - 0.01
η1	Basal Transcription Rate of OCT4	0 - 10
η2	Basal Degradation Rate of OCT4	10-7
η3	Basal Transcription Rate of SOX2	0.0001
η4	Basal Degradation Rate of SOX2	10-7
η5	Basal Transcription Rate of NANOG	0 - 10
η6	Basal Degradation Rate of NANOG	10-7
η7	Basal Transcription Rate of the Target Gene Product	10-4; 10-5
η8	Basal Degradation Rate of the Target Gene Product	10-7
η9	Basal Transcription Rate of the Target Gene Product for Self-Renewal	0.0001
η10	Basal Degradation Rate of the Target Gene Product for Self-Renewal	10-7
η11	Basal Transcription Rate of the Target Gene Product for Differentiation	1
η12	Basal Degradation Rate of the Target Gene Product for Differentiation	0.001
γ1	Degradation Constant of OCT4	0 - 5
γ2	Degradation Constant of SOX2	1
γ3	Degradation Constant of NANOG	0 - 2
γ4	Degradation Constant of the Target Gene Product	0.01
γ5	Degradation Constant of the Target Gene Product for Self-Renewal	0.05
γ6	Degradation Constant of the Target Gene Product for Dissociation	0.01

Results

This paper hypothesizes that the self-renewal or differentiation of a stem cell depends on the concentration of outside signals and that this behavior exhibits bistability. Further, by either changing the binding and repressive strength between OCT4-SOX2 and the transcription factors or instead increasing the transcription rate of NANOG, this will become an irreversible switch. The differential equations created to model this system were evaluated in Mathematica and used to produce a series of figures which support this hypothesis.

Figure 3

Bifurcation plots of OCT4-SOX2 and NANOG concentration as dependent on A+ concentration from Chickarmane et al. The arrows and "SN" refer to the two saddle node bifurcations.

Bifurcation plot of OCT4-SOX2 concentration from A+ concentration in Mathematica Bifurcation plot of NANOG concentration from A+ concentration in Mathematica

Figure 3 includes two steady-state plots, or bifurcation diagrams, of OCT4-SOX2 and NANOG concentration with varying concentrations of signal A⁺ using the parameters from Table 24. The concentration of signal B^- is set to 0.01 as specified by the paper. These are created by setting the four differential equations governing the concentrations of the transcription factors to 0 and solving for the concentrations of the transcription factors.

For A⁺ from about 0 to 87 (arbitrary units), the system is "off". When A⁺ is increased beyond about 87, the system turns "on", meaning the concentrations of the transcription factors OCT4, SOX2, and NANOG are high. The system remains "on" until the concentration of A⁺ is decreased below about 58. The figures created in Mathematica match the figure from the paper. This system shows an example of a hysteretic curve, where there are two saddle node bifurcations where the equilibrium changes from stable to unstable. This demonstrates that the system includes a bistable switch dependent on the concentration of A⁺, one key piece of the hypothesis.

Figure 4

Bifurcation plots of OCT4-SOX2 and NANOG concentration as dependent on B- concentration from Chickarmane et al.

Bifurcation plot of OCT4-SOX2 concentration from B- concentration in Mathematica Bifurcation plot of NANOG concentration from B- concentration in Mathematica

Figure 4 is the same as Figure 3 but dependent on the concentration of signal B^- instead, also using the parameters from Table 24. The concentration of A⁺ is set to 100 as specified by the paper.

For B^- from about 0 to 34, the system remains "on". When B^- is increased beyond 34, the system turns "off", so the concentrations of the transcription factors fall. As B^- is then decreased beyond 15, the system returns to its "on" state. This figure and these findings match those from the paper, and this shows that the bistability is also found through a negatively regulating signal, showing the second piece of the first part of the hypothesis.

Figure 6

Bifurcation plots of the concentration of NANOG and target genes as dependent on the concentration of OCT4-SOX2 with and without NANOG autoregulation from Chickarmane et al. Figure A shows NANOG concentration with and without NANOG autoregulation, and Figure B shows target gene concentration for the incoherent model. Figure C shows differentiation and self-renewal target gene concentration for the coherent model with NANOG autoregulation.

Bifurcation plot of NANOG concentration against OCT4-SOX2 concentration with and without NANOG autoregulation in Mathematica Bifurcation plot of target gene concentration from the incoherent model against OCT4-SOX2 concentration with and without NANOG autoregulation in Mathematica Bifurcation plot of differentiation and self-renewal target gene concentration from the incoherent model against OCT4-SOX2 concentration with autoregulation in Mathematica

Figure 6 includes three steady-state plots of the concentrations of NANOG or differentiation and self-renewal target genes as a function of the concentration of OCT4-SOX2. For this, the concentration of OCT4-SOX2 is assumed to be constant. Figure 6A was formed by setting the differential equations of OCT4, SOX2, and NANOG to 0 and solving. Figures 6B and 6C were formed by setting the differential equations of OCT4, SOX2, NANOG, and the incoherent target gene concentrations to 0 and solving. The concentration of B^- is held at 0.1.

Figure 6A shows the change in concentration of NANOG based on NANOG autoregulation being present or not, and Figure 6B shows the change in concentration of the differentiation target genes. When no NANOG autoregulation is present, the values of e₂ and f₂, the values of the binding and repression strength of NANOG to OCT4-SOX2-NANOG, are set to 0. The "FF" in the legend of the figure from the paper stands for feedforward, the type of loop studied in this example. The target gene concentrations are modeled using the incoherent model for this figure, with the differentiation target gene concentrations modeled with h₂ = 0.05 and the self-renewal target gene concentrations modeled with h₂ = 0.001 to achieve the different results. The value of f₃ is set to 0 to ignore the effects of B^-, which has the same effect as holding B^- constant. The parameters are found in Tables 24 and 25 with the following modifications:

Table 13

Parameters for All	Value
f3	0
6A No Autoregulation	Value
e2	0
f2	0
6B Autoregulation	Value
h1	0
h2	0.09
6B No Autoregulation	Value
e2	0
f2	0
h1	0
h2	0.09
Self-Renewal Target Genes	Value
h2	0.001

The lack of NANOG autoregulation in Figure 6A demonstrates that NANOG autoregulation is necessary for NANOG to increase when OCT4-SOX2 increases. Then, in Figure 6B, it is shown that without NANOG autoregulation, the differentiation genes continue to have a high concentration as OCT4-SOX2 is increased, so the autoregulation is necessary to achieve the switch between differentiation and self-renewal in stem cells. Finally, Figure 6C shows how the concentration of OCT4-SOX2 affects the concentration of differentiation or stem cell target genes, illustrating how as the concentration of OCT4-SOX2 increases, the concentration of the differentiation genes falls and the concentration of the self-renewal genes rises. The next figure will confirm that an increase in the signal that causes the concentration of OCT4-SOX2 to rise increases the concentration of the self-renewal genes. While this figure does not directly support the hypothesis, it demonstrates that the model was created correctly, since NANOG must autoregulate itself.

When using the specified parameter values for Figure 6B, there are discrepancies between the figure in the paper and the figure created in Mathematica. The curve found using the specified values has peaks that are higher and for larger concentrations of A⁺. The three parameters possibly responsible for this discrepancy are g₁, h₁, and h₂. The parameters g₁ and h₁, which are the binding and repression strengths between the target genes and OCT4-SOX2 in the incoherent model, both raise or lower the maximum for the two curves. By increasing h₂, which is the repression strength between the target genes and OCT4-SOX2-NANOG, the peaks of the curves moves to the left, so the curves peak for lower concentrations. By raising h₂ to 0.09, the correct location of the peaks are obtained, and then h₁ can be changed to 0 while leaving g₁ the same to obtain the correct curves. The differences between these two figures are shown here:

Bifurcation plot of target gene concentration using the specified parameters from the incoherent model against OCT4-SOX2 concentration with and without NANOG autoregulation in Mathematica

Bifurcation plot of target gene concentration using the modified parameters from the incoherent model against OCT4-SOX2 concentration with and without NANOG autoregulation in Mathematica

While these changes produce the correct figure, these modifications do not make sense. The value of h₁ should not equal 0, since this would mean that OCT4-SOX2 has no effect on the repression of the target genes, which is false based on the assumptions governing all other results. The value of h₂ is used to determine if the differentiation or self-renewal target genes are being shown, so h₂ should only equal 0.05 or 0.001. This suggests that the parameters specified in the paper are actually correct, though they produce a different figure

Figure 7

Bifurcation plots of the differentiation and self-renewal target gene concentrations from the incoherent model as dependent on A+ concentration from Chickarmane et al.

Bifurcation plot of the differentiation target gene concentration from the incoherent model as dependent on A+ concentration in Mathematica Bifurcation plot of the self-renewal target gene concentration from the incoherent model as dependent on A+ concentration in Mathematica

Figure 7 includes steady-state plots of the concentration of the differentiation and self-renewal genes using the incoherent model as a function of varying concentration of signal A⁺. These plots were made by setting the differential equations for the transcription factors and incoherent target gene equation to 0. The concentration of B^- is held at 20.

Using the parameters from Tables 1 and 2, the following modifications were made. The value of f₃ is set to 0 to ignore the effects of B^-, and the values of h₁ and η₇ were set to the values specified below to better shape the figure. The value of h₂ is set to 1 to represent the differentiation genes and to 0.001 to represent the self-renewal genes. The changes are summarized in the following table:

Table 14

Parameters for Both	Value
f3	0
h1	10-3
η7	10-4
Figure 7A	Value
h2	1
Figure 7B	Value
h2	0.001

These two figures show how changing the concentration of signal A⁺ leads to either increased concentration of the self-renewal target genes or the differentiation target genes. For this figure, when A⁺ is increased past about 88, the concentration of the differentiation target genes goes to close to zero while the concentration of the self-renewal genes remains high. Then, when A⁺ is decreased past about 67, the concentration of the self-renewal genes starts to fall and the concentration of the differentiation genes rises. This supports the hypothesis that an outside signal causes the change between differentiation and self-renewal target gene expression with a bistable switch for the incoherent model.

Figure 8

Bifurcation plots of the differentiation and self-renewal target gene concentrations from the coherent model as dependent on A+ concentration from Chickarmane et al.

Bifurcation plot of the self-renewal target gene concentration from the coherent model as dependent on A+ concentration in Mathematica Bifurcation plot of the differentiation target gene concentration from the coherent model as dependent on A+ concentration in Mathematica

Figure 8 includes steady-state plots of the concentration of the differentiation and self-renewal genes based on the coherent model instead of the incoherent model as a function of varying concentration of signal A⁺. The differential equations for the transcription factors and either the self-renewal or differentiation target gene concentration equation were set to 0, using the parameters from Tables 1 and 3. The value of f₃ is set to 0 to ignore the effects of B^-, as described above. The concentration of B^- is held at 20.

Table 15

Parameters	Value
f3	0

These two plots also demonstrate how the concentration of the self-renewal and differentiation genes changes with the concentration of signal A⁺. While these plots look qualitatively different from Figure 7, the two switches of the bistable system occur at the same concentrations of about 67 and 88. This shows that both the incoherent and coherent models fulfill the same bistable switch for the same concentrations of A⁺. Therefore, these figures cannot be used to decide if the incoherent or coherent model is correct, so further experimentation and modeling will be necessary. However, this figure still supports the hypothesis that outside signals within a system with a bistable switch determine the fate of the stem cell.

Bifurcation plots of the differentiation target gene concentration from the incoherent model as dependent on A+ concentration Bifurcation plots of the self-renewal target gene concentration from the incoherent model as dependent on A+ concentration

Bifurcation plots of the differentiation target gene concentration from the incoherent model as dependent on A+ concentration Bifurcation plots of the self-renewal target gene concentration from the incoherent model as dependent on A+ concentration

To best compare the effects of the incoherent and coherent models, the scale of A⁺ concentration was adjusted for Figure 8 and the figures were placed side-by-side. While the amplitude of the curves differs between models, the switches still occur at the same values of A⁺ equals about 58 and 88, demonstrating that both models show the same bistable effect based on signal concentration. This does not allow us to determine which model is correct.

Figure 9

Figure A is a bifurcation plot of the transcription factor concentrations as dependent on A+ concentration with weak NANOG feedback, and Figure B is a bifurcation plot of the differentiation target gene concentration from the incoherent model as dependent on A+ concentration with weak and strong NANOG feedback from Chickarmane et al.

Bifurcation plot of the transcription factor concentrations as dependent on A+ concentration with weak NANOG feedback in Mathematica Bifurcation plot of the differentiation target gene concentration from the incoherent model as dependent on A+ concentration with weak and strong NANOG feedback in Mathematica

Figure 9 includes steady-state plots showing and comparing the effects of weak NANOG feedback on the concentrations of the transcription factors and differentiation target gene concentration. The concentration of B^- is held at 0.1.

For both figures, changes to the parameters from Tables 1 and 2 were made to e₂ and f₂ to adjust the size of the graphs, and f₃ was set to 0 to ignore the effect of B^-. Figure 9A and the curve for weak NANOG feedback in Figure 9B were created by setting the differential equations for the transcription factors and target genes of the coherent model to 0. The values of a₃ and c₃ were set to 0.05, lowering the binding strengths between OCT4 and SOX2 to OCT4-SOX2-NANOG, and the values of b₃ and d₃ were set to 5.5*10^-4, lowering the repression strength of OCT4-SOX2-NANOG on OCT4 and SOX2. This gives the effect of lowering NANOG feedback to OCT4 and SOX2. For the curve for strong NANOG feedback in Figure 9B, the same equations were used, but the values of a₃ and c₃ were set to 1 and b₃ and d₃ to 1.5*10^-3, increasing the strength of NANOG feedback.

Table 16

Parameters for All	Value
e2	0.025
f2	9.25*10-4
f3	0
Figure 9A and Figure 9B Weak	Value
a3	0.05
b3	5.5*10-4
c3	0.05
d3	5.5*-4
Figure 9B Strong	Value
a3	1
b3	1.5*10-3
c3	1
d3	1.5*-3

Figure 9A shows that the bistable switch for the concentrations of the transcription factors is not present when NANOG feedback is lowered, as the concentrations steadily increase instead. This shows that NANOG feedback is a necessary part of the model to show the bistable behavior. Figure 9B compares the bistable switch that is present for strong NANOG feedback to the lack of hysteresis that is present for weak NANOG feedback. When NANOG feedback is weak, the concentration of the differentiation target genes lowers as the concentration of signal A⁺ rises, but that concentration never falls low enough that the differentiation genes are not expressed, further demonstrating that NANOG feedback is a significant part of this system.

Figure 10

Bifurcation plots of OCT4-SOX2 and NANOG concentration as dependent on A+ concentration with modified binding and repression strength from OCT4-SOX2-NANOG on OCT4, SOX2, and NANOG from Chickarmane et al.

Bifurcation plot of OCT4-SOX2 concentration from A+ concentration with modified binding and repression strength from OCT4-SOX2-NANOG in Mathematica Bifurcation plot of NANOG concentration from A+ concentration in with modified binding and repression strength from OCT4-SOX2-NANOG in Mathematica

Figure 10 is similar to Figure 3, but the binding and repression strengths from OCT4-SOX2-NANOG on the other transcription factors are changed, leading to an irreversible switch. The differential equations for the transcription factors are set to zero, and the concentration of B^- is held at 0.1.

The main set of parameters are found in Table 1. The values of a₃ and c₃ are set to 0.5, increasing the binding strength of OCT4-SOX2-NANOG and OCT4 and SOX2, and the values of b₃ and d₃ are set to 0.001, decreasing the repression strength of OCT4-SOX2-NANOG on OCT4 and SOX2. The value of e₁ is set to 0.01, increasing the binding strength between NANOG and OCT4-SOX2. The value of f₂ is set to 0.001, increasing the repression strength between NANOG and OCT4-SOX2-NANOG slightly, and the value of f₃ is set to 0.05, increasing the strength of repression of NANOG due to the signal B^-.

Table 17

Parameters for All	Value
a3	0.5
b3	0.001
c3	0.5
d3	0.001
e1	0.01
f2	0.001
f3	0.05

In Figure 10, when the concentration of A⁺ increases beyond about 81, the system is switched "on" and remains in that state, even when the signal A⁺ is removed. Therefore, by modifying the binding and repression strengths of OCT4-SOX2-NANOG to increase its binding efficiency, stem cells can be induced to differentiate with a short burst of signal A⁺ instead of continuous input. This supports the hypothesis that the bistable switch can be made irreversible.

Figure 11

Bifurcation plots of OCT4-SOX2 and NANOG concentration as dependent on B- concentration with modified binding and repression strength from OCT4-SOX2-NANOG from Chickarmane et al.

Bifurcation plot of OCT4-SOX2 concentration from B- concentration with modified binding and repression strength from OCT4-SOX2-NANOG in Mathematica Bifurcation plot of NANOG concentration from B- concentration with modified binding and repression strength from OCT4-SOX2-NANOG in Mathematica

Figure 11 is similar to Figure 4 but with the modifications explained above for Figure 10. The concentration of A⁺ is held at 100.

Table 18

Parameters for All	Value
a3	0.5
b3	0.001
c3	0.5
d3	0.001
e1	0.01
f2	0.001
f3	0.05

Figure 11 shows how the concentration of signal B^- overpowers the input concentration of signal A⁺, which is held constant at a level high enough to turn the system "on" with low signal B^- present. When the concentration of B^- is increased beyond about 42, the system switches "off" and cannot be turned "on" again by lowering the concentration of signal B^-. When used with Figure 11, this supports the part of the hypothesis that claims that by increasing the binding efficiency of OCT4-SOX2-NANOG to the transcription factors allows the fate of the stem cell to be chosen by directly modifying the concentration of signals for a short time.

While the steady-state line where the concentrations are close to 0 ends at B^- = 100 in the figure from the paper, that portion of the diagram should continue for all values of B^- as shown in the figures created in Mathematica, since the concentration will always be low for high values of B^-.

Figure 12

Bifurcation plots of OCT4-SOX2 and NANOG concentration as dependent on A+ concentration with high and low basal NANOG transcription rates from Chickarmane et al.

Bifurcation plot of OCT4-SOX2 and NANOG concentration as dependent on A+ concentration with a high basal NANOG transcription rate in Mathematica Bifurcation plot of OCT4-SOX2 and NANOG concentration as dependent on A+ concentration with a low basal NANOG transcription rate in Mathematica

Figure 12 includes the steady-state plots of the concentrations of OCT4-SOX2 and NANOG as a function of A⁺ concentration for high and low basal NANOG transcription rates. This required setting the transcription factor differential equations to zero and solving, and the main set of parameters used are found in Table 1. The concentration of B^- is held at 0.1.

For Figure 12A, the value of η₅ was set to 35, so the basal transcription rate of NANOG was increased, and the value of η₆ was set to 0.035, so basal degradation rate of NANOG was increased. Combined, these realistically simulate a higher basal transcription rate for NANOG. The value of f₃ was set to 0 as described above.

Table 19

Parameters for Both	Value
f3	0
Parameters for Figure 12A	Value
η5	35
η6	0.035

Similar to Figures 10 and 11, Figure 12A shows a hysteresis curve with only one saddle node present for positive concentrations of A⁺, meaning the system undergoes an irreversible switch to "on" when the concentration of A⁺ is increased beyond about 11 in a stem cell with a high basal transcription rate of NANOG. Figure 12B shows the lower basal transcription rate used for all other figures, with the two switches occurring at about 58 and 87 as above. By comparing these two figures, it is shown that by increasing the transcription rate of NANOG, it is also possible to acquire a stem cell that will differentiate on demand similar to the system created in Figures 10 and 11 by changing the binding efficiency of OCT4-SOX2-NANOG, as stated in the hypothesis.

Figure S1

Two-parameter bifurcation plots showing the concentration of A+ where saddle node bifurcations occur for varying values of a3, a2, e2, and e1, which are binding strengths of the transcription factors to the transcription factor complexes, from Chickarmane et al.

Two-parameter bifurcation plots showing the concentration of A+ where saddle node bifurcations occur for varying values of binding strengths of the transcription factors to the transcription factor complexes in Mathematica

Figure S1 is a two-parameter bifurcation plot that plots the values of a₃, a₂, e₂, and e₁ against the concentration of A⁺ where the switches occur. This was created by solving the system of transcription factor differential equations set to 0 without the specified parameter defined, using the parameters from Table 1 without other modifications, and then finding the value of A⁺ that gives the minimum of the function plotting the concentration of OCT4-SOX2 against A⁺ for varying values of the specified parameter.

The value of a₃ corresponds to the binding strength between OCT4 and OCT4-SOX2-NANOG, and a₂ corresponds to the binding strength between OCT4 and OCT4-SOX2. The value of e₂ corresponds to the binding strength between NANOG and OCT4-SOX2-NANOG, and e₁ corresponds to the binding strength between NANOG and OCT4-SOX2. Each of these is varied between 0 and 10. The concentration of B^- was set to 0.1.

Where there are two values of the concentration of signal A⁺ for a value of the parameter, a bistable switch is present. When only one value of A⁺ exists, an irreversible switch is present, and when no values of A⁺ exist, then there is no bifurcation. For a₃, a bistable switch is present for a₃ between about 0.1 and 0.3, with no switch below 0.1 and an irreversible switch about 0.3. This shows that increasing the binding strength between OCT4 and OCT4-SOX2-NANOG produces a stem cell where differentiation can be induced with short bursts of signals, as shown in Figures 10 and 11. For a₂, a bistable switch was found to be present for a₂ up to about 8.6, though the paper shows that this bistable switch is present for all values of a₂ shown. This discrepancy is likely the result of the paper not reporting the correct values of some parameters, but the figures are still very similar.

For e₂, a bistable switch is present for values of e₂ very close to 0.1, with no bifurcation below about 0.1 and an irreversible switch above about 0.1. This shows that raising the binding strength between NANOG and OCT4-SOX2-NANOG produces a stem cell that can easily be switched between self-renewal and differentiation, as shown in Figures 10 and 11. For e₁, there is a bistable switch for values between about 0.1 and 1.8 and an irreversible switch above 1.8. The figure from the paper only shows one line for the values of e₁, but by exploring the bifurcation plots of OCT4-SOX2 for A⁺ with varying values of e₁, it is clear that there should be two values of A⁺ where switches occur for e₁ between 0.1 and 1.8. The effects of varying e₁ are not mentioned in the paper, and all parameter values were checked, so it is unknown why this line was found using Mathematica but not present in the figure in the paper.

Figure S2

Two-parameter bifurcation plots showing the concentration of A+ where saddle node bifurcations occur for varying values of k1c, k2c, and k3c, which are the formation, dissociation, and degradation rates of OCT4-SOX2, from Chickarmane et al.

Two-parameter bifurcation plots showing the concentration of A+ where saddle node bifurcations occur for varying values of the formation, dissociation, and degradation rates of OCT4-SOX2 in Mathematica

Figure S2 is a two-parameter bifurcation plot that plots the values of k_1c, k_2c, and k_3c against the concentration of A⁺ where the switches occur. This was created using the same method as Figure S1 with the same parameter values. The values of k_1c, k_2c, and k_3c correspond to the formation, dissociation, and degradation rates of OCT4-SOX2. The concentration of B^- was set to 0.1.

There is a bistable switch present for values of k_1c from 0 to about 0.5, after which there is an irreversible switch. This shows that increasing the formation rate of the OCT4-SOX2 complex beyond about 0.5 makes the stem cell receptive to short bursts of signal, similarly to Figure 10 through 12. For k_2c and k_3c, there are bistable switches for all values of A⁺ shown by this figure, but the curves imply that beyond particular values of k_2c and k_3c the switch-like behavior disappears, implying that the dissociation and degradation rates of OCT4-SOX2 cannot be too large for the system to be able to change between high self-renewal and differentiation target gene concentration.

The plot for k_3c is not smooth at about A⁺ equal to 20. When investigated, it was found that Mathematica cut the solution slightly short in that area for low values of k_3c, making the switch appear to occur at higher values of A⁺.

Figure S3

Two-parameter bifurcation plots showing the concentration of A+ where saddle node bifurcations occur for varying values of γ1 and γ3, which are the degradation rates of OCT4 and NANOG, from Chickarmane et al.

Two-parameter bifurcation plots showing the concentration of A+ where saddle node bifurcations occur for varying values the degradation rates of OCT4 and NANOG in Mathematica

Figure S3 is a two-parameter bifurcation plot that plots the values of γ₁ and γ₃ against the concentration of A⁺ where the switches occur. This was created using the same method as Figure S1 with the same parameter values. The values of γ₁ and γ₃ correspond to the degradation rates of OCT4 and NANOG, and γ₁ can be considered as equivalent to γ₂, the degradation rate of SOX2. The concentration of B^- was set to 0.1.

Both curves in this figure were cut off near the upper tips, due to the limitations of the accuracy of Mathematica using the FindMinimum function where two minimums are very close to each other. For γ₁, a bistable switch is present from 0 to about 1.2 by the figure created in Mathematica, or 1.4 by the figure from the paper, after which no switch is present. Similarly, for γ₃, a bistable switch is present from 0 to about 4.5 or above 5 according to either the figure created in Mathematica or from the paper, after which no switch is present. These show that when the degradation rates of the transcription factors are increased too high, it is not possible for the stem cell to switch between self-renewal and differentiation and both sets of target genes will be expressed.

Figure S4

Two-parameter bifurcation plots showing the concentration of A+ where saddle node bifurcations occur for varying values of η1 and η5, which are the basal transcription rates of OCT4 and NANOG, from Chickarmane et al.

Two-parameter bifurcation plots showing the concentration of A+ where saddle node bifurcations occur for varying values the basal transcription rates of OCT4 and NANOG in Mathematica

Figure S4 is a two-parameter bifurcation plot that plots the values of η₁ and η₅ against the concentration of A⁺ where the switches occur. This was created using the same method as Figure S1 with the same parameter values. The values of η₁ and η₅ correspond to the basal transcription rates of OCT4 and NANOG, and η₁ can be considered as equivalent to η₃, the basal transcription rate of SOX2. The concentration of B^- was set to 0.1.

The values of η₁ and η₅ are typically 10^-4, so these bifurcation plots show values of η₁ and η₅ significantly larger than the typical values. For these large values, an irreversible switch is present, which shows that increasing the transcription rate of the transcription factors leads to the easily-controlled stem cell as described above, as found in Figure 12 for NANOG.

Figure S5

Bifurcation plots of target gene concentration as dependent on OCT4-SOX2 concentration with high and low values of e2, g1, and h2, which are the binding strength of OCT4-SOX2-NANOG to NANOG gene, the binding strength of OCT4-SOX2 to the target genes, and the repression of the target genes by OCT4-SOX2-NANOG, from Chickarmane et al.

Bifurcation plots of target gene concentration as dependent on OCT4-SOX2 concentration with high and low values of e2 in Mathematica Bifurcation plots of target gene concentration as dependent on OCT4-SOX2 concentration with high and low values of g1 in Mathematica Bifurcation plots of target gene concentration as dependent on OCT4-SOX2 concentration with high and low values of h2 in Mathematica

Figure S5 is an extension of Figure 6B, and is thus created in the same manner, but with varying values of e₂, g₁, and h₂. The value of e₂ represents the binding strength between OCT4-SOX2-NANOG and NANOG, which has a normal value of 0.1 and has values of 0.001 and 0.0001 in the figure. The value of g₁ represents the binding strength between OCT4-SOX2 and the target genes, which has a normal value of 0.1 and has values of 0.1, 0.85, and 1 in the figure. The value of h₂ represents the repression strength of OCT4-SOX2-NANOG on the target genes, which has a normal value of 0.05 and has values of 0.05 and 0.5 in the figure. The parameters used are found in Tables 1 and 2, with the value of f₃ set to 0 as described above and with the concentration of B^- set to 0.1.

Table 20

Parameters	Value
f3	0

In Figure S5A, the larger value of e₂ allows the concentration of the differentiation genes to be lowered with higher concentrations of OCT4-SOX2, while the smaller value has continuous expression of the differentiation genes, showing that e₂ must be sufficiently large for the system to switch between self-renewal and differentiation.

In Figure S5B, the concentration of the differentiation genes falls after about 65 for all values of g₁. The larger values of g₁ cause a significantly larger concentration of the differentiation target genes at all concentrations of OCT4-SOX2, to the extent that the concentration of differentiation target genes may still be too high after the switch. The peak concentration of the target genes for g₁ = 1 is higher in the figure in Mathematica than from the paper, which is unusual since the peak concentration appears to be a function of the value of g₁, which is specified by the paper. When g₁ is set to 0.85 instead, the correct curve is obtained.

In Figure S5C, the larger value of h₂ leads to a smaller peak concentration of target genes than the smaller value, but both values exhibit a switch between expression of the differentiation genes and no expression. This shows that a lower repression strength of OCT4-SOX2-NANOG on the target genes causes greater expression of the target genes, as expected.

Extension Figure 1

Time plot of NANOG concentration for A+ = 100 and B- = 0.1	Time plot of OCT4-SOX2 concentration for A+ = 100 and B- = 0.1
Time plot of NANOG concentration with noise for A+ = 100 and B- = 0.1	Time plot of OCT4-SOX2 concentration with noise for A+ = 100 and B- = 0.1
Time plot of NANOG concentration for A+ = 10 and B- = 100	Time plot of OCT4-SOX2 concentration for A+ = 10 and B- = 100
Time plot of NANOG concentration with noise for A+ = 10 and B- = 100	Time plot of OCT4-SOX2 concentration for A+ = 10 and B- = 100
Time plot of NANOG concentration for A+ = 100 and B- = 100	Time plot of OCT4-SOX2 concentration for A+ = 100 and B- = 100
Time plot of NANOG concentration with noise for A+ = 100 and B- = 100	Time plot of OCT4-SOX2 concentration with noise for A+ = 100 and B- = 100

These time plots show the concentrations of NANOG and OCT4-SOX2 over time for three combinations of A⁺ and B^-. The first set of figures uses a high concentration of A⁺ and a low concentration of B^-, using 100 and 0.1 respectively. The second set uses low A⁺ and high B^-, setting them to 10 and 100. The final set uses high A⁺ and high B^-, setting both to 100. The time plots without use the four original equations governing the concentration of the three transcription factors and OCT4-SOX2, and the time plots with noise use the modified equations which add stochastic noise to the transcription factor concentrations. The parameters used to create these plots are found in Table 24, and the initial conditions are found below. These initial conditions were chosen to prevent the solutions from falling below 0, which is not biologically possible.

Table 21

State Variable	Initial Condition
[OCT4]	5
[SOX2]	5
[NANOG]	5
[OCT4-SOX2]	0

The steady-state concentration of NANOG and OCT4-SOX2 correspond to the values found in Figure 3 as expected. Additionally, while there is no steady-state for the time plots with noise, these plots still follow the time plots without noise instead of wildly deviating, confirming that this model is a good simplification of the system.

Extension Figure 2

Time plot of differentiation target gene concentration using the incoherent model for A+ = 100 and B- = 0.1	Time plot of self-renewal target gene concentration using the incoherent model for A+ = 100 and B- = 0.1
Time plot of differentiation target gene concentration with noise using the incoherent model for A+ = 100 and B- = 0.1	Time plot of self-renewal target gene concentration with noise using the incoherent model for A+ = 100 and B- = 0.1
Time plot of differentiation target gene concentration using the incoherent model for A+ = 10 and B- = 100	Time plot of self-renewal target gene concentration using the incoherent model for A+ = 10 and B- = 100
Time plot of differentiation target gene concentration with noise using the incoherent model for A+ = 10 and B- = 100	Time plot of self-renewal target gene concentration with noise using the incoherent model for A+ = 10 and B- = 100
Time plot of differentiation target gene concentration using the incoherent model for A+ = 100 and B- = 100	Time plot of self-renewal target gene concentration using the incoherent model for A+ = 100 and B- = 100
Time plot of differentiation target gene concentration with noise using the incoherent model for A+ = 100 and B- = 100	Time plot of self-renewal target gene concentration with noise using the incoherent model for A+ = 100 and B- = 100

Similar to Extension Figure 1, these show the time plots for differentiation and self-renewal target gene concentrations from the incoherent model for varying combinations of A⁺ and B^- over time. The combinations of signal concentrations correspond to the values described above. These figures use the transcription factor equations and the incoherent target gene equation, which changes parameter values to show the self-renewal target genes. The figures with noise use the transcription factor equations with noise and the normal OCT4-SOX2 and target gene equations. The parameter f₃ is set to 0 to ignore the effects of B₊ as described above. The parameters used come from Tables 24 and 25, with the changes and initial conditions described in the table below.

Table 22

State Variable	Initial Condition
[OCT4]	10
[SOX2]	10
[NANOG]	10
[OCT4-SOX2]	0
[Target Genes]	0
Parameters for All	Value
f3	0
Parameters for Self-Renewal	Value
h2	0.005

These figures correspond to the transcription factor figures above. Similarly to above, when comparing the time plots to the time plots with noise, it is clear that the time plots are a good approximation of the time plots including noise.

Extension Figure 3

Time plot of differentiation target gene concentration using the coherent model for A+ = 100 and B- = 0.1	Time plot of self-renewal target gene concentration using the coherent model for A+ = 100 and B- = 0.1
Time plot of differentiation target gene concentration with noise using the coherent model for A+ = 100 and B- = 0.1	Time plot of self-renewal target gene concentration with noise using the coherent model for A+ = 100 and B- = 0.1
Time plot of differentiation target gene concentration using the coherent model for A+ = 10 and B- = 100	Time plot of self-renewal target gene concentration using the coherent model for A+ = 10 and B- = 100
Time plot of differentiation target gene concentration with noise using the coherent model for A+ = 10 and B- = 100	Time plot of self-renewal target gene concentration with noise using the coherent model for A+ = 10 and B- = 100
Time plot of differentiation target gene concentration using the coherent model for A+ = 100 and B- = 100	Time plot of self-renewal target gene concentration using the coherent model for A+ = 100 and B- = 100
Time plot of differentiation target gene concentration with noise using the coherent model for A+ = 100 and B- = 100	Time plot of self-renewal target gene concentration with noise using the coherent model for A+ = 100 and B- = 100

Extension Figure 3 is created in the same way as Extension Figure 2 but uses the coherent model instead. These are another set of time plots for differentiation and self-renewal target gene concentrations for varying combinations of A+ and B- over time, with the same combinations of signal concentrations. The equations used are the transcription factor and coherent target gene equations, with the extension transcription factor equations used for the noisy plots. The parameter f3 is again set to 0, and all other parameters are found in Tables 24 and 26. The change and initial conditions are summarized in the table below.

Table 23

State Variable	Initial Condition
[OCT4]	10
[SOX2]	10
[NANOG]	10
[OCT4-SOX2]	0
[Differentiation Target Genes]	0
[Self-Renewal Target Genes]	0
Parameters for All	Value
f3	0

The time plots again match the time plots with noise, and these also are qualitatively similar to the figures in Extension Figure 2. Overall, we can conclude that the equations used to model this system are a good simplification of the more realistic equations involving stochastic noise.

Mathematica Notebook

All necessary code with comments can be found in the following Mathematica notebook: Stem Cell Differentiation and Self-Renewal Notebook

Discussion

The hypothesis explored by this paper was that signal concentration causes stem cells to differentiate or self-renew and that this behavior exhibits bistability. Additionally, by changing the binding and repressive strength of OCT4-SOX2 or the transcription rate of NANOG, the stem cell can be forced to undergo irreversible switches between differentiation and self-renewal.

Figures 3 and 4 demonstrate that a bistable switch is present for varying values of signals A⁺ and B^- in the concentrations of the transcription factors OCT4, SOX2, NANOG, and OCT4-SOX2, and Figures 7 and 8 demonstrate this in the concentrations of the differentiation and self-renewal target gene products in both the incoherent and coherent models. Further, Figures S1 through S5 demonstrate this bistability is persistent over varying parameter values for binding strengths of OCT4-SOX2 and OCT4-SOX2-NANOG, kinetic constants of OCT4-SOX2, and degradation rates . Additionally, Figures 10 and 11 show that by increasing the binding and repression strengths of OCT4-SOX2 on other transcription factors would cause the bistable switch to become irreversible, making it easier for one to control if the differentiation or self-renewal genes are expressed. Figure 12 show that this is also true for increasing the basal transcription rate of NANOG. These directly support the hypothesis.

Limitations

The model simplifies the various outside signals the system could receive into the positive A⁺ signal and the negative B^- signal, which may be an oversimplification of the signals, and signals represented by A⁺ and B^-, such as BMP4 and WNT, have been shown to effect each other, which was not addressed in this model . The model assumes that the signals and transcription factors do not directly interact, which could change results. Additionally, the model fails to determine if the system is an incoherent or coherent network and also relies on nondimensionalized variables, so it is difficult to quantitatively apply this to experiments .

Further, this model ignores stochasticity in transcription, which is known to occur . However, this was added and addressed in the extension, showing that this assumption does not lead to incorrect conclusions.

While this model included the three main transcription factors and the complexes they form, many recent papers have explored the effects that other transcription factors have on this system. Three of the most significant include Klf4, which activates OCT4 and SOX2, PRDM14, which also activates OCT4 and SOX2, and FoxD3, which activates NANOG . In addition, at least thirteen other transcription factors have shown to be part of this network, as well as miRNA .

Regardless, this model still demonstrates biologically significant results that match experimental data.

Discrepancies

In Figure 6B and S5B, the maximum transcription of the differentiation target genes is too high when there is either no autoregulation in Figure 6B or the binding strength between OCT4-SOX2 and the target genes is high. These discrepancies were addressed by modifying the parameter values to give the correct figure and noting the changes. For Figure 6B, the changes were not reasonable, which suggest that the figure found with the specified parameter values is correct. In Figure S1, bistability disappears for a comparatively low value of a₂, but the figures are qualitatively the same. Additionally, bistability is present for a large range of values of e₁ in the created figure while the figure from the paper does not appear to have bistability for e₁ . Through the paper and recreated figures, it has been shown that bistability is present for e₁ equal to 0.0001 through 0.01, and additional figures created in Mathematica that have not been shown demonstrate that the two-parameter bifurcation diagram created in Mathematica is correct .

Outside Literature

Other work in this field has also shown that the transcription factors create a network that creates a bistable switch between self-renewal and differentiation . The overexpression of OCT4, SOX2, and NANOG, as well as KLf4, a transcription factor not included in this model, have been shown to induce pluripotency in cells that had already differentiated . The network modeled in this paper was also confirmed to include a large portion of the transcription factor network, though recent papers have discovered the importance of additional proteins . One of these papers has built off of the model presented here, including additional transcription factors, which still demonstrates the bistability found in this model . Further, the findings from the time plots with stochasticity reflect the findings of the papers that led us to include stochasticity as our extension .

Additionally, while the part of the hypothesis suggesting ways to control the pluripotency of stem cells using outside methods was explored through the model, no further research was done on how the findings could be applied to actual cells. This is likely due to the difficulties implementing strategies such as changing binding efficiencies or basal transcription rates. However, one extension of the model showed that an irreversible switch was present for high levels of the decay rate of the additional transcription factor REST, similar to findings from this model . Further, experiments have shown that by forcing the expression of NANOG, expression of differentiation target genes is reduced, which is parallel to increasing the basal transcription rate of NANOG .

Future Work

To further understand this system, steady-state plots for the concentration of target gene products due to varying concentration of signal A⁺ for several values of e₁ should be created to confirm the discrepancy in Figure S1. Also a figure similar to Figure S4 for values of η₁ and η₅ between 0 and 0.01, which is expected to show the change between a bistable and irreversible switch.

Other future goals may include determining if the system is incoherent or coherent and extending the model to involve more of the associated transcription factors discovered by other researchers . To explore the network that the coherent and incoherent systems are built on by using the model, one can create multiple steady-state plots to explore the results of changing NANOG binding efficiency to the target genes, which would either cause the incoherent system to have decreased self-renewal target gene expression while or the expression to increase for the coherent system .

Parameter Tables

Table 24: Normal Transcription Factors

Parameter	Value
k1c	0.05
k2c	0.001
k3c	5
a1, a2, a3	1, 0.01, 0.2
b1, b2, b3	0.0011, 0.001, 0.007
c1, c2, c3	1, 0.01, 0.2
d1, d2, d3	0.0011, 0.001, 0.0007
e1, e2	0.005, 0.1
f1, f2, f3	0.001, 9.95*10-4, 0.01
η1	10-4
η2	10-7
η3	10-4
η4	10-7
η5	10-4
η6	10-7
γ1	1
γ2	1
γ3	1

Table 25: Incoherent Normal Target Genes

Parameter	Value
e1, e2	10-4, 10-3
f1, f2	9.01*10-4, 10-3
g1	0.1
h1, h2	0.0019, 0.05
η5	10-5
η6	10-7
η7	10-5
η8	10-7
γ3	0.05
γ4	10-2

Table 26: Coherent Normal Target Genes

Parameter	Value
m1, m2	0.1, 0.1
n1, n2	10-3, 10-2
q1, q2	0.001, 0.01
η9	10-4
η10	10-7
η11	1
η12	0.001
γ5	0.05
γ6	10-2