Team:Peking/Model

HOME

TEAM

Team Members

Collaborations

PROJECT

Description

Design

Experiments

Notebook

Contribution

Results

Demonstrate

Improve

Attributions

PARTS

Parts Overview

Basic Parts

Composite Parts

Part Collection

SAFETY

HUMAN PRACTICES

Human Practices

Education & Engagement

AWARDS

Entrepreneurship

Hardware

Measurement

Model

Plant

Software

JUDGING FORM ⇗

• 
• Home
• Project

  • Description
  • Design
  • Demonstrate
  • Perspective
• Achievement

  • Parts
  • Judging Form
  • Improve
• Modeling
• Human Practices

  • Integrated Human Practices
  • Public Education
• Team

  • Team Members
  • Attributions
  • Collaborations
• Lab

  • Notebook
  • Safety
  • Acknowledgement
  • Measurement

• 
• 
• 

• Home
©2019 PEKING IGEM. All Rights Reserved.
• Home

• Regulatory
• dCas9 & Growth Rate
• Plasmid Copy Number
• Noise Control
• Productivity
• Quorum Sensing

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

$$v = \frac{v_{\max }\cdot[\mathrm{Ara}]^n}{[\mathrm{Ara}]^n + K^n}$$

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:

$$\begin{cases} \displaystyle \frac{\mathrm{d}}{\mathrm{d}t}m = \frac{[\mathrm{Ara}]^{n_m}\alpha_m}{[\mathrm{Ara}]^{n_m}+K_m^{n_m}} - \gamma_m m \\ \displaystyle \frac{\mathrm{d}}{\mathrm{d}t}c = \beta_c m - \gamma_c c + k_-b - k_+cr\\ \displaystyle \frac{\mathrm{d}}{\mathrm{d}t}r = \frac{[\mathrm{IPTG}]^{n_r}\alpha_r}{[\mathrm{IPTG}]^{n_r}+K_r^{n_r}} - \gamma_r r + k_- b - k_+ c r\\ \displaystyle \frac{\mathrm{d}}{\mathrm{d}t}b = -k_- b + k_+ c r - \gamma_b b \end{cases}$$

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

$$\begin{cases} \displaystyle \frac{[\mathrm{Ara}]^{n_m}\alpha_m}{[\mathrm{Ara}]^{n_m}+K_m^{n_m}} - \gamma_m m = 0 \\ \displaystyle \beta_c m - \gamma_c c + k_-b - k_+cr = 0\\ \displaystyle \frac{[\mathrm{IPTG}]^{n_r}\alpha_r}{[\mathrm{IPTG}]^{n_r}+K_r^{n_r}} - \gamma_r r + k_- b - k_+ c r = 0\\ \displaystyle -k_- b + k_+ c r - \gamma_b b = 0 \end{cases}$$

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

$$\frac{k_+ \gamma_b^2}{\gamma_c\gamma_r}b^2 - (k_- + \gamma_b + \frac{k_+ \gamma_b \alpha_r^*}{\gamma_c\gamma_r} + \frac{k_+ \gamma_b \alpha_r^*\beta_c}{\gamma_c\gamma_r\gamma_m})b +\frac{k_+ \alpha_m^* \alpha_r^*\beta_c}{\gamma_c\gamma_r\gamma_m} = 0$$

where

$$\begin{cases} \displaystyle \alpha_r^* = \frac{[\mathrm{IPTG}]^{n_r}\alpha_r}{[\mathrm{IPTG}]^{n_r}+K_r^{n_r}} \\ \displaystyle \alpha_m^* = \frac{[\mathrm{Ara}]^{n_m}\alpha_m}{[\mathrm{Ara}]^{n_m}+K_m^{n_m}} \end{cases}$$

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

$$b = \frac{B - \sqrt{B^2 - 4AC}}{2A}$$

where

$$\begin{cases} \displaystyle A = \frac{k_+ \gamma_b^2}{\gamma_c\gamma_r}b^2 \\ \displaystyle B = k_- + \gamma_b + \frac{k_+ \gamma_b}{\gamma_c\gamma_r}\frac{[\mathrm{IPTG}]^{n_r}\alpha_r}{[\mathrm{IPTG}]^{n_r}+K_r^{n_r}} + \frac{k_+ \gamma_b \beta_c}{\gamma_c\gamma_r\gamma_m}\frac{[\mathrm{IPTG}]^{n_r}\alpha_r}{[\mathrm{IPTG}]^{n_r}+K_r^{n_r}} \\ \displaystyle C = \frac{k_+ \beta_c}{\gamma_c\gamma_r\gamma_m}\frac{[\mathrm{IPTG}]^{n_r}\alpha_r}{[\mathrm{IPTG}]^{n_r}+K_r^{n_r}}\frac{[\mathrm{Ara}]^{n_m}\alpha_m}{[\mathrm{Ara}]^{n_m}+K_m^{n_m}} \end{cases}$$

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. https://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:

