Modeling

Our team attach great importance to quantitative description of our biological systems and experimental phenomena. We applied modeling and quantitative methods in throughout our project, from macroscopic to microcosmic, from mechanics to applications.

To make the DNA replication control system precise and versatile, we focused on bacterial behaviors related to DNA replication, such as cell growth, etc. We built a single-cell stochastic model to describe that DNA replication would be delayed due to the expression of dCas9-sgRNA. Another model predicted the system could improve the metabolic efficiency, like to produce more high value bio-products by interfering DNA replication. The prediction was proved by our experiments. We also used models to simulate the behaviors of plasmid replication control system and quorum sensing-based DNA replication control system. These phenomena were validated by subsequent experiments. Last but not least, with the modelling, we have explored our system enabled to control the gene expression noise due to the variation of gene copy number. All the theoretical models helped us along the way of fully exploit potentials of "Dr. Control".

Regulatory Part

Introduction

In this section, we apply a deterministic model to predict the behavior of the regulatory part in our project.

Originally, we did not include the regulation of sgRNA. In this case, the sgRNA is constantly transcripted (By a J23119 promoter, see Project Description) However, we found some works reporting that dCas9 itself can influence the cell's growth rate, thus we have to explain that our system works in our exprected ways. Aditionally, we know that the promoters, even if regulated by an inducer or an inhibitor, may perfrom a "leakage expression", hence we have to apply some other ways to gain a better control of the system. Both of these motivated us to design a second switch for sgRNA.

Our regulatory part includes an IPTG-activated pLac-promoted sgRNA and an Arabinose-activated pBAD-promoted dCas9. The mechanics are described in Pic 1. The formation of the effective dCas9-sgRNA complex depends on both IPTG and Arabinose inducers.

Mono-regulation

We started from the simpler case: when there is only one regulated variable, e.g. dCas9 regulated by Arabinose. Since the expression of dCas9 is promoted by pBAD promoter and regulated by Arabinose, it is reasonable to describe this process by Hill's equation

in which $v$ is the producing rate of dCas9 mRNA, $n$ is the Hill's coefficient of pBAD promoter indicating the sharpness of induction rate and $K$ is the coefficient indicating the IPTG concentration at which $v = v_{\max }/2$ .

Bi-regulation

We introduced the regulation of expresion of sgRNA by changing the promoter of sgRNA from J23119 to T7. This T7 promoter's activation is regulated by a pTAC promoter. Therefore, we can further get another Hill's function for this promoter. Besides the transcription of sgRNA and dCas9 mRNA, the translation of dCas9, the binding of dCas9 and sgRNA and the degradation of all the mentioned matters are to be considered. Here, we describe this system by an ordinary differential equation system:

Where $m$ represents the concentration of dCas9 mRNA, $c$ represents the concentration of dCas9, $r$ represents the concentration of sgRNA and $b$ represents the concentration of dCas9-sgRNA concentration. The parameters are determined from both others' works and our experiments and listed as follows.

parameter	description	value
$n_m$	Hill's coefficient of pBAD promoter[1]	1.427
$K_m$	Half-activate Arabinose concentration of pBAD promoter[1]	0.5736%
$\alpha_m$	Maximum dCas9 mRNA expression rate[2]	0.0011 nM·min-1
$n_r$	Hill's coefficient of pTac promoter[1]	3.859
$K_r$	Half-activate Arabinose concentration of pTac promoter[1]	0.05137
$\alpha_r$	Maximum sgRNA expression rate[2]	0.0011 nM·min-1
$\beta_c$	dCas9 protein synthesis rate from dCas9 mRNA[2]	0.0057 protein·transcript-1·min-1
$\gamma_c$	dCas9 degradation rate[2]	5.6408×10-4min-1
$k_+$	Rate of dimerization of dCas9 and sgRNA[2]	1.0

By numerically solving the ODE system under different arabinose and IPTG concentration, we can show that the expression output follows the AND gate logic: the concentration of dCas9-sgRNA is only considerable when the concentration of both IPTG and Arabinose are high enough.

Figure 1: the mechanism of bi-regulation. Both the expression of sgRNA and dCas9 are regulated by the inducers, and they bind to form effective dCas9-sgRNA complex

