Exemplary Results Draft 4
Results
This paper hypothesizes that changes in concentration of outside signals causes stem cells to switch from maintaining its state as a stem cell or differentiating, and that this behavior exhibits bistability. Further, it is proposed that a stem cell could be modified to become self-renewing in two ways: either by changing the binding efficiency between OCT4-SOX2-NANOG, a transcription factor complex, and OCT4, SOX2, and NANOG, the transcription factors, or increasing the transcription rate of NANOG, one of the transcription factors. 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
| |
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 1. 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
| |
Figure 4 is the same as Figure 3 but dependent on the concentration of signal B- instead, also using the parameters from Table 1. The concentration of A+ is set to 100 as specified by the paper.
For B- from about 0 to 34, the system remains "on". The paper claimed that this switch occurred at B- = 36, which is very close to the value that was found using Mathematica. Since this value is almost identical to the value found, no investigation was made past rechecking all parameter values, which were correct. 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 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 the original paper. Figure A shows NANOG concentration with and without NANOG autoregulation, and Figure B shows target gene concentration for the integrated model. Figure C shows differentiation and self-renewal target gene concentration for the incoherent model with NANOG autoregulation. |
|
Figure 6 includes three steady-state plots of the concentrations of NANOG or differentiation or 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 integrated 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 e2 and f2, 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 integrated model for this figure, with the differentiation target gene concentrations modeled with h2 = 0.05 and the self-renewal target gene concentrations modeled with h2 = 0.001 to achieve the different results. The value of f3 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 1 and 2 with the following modifications:
| Parameters for All | Value |
|---|---|
| f3 | 0 |
| No Autoregulation | Value |
| e2 | 0 |
| f2 | 0 |
| Self-Renewal Target Genes | Value |
| h2 | 0.001 |
NANOG autoregulation 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.
The curve for the concentration of the target genes without NANOG autoregulation in Figure 6B in Mathematica curves higher with a peak that requires a larger concentration of OCT4-SOX2 than the curve from the paper. All parameter values were checked, and relevant parameters, including the binding strength between target genes and OCT4-SOX2 and the repression strength between target genes and OCT4-SOX2 and OCT4-SOX2-NANOG, were varied in order to see if a mistake was made in parameter values. An explanation for this discrepancy was not found, but the two graphs are qualitatively the same, so a simple change in parameters is likely the cause.
Figure 7
| |
Figure 7 includes steady state plots of the concentration of the differentiation and self-renewal genes using the integrated model as a function of varying concentration of signal A+. These plots were made by setting the differential equations for the transcription factors and integrated 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 f3 is set to 0 to ignore the effects of B-, and the values of h1 and η7 were set to the values specified below to better shape the figure. The value of h2 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:
| 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 integrated model.
Figure 8
| |
Figure 8 includes steady state plots of the concentration of the differentiation and self-renewal genes based on the incoherent model instead of the integrated 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 f3 is set to 0 to ignore the effects of B-, as described above. The concentration of B- is held at 20.
| 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 integrated, or coherent, and incoherent models fulfill the same bistable switch for the same concentrations of A+. Therefore, these figures cannot be used to decide if the integrated or incoherent 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.
Figure 9
 Figure A is a bifurcation plot of the transcription factor concentrations as dependent on A+ concentration with weak NANOG feedback from the original paper, and Figure B is a bifurcation plot of the differentiation target gene concentration from the integrated model as dependent on A+ concentration with weak and strong NANOG feedback |
|
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 e2 and f2 to adjust the size of the graphs, and f3 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 incoherent model to 0. The values of a3 and c3 were set to 0.05, lowering the binding strengths between OCT4 and SOX2 to OCT4-SOX2-NANOG, and the values of b3 and d3 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 a3 and c3 were set to 1 and b3 and d3 to 1.5*10-3, increasing the strength of NANOG feedback.
| 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
| |
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 a3 and c3 are set to 0.5, increasing the binding strength of OCT4-SOX2-NANOG and OCT4 and SOX2, and the values of b3 and d3 are set to 0.001, decreasing the repression strength of OCT4-SOX2-NANOG on OCT4 and SOX2. The value of e1 is set to 0.01, increasing the binding strength between NANOG and OCT4-SOX2. The value of f2 is set to 0.001, increasing the repression strength between NANOG and OCT4-SOX2-NANOG slightly, and the value of f3 is set to 0.05, increasing the strength of repression of NANOG due to the signal B-.
| 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.
Figure 11
| |
Figure 11 is similar to Figure 4 but with the modifications explained above for Figure 10. The concentration of A+ is held at 100.
| 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.
Figure 12
| |
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 η5 was set to 35, so the basal transcription rate of NANOG was increased, and the value of η6 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 f3 was set to 0 as described above.
| 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
Figure S1 is a two parameter bifurcation plot that plots the values of a3, a2, e2, and e1 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 a3 corresponds to the binding strength between OCT4 and OCT4-SOX2-NANOG, and a2 corresponds to the binding strength between OCT4 and OCT4-SOX2. The value of e2 corresponds to the binding strength between NANOG and OCT4-SOX2-NANOG, and e1 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 a3, a bistable switch is present for a3 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 a2, a bistable switch was found to be present for a2 up to about 8.6, though the paper shows that this bistable switch is present for all values of a2 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 e2, a bistable switch is present for values of e2 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 e1, 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 e1, but by exploring the bifurcation plots of OCT4-SOX2 for A+ with varying values of e1, it is clear that there should be two values of A+ where switches occur for e1 between 0.1 and 1.8. The effects of varying e1 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
Figure S2 is a two parameter bifurcation plot that plots the values of k1c, k2c, and k3c 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 k1c, k2c, and k3c 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 k1c 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 k2c and k3c, there are bistable switches for all values of A+ shown by this figure, but the curves imply that beyond particular values of k2c and k3c 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 k3c 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 k3c, making the switch appear to occur at higher values of A+.
Figure S3
Figure S3 is a two parameter bifurcation plot that plots the values of γ1 and γ3 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 γ1 and γ3 correspond to the degradation rates of OCT4 and NANOG, and γ1 can be considered as equivalent to γ2, 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 γ1, 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 γ3, 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
Figure S4 is a two parameter bifurcation plot that plots the values of η1 and η5 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 η1 and η5 correspond to the basal transcription rates of OCT4 and NANOG, and η1 can be considered as equivalent to η3, the basal transcription rate of SOX2. The concentration of B- was set to 0.1.
The values of η1 and η5 are typically 10-4, so these bifurcation plots show values of η1 and η5 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 the original paper |
|
Figure S5 is an extension of Figure 6B, and is thus created in the same manner, but with varying values of e2, g1, and h2. The value of e2 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 g1 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 and 1 in the figure. The value of h2 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 f3 set to 0 as described above and with the concentration of B- set to 0.1.
| Parameters | Value |
|---|---|
| f3 | 0 |
In Figure S5A, the larger value of e2 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 e2 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 both values of g1, but the larger value of g1 causes 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 the larger value of g1 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 g1, which is specified by the paper. It is believed that a different value of g1 was actually used instead of 1.
In Figure S5C, the larger value of h2 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.
Mathematica Notebook
All necessary code with comments can be found in the following Mathematica notebook: Stem Cell Differentiation and Self-Renewal Notebook
Parameter Tables
Table 1: 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 2: Integrated 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 3: Normal Differentiation and Self-Renewal 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 |
