0 Overview

In order to further explore the intrinsic mechanism of the genetic circuits designed in this project, to provide a reasonable analysis and prediction of the complex response system, and to be able to provide some feedback and suggestions to the experiments, we have built a mathematical model and resorted to a number of software tools to guide the experiments in obtaining the data used to fit the model and to optimize the model further. This can provide more inspirations for our project.

In this module, our ultimate goal is to verify the feasibility of the DNA damage effect evaluation system and give the optimization direction of the genetic circuit from the perspective of gene components. In order to achieve this purpose, our main work is divided into two parts: (1) In the overall system, based on Hill equation to establish the sensor and amplifier model of the DNA damage detection system, to obtain the amplification and dynamic range of the system, and explore the optimization direction from the values of the parameters obtained from the fitting; (2) In the design of the circuit parts, combining with the results in (1), we analyze the sequences of the promoter RecA and RBS , with the help of a web tool. And then verify that the optimization of the sequences can improve the expression and amplification coefficients of the system in experiments.

Fig.1 Modeling Flowchart

1 System: The Analysis of DNA Damage Effect Evaluation System

1.1 Purpose

Hill equation, which is generally used to describe the ratio of macromolecules saturated by ligand molecules, can be used to determine the degree of receptors binding to an enzyme or receptor, and can well describe the process in which damaged DNA fragments bind to a protein to activate the expression of downstream genes in the whole-cell DNA damage detection system designed in our project. We quantitatively analyzes the system by using NA to damage DNA, measuring the expression intensity data of downstream fluorescent reporter genes, and establishing the relationship between the inputs and fluorescence outputs of different levels of damage.

The modeling process includes: (1) How the binding of LexA inhibits the transcription; (2) How the damaged DNA activate the inhibited promoters after adding the inducer; (3) Establish the relationship between the input and output of the sensor; (4) The process of amplification.

1.2 Method and Model

The modeling process of the sensors and amplifiers of this system can be divided into the following four processes[1]：

(1) Binding of LexA Inhibit Transcription

Whether or not damage inducers are added, repressor protein LexA ( denoted as $L$ ) recognizes the promoter-specific site $D$ and binds to DNA, forming complex $LD$ to inhibit the transcription of downstream genes. The total concentration of the specific site $D$ in the host cell is denoted as $[D_{total}]$, and the concentration of the free site is recorded as $[D^{\ast }]$. According to the law of conservation of matter:

$$\begin{array}{}\text{(1)}& [D^{\ast }]+[LD]=D_{total}\end{array}$$

Set the rate constants for the collision binding of the repressor protein and the site $D$ and for the breakdown process of the complex $LD$ as $k_{on1}$ and $k_{off1}$ respectively. Then the change in the concentration of the complex $LD$ can be expressed by the following differential equation:

$$\begin{array}{}\text{(2)}& \frac{d[LD]}{dt}=k_{on}[L][D^{\ast }]-k_{off}[LD]\end{array}$$

When the process of LexA binding is at steady state, $[LD]$ will no longer change, that is $d[LD]/dt=0$. According to the two equations $(1)(2)$, there is：

$$\begin{array}{}\text{(3)}& \frac{[D^{\ast }]}{[D_{total}]}=\frac{K_d}{K_d+[L]}\end{array}$$

In the equation above, $K_d$ is the apparent dissociation constant of complex $LD$, that is $K_d=[L][D^{\ast }]/[LD]$. This equation describes the probability of site $D$ free in the presence of the repressor protein LexA. In this system, we define $\alpha$ as the basal expression rate of the promoter in the absence of the inducer. The value of $\alpha$ is due to the equilibrium of $LD$ binding and dissociation. According to the equation $(3)$, $0<\alpha <1$ .

(2) The Damaged DNA Binds to RecA to Activate the Inhibited Promoter

Add inducers to cause DNA damage. The broken DNA fragments are recorded as $S$, and several pieces combined with protein RecA $R$ can form the complex $[(S{)}_{n}R]$, which can send a signal to the corresponding LexA to release the inhibitory effect, thereby activating the transcription of downstream genes and reporting expression intensity.