$${\bf A} = \begin{pmatrix} -k_1-\kappa & k_1 & \kappa \\ k_2 & -k_2 & 0 \\ 0 & 0 & 0\end{pmatrix}$$

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_1 = k_1^*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

$${\bf F}(t) = \begin{pmatrix} 1 & 0 & 0 \end{pmatrix} \times {\rm e}^{{\bf A}t}$$

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)$ :

$$F(t) = 1 - \frac{\lambda_2(\lambda_1 + k_2){\rm e}^{\lambda_1t} - \lambda_1(\lambda_2 + k_2){\rm e}^{\lambda_2t}}{k_2(\lambda_2 - \lambda_1)}$$

where

$$\begin{align*} \lambda_1 = \frac{-\kappa -k_1 - k_2 - \sqrt{(\kappa + k_1 + k_2)^2 - 4\kappa k_2}}2\\ \lambda_2 = \frac{-\kappa -k_1 - k_2 + \sqrt{(\kappa + k_1 + k_2)^2 - 4\kappa k_2}}2 \end{align*}$$

The average waiting time is:

$${\mathbb E}(T) = \int_0^{+\infty}tF\ '(t){\rm d}t = \frac{k_1 + k_2}{k_2\kappa} = \frac{k_1^*c + k_2}{k_2\kappa}$$

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 (Fig.6). 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 (Fig.7). 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

$$K(n_N, n_C) = \kappa\cdot n_N\cdot\left(1 - \frac{n_N + n_C}{n_\max}\right)$$

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.

$$\frac{\mathrm{d}}{\mathrm{d}t}p_{N, C} = k_C\cdot p_{N+1, C-1} + k_N\cdot p_{N-1, C+1} + K(n_{N-1}, n_C)\cdot p_{N-1, C}-(k_C + k_N + K(n_N, n_C))\cdot p_{N, C}$$

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

$$\begin{cases} n_N \rightarrow n_N'\\ n_C \rightarrow n_C' \end{cases}$$

where

$$n_N' \sim \mathcal{B}(n_N, \frac{1}{2}), \ \ \ n_C' \sim \mathcal{B}(n_C, \frac{1}{2})$$

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 Fig 11, 12)

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:

$$g_{q}(\lambda)=g_{0, q} \exp \left(\lambda\left[\left(1-x_{q}\right) C+D\right]\right)$$

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.

Figure 14: dCas9 controls the copy number and the copy number fluctuation by blocking the replication fork formation

Hypothesis

1. Each cell contains only one set of genome, which may contain multiple replication forks so that each gene's copy number in the genome may be different. The cell divide immediately after genome replication. (Our experiment results show that some cells actually contains more than one set of genome. These cells are very long and in each cell genomes are distant from each other. We treat this kind of cell as a chain of multiple cells)
2. The cell containing the studied genome has been replicating exponentially for several generations in a stable environment, thus the replication fork's distribution on the genome is steady (i.e. sampled from a fixed distribution)
3. The replication fork only forms at OriC and the formations of replication forks are independent with each other.
4. The frequency that of a single OriC site forming a new replication fork $k$ and the time it takes for the replication complex to replicate the whole genome $T$ are constant.
5. Our system works by decreasing $k$ . This process is a fast grade 1 process.

Gene Copy Number Fluctuation and its Control

To study the copy number of each in a given branchy genome, we firstly need to know how this genome and its replication forks are formed. Considering a newly formed gene. $T$ ago, the replication fork forming the current genome was just newly formed and the OriC number of this branch is 1. In the next $T$ period, DNA replication complex bind to OriC to form new replication fork, creating a new OriC for this branch. According to hypothesis 3 and 4, the copy increasing process of OriC through the time period $T$ is a Yule-Furry process (a stochastic counterpart to the deterministic exponential growth model) with parameter $k$ . According to the theory of Yule-Furry process, the probabilty that the number of OriC is $m$ at time $t\ (0\leq t \leq T)$ is:

