Delay System Simulation

GeneticCircuit iGEM DelaySystem

---

---

Introduction

One of the most novel designs in our project is the delay system of the cordycepin expression. Cordycepin, which is also called 3'-deoxyadenosine will be easily converted to functionless substance 3'-deoxyinosine (3'-dI) by adenosine deaminase (ADA) in cell when its concentration is high. Pentostatin in our project, acts as the protector against ADA. In this system, pentostatin needs to be expressed before the expression of cordycepin to protect it. Based on the background, we designed a delay system using a pMET3 promoter, which is suppressed by the presence of methionine. More details about delay system.

As the natural consumption of methionine, the expression of cns1-cns2 fusion protein will turn on when methionine is depleted. Therefore, we realize delay production of cordycepin. Additionally, we can adjust the delayed time by changing the initial concentration of added methionine.

Prerequisite Knowledge

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 expression:

$$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}]$$

Analysis of delay system

The interaction of inducer, transcription factors and promoters are illustrated below:

Figure. Illustration of the details in the expression of cns1A2. The repressor of pMET3 promoter is shown in dimer but in reality, the internal mechanism of transcription regulation of pMET3 is still unknown. On the contrary, GAL4 activator is reported to work as a dimer. In the figure, when methionine is depleted, the expression of cns1A2 will switch on.

One of the most important parts of the delay system simulation is the methionine consumption. The consumption of methionine controls the switching time of the system, so the model used in this session needs to be carefully evaluated.

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, the rate of increase of the population 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 th 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. 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)}$$

Methionince 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 quikly 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.

Genetic Circuit Simulation of Delay System

In our project, cns3 is constitutively expressed so it will not be affected by the concentration of the transcription factor.

In S. cerevisiae, the repression of MET3 by methionine is unusually stringent, with absence of leaky expression. [Care R S et al. (1999)]. Thus, in this modelling project, we do not consider the leaky expression in the equations.

The expression of the whole system is represented by the following equations:

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

$$\frac{d[\text{mRNA}_\text{GAL4}]}{dt} = \frac{k_\text{TC-GAL4}}{1 + \left(\frac{[\text{met}]}{K_\text{A-MET3}}\right)^n} - \delta_\text{GAL4}[\text{mRNA}_\text{GAL4}]$$

$$\frac{d[\text{GAL4}]}{dt} = k_\text{TL-GAL4}[\text{mRNA}_\text{GAL4}] - \beta_\text{GAL4} [\text{GAL4}]$$

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

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

$$\frac{d[\text{mRNA}_\text{cns3}]}{dt} = k_\text{TC-cns3} - \delta_\text{cns3}[\text{mRNA}_\text{cns3}]$$

$$\frac{d[\text{cns3}]}{dt} = k_\text{TL-cns3}[\text{mRNA}_\text{cns3}] - \beta_\text{cns3} [\text{cns3}]$$

Parameter Determination

Methionine Consumption

Total onsumption 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 strain growth. More details about YM4271 growth measurement.

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

Figure. 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. reseach, there is a conversion of unit between $mM$ and $\text{molecules/cell}$, it assumes the cell vlolume 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 avaible 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 eqaution 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);

Result

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:

Figure. Illustration of expression level of cns1A2 and HisG protein under 1.0mM initial methionine when A = 1.900786, e = 28.689490, b = 0.270983, kmet = 1.15E-4, kTC1 = 2.37E-7, kA1 = 4E-8, D1 = 2.60, kTL1 = 23117, B1 = 0.97, kTC2 = 1.48E-6, kA2 = 4.9E-9, D2 = 3.78, kTL2 = 33502, B2 = 0.97, kTC3 = 6.83E-7, D3 = 3.78, kTL3 = 68774, B3 = 0.97, N = 2.

The figure above demonstrates our design of the delay system as methionine successfully delays the expression time of cns1A2.

Our goal of designing delay system is flexible manipulation of expression time of our target protein. We would like to investigate the connection of initial methionine concentration and its corresponding delayed time. We tested the methionine concentration between 0.1mM and 0.5mM and plotted the result in the following figure.

Figure. Illustration of concentration of cns1A2 under various initial methionine when A = 1.900786, e = 28.689490, b = 0.270983, kmet = 1.15E-4, kTC1 = 2.37E-7, kA1 = 4E-8, D1 = 2.60, kTL1 = 23117, B1 = 0.97, kTC2 = 1.48E-6, kA2 = 4.9E-9, D2 = 3.78, kTL2 = 33502, B2 = 0.97, kTC3 = 6.83E-7, D3 = 3.78, kTL3 = 68774, B3 = 0.97, N = 2.

