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

---