$$p_m(t) = \mathrm{Pr}[N(t) = m] = \mathrm{e}^{-k t}(1 - \mathrm{e}^{-k t})^{m - 1}$$

(Note: according to hypothesis 1, the initial number of OriC is 1, therefore we are using the formula under the condition $n_0=1$ )

Furthermore, the DNA replication complex is moving from OriC to Ter in a constant speed after forming the replication fork at OriC. Suppose we have known that the number of OriC at time $t$ is $m$ , then after $T - t$ forms the genome we are studying. In this $T - t$ period, all the DNA replication complex move forward a distance $(T - t)/T$ relative to the genome. Therefore, after $T$ the copy number of gene at $(T - t)/T$ relative to OriC, or $t/T$ relative to Ter, equals to $m$ , the number of OriC at time $t$ .

Figure 15: Because all the replication fork moves in a constant speed, the copy number of the gene also "moves". This allows us to estimate the gene copy number of every site at every time.

We introduce $x=t/T$ as a site's relative distance to Ter in the genome ( $0\leq x \leq 1$ , $x = 0$ corresponding to Ter, $x = 1$ corresponding to OriC), then we get

$$p_m(x) = \mathrm{e}^{-k T x}(1 - \mathrm{e}^{-k T x})^{m - 1}$$

This equation determines the copy number distribution of all the genes at all the sites. We can further deduce

$$\mathbb{E}[N(x)] = \mathrm{e}^{k Tx}$$

and

$$\mathbb{Var}[N(x)] = \mathrm{e}^{2k Tx} - \mathrm{e}^{k Tx}$$

The nearer the gene is to OriC, the greater $x$ is, the greater both $\mathbb{E}[N(x)]$ and $\mathbb{Var}[N(x)]$ are. $\mathbb{Var}[N(0)]=0$ and $\mathbb{E}[N(0)] = 1$ corresponds to the model hypothesis that the copy number of Ter is invariantly 1.

Variance itself is not sufficient enough to indicate the intensity of fluctuation because the average expression is also increasing. We apply CV value, the ratio of standard deviation to the average, to more precisely describe the intensity of fluctuation:

$$\mathrm{CV}[N(x)] = \frac{\sqrt{\mathbb{Var}N(x)}}{\mathbb{E}N(x)} = \sqrt{1-\mathrm{e}^{-k T x}}$$

$\mathrm{CV}[N(0)]=0$ , indicating that that the copy number of Ter is invariant. $\mathrm{CV}[N(x)]$ increases with $x$ , indicating that the copy number of genes near OriC in sequence tend to vary more than that of genes far from OriC in sequence.

It is obvious that the CV value decreases when the value of $k$ and $x$ becomes smaller. To further demonstrate this relation, we numerically calculate the analytical results and visualize them. The parameters' values are given as follows:

parameter	description	value
$T$	Bacteria replication time[4]	41 min
$k$	binding frequency of dCas9 to the genome	0~0.04 min-1

$T$ is relatively stable^[2]. However, $k$ , as stated before, can be hacked by our system. Therefore we are extremely interested in figuring out how the noise changes with $k$ . Fig 16 shows resulting noise amplitude at different $k$ value. It turns out that the noise amplitude is significantly greater when $k$ is bigger. Fig.17 shows the copy number noise of gene at the middle of the entire genome varying with $k$ .

Figure 16: The gene copy number noise increases when k is small and when the gene is located near OriC (or far from Ter)

Figure 17 left: Copy number CV for different k and genes of different sites; right: CV for the copy number of gene at the middle of the genome at different k.

In this section, we explore the potential of our system as a bacterial noise control component. The change of gene copy number is an important part of expression noise. By controlling the replication of bacterial genome, our system can reduce the fluctuation of gene copy number to control the noise of bacterial transcription, so as to increase the stability of various synthetic biological gene circuits.

References