The delayed time under 0.1mM methionine in our model is about 1 or 2 hours, comparing to around 70 hours under 50mM methionine. We noticed that the upper bound of cns1A2 protein under various methionine concentrations is constant, suggesting that methionine only affects the delayed time.

This result testifies the flexibility of the adjustment on the delayed time by proper determination of methionine concentration. However, 70 hours is not the limit of the system. The solubility of methionine is about 379mM. Therefore, we further investigated the delayed time when increasing methionine concentration.

Figure. Illustration of concentration of cns1A2 under initial methionine concentration of 100mM, 200mM, 300mM and 379mM when A = 1.900786, e = 28.689490, b = 0.270983, kmet = 1.15E-4, kTC1 = 2.37E-7, kA1 = 4E-8, D1 = 2.60, kTL1 = 23117, B1 = 0.97, kTC2 = 1.48E-6, kA2 = 4.9E-9, D2 = 3.78, kTL2 = 33502, B2 = 0.97, kTC3 = 6.83E-7, D3 = 3.78, kTL3 = 68774, B3 = 0.97, N = 2.

Initial concentration from 100mM to 379mM causes an extension of delayed time from around 100 hours to 360 hours. So that we can theoretically delay the expression of cns1A2 to maximal 360 hours, which is 15 days. Definitely yeast cells are not likely to live more than 15 days in normal cases, so in theory, our design of this delay system can satisfy any requirement of the delayed time.

Conclusion

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 yeas 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.

Appendix

Symbol list

Symbol	Description
[ ]	The concentration symbol
$N$	population size of the cells in the culture
$\text{met}$	methionine
$Met$	initial concentration of methionine
$\text{mRNA}_{\text{GAL4}}$	mRNA transcribed by gal4 gene
$\text{GAL4}$	GAL4 protein, which is translated by $mRNA_{GAL4}$
$\text{mRNA}_{\text{cns1A2}}$	mRNA transcribed by cns1-cns2 fusion protein
$\text{cns1A2}$	fusion protein of cns1-cns2
$\text{mRNA}_{\text{cns3}}$	keto form of the product of cns2, at location of enzyme cns1
$\text{cns3}$	cns3 protein, which is tranlated by $\text{mRNA}_{\text{cns3}}$

Constant list

Constant	Description	Value	Reference
$A$	Carrying capacity (in OD value) of the system	1.900786	nonlinear regression
$b$	rate of populational increase	0.270983	nonlinear regression
$e$	x-axis offset of the inflection point	28.689490	nonlinear regression
$n$	Hill constant	2	Traven et al. (2006)
$k_{\text{met}}$	Consumption rate of methionine	$1.15\times 10^{-4} M/10^\text{OD}$	British_Columbia 2012
$k_{\text{TC-GAL4}}$	Transcription strength of binded promoter of gal4 gene	$2.37\times 10^{-7} M/h$	Miller C etal. (2014)
$k_{\text{A-MET}}$	ligand concentration producing half occupation of pMET3 promoter	$4\times 10^{-5}M$	Mao et al. (2002)
$\delta_{\text{GAL4}}$	Degradation rate of $\text{mRNA}_{\text{GAL4}}$	$2.60\quad {h}^{-1}$	Mitre et al. (2016)
$k_\text{TL-GAL4}$	Translation rate of $\text{mRNA}_{\text{GAL4}}$	$23117 \quad h^{-1}$	Arava et al. (2003)
$\beta_\text{GAL4}$	Degradation rate of GAL4 protein	$0.97h^{-1}$	Belle et al. (2006)
$k_{\text{TC-cns1A2}}$	Transcription strength of binded promoter of cns1A2 gene	$1.48 \times 10^{-6} M/h$	Mitre et al. (2016)
$k_{\text{A-GAL4}}$	ligand (GAL4) concentration producing half occupation of pGAL1 promoter	$4.9\times 10^{-9} \text{M}$	Hong et al. (2008)
$\delta_{\text{cns1A2}}$	Degradation rate of $\text{mRNA}_{\text{cns1A2}}$	$3.78\quad h^{-1}$	Miller C et al. (2004)
$k_\text{TL-cns1A2}$	Translation rate of $\text{mRNA}_{\text{cns1A2}}$	$33502 \quad h^{-1}$	Arava et al. (2003)
$\beta_\text{cns1A2}$	Degradation rate of cns1A2 protein	$0.97\quad h^{-1}$	Belle et al. (2006)
$k_{\text{TC-cns3}}$	Transcription strength of binded promoter of cns3 gene	$6.83\times 10^{-7} M/h$	Miller C etal. (2014)
$\delta_{\text{cns3}}$	Degradation rate of $\text{mRNA}_{\text{cns3}}$	$3.78\quad h^{-1}$	Miller C et al. (2004)
$k_\text{TL-cns3}$	Translation rate of $\text{mRNA}_{\text{cns3}}$	$68774 \quad h^{-1}$	Arava et al. (2003)
$\beta_\text{cns3}$	Degradation rate of cns3 protein	$0.97h^{-1}$	Belle et al. (2006)

