General Description

Considering the gene expression process, we can divide the process into two major sections, protein synthesis and protein degradation.

Assumption 1: Protein synthesis and degradation are deterministic processes in the system.

Therefore, we can derive following equation:

$\frac{d[P]}{dt} = v_{ProteinSyn} - \beta P$

where in this equation, $[P]$ represents the concentration of the protein, while the $v_{synthesis}$ represents the velocity of the synthesis and $\beta$ represents the degradation rate of the protein.

While focusing on the synthesis of the protein, it can be divided into transcription, translation and degradation of mRNA.

Assumption 2: Protein Synthesis is considered to be the process of transcription, translation and mRNA degradation.

Thus, we can get these expressions:

$v_{ProteinSyn} = k_{TL}[mRNA]$

$\frac{d[mRNA]}{dt} = v_{mRNASyn} - \delta[mRNA]$

If transcription factor (TF) is an activator that increases the production of the protein, $P$, then the rate of production of $P$ will be proportional to the fraction of the factors that are bound, $\theta$. If TF is a repressor, then the rate will be
proportional to $(1 − θ)$. $\delta$ is the degradation constant of mRNA. (Chandran D et al., 2008)

$v_{mRNASyn}$ can be represented by the transcriptional constant $k_{TC}$ and the fraction of TF that are bound:

$v_{mRNASyn} = k_{TC}\theta$

if the TF is the activator, while

$v_{mRNASyn} = k_{TC}(1-\theta)$

if the TF is the repressor.

As a whole, the kinetic equations of a protein is:

$\frac{d[P]}{dt} = k_{TL}[mRNA] - \beta [P]$

$\frac{d[mRNA]}{dt} = k_{TC}\mathcal{F^{+-}} - \delta[mRNA]$