Fig.2 DNA Damage Induces Self-repair Response

The unbound RecA protein is denoted as $R^{\ast }$ and the total concentration is denoted as $[R_{total}]$. According to the law of conservation of matter:

$$\begin{array}{}\text{(4)}& [(S{)}_{n}R]+[R^{\ast }]=[R_{total}]\end{array}$$$$\begin{array}{}\text{(5)}& \frac{d[(S{)}_{n}R]}{dt}=k_{on1}[R^{\ast }][S{]}^n-k_{off1}[(S{)}_{n}R]\end{array}$$

In the equation above, $k_{on1}$ and $k_{off1}$ represent the rate constants for the synthesis and decomposition of complex $(S{)}_{n}R$, respectively.

When the reaction reaches a steady state, $[(S{)}_{n}R]$ no longer changes. Set $K_{d1}=k_{off1}/k_{on1}$, or the apparent dissociation constant of the complex. According to the two equations $(4)(5)$:

$$\begin{array}{}\text{(6)}& \frac{[(S{)}_{n}R]}{[R_{total}]}=\frac{[S{]}^n}{K_{d1}+[S{]}^n}\end{array}$$

This is the Hill equation that describes how damaged DNA binds to the protein RecA. Combined with process (1), the above equation can be understood as the proportion of activated RecA proteins in our system.

(3) The Relationship Between the Input and Output of the Sensor

The activated RecA protein can cause the repressor protein LexA to self-decomposition, thereby eliminating the inhibition of SOS self-repair responses. The equation can represent the activation rate of the promoters whose transcription is inhibited by LexA after the addition of an injuring agent. In order to quantify the expression of the reporter gene, assume the promoter have a maximum transcription rate of $k$, and set the degradation rate of the transcripts of reporter genes as $b$, so the expression rate of the downstream reporter gene $G$ can be expressed by the following equation[2]:

$$\begin{array}{}\text{(7)}& \frac{d[G]}{dt}=k\ast (\alpha +\frac{[S{]}^n}{K_{d1}+[S{]}^n})-b\ast [G]\end{array}$$

$\alpha$ is the proportion of the promoters that do not bind the repressor protein LexA even if no inducer is added, mentioned in the process (1). In our experiment, we measure the fluorescence values at steady state, meaning that $[G]$ no longer changes. So the steady state solution obtained from the equation $(7)$:

$$\begin{array}{}\text{(8)}& [G]=k_1\ast (\alpha +\frac{[S{]}^n}{K_{d1}+[S{]}^n})\end{array}$$

In the equation above, $k_1=\frac{k}{b}$。We describe the relationship between DNA damage and downstream reporter gene expression by this equation.

(4) The HrpRS Complex Binds HrpL to Initiate Transcription and Implement Signal Amplification

Fig.3 Hrp Amplification Circuit

​ The HrpRS complex $RS$ binds to the upstream activation sequence $H$ of the HrpL promoter, which remodels the conformation of the σ54-RNAP-HrpL tight promoter complex into an open conformation for transcriptional activation, resulting in an amplification of the signals detected by the sensor. This process can be described by Hill equation as:

$$\begin{array}{}\text{(9)}& \frac{d[G^{{}^{\prime }}]}{dt}=k_L\ast \frac{[RS{]}^{n_{RS}}}{K_{d2}+[RS{]}^{n_{RS}}}-b_L[G^{{}^{\prime }}]\end{array}$$

In the equation above, $[RS]$ is the concentration of HrpRS complex protein at steady state, $k_L$ and $b_L$ respectively correspond to the maximum expression rate of the HrpL promoter and the degradation rate of the downstream reporter gene $G^{{}^{\prime }}$ expression product, and $K_{d2}$ corresponds to the apparent dissociation constant of the process of the HrpRS complex protein binding to HrpL. Similarly, we measure the fluorescence intensity at steady state, so the input-output relationship of the amplifier is:

$$\begin{array}{}\text{(10)}& [G^{{}^{\prime }}]=k_2\ast \frac{[RS{]}^{n_{RS}}}{K_{d2}+[RS{]}^{n_{RS}}}\end{array}$$

