Model Description

Model Description

Introduction

Our model consists of four parts. The first is our basic but pillar design of three delicate circuits. They are, separately, the Self-renewal design, which depends on the Self-renewal LuxR in our system to ensure reusability; the threshold guard, which turns on only when the inducer strength exceeds a certain threshold; and the suicide design, depends on IPTG tag. It deserves special mention that we modeled and found that the threshold may be between 0.6 and 0.8 relative to unit strength, so we designed an experiment for our wet lab, surprisingly found the real threshold, and determined its strength to be 680 $\mu mol/L$ .

The second is our visualization of the diffusion process. In this part, we use basic physical rules innovatively combined with computer vision and convolution on two-dimensional images. Finally, we designed a diffusion APP interface with custom interaction and visualization.

The third part is the model of a calibration curve, finely determining the relation between our inducer and reporter. This curve is expected to be well fitted to the experimental data from our wet lab, so it can be more practical. In this part, we fundamentally established a series of rules referred to classical automatic control theory. Each type of transfer function in the frequency domain represents a certain type of rule in the biological context.

The forth part is the interaction with Hardware. We use COMSOL to simulate the streamline in different sizes of containers, and finally get the suitable size for our hardware.

The last part is the data exchanged with our wet lab. We fit the threshold of the threshold guard, calculate the inhibition coefficient for further modeling. We fit the attenuation coefficient of fluorescence, and according to the relationship between fluorescence and product, the product decay coefficient in the model can be obtained. Details are as follows. Codes and references can be avilable via our link.

a. Three Basic Models

Self-renewal Design

Given the balance between efficiency and economy, the advantages and disadvantages of the two designs were compared before the experiment for the experimental group to choose. One is the completely positive feedback design, and the other is the self-renewal design with $LuxR$ as the restrictions. We set cl in the downstream of lucl and Cl can inhibit the promoter of luxR, so the total amount of $LuxR$ can be seen as a constant after trigger. The completely positive design is built by Simbiology, and for comparation the renewal design is built too. Their comparation are summaried as follows.

We use GFP as our reporter here.

a. Completely Positive Feedback Design

The Amount of LuxR and GFP vs Time after Trigger in Positive Feedback Design

Using the logarithmic time axis, it shows that LuxR, RA, and GFP rose from front to back with a certain lag after adding the trigger. But in logarithmic coordinates, we can see that the lag action is 10s-20s. This is the timeliness of the positive feedback of the first gene line.

b. Self-renewal Design

The Amount of LuxR and GFP vs Time after Trigger in Self-renewal Design

The logarithmic coordinate of substance concentration shows that LuxR will have an exponential decline after trigger, and GFP will also have an exponential decline after a lag period, and finally will be below 0.1.

We hope to reduce the cost of using the recycling system in the future.

Overall circuit

Overall circuit of Self-renewal Design

To deduct the transition model to demonstrate the relationship between the initial total amount of LuxR and our production, that is, P, we derive a series of equations and then come to a final function as $P(R, L)$ using equations and physical formulas.

Where R is the Free LuxR in the system and L is the free AHL in the system. They are, separately, controller and responsor of our gene circuit.

Notation

Note that the mRNA and DNA of each component are represented by the combination of two parts, the above substances as subscripts and the symbols as different rates. Symbols are as follows. We use $\mu$ to denote the transcription rate of genes, $\sigma$ to denote the degradation of genes, $k$ to denote the translation rate of RNA and $d$ to denote the degradation of proteins. Additionally, we use $f$ to represent the expression rate of genes, $K$ to represent the forward reaction rate, and $\mu_0$ to express the basal rate of each reaction.

Assumptions

1. The number of genes is constant.
2. A set of ordinary differential equations (ODEs) describes the average behavior of a cell population in a deterministic approach.
3. Transfer functions are measured at the steady state.
4. The total amount of R is constant and controllable.
5. All chemical combinations reach equilibrium quickly.

Reactions and substances

The responsor AHL

We first ignore the binding effect of AHL and LuxR to test the reponse effect. The detected small molecules are enzymatically disassembled and enter the system to bind to XylS, which promotes the expression of the downstream gene luxI from the Pm promoter, activating the synthesis of AHL. The process can be divided into the following sections.

1. Receptor-ligand binding.
2. LuxI transition.
3. AHL synthesis.

For the first part, we wish to use the diffusion model to estimate the delay time. For the second and third processes we use ODEs and assume that the transcription of luxI reaches equilibrium quickly, thus the mRNA of AHL in the system is constant.

We denote the activator as A, the activated gene luxI as [ALuxI], and the magnification effect is represented by $M_l$. Our equations are summarized as follows:

As stated, the second reaction reaches equilibrium quickly, and the basal rates are ignored, we have a stable concentration of AHL