$\mathcal{F^{+-}}=\left\{ \begin{aligned} & \theta \quad\quad\quad, \text{activator} \\ & (1-\theta) \quad, \text{repressor} \\ \end{aligned} \right.$

Hill-Langmuir Equation

Hill function is widely used in describing the interaction between promoter and its transcription factor. (Goldenfeld N et al., 2007)

A useful function that describes many real gene input functions is called the Hill function. The Hill function can be derived from considering the equilibrium binding of the transcription factor to its site on the promoter.

Likewise, the description of Hill equation in Wikipedia is:

In biochemistry and pharmacology, the Hill equation refers to two closely related equations that reflect the binding of ligands to macromolecules, as a function of the ligand concentration.

According to Hill-Langmuir equation, the fraction of the receptor protein concentration that is bound by the ligand $\theta$ is represented by the concentration of the ligand $[L]$, the dissociation constant $K_d$ and the Hill coefficient $n$.

$\theta = \frac{[L]^n}{K_A^n+[L]^n}$

In this equation, $K_A$ is the ligand concentration producing half occupation.

Kinetic Simulation of a protein

As illustrated before, the kinetic equations of a protein is:

$\frac{d[P]}{dt} = k_{TL}[\text{mRNA}] - \beta P$

$\frac{d[\text{mRNA}]}{dt} = v_{\text{mRNASyn}} - \delta[\text{mRNA}]$

If the TF is an activator, the $v_{\text{mRNASyn}}$ is represented as:

$v_{\text{mRNASyn}} = k_{TC}\frac{[\text{acti}]^n}{K_A^n+[\text{acti}]^n} = \frac{k_{TC}}{1 + (\frac{K_A}{[\text{acti}]})^n}$

While if the TF is a repressor, the $v_\text{mRNASyn}$ is represented as:

$v_\text{mRNASyn} = k_{TC}(1-\frac{[\text{repr}]^n}{K_A^n+[\text{repr}]^n}) = k_{TC}\frac{K_A^n}{K_A^2+[\text{repr}]^n} = \frac{k_{TC}}{1 + (\frac{[\text{repr}]}{K_A})^n}$

Therefore, we can get that, if the TF is an activator:

$\frac{d[P]}{dt} = k_{TL}[\text{mRNA}] - \beta P$

$\frac{d[\text{mRNA}]}{dt} = \frac{k_{TC}}{1 + (\frac{K_A}{[\text{acti}]})^n} - \delta[\text{mRNA}]$

If the TF us a repressor:

$\frac{d[P]}{dt} = k_{TL}[\text{mRNA}] - \beta P$

$\frac{d[\text{mRNA}]}{dt} = \frac{k_{TC}}{1 + \left(\frac{[\text{repr}]}{K_A}\right)^n} - \delta[\text{mRNA}]$

Methionine Consumption

Since we use the natural consumption as the time delay to postpone the expression of fusion protein cns1A2, we'd like to simulate the consumption of the methionine.

Assumption 3: The consumption of the methionine of a yeast cell per unit time is fixed, and the total consumption is dependent on the population of the yeast cells and the time.

Logistic Population Model

The culture of engineered yeast cells is limited owing to the restricted space and nutrition. One of the methods to describe the population in symmetric way is logistic population model.

According to the population dynamics, population increasing rate is proportional to the population size.

$\frac{dN}{dt} = r \cdot N$

where $N$ is the population size, and r is the intrinsic rate of increase.

In logistic population model, the carrying capacity of the system being studied is defined as $K$. Carrying capacity is the population level at which the birth and death rates of a species precisely match, resulting in a stable population over time.

$r(t) = r_0(1-\frac{N}{K})$

Therefore, the population size can be represented by the initial intrinsic rate of increase $r_0$, the carrying capacity $K$ and the population size itself.

$\frac{dN}{dt} = r_0(1-\frac{N}{K})\cdot N$

Gompertz Curve

In bacteria growth researches, it is found that a method called Gompertz Curve has a better performance in fitting with the experimental data (Zwietering M et al., 1990).

Gompertz curve is raised originally to describe a person's resistance to death decreases as his age increases. (Winsor et al., 1932)

Used in Biology to describe growth curve, it is written in:

$N(t) = ae^{-e^{b-ct}}$

where N is the population size which is dependent on the time, asymptote $a$, the x-axis skewing of the inflection point $b$ and the growth rate $c$.

Considering the differenciation of the Gompertz curve:

$\frac{dN}{dt} = cNe^{b-ct}$

Reparameterization of Growth Models

Since the parameters in the equation above (a, b, c) do not have biological meaning, it is difficult to estimate the values of the parameters. Therefore, all the growth models were rewritten to substitute the mathematical parameters with $A$, $\mu_m$, and $\lambda$. (Zwietering M et al., 1990)

Figure. 2 Biological meaning of parameter reparameterization of the growth model. $A$ is the asymptote, $\lambda$ is the x-axis skewing and $\mu_m$ is the max growth rate.

According to the article, the logistic equation can be reparameterized to:

$N = \frac{A}{\{1+exp[\frac{4\mu_m}{A}(\lambda-t) + 2]\}}$

and

$\frac{dN}{dt} = \frac{4\mu_m}{A}(1-\frac{N}{A}) N$

Likewise, Gompertz curve can be reparameterized in the following form:

$N = A\cdot exp\{-e^{\left[\frac{\mu_m\cdot e}{A}(\lambda-t)+1\right]}\}$

and

$\frac{dN}{dt} = \frac{\mu_m\cdot e}{A}N\cdot e^{\left(\frac{\mu_m\cdot e}{A}(\lambda-t)+1\right)}$

Methionine Comsumption Calculation

In the case of our project, the yeast strain we use is able to produce methionine intrinsically but it will not produce methionine until the externally added methionine is depleted. The consumption of methionine is dependent on the population of the yeast cell and time.

Assumption 4: The consumption of methionine will not affect the protein synthesis in the cell, since the yeast strain can intrinsically produce methionine, and the strength of the promoter will not be weakened because of it.

Since the cells consume methionine more quickly when it is abundant and slowly when it is nearly depleted， the change of concentration of methionine can be represented as:

$\frac{d[met]}{dt} = -k_{met}N\frac{[met]}{Met}$

In the equation above, $Met$ is the initial concentration of methionine, using which we eliminate the influence of the addition of $[met]$ to the equation. This equation shows that the cell consumes methionine more slowly as the gradual depletion, which is reasonable in the real world.

In our simulation for delay time system, we use logistic model to represent the cell growth and the equation above becomes:

$\frac{d[met]}{dt} = -k_\text{met}\cdot 10^{\frac{A}{1+e^{-b\cdot t+b\cdot e}}}\frac{[met]}{Met}$

In the equation above, $\frac{A}{1+e^{-b\cdot t+b\cdot e}}$ is the OD value of cell growth, and we simply evaluate the cell number based on this OD value. Therefore, $k_{met}$ is set optimally based on the result from British_Columbia 2012. $A$ is the asymptote, $e$ is the x-axis skewing of the inflection point and $b$ is the growth rate.

Simulation of cns1A2 Expression

Putting all the equations together, we can derive the following ordinary differential equation (ODE) system:

$\frac{d[met]}{dt} = -k_\text{met}\cdot 10^{\frac{A}{1+e^{-b\cdot t+b\cdot e}}}\frac{[met]}{Met}$

$\frac{d[mRNA_{GAL4}]}{dt} = \frac{k_{TC-GAL4}}{1 + (\frac{[met]}{K_{A-GAL4}})^n} - \delta_{GAL4}[mRNA_{GAL4}]$

$\frac{d[GAL4]}{dt} = k_{TL-GAL4}[mRNA_{GAL4}] - \beta_{GAL4} [GAL4]$

$\frac{d[mRNA_{cns1A2}]}{dt} = \frac{k_{TC-cns1A2}}{1 + (\frac{K_{A-cns1A2}}{[GAL4]})^n} - \delta_{cns1A2}[mRNA_{cns1A2}]$

$\frac{d[cns1A2]}{dt} = k_{TL-GAL4}[mRNA_{cns1A2}] - \beta_{cns1A2} [cns1A2]$

In the equations, since the ordinary differential equation of methionine concentration causes minus concentration after a certain period, we did an alteration, multiplying $\frac{[met]}{Met}$, to keep the concentration around zero when it is depleted. The equations above describes the genetic circuit involved in delay system, leaving alone the expression of cns3-HisG domain.

Parameter Determination

All the parameters used in delay system model are based on experimental data or derived from research articles, which make the quantitative model reliable to use.

Methionine Consumption

Total consumption rate of methionine is related to the amount of yeast cell in culture, and it is evaluated by logistic model. The parameters (A, b, e) are from the fitting of logistic model to the experimental data of S.cerevisiae YM4271 straingrowth. More details about YM4271 growth measurement.

We estimated the parameters using nonlinear regression of logistic model in R and get the result below.

Figure. 3 Data fitting of experimental data and logistic model and parameter estimates of the model. In the figure, b is the opposite number of growth rate, d is the asymptote and e is the x-axis skewing.

From the p-value of the parameter estimates, we could regard it a reliable simulation. Therefore, we chose the parameters of yeast growth based on this result.

mRNA Synthesis

The transcriptional activator GAL4 binds as a dimer, so the hill constant of it is $2$. Originally, there is an interaction between GAL4 and GAl80, which has more interaction with other proteins like GAL3. However, in our project, after consulting Prof. Liu for genetic engineering strategies in yeast, we change our yeast strain from AH109 to YM4271, which is GAL4 and GAL80 deficient (GAL4$\Delta$, GAL80$\Delta$). Therefore, the GAL4 protein acts as the pure activator in this case. (Traven A et al., 2006)

Hong et al. discovered that Gal4 binds DNA with an apparent dissociation constant ($K_d$) of 24 nM. Since $K_A$ is defined as $K_A^n = K_d$ and $n = 2$, $K_\text{A-GAL4}$ equals 4.9 nM. (Hong M et al., 2008)

In the equations, parameters like mRNA synthesis rate $k_{TC}$, mRNA translation rate $k_{TL}$, mRNA decay rate $\delta$ and protein decay rate $\beta$ are of vital importance in the simulation. According to Miller et al, mRNA decay and synthesis are functionally independent during normal growth.
(Miller C et al., 2011)

In dynamic transcriptome analysis of yeast, the average synthesis rate of transcription factor (TF) is around $25$ mRNA/cell/cellCycle. A cell cycle is defined as 150 min. Here, we specify the transcription rate of GAL4 controlled by pMET3 to be this average value. (Miller C et al., 2011)

Referring to Mitre et al. research, there is a conversion of unit between $mM$ and $\text{molecules/cell}$, it assumes the cell volume of yeast is $70 \mu m^3$: (Mitre T et al., 2016)

$1mM = \frac{10^{-3}\text{mol}}{1 L} = \frac{6.02214129 \times 10^{23}\times 10^{-3} \text{molecules}}{(1 \text{dm})^{3}} =$

$=\frac{6.02214129 \times 10^{20} \text{molecules}}{(10^{-5} \mu m)^{3}} \times \frac{70 \mu m^3}{\text{cell}} = 4.2154989\times 10^{7} \text{molecules/cell}$

Using this calculation in our project,

$k_\text{TC-GAL4} = \frac{25}{4.2154989\times 10^{7} \times 150/60} = 2.37\times 10^{-7}$

Likely, paper from Sun et al. reveals the mean mRNA synthesis rate, $72$ mRNA/cell/cellCycle. Therefore, we define the transcription rate of our uncharacterised protein cns3 to this value. (Sun M et al., 2012)

$k_\text{TC-cns3} = \frac{72}{4.2154989\times 10^{7} \times 150/60} = 6.83\times 10^{-7}$

In another paper, the transcription rate (synthesis rate) of GAL1 gene after pGAL1 promoter is 1.042 molecules/(cell × min) (Mitre T et al., 2016), which can be converted to M/h:

$k_{TC-cns1A2} = \frac{1.04\times 60}{4.2154989\times 10^{7}} = 1.48 \times 10^{-6}$

The unit of these constants are $M/h$.

The mechanism and kinetic of repression by methionine to MET3 promoter is is poorly known, resulting in no data for Hill coefficient of methionine and dissociation constant binding to its target.

From the only paper found to describe the $\beta$-galactosidase activity driven by MET3 promoter, we estimate the $K_{A-MET3}$, which is the ligand concentration producing half occupation, to be $40 nM$, which is standardized to $4\times 10^{-8}M$. (Mao X et al., (2002)

mRNA Degradation

Considering the decay rate of mRNAs, we can derive the average decay rate of the transcription factor which substitutes the specific decay rate for each mRNA, since little data are available for the decay rate of $\text{mRNA}_{\text{GAL4}}$.

The decay rate is derived from the half-life of mRNA.
(Pérez‐Ortín J et al., 2013)
$\text{HL} = \frac{\text{ln}2}{k_d}$

In which $\text{HL}$ is the mRNA half-life and $k_d$ represents the degradation constant.

In dynamic transcriptome analysis, the median half-lives of transcription factor mRNAs calculated in yeast cell is around 11 min (Miller C et al., 2011), from which we can get the median decay rate in unit of $h^{-1}$.

$k_d = \frac{\text{ln}2}{\text{HL}} = \frac{\text{ln}2}{\frac{11}{60}} = 3.78$

Since cns3 and cns1A2 protein are poorly characterised, we define the decay rate of the mRNAs of these two protein the median decay rate of yeast mRNAs, which are calculated above, $3.78 h^{-1}$.

However, the decay rate of GAL4 mRNA is described in the literature, which is around $43.32\times 10^{-3} \text{min}^{-1}$ which is derived from the half-life of 16 min (Mitre T et al., 2016).

From the data above, we can get the decay rate of the GAL4 mRNA.

$\delta_{GAL4} = 43.32 \times 10^{-3} \text{min}^{-1} = 2.60 \quad h^{-1}$

Protein Synthesis

Considering the translation rate of the mRNAs, Mitre et al. estimated the parameters based on the equation below:

$k_{TL} = \frac{\text{The fraction of translated mRNA} \times \text{elongation rate}}{protein lenght} \times \frac{\text{numbers of proteins}}{mRNA}$

In, Arava et al. (2003), "fraction of translated mRNA" is defined to be relative translation rate, which is unitless:

$\text{relative translation rate} = \text{ribosome occupancy}\times \text{ribosome density}$

Since there is no data available for the relative translation rate of GAL4, we use the median relative translation rates of all associated GAL mRNAs reported, which is given by 0.145. (Arava Y et al., 2003)

As the ORF length increased, the ribosome density dropped from an average of 1.2 ribosomes per 100 nts for ORFs shorter than 400 nts to an average of 0.14 ribosomes per 100 nts for ORFs longer than 3,600 nts.

The length of cns3 protein is around $800$ aa, while the length of fusion protein of cns1 and cns2, cns1A2, is around $1200$ aa. According to the relation illustrated in the literature, the ribosome density is $0.6$ for cns3 protein and $0.4$ for cns1A2.

With the data given that the average ribosome occupancy is $71\%$, we can calculate the relative translation rate of cns3 is $0.426$ while that of cns1A2 is $0.284$.

Additionally, from the article of Bonven et al. (1979), the mRNA elongation rate if around 9.3 a.a per second. (Bonven B et al., 1979)

Considering the paper of Ideker et al., ratios of proteins to mRNA lies between 4200 and 4800. (Ideker T et al., 2001)

Integrating the constant about translation above, we can calculate the translation rate of the proteins used in our project.

$k_\text{TL-GAL4} = \frac{0.145\times 9.3 \times\ 3600\quad \text{a.a/h}}{882 \quad a.a}\times 4200\frac{\text{proteins}}{\text{mRNA}} = 23117\frac{\text{proteins}}{ \text{mRNA}\times h}$

$k_\text{TL-cns1A2} = \frac{0.284\times 9.3\times\ 3600\quad \text{a.a/h}}{1192 \quad a.a}\times 4200\frac{\text{proteins}}{\text{mRNA}} = 33502\frac{\text{proteins}}{ \text{mRNA}\times h}$

$k_\text{TL-cns3} = \frac{0.426\times 9.3\times\ 3600\quad \text{a.a/h}}{871 \quad a.a}\times 4200\frac{\text{proteins}}{\text{mRNA}} = 68774\frac{\text{proteins}}{ \text{mRNA}\times h}$

Protein Degradation

Using an epitope-tagged strain collection, researchers find that the mean half-life of proteins in yeast is about $43 min$ (Belle A et al., 2006). From the equation in mRNA degradation, we can calculate the degradation rate, which is $0.97h^{-1}$.

Since we can no direct data for the degradation rate of GAL4, cns1A2 and cns3, we use this mean degradation rate.

MATLAB code

MATLAB code to solve the equation system:

function [f] = dXdT(t, x)

%constant
A = 1.900786;
e = 28.689490;
b = 0.270983;
kmet = 1.15*10**(-4);
kTC1 = 2.37*10**(-7);
kA1 = 4*10**(-8);
D1 = 2.60;
kTL1 = 23117;
B1 = 0.97;
kTC2 = 1.48*10**(-6);
kA2 = 4.9*10**(-9);
D2 = 3.78;
kTL2 = 33502;
B2 = 0.97;
kTC3 = 6.83*10**(-7);
D3 = 3.78;
kTL3 = 68774;
B3 = 0.97;
N = 2;
MET=0.0001; %Same with the initial concentration of methionine

%parameter
met=x(1);
mgal=x(2);
gal=x(3);
m1A2=x(4);
p1A2=x(5);
mc3=x(6);
c3=x(7);

%ODE
dmetdt = -kmet*(met/MET)*10**(A/(1+exp(-b*t+b*e)));
dmgaldt = kTC1/(1+(met/kA1)**N)-D1*mgal;
dgaldt = kTL1*mgal-B1*gal;
dm1A2dt = kTC2/(1+(kA2/(gal+10**(-10)))**N)-D2*m1A2;
dp1A2dt = kTL2*m1A2-B2*p1A2;
dmc3dt = kTC3-D3*mc3;
dc3dt = kTL3*mc3-B3*c3;
f = [dmetdt; dmgaldt; dgaldt; dm1A2dt; dp1A2dt; dmc3dt; dc3dt];

end

tspan = 0:0.001:60;
initial = [0.0001,0,0,0,0,0,0];

[t,x] = ode45( @dXdT, tspan, initial);

The result of the ODE solution is plotted below. We compared the expression of cns1A2 and HisG in our delay system and concluded that we successfully delayed the expression of cns1A2. Since the typical concentration of methionine used as repressing signal is 1mM, the delayed time under 1mM methionine is around 20 hours, which means that cns1A2 starts expression 20 hours after the addition of methionine:

In the modeling part of the delay system, we used the Hill equation to describe the binding between transcription factor and promoter. We took transcription, translation, RNA degradation, and protein degradation into consideration. We ignore the leaking expression of GAL4 protein because pMET3 is reported as a tight promoter. We used the logistic model as the population simulation of yeast cells and calculate the methionine consumption based on this model. All the parameters are chosen according to corresponding research articles. The model shows a range of delayed time from 1 hour to maximal 360 hours in theory.

The result of modeling in the delay system will give not only a profile of our project but also a guideline for the experiment. According to the details of our demand, we can alter the concentration of methionine and then adjust the delayed time of cordycepin production. In conclusion, we successfully realized flexible control of delayed time using our delay time system.