In the equation above, $k_2=\frac{k_L}{b_L}$.

In summary, we describe the input-output relationship of the DNA damage detection sensor by Eq. $(8)$ and establish the input-output relationships of the signal amplification section by Eq.$(10)$:

$$\begin{gathered} [G]=k_1\ast (\alpha +\frac{[S{]}^n}{K_{d1}+[S{]}^n}) \\ [G^{{}^{\prime }}]=k_2\ast \frac{[RS{]}^{n_{RS}}}{K_{d2}+[RS{]}^{n_{RS}}} \end{gathered}$$

The established assumptions include:

1. The measured data as well as the analysis are carried out in the steady state of the system;

2. Generally the activity can only be fully stimulated when the number of ligand-protein binding reaches the corresponding $n$ or $n_{RS}$, so the treatment is simplified by disregarding the case where the binding process does not reach the corresponding $n$ or $n_{RS}$.

Table 1. Notations

1.3 Results

1.3.1 Experimental Data Acquisition

We choose NA as the DNA damage inducer to characterize the effectiveness of the DNA damage detection system. Add a series of DNA damage agents in concentration gradients of 0.3, 0.625, 1.25, 1.5, 5, and 10 μ mol/L to the cultured bacterial fluids, and we measure the absorbance OD values of the samples after full reaction. The experiment data are processed using the following formula:

$$SFU=\frac{RFU}{O{D}_{600}}$$

Where $RFU$ is the relative fluorescence intensity unit determined by the enzyme marker, and $O{D}_{600}$ is the OD value of the same sample, which represents the fluorescence intensity per unit OD. The $SFU$ reflects the expression intensity of fluorescent protein of the bacterium, and the analysis with this can avoid the fluorescence intensity bias caused by the difference in the number of bacteria between different samples.

1.3.2 Fitting and Analysis of the System

We measure the expression of wild strains (without Hrp amplification lines) induced by different concentrations of NA. And Eq. (8) is fitted by nonlinear fitting to obtain the parameters of the sensor input-output model:

Fig.4 The Expression of Wild Strains Induced by Different Concentration of NA

This results in a DNA damage detection sensor input-output model. The model also quantitatively characterizes the RecA promoter transcription level in the designed circuit. Combined with the previous analysis, the parameter $k_1$ is the maximum transcript level of this promoter in steady state, and we consider the output of this model as the transcript input of the amplifier when analyzing the amplifier.

The expression of strains containing Hrp amplification lines at the same concentration gradient are measured as the output of the amplifier. Parameters are obtained after nonlinear fitting:

Fig.5 Responses of the Amplifier

According to the relational equation of the amplifier input and output expressions obtained from the analysis of process (4) in 1.2, with an $adjust-R^2$ of 0.9844. it is obtained that the amplification of the detected signals by this system is 2.809, and the dynamic range of its detection is in the range of 2200~31000 a.u.. The $R^2$ is 0.9775.

The amplifier designed for our system has some discrepancy with the results in Ref. By comparing the values of the parameters fitted in the equations on both sides, it is noted that there is a gap in the value of the parameter $k_2$, which characterizes the maximum expression level in the promoter steady state. In addition, the RBS sequence before the HrpR protein coding sequence that we used is also relatively short. The implications of this part will be discussed in depth in 2.3.

Fig.6 The Comparison between Amplifier Model Prediction and Experimental Data

To further verify the reliability of the model, we used the above model to calculate the output of the amplifier at different times and at different inducer concentrations then compare them with the results of the experiments. The Spearman correlation coefficient of between the model predicted and experimentally characterized responses is 0.88333, indicating a strong correlation between the two with a $p-value$ of 0.00159.

1.3.3 Continuous Monitoring of DNA Damage Effect Evaluation System

In addition, in order to further explore the output of this system over time, we continuously monitor the wild and Hrp strains after the addition of the inducer for 3 hours, and the OD values and fluorescence are obtained every half hour. The concentrations of the inducer are 1, 2, 3, 4, and 5 μmol/L, and three sets of parallel experiments are performed at each concentration to take the average value.