Figure 2: Comparison of dCas9 expression at different dose of inducers. The concentration of dCas9-sgRNA is high only when both Arabinose and IPTG are sufficient

Steady State of the Bi-regulation Model

To study the steady state of the system more carefully, we solve the equation when all particles' concentrations do not vary with time

Since the effective molecule in our system is the dCas9-sgRNA complex, we are interested in the concentration of dCas9-sgRNA complex (e.g. the value of $b$ ).

from the system above, we can yield

where

This quadratic equation has two possitive solutions. However,only the smaller solution is reasonable (consider when $\gamma_m \rightarrow +\infty$ , the smaller solution is 0 and the bigger one is a positive real number $(\gamma_c \gamma_r k_- + \gamma_c \gamma_b \gamma_r + k_+ \gamma_b \alpha_r)/(k_+ \gamma_b^2)$ In a molecular biology context, when $\gamma_m$ is extremely big, mRNA degrades so frequently that there is almost no dCas9 mRNA, hence no dCas9, and subsequentially no dCas9-sgRNA complex, so the zero solution is the one corresponding to the real biological situation)

Thus, we get the final solution

where

This equation gives the precise expression of the relation between the concentration of dCas9-sgRNA and the amount of inducer we add. Applying the given parameters, we can demonstrate this relationship by a heat plot (Fig 3).

Figure 3: The steady-state concentration of dCas9-sgRNA at different inducer dose.

Here, the AND gate logic picture is more clear. We get a dCas9-sgRNA concentration over 10^-3nM only when the concentrations of both IPTG and sgRNA are high enough. It can also be seen that the Arabinose regulation has a wider dynamic range. This guided us to use more Arabinose regulation in our actual experiments.

Conclusion

This model theoretically explains the mechanisms of our system's regulation parts. Regulation parts in our system, and also in all synthetic systems, are essential both for experiments and for practical usage. By modelling the regulation, we gained an insight of the ways to improve the performance of the "switch" in our system.

References and Comments

1. Fitted from our CRISPRi experiment data
2. http://2014.igem.org/Team:Waterloo/Math_Book/CRISPRi
3. Bakshi S, Siryaporn A, Goulian M, Weisshaar JC. Superresolution imaging of ribosomes and RNA polymerase in live Escherichia coli cells. Mol Microbiol. 2012 Jul85(1): 21-38. doi: 10.1111/j.1365-2958.2012.08081.x. p.22 right column top paragraph

The Relationship between the Expression of dCas9 and the Cell's Growth Rate

Introduction

In this section we show that our system hacks the bacteria's replication rate. Specifically, the average time of a bacteria's cell cycle depends linearly on the concentration of the dCas9-sgRNA molecule copy number. This is deduced from a coarse-grained model including the process of dCas9 and DnaA binding to the OriC site and the cell's replication.

The bacteria's replication is initiated by an accummulation of DnaA on the genome's OriC. To simplify the cases but still remaining the essential mechanics, we regard the dCas9's replication hacking is because of competing binding to OriC with DnaA. Once OriC is occupied by dCas9, it cannot bind DnaA-ATP, and hence cannot replicate.

Figure 4: Three probable states of a bacterium - dCas9 unbound to OriC, dCas9 bound to OriC and replication initiated.

List of hypotheses

• The bingding process of both dCas9 and DnaA are reversible
• Genome DNA replication can start only when the OriC is bound to DnaA

Stochastic model of dCas9 binding and replication initiation

A Markov chain on continuous time is built, with the hypothesis that the replication time is subjected to a Poisson distribution. A cell has three possible states at a certain time

1. OriC is bound to a dCas9 molecule
2. OriC is naked
3. The bacteria start replication

A cell can transform from state 1 to state 2 ( $k_1$ ), from state 2 to state 1( $k_2$ ), and from state 2 to state 3 ( $\kappa$ ). A infinitesimal transition matrix is written to describe the transition among these states:

where

• $k_1$ represents the frequency that an exposed OriC is bound by a dCas9-sgRNA molecule (state 1 $\rightarrow$ state 2). This value is proportion to the concentration of dCas9-sgRNA (denoted as $c$ ):