as well as the dynamic concentration of AHL:

The controller LuxR

As to LuxR, we design a renewable system based on negative feedback. The expression of cl in the downstream of luxI can inhibit the production of LuxR. Given the high amount of initial LuxR and the great deterrent effect of Cl, we assume that when reactions are triggered, our LuxR will not be produced but only be consumed. Hence we may reach a similar assumption that the total amount of LuxR after triggered is constant.

Firstly, we rigorously derive the LuxR content as follows.

Note that $f_R$ is the expression rate of [DNA_R], and it can be inhibited by Cl, so its size can be estimated as:

When n is 0, $f_R = f_{Rm}/2$, and when Cl is 0, $f_R = f_{Rm}$.

Therefore we have the differential equation of R:

If equating Cl with AHL, as it is in equation (2-1), we further have:

The reporter P

Considering the mRNA of our product P, note as $mRNA_P$, is transcribed by the combination of DNA and $R_{2}L$, that is $Q_P$, hence

Where $\sigma_P$ represents the degradation rate of mRNA, $\mu_P$ is the transcription rate of gene p regulated by the $Q_p$, and $\mu_0$ is the transcription rate of gene in the absence of the transcription factor. Since we use an orthogonal system to avoid leakage, $\mu_0$ can be ignored.

The chemical combinations are as follows.

The second equation represents the dimerization of R, and the third represents the combination of AHL and the dimerization of R. Based on assumption 5), at equilibrium for each of these transformations, we can get their equilibrium concentration:

Where $\rho,K_1,K_2$ represent the forward reaction rate of each formulation aligned in (5).

The production of P is based on the concentration of $mRNA_p$, with the same form of ODE, we have

Where $k_P$ is the translation rate of $mRNA_P$, $d_P$ is the degradation rate of P.

With the equation above, and $L$ represents $AHL$ deducted in responsor part, we can obtain the concentration of P as:

Since gene $Q$ is seen as a constant, $P$ is a function of $R, L$, which is our controller. Note that they are theoritically both the free amount of each component. And the combination effect should be taken into consideration so the free $L$ and $R$ can be thought like equation (4-2) and (2-3). But we will use a different method to explain these..

Differently, AHL can be complemented by the expression of luxI but LuxR and will delete more quickly since its DNA is inhibited by Cl. So we let LuxR represent the total R, and R is the free LuxR.

Conclusion

We take the chemical reaction as the second process, and the production as the first process. And the chemical reactions take place after the first process reaches equilibrium, that is to say, we break the (4-2) and (2-3) into two parts.

Firstly, we calculate the total amount of $AHL$ and $LuxR$ with the three equations (2-2) (4-1), (8)

Free AHL will be modified according to (11-1) and (5). Free LuxR can be calculated according to (11-2) and (5). And combined with the equation (12), that is

The substitution is as follows:

Where $M_p$ can be translated Into the total magnifying effect of our system, and other $\tau,K$ is

the identification parameters of our transition system, which we will simulate using the classic controlling theory. We will simulate both the dynamic and the stable process of our reporter and controller in following parts.

Results

We will simulate this model in Part c and get results there.

Threshold guard

Overall circuit

Overall circuit of Threshold Guard

Notation

We decompose the function into three parts, the activation of our gene, the transcription of both gene( $G$) and activated gene( $G_a$), and the translation of the mRNA( $M$) as the basic part of this design. And to ensure the effectiveness and orthogonality of the switch, we hope the activation effect, which can be estimated as $k_1/k_{-1}$ to be higher, as well as the expression rate of the activated gene.

Reactions and subatances

The basal part as described as follows:

To measure the negative effect of the sRNA(mS) and repressor(R), we use two sets of principles. First, the transcriptional repressor R will inhibit the translation of the gene, and we used decreased gene transcriptional activity to characterize these effects. Additionally, we use the Hill equation to model the negative effect, that is

Since the repression effect of sRNA is not identical, we use another method to model the repressive effect of mS on RBS. In our design, the structure of sRNA is specific and can recruit the RNA chaperone $Hfq$, which is known to facilitate the hybridization of sRNA and target mRNA as well as RNA degradation. We assume the number of $Hfp$ is enough when sRNA is enough, and then the effect can be measured as both combination of sRNA and mRNA and the facilitation of mRNA degradation.

The combination reaction is as

We assumed $k_{−4} = 0$ to simplify the model.

Results

a. We used the discretization method to simulate and get the concentration of two repressors mS and R at the equilibrium state. We can see in the left figure that with the inducer, when the repression strength is stronger, the number of R will be higher when equilibrium. On the contrary, without the inducer, when the repression strength is stronger, the number of mS will be higher when equilibrium. It indicates that the two repressors will work in different states in our orthogonal design.