Fig.7 SFU of Wild Wtrains and Hrp Strains in Continuous Monitoring

This shows that the DNA damage detection system reaches relatively stable results after 2-2.5h. And it is very obvious that Hrp strain (i.e., the strain containing the complete DNA damage effect evaluation system) are able to achieve significant results in a much shorter period of time.

2 Parts: Optimization of Genetic Circuit Components

2.1 Overview

In order to further optimize the amplification effect of the system, we analyzed the RecA promoter, which is in the SOS response family of DNA damage, through web tools from the perspective of promoter engineering and RBS transformation, and found HrpR and HrpS in the igem.parts library. The RBS class related to HrpR and HrpS is screened according to its translation initiation rate to achieve the effect of reducing the background noise and improving the sensitivity.

The web tool we used for filtering is RBS calculator from https://salislab.net/software/.

2.2 Analysis of RecA Promoter

2.2.1 Theorem

Stress injury response ( SOS Response ) is a stress response of the organism itself. When bacterial DNA is damaged to a certain extent, it will repair the activity of related genes by activating SOS response. In the genome of E.coli, many genes were found to be involved in the damage repair process of SOS response. Under normal circumstances, the activity of DNA damage repair genes is inhibited by LexA repressor protein. When SOS response occurs, the activity of RecA protein is activated, which triggers the self-cleavage of LexA protein, and then initiates the transcriptional activity of downstream genes. Therefore, the RecA promoter plays an important role in the initiation of biosensors.

In order to further analyze the RecA promoter, we analyzed the RecA promoter based on the Promoter Caculator in online software https://www.denovodna.com/software/.

Promoter Caculator is an online software for predicting transcription rate, which performed large-scale parallel in vitro experiments on the designed promoter sequence with the designed barcode to systematically measure the interaction of control site-specific transcription on the σ70 promoter. Based on these data, Promoter Caculator developed a statistical thermodynamic model to calculate how RNAP / σ70 interacts with any DNA to predict the transcription initiation rate at each location. The model has only 346 interaction energy parameters, but can accurately predict the transcription rate of 22,132 bacterial promoters with different sequences. [3]

The Promoter Caculator established a free energy model to calculate $\mathrm{\Delta }{G}_{total}$ from sequence information. The model considered the interactions associated with each sequence, and then aggregates them into linear, additive contributions.

The free energy model is :

$$\mathrm{\Delta }{G}_{total}=\mathrm{\Delta }{G}_{-10}+\mathrm{\Delta }{G}_{-35}+\mathrm{\Delta }{G}_{up}+\mathrm{\Delta }{G}_{-10ext}+\mathrm{\Delta }{G}_{spacer}+\mathrm{\Delta }{G}_{DISC}+\mathrm{\Delta }{G}_{ITR}$$

2.2.2 Design

We analyzed the RecA promoter fragment used in the experiment, and predicted the different speed records of each site to start transcription, so as to better understand the RecA promoter.

Fig.8 Tx Rate at each position

2.2.3 Discussion

In the 342-346 position of the promoter sequence, we can see a peak of transcription rate, with an average value of about 13442 ( au ), and we further look at the free energy of this part. It can also be seen that the $\mathrm{\Delta }{G}_{total}$ of this part is lower than that of other positions.

Table 2. Free Energy from 342 to 346

2.3 RBS Screening Based on Translation Initiation Rate

2.3.1 Theorem

In the four stages of prokaryote's translation, translation initiation is the rate-limiting step in most cases. Among them, the translation initiation rate is mainly affected by factors such as the RBS sequence on the mRNA, the initiation codon, and the secondary structure before translation initiation, which has an important influence on protein expression.