• $k_2$ represents the frequency that a dCas9-sgRNA molecule decouple with OriC (state 2 $\rightarrow$ state 1).
• $\kappa$ represents the frequency that a exposed OriC is bound by replication initiators and the replication starts (state 1 $\rightarrow$ state 3)

suppose initially, the bacterium's OriC not bound with dCas9-sgRNA (i.e. the initial state distribution is $(1,0, 0)$ ). We calculate the distribution function $F$ of the waiting time before replication. To do this, we only have to calculate

thus the third element of ${\bf F}(t)$ will be the probability $F(t)$ we want.

It is not hard to get an explicit expression of $F(t)$ :

where

The average waiting time is:

linearly related to the concentration of dCas9-sgRNA.

Conclusion

This model provides the core mechanism of our system - that the dCas9 hacks the replication by binding to OriC, keeping the genome from replication initiation.

Plasmid Copy Number Hacking

Introduction

In this section we discuss how our system can control the plasmid copy number (see Project Description). Typically the plasmids replicate themselves during the cell cycle^[1]. After a cell's division, the plasmids are equivalently distributed to both of the children cells. Here we discuss a cell line's behavior over time.

Figure 5: A "Cell Line", analogous to a "path" in a tree structure"

Figure 6: An illustration of plasmid dilution. In a cell division, every plasmid in the mother cell have 1/2 probability to be atttributed to both of the children cells."

Figure 7: The plasmid replication and dilution process. The stability of plasmid copy number relies on a balance of plasmid replication and dilution.

Figure 8: The plasmid clearacne process. When the plasmids are hardly replicating, they will be diluted until finally cleared in the cell line.

In an unbiased cell division, the plasmids in a cell are distributed equivalently to both of the children cells. Therefore, while tracing a cell line, we see that a cell lose half of its plasmids after each division (pic plasmid_dilution.png). It is reasonable to assume this process as a random process in which each plasmid has 1/2 probability to "disappear" (actually entering its sibling cell), and 1/2 probability to remain in the cell line we are interested in. This corresponds to "dilution" on a macro scale.

Besides, the plasmids may also replicate itself and bind or unbind to dCas9, These process, together with dilution, make up all the plasmids' behaviors we are interested in (pic plasmid_transition ). When bound to dCas9, the plasmid cannot replicate. Here we introduce two parameters $k_N$ , the frequency that a plasmid bind to a dCas9 molecule, and $k_C$ , the frequency that a plasmid and dCas9 decouple. shows all the probable transition from one state to another.

Figure 9: Probable transitions a cell can perform.

Figure 10: An illustration of the modified Gillespie algorithm.

Another factor taken into account is that the replication of cells are also regulated by the intracellular environment to ensure that the plasmid number do not increase uncontrolled. Hence we assume that during a cell's cycle, the growth of the cell's plasmid number is a logistic process. Precisely, The frequency of plasmid replication is

where $n_N$ and $n_C$ are the number of plasmids unbound or bound to dCas9, $n_\max$ is the maximum number of the plasmid allowed in a single cell, and $\kappa$ is a constant representing the maximal replication rate.

Since the cell cycle length is relatively invariant in a stable environment, we use a fixed cell cycle duration $T$ , and $\lambda = T^{-1}$ as the cell's replication frequency. At each time $nT\ (n\in \mathbb{N^*})$ , the cell's plasmid number undergo a rappid change

where