Reference

[Chandran D et al. (2008)]
Chandran, D., Copeland, W. B., Sleight, S. C., & Sauro, H. M. (2008). Mathematical modeling and synthetic biology. Drug Discovery Today: Disease Models, 5(4), 299-309.

[Goldenfeld N et al. (2007)]
Goldenfeld, N. (2007). An introduction to systems biology: design principles of biological circuits; systems biology: properties of reconstructed networks. Physics Today, 60(6), 63-64.

[Zwietering M et al. (1990)]
Zwietering, M. H. . (1990). Modeling of the bacterial growth curve. Appl. Environ. Microbiol., 56(6), 1875-1881.

[Winsor et al. (1932)]
Winsor, & C., P. . (1932). The gompertz curve as a growth curve. Proceedings of the National Academy of Sciences, 18(1), 1-8.

[Care R S et al. (1999)]
Care, R. S. , Trevethick, J. , Binley, K. M. , & Sudbery, P. E. . (1999). The met3 promoter: a new tool for candida albicans molecular genetics. Molecular Microbiology, 34.

[Traven A et al. (2006)]
Traven, A. , Jelicic, B. , & Sopta, M. . (2006). Yeast gal4: a transcriptional paradigm revisited. Embo Journal, 7(5), 496-499.

[Hong M et al. (2008)]
Hong, M. , Fitzgerald, M. X. , Harper, S. , Luo, C. , Speicher, D. W. , & Marmorstein, R. . (2008). Structural basis for dimerization in dna recognition by gal4. Structure, 16(7), 1019-1026.

[Miller C et al. (2011)]
Miller, C., Schwalb, B., Maier, K., Schulz, D., Dümcke, S., Zacher, B., ... & Martin, D. E. (2011). Dynamic transcriptome analysis measures rates of mRNA synthesis and decay in yeast. Molecular systems biology, 7(1).

[Mitre T et al. (2016)]
Mitre, T. M., Mackey, M. C., & Khadra, A. (2016). Mathematical model of galactose regulation and metabolic consumption in yeast. Journal of theoretical biology, 407, 238-258.

[Sun M et al. (2012)]
Sun, M., Schwalb, B., Schulz, D., Pirkl, N., Etzold, S., Larivière, L., ... & Cramer, P. (2012). Comparative dynamic transcriptome analysis (cDTA) reveals mutual feedback between mRNA synthesis and degradation. Genome research, 22(7), 1350-1359.

[Mao X et al. (2002)]
Mao, X., Hu, Y., Liang, C., & Lu, C. (2002). MET3 promoter: a tightly regulated promoter and its application in construction of conditional lethal strain. Current microbiology, 45(1), 37-40.

[Pérez‐Ortín J et al. (2013)]
Pérez‐Ortín, J. E., Medina, D. A., Chávez, S., & Moreno, J. (2013). What do you mean by transcription rate? The conceptual difference between nascent transcription rate and mRNA synthesis rate is essential for the proper understanding of transcriptomic analyses. BioEssays, 35(12), 1056-1062.

[Arava Y et al. (2003)]
Arava, Y., Wang, Y., Storey, J. D., Liu, C. L., Brown, P. O., & Herschlag, D. (2003). Genome-wide analysis of mRNA translation profiles in Saccharomyces cerevisiae. Proceedings of the National Academy of Sciences, 100(7), 3889-3894.

[Bonven B et al. (1979)]
Bonven, B., & Gulløv, K. (1979). Peptide chain elongation rate and ribosomal activity in Saccharomyces cerevisiae as a function of the growth rate. Molecular and General Genetics MGG, 170(2), 225-230.

[Ideker T et al. (2001)]
Ideker, T., Thorsson, V., Ranish, J. A., Christmas, R., Buhler, J., Eng, J. K., ... & Hood, L. (2001). Integrated genomic and proteomic analyses of a systematically perturbed metabolic network. Science, 292(5518), 929-934.

[Belle A et al. (2006)]
Belle, A., Tanay, A., Bitincka, L., Shamir, R., & O’Shea, E. K. (2006). Quantification of protein half-lives in the budding yeast proteome. Proceedings of the National Academy of Sciences, 103(35), 13004-13009.