The Amount of R vs mS in Stable State with Different Trigger

b. To test the sensitivity of the threshold guard, and for the smaller consumption of calculating cost, we use the result in a., and we accordingly set repression strength between 5 and 10. Then we will control a set of inducer levels to see how the amount of R will change with time. The results are as follows.

The Amount of R vs Time with Different Inducer Strength

left(Repressor strength 10), right(Repressor strength 5)

From result b., we can see again that, on one hand, when the repression strength is stronger, the threshold will be lower, and the response curve is steeper, which means the switch will be more sensitive. On the other hand, when the repression strength is weak, the threshold will be higher, and the tolerance to noise will increase. We need to choose the ideal set of threshold and repression strength to ensure both sensitivity and anti-noise ability. We get the suitable level of inducer from b., and set the inducer level between 0.56 and 1. The amount of R at equilibrium state are measured. The R and A are the amount of $\log(·)$.

The Amount of R vs A in Stable State with Different Repressor Strength

We can expect the threshold at about 70% of the max inducer amount. We designed an experiment for our wet lab, surprisingly found the real threshold, and determined its strength to be 680 $\mu mol/L$ , and measured the rough strength of repressor to be exceeding 10.

Suicide Design

b. Diffusion Model

Principles

Thermal Diffusion

1-Dim Discrete Form

2-Dim Discrete Form

Translation

Formula All-in-one

Display Logarithmically

The Demo of Diffusion APP

c. Calibration Model

Basic transformation rules

We build this model using cybernetics. To create the one-to-one relation of the frequency-domain functions and the time-domain differential equations, which are usually used in the biological context, we fundamentally established a series of rules referred to classical automatic control theory. Each type of transfer function in the frequency domain represents a certain type of rule in a biological context, such as transcription or degradation. We will take some examples here.

Examples

a. Transfer functions

• Transcription & translation

Gene Transcription

mRNA Translation

So the transfer function is

• Degradation

Degradation

Here we use the transformation rule as

Here $R(s) = \dfrac{1}{s}$. So the transfer function is

• Combination of Transcription and Degradation

So the transfer function is

b. Sources

• Constant

Constant Input

If the source is constant, then the transfer function is

The K is the value of a constant.

• Periodic

Periodic Input

If the source is periodic, it can take a variety of forms, such as sine wave and sawtooth wave. For example, if we use sine wave, we have its transfer function as

The $\omega$ is the frequency of source.

Simulink Simulation

Basic blocks and linealization

Using our previous work in Part a.1 self-renewal design, we can change the system can into following blocks, the signal will flow to final reporter P, and the simulus would be A, the controller will be free R.

The Block Diagram of LuxR

The Block Diagram of AHL

And we can Linearlize the equation by apply the Taylor's formula:

If we work at the stable work point, we can substitute two magnifying parameters to build a Linear System. However, with the aid of Simulink, we can product two signals and need linearization. Notwithstanding the simplicity of product them, we propose that with the Linearized Model, we can apply more classic control theories such as the Bode graph.

And the whole graph will be like Figure 17.

The BLock Diagram of Activated System

The Transfer Function

Then we should calculate the transfer function $G_1$ and $G_2$, and note that there is a problem with the function $G_1$. As we stated before, if the system is a linear model, it would be simple to rectify it. But as shown in equation (11-2), the transfer function between A and the out free R is a polynomial type. It comes two solutions, one is the Tayor formula applied in the stable work point, and the other is the real model using the initial function. Here we apply the first one, for we just have two work points, that is, A is 0 and A is 1. The transfer function will be:

Then we substitude $A_0$ with constant 1 and 0, and we have:

So we adjust our model shown by Figure 17. The absense and presence of A channels will be denote as two signals with different transfer functions. As is shown by Figure 18.

The Block Diagram of The Whole System

Through the correspondence between the time domain and frequency domain in automatic control theory, and more details and code will be available in our Codes and References, we obtain that

Finally, we built our SISO system, which can measure the transfer of input and output both statically and dynamically, giving help to us to analyze the effectiveness of our design. Our Wet Lab gives us hints of parameters and some measurements, then we will fit our model to our data in the final session.

SimuLink Sitimulation of SISO System

After deriving the system transfer functions, we use the Control system theory to build a complete block of our design, and we first apply the parameters naively, then adjust them to fit our results, we will further change the peak input of Activtaor, as shown in the left of Figure 19. , to measure the output of the P, as shown in the right side.

The System Build by Simulink

Results

Product vs Inducer in Stable State

The Naive parameters work well with our assumptions, and we can see that the LuxR can be replenished after the activator functions, and P can reach a certain amount after the trigger.

The Product vs Time under Different Step Inputs

d. Hardware Simulation with COMSOL