1. Elowitz, Michael B., et al. Science 297.5584 (2002): 1183-1186.
2. Chen Liao, Andrew E. Blanchard and Ting Lu. An integrative circuit–host modelling framework for predicting synthetic gene network behaviours. Nature Microbiology volume 2, pages1658–1666 (2017)
3. Cooper S, Helmstetter CE. Chromosome replication and the division cycle of Escherichia coli B/r. J Mol Biol. 1968 Feb 14;31(3):519-40.
4. Helmstetter CE. DNA synthesis during the division cycle of rapidly growing Escherichia coli B/r. J Mol Biol. 1968 Feb 14 31(3):507-18. P.518 top paragraph
5. Blattner FR1, Plunkett G 3rd, Bloch CA, Perna NT, Burland V, Riley M, Collado-Vides J, Glasner JD, Rode CK, Mayhew GF, Gregor J, Davis NW, Kirkpatrick HA, Goeden MA, Rose DJ, Mau B, Shao Y. The complete genome sequence of Escherichia coli K-12. Science. 1997 Sep 5;277(5331):1453-62.
6. Méchali M. Eukaryotic DNA replication origins: many choices for appropriate answers. Nat Rev Mol Cell Biol. 2010 Oct11(10):728-38. doi: 10.1038/nrm2976. p.728 left column 2nd paragraph

Productivity

In this section we provide a model to explain why our system is able to increase the cells' productivity to multiple types of products like GFP and indigo. It is somehow counterintuitive to think that a decline in microbiomes' growth rate may increase the microbial productivity. Therefore, it is necessary to carefully test this issue in a quantitive method.

The undesired productivity of microbiomes is thwarting biological production technology from being more widely applicable. Therefore, it is always concerned how to improve the productivity of microbiomes, or how to enforce these wild or engineered cells to turn more of the substrates we feed them into the bio-products we want.

A non-negligible factor thwarding the bicrobiomes from producing our bioproduct is the growth of the bicrobiomes themselve. All the cell products we want to yield from the microbiomes are products of the microbiomes' own metabolism forming the microbiomes' own biomass. Therefore, their production is under strict regulation of the microbiomes. In most engineered microbiomes, the genes we introduce to the microbiomes are foreign genes unrelated to the cells' growth and normal metabolism, and the production of these genes are undesired for the microbiomes themselves. The engineered microbiomes, stressed by these products, will activate all the possible pathways to change its gene expression pattern and to reallocate its nutrients and enzymes to produce more "necessities" for themselves. This is particularly true for bicrobiomes under steady growth rate, whose nutrients are mainly used for growth.

Actually, many studies concerns the relationship between microbiomes' growth and the way they allocate their resources. Terence Hwa et.al^[1] studied E. coli's metabolism under different nutrient and antibiotic condition by dividing the cells' proteins into different sectors. Particularly in this work, a sector defined as "Unnecessary Expression"(corresponding to the bio-products we want) is found to be negatively related to the growth rate. The authors described this relation semi-quantitively as

$$\lambda(\phi_U) = \lambda(\phi_U = 0)\left(1 - \frac{\phi_U}{\phi_c}\right)$$

were $\phi_U$ is the ratio of "Unnecessary" protein mass to the total protein mass of a cell, and $\lambda$ is the growth rate corresponding to the ratio $\phi_U$ . Ting Lu et. al^[2] built a coarse-grained whole-cell model including the productions and functions of proteins of different sectors and further validated Hwa's model. In this model, the regulation performed by ppGpp is included and this provides an explanation of the phenomenological result in [1]. ppGpp is an important gene expression regulation molecule responding to various types of environmental stress and cell's abnormal pnysiology. One of its function is to down-regulate the expression of ribosomal RNA, protein and expression affiliated proteins while up-regulate the expression of some enzymes relating to the cell's core metabolism, adapting the cell from a fast-growing state to a slow-growing state. Some optimization models not explicitly including ppGpp also discovers a similarity between cells' optimal resource allocation strategies and the ppGpp regulation strategies.

All the works mentioned above regard the growth rate as the result or equivalence of the microbiomes' protein accumulation. In our system, however, the growth rate is hacked. Because of this, the causal relationship between the cell physiology (growth rate) and the nutrient reallocation is different from these works. A better explanation of the mechanics of our system is that the hacked growth rate freed the cells from producing too much growth-necessary proteins and enabled them to produce the bio-product we want. We emphasize that our model is still a phenominological model in which the resource allocation regulation is finished in a "black box", not explicitly related to the regulation of any single pathway.

Figure 18: The nutrient flow of a bacterium. A bacterium uptakes some nutrients and use them for both replication and production. It allocates the nutrients according its growth rate.