Putting all these together, we can perform a modified Gillespie simulation over time (similar to traditional Gillespie algorithm except that the dilution happens exactly at $T, 2T, 3T \cdots$ . The parameters are given as follows

parameter	description	value
$T$	time of cell cycle[2]	100 min
$k_N$	dCas9 unbinding frequency to the plasmids	0.067 min-1
$k_C$	dCas9 binding frequency to the plasmids	0min ~ 0.05 min-1
$\kappa$	the plasmid replication rate	0.02857 min-1
$n_\max$	maximal cell copy number in the c[3]	20

We are especially interested in the parameter $k_C$ , which varies with both the binding box's affinity and the concentration of dCas9. Therefore, we run the simulation with different $k_C$ values while other parameters are unchanged . The simulation result shows that when $k_C$ increase, the ratio of dCas9 bound to sgRNA tend to be greater and the cell line run out of the plasmid in less generations. (see pic plasmid_time_series, plasmid_clearance_generation.png)

Figure 11: A typical simulation result. As the value of kc increases, the number of plasmids significantly decreases and reach zero at an earlier time.

Figure 12: Generations of clearance for different kc

It can be seen from Figure 11 that when dcas9 is not expressed ( $k_C=0$ ), the number of plasmid copies remains around 20. However, after the introduction of dcas9 ( $k_C>0$ ), the number of plasmid copies begins to decrease and the proportion of plasmid combined with dcas9 increases. With the increase of $k_C$ , dcas9 can clear plasmids earlier.

To more accurately describe this process, we have simulated for a lot of times (Fig. 12). It can be seen that with the increase of $k_C$ , the number of generations spent to clear plasmids in cell lines will decrease.

It should be noted that in our system, the elimination of the plasmid by one cell line does not mean that the plasmid is completely eliminated in the whole bacterial culture system. Because we did not introduce the degradation and loss of plasmids into our model, the only factor to reduce the concentration of plasmids was dilution due to division, so the number of plasmids in the whole bacterial culture system has been increasing. But its increasing speed is not as fast as that of bacteria. The plasmids are always diluted until the concentration is very close to 0, so we can think that the plasmids are eliminated. The significance of our model is that when the average algebraic value of a cell line to clear plasmids is very small, it can be predicted that plasmids will be diluted to a very thin concentration in a few algebras. In general, under the condition of limited nutrition, cells will inevitably have a certain mortality rate, and plasmids will also have a certain probability of loss. Taking these factors into account, the removal effect will be significant.

References

1. ZHAO Yue-e, ZHU Shun-ya, MA Yong- ping. Control Mechanism of Bacterial Plasmid Replication. LETTERS IN BIOTECHNOLOGY Vol.18 No.3 May, 2007
2. Michelsen O, Teixeira de Mattos MJ, Jensen PR, Hansen FG. Precise determinations of C and D periods by flow cytometry in Escherichia coli K-12 and B/r. Microbiology. 2003 Apr149(Pt 4):1001-10.
3. https://www.qiagen.com/cn/service-and-support/learning-hub/technologies-and-research-topics/plasmid-resource-center/growth-of-bacterial-cultures/#tab2

Gene Expression Noise Control

In this section we model to illustrate how our system can control the expression noise in a cell. Gene expression noise is explained as fluctuation of "very low copy numbers of many components" leading to "large amounts of cell-cell variation observed in isogenic populations"^[1]. This noise can be either intrinsic or extrinsic. Extrinsic noise include most of the environmental factors like the nutrient and antibiotics, gene expression regulation by inhibitor or enhancer, and also the gene copy number in the cell.

In a fast growing cell, the copy number of OriC can exceed ten. In these cases, genes near the OriC may express more than genes far from OriC. Ting Lu^[2] deduced from Helmstetter-Copper^[3] model a relation between the gene's relative location to OriC and gene's copy number:

where $g_q$ is the gene's copy number, $\lambda$ is the cell's growth rate and $x_q$ is the gene's relative location to OriC, and ohter parameters are constants. However, both Helmstetter and Cooper's model and Ting Lu's model are deterministic and unable to be applied to analyze the random factors. Daniel L. Jones et.al included gene copy number variance as a factor of noise, but the model was coarse-grained and only genes with only one replication fork are considered. To fully expose the "noisy" nature of intracellular gene copy number variance, we describe the genome's replication as a stochastic process. We introduce a parameter $k$ representing the frequency that of a single OriC site forming a new replication fork, and $T$ as the time it takes for the replication complex to replicate the whole genome. We deduce that this process is a Yule-Furry process to gene's relative location to OriC with parameter $kT$ .

Figure 13: The formation of replication forks.

In our system, the genome or plasmid DNA replication is blocked by dCas9. This can prevent the genome from forming new replication forks. Specifically, our system decreases $k$ , the frequency of replication fork formation.