We collaborated with the hardware group to design the reaction chamber: By using COMSOL Multiphysics to model the diffusion and fluorescence of the analyte introduced into the reaction chamber, we evaluated the effects of different chamber shapes on the fluorescence emission.

In this model, to simplify the process, the analyte diffusion was modeled as a fluid with streamlines and flux tracking the diffusion and distribution over time. The light intensity on the photoresistor was then evaluated at specific time points based on the analyte distribution, assessing the detector design.

The fluorescence in the different size of containers are shown in the gif, and we finally choose the 0.2 * 0.2 * 0.4 container. We can see that as time goes by, the fluorescence strength is growing gradually, and the streamline is relatively stable.

a. Different sizes of containers

Different Sizes of Containers

We can see when the container is smaller, the streamline is sparser. But the cost will be higher when the size is bigger. To balance it we choose the middle one.

b. 3D point diffusion

After choosing the size of the container, we further explore the real situation of our diffusion -point diffusion. We take a small square on the surface of the container as the diffusion starting point, and the concentration in the container grows like below.

The 3D Point Diffusion

c. 2D point diffusion

To further simulate the real point diffusion model, we use a two-dimensional cross-section to simulate. We take the container's width and height as the rectangle's length and width. To simulate point diffusion, we set a small semicircle as the starting point, set the flux to 0 on the container wall, and set the initial diffusion concentration inside the diffusion starting semicircle. Then we can see from the simple simulation that at the beginning the concentration in the container is low, and in the end, the concentration in the container is high and the concentration in the semicircle becomes low, represented by the color referred to the color bar.

The 2D Point Diffusion

e. Data fit

Threshold of threshold guard and Strength of repressor

Since we model the threshold guard with computer discretization method and find a threshold between 0.6 and 0.8 of the whole strength, we direct our wet lab to conduct a series of experiments to find out the real threshold, and to determine the growth coefficient of the product in terms of the inducer.

Our wet lab design the concentration gradient between 100 and 20000, that is the amount of A. In order to get the amount of R, we use the GFP, which is closely related to R downstream, to show the final product. Finally we get the double log scatter plot, as shown in Figure 25. Each point is the average of three sets of repeated experiments.

Fluorescence Intensity vs Inducer Strength of Wet Lab

We use the logarithm of the product and inducer to represent the strength of the them. And in order to keep consistent of the simulation, we take logarithms of the product strength again to find the threshold. From the Figure 26. we calculate that the switch turns on at about $\mu mol/L$, and reach the intensity to turns on at about $10^{2.8} = 680\mu mol/L$ , which is exactly the same as our simulations. Note that we calculate the threshold at which the desired concentration reaches the turn-on threshold, which is 1 $\mu mol/L$. While at the point where the slope is greatest, that is, 100-150 $\mu mol/L$, which is the threshold our wet lab uses.

The Fit Curve of Fluorescence Intensity vs Inducer Strength of Wet Lab

We then fit the scatter plot and find that the inhibition coefficient is greater than ten times the unit strength strength, which proves the effeciency of our switch.

The Combination of Wet Lab Data and Model Simulation

The Decay of Relative Fluorescence Intensity

We conduct different treatments for the four groups of experiments, and set up repeated experiments, in which 1,2,3,4 represent the four groups of experiments, and A and B represent different treatment methods. We measure the decay of the relative fluorescence intensity in these 8 cases, and use R for linear fitting to obtain the relative linear coefficient of fluorescence intensity decay, and $R^2=-0.92, ~ p<2.2e-16$ showed significant linear correlation. The fitting equation was shown in the Figure 28. According to the relationship between fluorescence and product, the product decay coefficient in the model can be obtained.

The Fit curve of Decay of Relative Flurescence Intensity

Each Fit curve of Decay

Supplementary Model

Experiments of Two Threshold (Best Composite Part)

Experiment Part

https://2023.igem.wiki/ucas-china/parts

Modeling Part

Raw Data Visualization

We drew our experimental data and showed the error bar, and from the bar chart, we could see that there was indeed a threshold.

Bar Graph of PobR

Bar Graph of XylS

Data Fit

We perform logarithmic processing on the experimental data, and perform exponential fitting, so that the properties of the switch can be estimated at each concentration, which can be used to predict and establish the switch characteristic curve.

The Switch curve of PobR

The Switch curve of XylS

T-test and Get the Threshold with the lowest P-value

We conducted the paired T-test on the two adjacent groups respectively and obtained the P value, which was used as a representation of adjacency difference. By detecting the lowest P-value, we determined that our thresholds are located at 0.20-0.21 $mM$ and 100 $\mu M$, respectively.

The P Value of Adjacement groups of PobR

The P Value of Adjacement groups of XylS

Time Estimation of the Whole Cascade

The reaction rules

LuxI Synthesis