In our model, there is a single-source nutrient. An engineered cell uptakes the nutrient and uses it both for its own growth and for the production of the bio-product. The ratio of nutrient used for cell's own growth to the total nutrient uptaken by the microbiomes is $\alpha$ , and the ratio of nutrient used for bio-product production to the total nutrient uptaken by the microbiomes is $1 - \alpha$ . Thus, $\alpha=0$ corresponds to a "production-only" state and $\alpha=1$ corresponds to a "growth-only" state. This process is written as an ODE system

$$\begin{cases} \displaystyle \frac{\mathrm{d}}{\mathrm{d}t}n = -\frac{k_n n}{n+K_n}m\\ \displaystyle \frac{\mathrm{d}}{\mathrm{d}t}m = \frac{k_m n}{n+K_n}m\alpha\\ \displaystyle \frac{\mathrm{d}}{\mathrm{d}t}p = \frac{k_p n}{n+K_n}m(1-\alpha)\\ \end{cases}$$

where $n$ is the mass of nutrient, $m$ is the biomass, $p$ is the mass of the bioproduct, $k_n$ is the cell's nutrient uptake rate, $k_m$ is the maximum growth rate when the nutrient is sufficient and all the uptaken nutrient are used for cell's own growth, and $k_p$ is maximum production rate when the nutrient is sufficient and all the uptaken nutrient are used for the bio-product's production, $K_n$ is the nutrient corresponding to half-maximum nutrient uptaking rate, and $\alpha$ is the above-mentioned allocation ratio.

Noticing that for a given $\alpha$ , when the nutrient is sufficient, the second ODE can be re-written as

$$\frac{\mathrm{d}}{\mathrm{d}t}m = k_m\alpha\cdot m$$

This corresponds to an exponential growth with a growth rate $\lambda = k_m\alpha$ .

$$\begin{cases} \displaystyle \frac{\mathrm{d}}{\mathrm{d}t}n = -\frac{k_n n}{n+K_n}m\\ \displaystyle \frac{\mathrm{d}}{\mathrm{d}t}m = \frac{n\lambda}{n+K_n}m\\ \displaystyle \frac{\mathrm{d}}{\mathrm{d}t}p = \frac{k_p n}{n+K_n}\left(1-\frac{\lambda}{k_m}\right)m\\ \end{cases}$$

the equation explicitly includes the growth rate. The parameter values are given as follows

parameter	description	value
$k_m$	the theoretical maximal growth rate of E. coli[1]	2.85 h-1
$k_n$	the maximal nutrient uptake rate of E. coli [3]	5.70 h-1
$k_p$	the maximal producting rate of the bio-product	0.1425 h-1
$K_n$	the nutrient concentration corresponding to half-maximal growth rate	4.0 kg/ml
$\lambda$	the growth rate	0~ $k_m$

Figure 19: The dynamic of the amount of bacteria, nutrient and product at different growth rate.

The ODE is difficult to solve but we can analyze the steady states of it.

suppose initially $m(0) = m_0, n(0) = n_0, p(0) = p_0$ , noticing that

$$\frac{1}{k_n}\frac{\mathrm{d}}{\mathrm{d}t}n + \frac{1}{\lambda}\frac{\mathrm{d}}{\mathrm{d}t}m = \frac{1}{k_n}\frac{\mathrm{d}}{\mathrm{d}t}n + \frac{k_m}{k_p(k_m-\lambda)}\frac{\mathrm{d}}{\mathrm{d}t}p = 0$$

thus we know that $\displaystyle \frac{n}{k_n} + \frac{m}{\lambda}$ and $\displaystyle \frac{n}{k_n} + \frac{k_m p}{k_p(k_m-\lambda)}$ remain constant in the entire dynamic process. Provided the initial condition, we can deduce

$$\begin{cases} \displaystyle \frac{n}{k_n} + \frac{m}{\lambda} = \frac{n_0}{k_n} + \frac{m_0}{\lambda} \\ \displaystyle \frac{n}{k_n} + \frac{k_m p}{k_p(k_m-\lambda)} = \frac{n_0}{k_n} + \frac{k_m p_0}{k_p(k_m-\lambda)} \end{cases}$$

holds throughout the growth-and-production process. Finally, the nutrients are exhausted and the microbiomes stop both growing and producing. Denoting the concentration of all the materials at the final state as $n^{\mathrm{ss}}, m^{\mathrm{ss}}$ and $p^{\mathrm{ss}}$ , the above-mentioned equation still holds:

$$\begin{cases} \displaystyle \frac{n^{\mathrm{ss}}}{k_n} + \frac{m^{\mathrm{ss}}}{\lambda} = \frac{n_0}{k_n} + \frac{m_0}{\lambda} \\ \displaystyle \frac{n^{\mathrm{ss}}}{k_n} + \frac{k_m p^{\mathrm{ss}}}{k_p(k_m-\lambda)} = \frac{n_0}{k_n} + \frac{k_m p_0}{k_p(k_m-\lambda)} \end{cases}$$

moreover, at the final stage, the nutrients are exhausted so $n^{\mathrm{ss}}=0$ . Then we can solve

$$\begin{cases} \displaystyle m^{\mathrm{ss}} = m_0 + \frac{\lambda n_0}{k_n} \\ \displaystyle p^{\mathrm{ss}} = p_0 + \frac{n_0k_p}{k_n}(1 - \frac{\lambda}{k_m}) \end{cases}$$

It's worth noting that $p^{\mathrm{ss}}$ is negatively and linearly related to the growth rate. This result is in accordance with the empirical formula provided in [1].

In this section, we discussed the possibility of our system being used to increase production of biological products. Our experiments also confirmed that in our system, the yield of GFP and indigo would increase with the increase of the concentration of inducer. It is consistent with the prediction of our model. In general, we can control the growth rate of bacteria by reducing nutrition, and our system has successfully controlled the growth rate of bacteria under the condition of sufficient nutrition. At this time, the energy of bacteria will be allocated to more other metabolic activities, specifically, those metabolic activities that are not so closely related to their core life activities. This is the core mechanism that our model will explain.

References

1. Matthew Scott, Carl W. Gunderson, Eduard M. Mateescu, Zhongge Zhang, Terence Hwa. Interdependence of Cell Growth and Gene Expression: Origins and Consequences. Science 330, 1099 (2010)
2. Chen Liao, Andrew E. Blanchard and Ting Lu. An integrative circuit–host modelling framework for predicting synthetic gene network behaviours. Nature Microbiology volume 2, pages1658–1666 (2017)
3. Smil, 1998, Energies: An illustrated guide to the biosphere and civilization, MIT Press

Quorum Sensing

In this section we model to illustrate how our system works coupling with a quorum sensing system. We want out system to realize the population's auto-regulation - the cells stop growing fast when they sense a lot of other cells crowding around it. We realize this by applying LuxI-LuxR system, in which a kind of small molecule call AHL. This system was firstly discovered in V. fischeri as a means of intercellular communication. The cells produce exceeding amount of AHL which can move either into or out of the cells. When the population is large or dense in a region, the local AHL concentration increases and the cells sense this high concentration and respond to this by up- or down-regulating some genes' expression. In our system, the sensing of AHL results in an increment of dCas9, and consequently a down regulation in growth or replication rate (see design for the design and experiments)

A Simulation of Donor-Receptor Experiment

Our first series of experiments involve the testing of our fore-mentioned logic. We separately introduced the AHL producing parts and the AHL sensing parts to two strains of E. coli. The former is called "donor" and the latter is called "recipient". The recipient cells are evenly coated onto the solid medium and the donor cells were dropped at the center of the medium. It is expected that the donors AHL, the AHL difusses around and inhibits the growth of the receptors around it (see design)

Here we use a simple diffusion model and visualized simulation to describe this process. We suppose that the solid media is a 2D plane, the recipients are uniformly distributed on the plane, and the donor is dense at the center $(0, 0)$ of the media (subjected to a normal distribution), thus the AHL production is subjected to a normal distribution centered at $(0, 0)$ . A partial differential equation (PDE) can be derived from these hypothethes:

$$\frac{\partial u}{\partial t} = D(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2})u + C\mathrm{e}^{-\frac{x^2 + y^2}{\sigma^2}}$$

where $x, y$ are coordinates, $\sigma$ represents the radius of the donor colony, $C$ and $D$ denotes the diffusion constant. The randomly distributed recipients stop growing when the local concentration reach 10^-6μM/L.

Here we numerically solve the equation with finite difference method. We draw an animation visualize this dynamic process. Furthermore, we randomly place 400 receptors on the plane to show their colony formation. It can be seen that the receptors distant from the center grow into colonies, while most of the receptors nearby the center grow into relatively small colonies or cannot grow into colonies.

Figure 20: The simulation result.

Reference

1. Lupp C1, Ruby EG. Vibrio fischeri uses two quorum-sensing systems for the regulation of early and late colonization factors. J Bacteriol. 2005 Jun;187(11):3620-9.