Team:uOttawa/Project
Simulating Stem Cells for Tomorrow's Treatments
Stem cells are a hot subject of medical research in recent years. Perhaps the most interesting property that these cells possess is their ability to differentiate from pluripotent cells into specialists.
We decided to engineer a network that mimics this phenomonon of differentiation.
In 2014, iGEMuOttawa undertook a similar project, where we implemented a network called a tri-stable switch. Based on Sui Huang's 2009 paper this network allows the expression of one of two genes, A and B, that represent two differentiated states. But there is also a third state, where both A and B are expressed. This third state represents the pluripotent cell.
An example of a tri-stable switch in nature. SA, SB and SC stand for state A, B and C, where state C is the undifferentiated AB state. Figure from Huang et al. 2007.
However, our mathematical models have revealed that the dynamics of the genetic network we built cannot yield the third, pluripotent state. Thus, with new data from our model, we have completely re-designed the network with new dynamics to hopefully achieve tristability.
Our new genetic network.
Although our network this year still did not achieve tristability, we have discovered new methods of gene regulation by modifying promoters with transcription factor binding sites. With time, it will eventually be possible to costruct a genetic network that implements the tri-stable switch. This research will shed light on the inner workings of natural stem cells, and perhaps allow synthetic stem cells to be engineered in next-generation therapies.
The Tri-stable Switch
Tri-stability is predicted from a bi-stable "toggle" switch with self-activation (Huang 2009). The self-activation enables there to exist a "sweet spot" where the expression of both genes is equal and stable in a random and stochastic system like a cell.
Simple schematic for a tri-stable switch. Notice the "toggle" switch part with mutual repression, and the self-activation.
Huang (2009) suggests that this kind of switch is used in the process of differentiation between tropho-ectoderm and inner cell mass. Cell differentiation is a very complex process, and it is likely that cells use many of these switches in a cascade to achieve a complex "hierarchy" of specialization.
An example of a tri-stable switch in nature. SA, SB and SC stand for state A, B and C, where state C is the undifferentiated AB state. Figure from Huang et al. 2007.
In the top figure, it is evident that genes A and B must be both activated and repressed at the same time. More interestingly, the products of A and B must function as both activators and repressors at the same time. In 2014, iGEMuOttawa designed a set of dual-input promoters that are capable of being both activated and repressed. However, since activating and repressing transcription factors work on the same promoter, the exhibit cooperativity and lead to multiplicative behaviour. That is, the amount of repression is dependent on the amount of activation, and vice-versa.
Our 2014 network. Notice how each promoter is being induced and repressed at the same time.
We discovered that such multiplicative behaviour will not produce tri-stability. Instead, we needed additive behaviour. In additivity, the amount of repression is not dependent on the amount of activation. Adding x amount of repressor will always decrease the final protein concentration by a fixed amount, regardless of the quantity of the activator.
We re-designed our network with this in mind. We split each gene into its activating and repressing components, and created a new set of constitutive, repressible promoters.
Our new genetic network.
Promoters
It is evident from the network diagram that we need four distinct promoters: two inducible ones, and two repressible ones. The catch is that the same transcription factor must act as both an inducer and repressor. Based on promoters created by Tom Ellis (2009) and our own iGEM team in 2014. By placing transcription factor binding sites about 10 bp away from the TATA box, we can repress the gene when the transcription factor hinders the binding of the polymerase at the TATA box.
Our genetic network.
The transcription factors we decided to use were GEV, derived from the gal4 binding protein, and rtTA, from the tet-on/tet-off systems. Thus, we needed to have two promoters that are induced by GEV and rtTA, and two that have basal expression but are repressed. We had the former, and we needed to build the latter.
We decided to re-engineer two constitutive promoters: pTDH3 and pPGK1. Both of these promoters come from genes coding glycolytic enzymes, and so have extremely high basal expression rates. To make them repressible, we had two approaches:
- We took the upstream part of each promoter containing the UAS and combined it with the downstream part of the pGal promoter, which we have successfully engineered in GEV-repressible and rtTA-repressible variants;
- We added two binding sites close to what appeared to be the TATA box on each of the two promoters. This proved rather difficult, as pTDH3 and pPGK1 are TATA-less, don't have a well-defined transcription start site, and are very AT-rich in the region where the TATA box should exist.
This left us with 8 new promoters to test.
Schematics of all new promoters that we attempted to construct. Gx indicates GEV-repressibility, and Tx indicates rtTA-repressibility. Both pTDH3 and pPGK1 are strong constitutive promoters.
We ordered each promoter as gBlocks and assembled them into a chassis suitable for testing.
Results
We tested each promoter by transforming them into a suitable yeast chassis and using the promoter to drive GFP. Unfortunately, only one promoter (pPGK1Gx) confirmed to have transformed correctly, so we do not have as much data as we would like.
Since pPGK1Gx is repressed by GEV, we used beta-estradiol to modulate GEV activity. Similarly, we would have used doxycycline to modulate rtTA activity. Beta-estradiol binds to an estrogen receptor domain (the 'E' in GEV) and reveals a nuclear localization tag, causing GEV to be transported into the nucleus, where it binds to genomic DNA. Higher concentrations of estradiol upregulates GEV, which would in this case increase repression.
Since GEV is a transcriptional activator, it does to some degree recruit RNA polymerase. We did not expect pPGK1Gx to exhibit sigmoidal behaviour. Rather, we expected there to be repression at low concentrations of estradiol, and activation with high concentrations, leading to a U-shaped curve.
Expression rate of GFP driven by our new pPGK1Gx promoter. The x-axis has concentration of beta-estradiol in nM, and the y-axis has fluorescence in arbitrary units. Notice the U-shaped curve, and how the expression never reaches basal levels, even with extremely high concentrations of estradiol. The two data series represent two biological replicates.
Fluorescence was measured using a flow cytometer.
See this promoter in the registry.
Modelling
Equations
Here we define the equations that describe the tri-stable switch genetic network that we designed. We simplify our model from four ordinary differential equations (ODEs) to two ODEs. We make the assumption that due to the considerably faster degradation rate of mRNA compared to protein[citation], it reaches a steady state very quickly, such that \[ \frac{d[m_{1}]}{dt}=0\]
In addition to this assumption, we also simplify the model by not taking into account the inducible nature of our system (ie. The binding event between inducer and the protein activator/repressor). This means that in our equations, GEV and rtTA represent the induced population.
$$\frac{d[m_{1}]}{dt}=a_1\frac{[GEV]^n}{K^n_{a_1}+[GEV]^n}+b_1\frac{K^n_{b_1}}{K^n_{b_1}+[rtTA]^n}-\gamma_{m1} [m_1] \label{eq1}$$
Equation (1) describes the rate of change in the concentration of the mRNA of GEV. It consists of two terms describing the transcription of the GEV gene from two different promoters, one of which is auto-activated (the first term), and another which is repressed by rtTA (the second term). The final term describes the degradation of the mRNA at a constant rate.
| Parameter | Description |
| $$ a_1 $$ | Maximal transcription rate of GEV mRNA from the self-activated promoter |
| $$a_2$$ | Maximal transcription rate of GEV mRNA from the repressed promoter |
| $$K_{a1}$$ | Dissociation constant of GEV/promoter complex (range of [0,1]) |
| $$K_{b1}$$ | Dissociation constant of rtTAV/promoter complex (range of [0,1]) |
| $$\gamma_{m1}$$ | GEV mRNA degradation rate |
| $$n$$ | Hill coefficient (assumed to be the same for every promoter/protein complex) |
$$\frac{d[GEV]}{dt}=k_1[m_{1}]-\gamma_{p1}[GEV]$$
Equation (2) describes the rate of change in the concentration of GEV protein. Translation occurs at a constant rate and degradation of GEV protein occurs at a constant rate as well.
| Parameter | Description |
| $$ k_1 $$ | Maximal translation rate of GEV |
| $$\gamma_{p1}$$ | GEV protein degradation rate |
$$\frac{d[m_{2}]}{dt}=a_2\frac{[rtTA]^n}{K^n_{a_2}+[rtTA]^n}+b_2\frac{K^n_{b_2}}{K^n_{b_2}+[GEV]^n}-\gamma_{m2} [m_2]$$
Equation (3) describes the rate of change in the concentration of the mRNA of rtTA. It consists of two terms describing the transcription of the rtTA gene from two different promoters, one of which is auto-activated (the first term), and another which is repressed by rtTA (the second term). The final term describes the degradation of the mRNA at a constant rate.
| Parameter | Description |
| $$ a_2 $$ | Maximal transcription rate of rtTA mRNA from the self-activated promoter |
| $$b_2$$ | Maximal transcription rate of rtTA mRNA from the repressed promoter |
| $$ K_{a2}$$ | Dissociation constant of rtTA/promoter complex (range of [0,1]) |
| $$K_{b2}$$ | Dissociation constant of GEV/promoter complex (range of [0,1]) |
| $$\gamma_{m2} $$ | rtTA mRNA degradation rate |
| $$n$$ | Hill coefficient (assumed to be the same for every promoter/protein complex) |
$$\frac{d[rtTa]}{dt}=k_1[m_{2}]-\gamma_{p2}[rtTA] $$
Equation (4) describes the rate of change in the concentration of rtTA protein. Translation occurs at a constant rate and degradation of rtTA protein occurs at a constant rate as well.
| Parameter | Description |
| $$ k_2 $$ | Maximal translation rate of rtTA |
| $$\gamma_{p2}$$ | rtTA protein degradation rate |
We reduce the model by classifying [m1] and [m2] as fast variables which means the fast nature of mRNA dynamics lead to a quasi-steady-state in mRNA concentration. Mathematically, this assumption is represented by setting the rates of change of both mRNA variables to zero. We will only show how equations 1 and 2 can be reduced to just one equation since the same process is used for equations 3 and 4.
We start by setting [m1] to zero: \[ \frac{d[m_{1}]}{dt}=0.\]
This means that the following expression (Equation (5)) is now equal to zero:
$$0=a_1\frac{[GEV]^n}{K^n_{a_1}+[GEV]^n}+b_1\frac{K^n_{b_1}}{K^n_{b_1}+[rtTA]^n}-\gamma_{m1} [m_1]$$
To isolate [m1], we rearrange equation 5 to obtain equation 6,
$$[m_{1}]=\frac{a_1\frac{[GEV]^n}{K^n_{a_1}+[GEV]^n}+b_1\frac{K^n_{b_1}}{K^n_{b_1}+[rtTA]^n}}{\gamma_{m1}}$$
Equation 6 can now be substituted into equation 2 for the [m1]. We obtain equation 7:
$$\frac{d[GEV]}{dt}=k_1\frac{a_1\frac{[GEV]^n}{K^n_{a_1}+[GEV]^n}+b_1\frac{K^n_{b_1}}{K^n_{b_1}+[rtTA]^n}}{\gamma_{m1}}-\gamma_{p1}[GEV] \label{eq7}$$
The parameters in equation 7 can be grouped up to yield equation 8:
$$\frac{d[GEV]}{dt}=\frac{k_1a_1}{\gamma_{m1}}\frac{[GEV]^n}{K^n_{a_1}+[GEV]^n}+\frac{k_1b_1}{\gamma_{m1}}\frac{K^n_{b_1}}{K^n_{b_1}+[rtTA]^n}-\gamma_{p1}[GEV]$$
We can combine the parameters into 2 new parameters to further simplify the equation. By combining the following parameters we get equation 9:
$$a_1=\frac{k_1a_1}{\gamma_{m1}} \\ b_1=\frac{k_1b_1}{\gamma_{m1}}$$
$$\frac{d[GEV]}{dt}=a_1\frac{[GEV]^n}{K^n_{a_1}+[GEV]^n}+b_1\frac{K^n_{b_1}}{K^n_{b_1}+[rtTA]^n}-\gamma_1 G $$
By performing the same reduction to equations 3 and 4 we get the rate equation 10 for [rtTA]
$$\frac{d[rtTA]}{dt}=a_2\frac{[rtTA]^n}{K^n_{a_2}+[rtTA]^n}+b_2\frac{K^n_{b_2}}{K^n_{b_2}+[GEV]^n}-\gamma_2 G \label{eq10}$$
Our reduced model is following two equations:
$$\frac{d[GEV]}{dt}=a_1\frac{[GEV]^n}{K^n_{a_1}+[GEV]^n}+b_1\frac{K^n_{b_1}}{K^n_{b_1}+[rtTA]^n}-\gamma_1 G$$ $$\frac{d[rtTA]}{dt}=a_2\frac{[rtTA]^n}{K^n_{a_2}+[rtTA]^n}+b_2\frac{K^n_{b_2}}{K^n_{b_2}+[GEV]^n}-\gamma_2 R $$
Where,
$$\begin{align*} a_1=&\frac{k_1a_1}{\gamma_{m1}} & b_1=&\frac{k_1b_1}{\gamma_{m1}} \\ a_2=&\frac{k_2a_2}{\gamma_{m2}} & b_2=&\frac{k_2b_2}{\gamma_{m2}} \\ \end{align*}$$
| Parameter | Description | Sample Value |
| $$ a_1 $$ | Maximal transcription rate of GEV mRNA from the self-activated promoter | 1 |
| $$a_2$$ | Maximal transcription rate of rtTA mRNA from the repressed promoter | 1 |
| $$ b_1 $$ | Maximal transcription rate of GEV mRNA from the self-activated promoter | 1 |
| $$b_2$$ | Maximal transcription rate of rtTA mRNA from the repressed promoter | 1 |
| $$K_{a1}$$ | Dissociation constant of GEV/promoter complex (GEV gene) (range of [0,1]) | 0.5 |
| $$K_{a2}$$ | Dissociation constant of rtTA/promoter complex (rtTA gene) (range of [0,1]) | 0.5 |
| $$K_{b1}$$ | Dissociation constant of rtTA/promoter complex (GEV gene) (range of [0,1]) | 0.5 |
| $$K_{b2}$$ | Dissociation constant of GEV/promoter complex (rtTA gene) (range of [0,1]) | 0.5 |
| $$\gamma_{m1}$$ | GEV degradation rate | 1 |
| $$\gamma_{m1}$$ | rtTA degradation rate | 1 |
| $$n$$ | Hill coefficient (assumed to be the same for every promoter/protein complex) | 4 |
Equilibria and Nullclines
A fixed point or equilibrium point is a set of initial conditions which are equal to the output of a function. If we look at a hypothetical function f(x), c is a fixed point if, \[f(c)=c.\]
Nullclines help us determine the equilibria of a system. A simple system of equations like
$$\begin{align*} \frac{dx}{dt}=&f(x,y)\\ \frac{dy}{dt}=&g(x,y) \end{align*}$$
will have two nullclines: an x-nullcline and a y-nullcline. The x-nullcline is a solution curve that satisfies \[\frac{dx}{dt}=0.\]
Therefore the x-nullcline would be the curve that satisfies \[0=f(x,y)\]
An equilibrium point of our system would be a set of initial concentrations of GEV and rtTA that would not change upon iteration of the system. Our system is of the form
$$\begin{align*} \frac{d[GEV]}{dt}=&f([GEV],[rtTA])\\ \frac{d[rtTA]}{dt}=&g([GEV],[rtTA]). \end{align*}$$
This means that the GEV nullcline and rtTA nullcline will satisfy
$$\begin{align*} \frac{d[GEV]}{dt}=&0\\ \frac{d[rtTA]}{dt}=&0 \end{align*}$$
or $$0=a_1\frac{[GEV]^n}{K^n_{a_1}+[GEV]^n}+b_1\frac{K^n_{b_1}}{K^n_{b_1}+[rtTA]^n}-\gamma_1 G \\ 0=a_2\frac{[rtTA]^n}{K^n_{a_2}+[rtTA]^n}+b_2\frac{K^n_{b_2}}{K^n_{b_2}+[GEV]^n}-\gamma_2 R $$
To determine the equilibrium points of the system we must determine what concentrations of GEV and rtTA will not change over time. This occurs when the two nullclines intersect since at that point both rate equations are equal to zero. Using the software XPP Auto, we are able to plot the nullclines and find the equilibrium points.
The nullclines figure is a phase plane that has both of the nullclines plotted on it. We then determine that the equilibrium points are the following five points: (1,1),(0.0039216,1.9961),(1.9961, 0.0039216),(1.5116, 0.48842),(0.48842, 1.5116). From figure Flow we can interpret the behaviour of the system using the vector field. We see that three of the equilibrium points attract solutions: (1,1),(0.0039216,1.9961) and (1.9961, 0.0039216). These points represent the three stable states of our genetic switch. Meanwhile separating these attractors are two unstable equilibrium points: (1.5116, 0.48842),(0.48842, 1.5116). These points represent the threshold that must be trespassed to initiate switching between different states.
Nullclines
Phase plane of the model with both the [GEV]-nullcline (red) and [rtTA]-nullcline (green)
Fixed Points
Phase plane of the model with both the [GEV]-nullcline (red) and [rtTA]-nullcline (green). Fixed points are labelled with ovals for stable fixed points and diamonds for unstable fixed points.
Flow
Phase plane of the model with both the [GEV]-nullcline (red) and [rtTA]-nullcline (green). Fixed points are labelled with ovals for stable fixed points and diamonds for unstable fixed points. A vector field is overlayed showing rates of change of [GEV] and [rtTA] with different initial concentrations
Equilibrium Point Stability
The stability of equilibrium points in nonlinear systems is determined using the jacobian matrix. The jacobian matrix is a matrix of partial derivatives of every equation in the system relative to every variable in this system.
$$\frac{dx}{dt}=f(x,y)\\ \frac{dy}{dt}=g(x,y)$$
In a system of this form, the jacobian matrix would looks like: \[ J = \left| \begin{array}{ccc} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} \\ \frac{\partial g}{\partial x} & \frac{\partial g}{\partial y} \end{array} \right|.\]
If we evaluate the jacobian matrix at an equilibrium point and determine the eigenvalues of the system, we can make conclusions on the stability of that particular equilibrium point. Using the eigenvalue properties in the following table we can say whether an equilibrium point is stable or unstable. It is important to note that for an equilibrium point to be stable, both eigenvalues must be negative.
| Eigenvalue type | Stability |
| Positive real | Unstable node |
| Negative Real | Stable node |
| Complex with positive real | Unstable spiral |
| Complex with negative real | Stable spiral |
| Real, opposite signs | Saddle node |
Another way of visualizing the different types of equilibrium points is through a chart:
Chart of several different kinds of equilibrium points, and criteria for what causes them.
Using XPP Auto we are able to determine the eigenvalues of the jacobian matrix for each equilibrium point. The results in the table correlate with the vector field phase plane. We see the three stable steady states which represent the three steady states of our genetic network. We also see the two unstable equilibrium points. Therefore it is possible, with these sample parameter values, for a tri-stable genetic network with our design to work. Further testing with parametrization of our assembled network is necessary for us to mathematically validate the assembled system.
| Equilibrium point | Eigenvalues | Stability |
| (1,1) | $$\lambda_1=-1\\\lambda_2=-0.558095$$ | Stable Node |
| (0.0039216,1.9961) | $$\lambda_1=-1\\\lambda_2=-0.992188$$ | Stable node |
| (1.9961, 0.0039216) | $$\lambda_1=-1\\\lambda_2=-0.992188$$ | Stable Node |
| (1.5116, 0.48842) | $$\lambda_1=-1\\\lambda_2=1.072939$$ | Saddle Node |
| (0.48842, 1.5116) | $$\lambda_1=-1\\\lambda_2=1.072939$$ | Saddle Node |
References
Ellis, T., Wang, X., & Collins, J. (2009). Diversity-based, model-guided construction of synthetic gene networks with predicted functions. Nature Biotechnology, 27(5): 465-471.
Wang, Y., Huang, C., Tung, S., & Lin, Y. (2000). Competition with TATA Box-Binding Protein for Binding to the TATA Box Implicated in Human Cytomegalovirus IE2-Mediated Transcriptional Repression of Cellular Promoters. DNA and Cell Biology, 19(10): 613-619.
Brachman, C., Davies, A., Cost, G., Caputo, E., Li, J., Hieter, P., & Boeke, J. (1998). Designer Deletion Strains derived from Saccharomyces cerevisiae S288C: A Useful set of Strains and Plasmids for PCR-mediated Gene Disruption and Other Applications. Yeast, 14: 115-132.
Balázsi, G., Van Oudenaarden, A., & Collins, J. (2011). Cellular Decision Making and Biological Noise: From Microbes to Mammals. Cel,l 144(6): 910-925.
Huang, S. (2009). Reprogramming cell fates: Reconciling rarity with robustness. Bioessays, 31(5): 546-560.
Way, J., Collins, J., Keasling, J., & Silver, P. (2014). Integrating Biological Redesign: Where Synthetic Biology Came From and Where It Needs to Go. Cell, 157(1): 151-161.
Dueber, J., Mirsky, E., & Lim, W. (2007). Engineering synthetic signaling proteins with ultrasensitive input/output control. Nature Biotechnology, 25(6): 660-662.
Huang S, Guo YP, May G, Enver T. “Bifurcation dynamics in lineage-commitment in bipotent progenitor cells.” Developmental Biology (2007). 305:695-713.