In order to further improve the expression level of HrpR protein and HrpS protein in the line amplifier, so as to improve the line response sensitivity, we predict the translation initiation rate based on the online software RBS calculator ( https://www.denovodna.com/software/ ).

RBS Calculator is a nucleic acid analysis and design software designed and developed by the University of Pennsylvania. It has the ability to predict the translation initiation rate ( TIR ) of protein coding sequence, which can be directly run online on the website.

Fig.9 Changes in Free Energy of Ribosomal Binding mRNA

RBS Calculator has a powerful analysis and design function: after inputting the protein coding sequence and the host type, the total free energy $\mathrm{\Delta }{G}_{tot}$ changes of several reactions in the translation initiation process are calculated by the established thermodynamic model. These include the free energy $\mathrm{\Delta }{G}_{mRNA:rRNA}$ released when the ribosomal 16SrRNA binds to a specific site of the post-transcriptional sequence, the energy $\mathrm{\Delta }{G}_{start}$ released by the binding of the initiation codon to the tRNA, the energy $\mathrm{\Delta }{G}_{standby}$ released when the alternative site is folded, the free energy compensation $\mathrm{\Delta }{G}_{spacing}$ caused by the distance between the ribosomal 16SrRNA binding site and the initiation codon, and the work $\mathrm{\Delta }{G}_{mRNA}$ required to open the stable secondary structure of the pre-translational mRNA partial interval.[4]

$$\begin{array}{rl}\mathrm{\Delta }{G}_{tot}& =\mathrm{\Delta }{G}_{final}-\mathrm{\Delta }{G}_{initial}\\ & =(\mathrm{\Delta }{G}_{mRNA:rRNA}+\mathrm{\Delta }{G}_{start}+\mathrm{\Delta }{G}_{spacing}-\mathrm{\Delta }{G}_{standby})-\mathrm{\Delta }{G}_{mRNA}\end{array}$$

In addition, we note that ribosomes and mRNA transcripts are in a state of dynamic equilibrium during the exponential growth phase of the host cell. Based on this, we further assume that the translation initiation rate r is proportional to the number of ribosome-mRNA complexes, and derive:

$$\begin{array}{r}r\propto exp(-\beta \mathrm{\Delta }{G}_{tot})\end{array}$$

2.3.2 Design

We obtained a series of RBS sequences Community RBS Collection in the iGEM library, which are suitable for general protein expression in E.coli or other prokaryotes ( including HrpRS protein used in DNA damage detection system ), and cover a series of translation initiation rates. We connected the above sequences with the protein coding sequences of HrpR and HrpS, and predicted the corresponding transcription initiation rate through the RBS Calculator, so as to screen out the RBS sequences of HrpR and HrpS proteins that are most suitable for our application for experiments.

As shown in Table 3., when the encoded protein is HrpR, the starting rate corresponding to each RBS is as follows:

Table 3. The Predicted Transcription Rate of HrpR

When the encoded protein is HrpS, the starting rate of each RBS is as follows:

Table 4. The predicted transcription rate of HrpS

2.3.3 Discussion

In order to verify the reliability of the model for scoring the effect of RBS transformation, we perform experimental characterization of the above RBS sequences. Experiments are conducted with HrpR protein as the target, and four gene lines after RBS modification are constructed preferentially, and transformation and characterization experiments are conducted with the expectation that the change in transcription rate due to RBS modification will be reflected by the change in the average fluorescent protein expression over the same period of time. Currently, we are conducting characterization experiments.

3 References

[1] Alon, U. An Introduction To Systems Biology: Design Principles Of Biological Circuits.

[2] Wang B, Barahona M, Buck M. Engineering modular and tunable genetic amplifiers for scaling transcriptional signals in cascaded gene networks. Nucleic Acids Res. 2014 Aug;42(14):9484-92. doi: 10.1093/nar/gku593. Epub 2014 Jul 16. PMID: 25030903; PMCID: PMC4132719.

[3] LaFleur, T.L., Hossain, A. & Salis, H.M. Automated model-predictive design of synthetic promoters to control transcriptional profiles in bacteria. Nat Commun 13, 5159 (2022). https://doi.org/10.1038/s41467-022-32829-5

[4] Salis, H., Mirsky, E. & Voigt, C. Automated design of synthetic ribosome binding sites to control protein expression. Nat Biotechnol 27, 946–950 (2009). https://doi.org/10.1038/nbt.1568