Team:OUC-China/Model
Overview
Our model is the conjunction with experiments tightly. First, in order to choose the 3WJ repressors which have higher GFP fold reduction, we developed a score function to rank potential 3WJ repressors. Through calculating by using the MATLAB, the function displays the best parameters linear regression, which provides a coefficient of determination R2 of 0.6393. Second, based on the experiment data, we established ODE model to simulate the whole process of RNA switch. Our experimental results and modeling results are of the same order of magnitude. And according to sensitivities analysis, we found the degradation rate plays an important role in the ON/OFF ratio. So, we designed 5' end hairpin structure to improve the stability of RNA. What’s more, we designed the missing type of RNA logic gates—IMPLY and XOR. We used NUPACK to analyze their performance and chose the well-characterized sequences. Finally, they were verified by experiments, which were basically consistent with the modeling results.
3WJ Scoring Function
In order to choose the 3WJ repressors which have higher GFP fold reduction, we first developed a score function to rank potential 3WJ repressors.
The thermodynamic parameters were calculated using a local implementation of the NUPACK 3.0 functions mfe, pfunc, and complexes with Mathews et al., 1999 energy parameters1. A series of different linear regressions were performed using the thermodynamic parameters and log10(GFP fold reduction) obtained from the repressor library. These regressions were calculated using the Matlab regress function, which implemented a multiple linear regression algorithm using least squares.
Then we learned several parameters showing stronger correlations with the experimental data, so we used them to develop the score function Sk. And the function displays the best parameters linear regression obtained from experimental measurements of GFP fold reduction, which provided a coefficient of determination R2 of 0.6393.
Figure 1. Linear regression between Sk and log（GFP fold reduction）
Using thermodynamic parameters and experimental data of 3WJ calculated the scoring function Sk (more details shows following). Sk=0.0856*deltaG_min_hpin - 0.0717*deltaG_min_comp+ 0.0419*deltaG_hpin + 0.0334* deltaG_comp- 0.0397* net_deltaG_comp + 0.0908* recon_dup2link_deltaG. The function displays the best parameters linear regression of log10(GFP fold reduction), which provided a coefficient of determination R2 of 0.6393.
A 3WJ repressor library was get from a previous study[1]. A set of easy-to-calculate thermodynamic parameters were defined as shown in Table 1. The parameters were calculated using the subsequences, or sequence windows, also specified in Table 1. These parameters can be broadly classified into four different categories.
1. MFE of RNA strands and critical RNA subsequences: The minimum free energy (MFE) of the RNA strands and subsequences was calculated to assess the strength of the repressor secondary structures and the binding between the trigger and switch RNAs.
2. Measures of binding for critical RNA domains: The free energies of duplex structures formed between binding domains and 45-nt minimal trigger sequences were calculated to assess the binding strengths of these interactions. The duplexes were simulated as two separate molecules that hybridize with one another.
3. Net reaction free energies: The net change in free energy starting from the separate trigger and switch RNA strands and ending with the trigger-switch complex.
4. Influence of coding sequence secondary structure on translation in the inactive state: The MFE for the switch RNA subsequence immediately after the binding site for the trigger RNA was used to assess the degree of translation inhibition caused by sequestering the RBS and start codon in a stem-loop. These parameters were also calculated over several switch RNA subsequences.
The thermodynamic parameters were calculated for 36 3WJ repressors from the library using a local implementation of the NUPACK 3.0 functions mfe, pfunc, and complexes with Mathews et al., 1999 energy parameters1. A series of different linear regressions were then performed using the thermodynamic parameters and log10(GFP fold reduction) obtained from the repressor library. These regressions were calculated using the Matlab regress function, which implemented a multiple linear regression algorithm using least squares. This analysis identified several parameters showing stronger correlations with the experimental data and suggested an extension of the approach to three parameter regressions.
To develop a scoring function to rank potential 3WJ repressors, all possible combinations of up to three different thermodynamic parameters were first computed. This search was limited to three parameters to evaluate a number of parameter combinations that could be computed reasonably quickly since the number of combinations increases exponentially with the number of fitting parameters. Regression coefficients from the top 7 linear regressions were extracted to generate an overall scoring function. The combined regression coefficients were generated using equation(1)：
where Bi is the combined regression coefficient for parameter i，Cij is the regression coefficient for parameter i in regression j, and N= 7 is the number of regressions combined. In cases where parameter i was not used in the linear regression, the regression coefficient was equal to zero. The scoring function for ranking designs was then of the following form(2):
Where Sk is the score computed for 3WJ repressor design, k and ΔGik is the free energy calculated for parameter i for 3WJ repressor design k. Repressors predicted to offer higher performance were thus given higher scores. The final equation for the scoring function, using the parameter names from Table 1, is:
This function displays the best parameters linear regression obtained from experimental measurements of GFP fold reduction, which provided a coefficient of determination R2 of 0.6393.
Table 1. The parameters used in calculating the scoring function
Reference
[1] Kim, J., Zhou, Y., Carlson, P.D. et al. De novo-designed translation-repressing riboregulators for multi-input cellular logic. Nat Chem Biol 15, 1173–1182 (2019). https://doi.org/10.1038/s41589-019-0388-1
ODE
Rate Of Activated Transcription
The rate of activated transcription from a regular gene depends on the promoter occupancy. P is an activator, then the rate is proportional to the occupancy:
[P] is the concentration of P, α is the maximal transcription rate and K is the dissociation constant of the binding event, is the half-saturating concentration for the transcription factor P.
The rate of transcription of the three promoters (pT7,pTet,pLux) can be calculated by substituting the experimental data. Here shows the result.
Table 1. the α and K using in the equation
|
|
Figure 2. The relation between transcription rate of three promoters and the concentration of inducer simulated by MATLAB
(A) The maximum transcription rate of pT7 is close to 3h-1. And when the inducer concentration was lower than 0.1mol/L, the transcription rate increased rapidly. (B) The maximum transcription rate of pTet is close to 2.8h-1 (shows in color blue). The maximum transcription rate of pLux is close to 2.6h-1(shows in color orange). And when the inducer concentration was lower than 0.1mg/mL, the transcription rate of both pTet and pLux increased rapidly.
ODE Of The Switch
In our work, our aim is to have higher ratio of ON/OFF so that the logic gate 0/1 is more distinguishable. Therefore, in order to understand which parameter plays a key role in our switch, we used the ODE model to simulate the whole process of it.
Reaction
We can describe the process of 3wj as follow:
① null -> 3WJon
② null -> Trigger
③ 3WJon + Trigger <-> 3WJoff
④ 3WJon -> GOI + 3WJon
⑤ Trigger -> null
⑥ GOI -> null
⑦ 3WJon -> null
⑧ 3WJoff -> null
ODES:
To simulate the dynamics of sfGFP, we use ordinary differential equations to model the reactions above. And ODEs are given as follows:
For readability, the complex symbol is simplified as:
Species, Symbols, and Parameters
Table 2. Species, symbols and parameters.
The parament we used in the ODEs is listed in the following table:
Table 3. The parament we used in the ODEs.
Simulation:
With the help of the Simbiology toolbox in Matlab，we simulated our 3wj for 12h, the result can be seen in Figure 3. The final ON/OFF ratio is about 4. The results of our experiment are of the same order of magnitude.
Figure 3. The dynamics of sfGFP by model prediction
Simulated the output of 3WJ for 12h by using MATLAB Simbiology toolbox. A) shows the state of ON，B) shows the state of OFF. The final ON/OFF ratio is about 4. The results of our experiment are of the same order of magnitude.
Sensitivities Analysis
In order to explore which parameter plays a crucial role in the translation level of our switch system, we carried out sensitivity analysis on the parameters obtained by the ordinary differential equations.
Figure 4. Sensitivity analysis of switch
The first one is the sensitivity analysis of 3wj repressor and the second one is the sensitivity analysis of toehold switch. In both results, the rate of transcription, the rate of translation and the rate of degradation play an important role in the output of our system.
From figure 4. we can see that except for the rate of transcription, the rate of degradation also plays an important role in the output of our system. So we decided to improve the stability of the RNA to promote its ON/OFF ratio. According to the literature, we have read, we know that we can add a hairpin structure at the 5’ end of RNA to achieve our goal. So we designed a few triggers with 5’ hairpin. However, due to the impact of coronavirus, we didn’t have enough time to test this design. We will continue to explore its effects in the future.
References
[1] Meyer, S., Chappell, J., Sankar, S., Chew, R. and Lucks, J.B. (2016), Improving fold activation of small transcription activating RNAs (STARs) with rational RNA engineering strategies. Biotechnol. Bioeng., 113: 216-225. doi:10.1002/bit.25693
[2]Carrier, T.A. and Keasling, J.D. (1999), Library of Synthetic 5′ Secondary Structures To Manipulate mRNA Stability in Escherichia coli. Biotechnol Progress, 15: 58-64. doi:10.1021/bp9801143
Sequence design
Improve Stability Of The Trigger
In order to improve the ON/OFF ratio, we tried to increase the stability of triggers. By referring to the literature, we learned that adding hairpin structure at the 5 'end of RNA could prevent the degradation of RNA, thus improving the stability of RNA[1]. Therefore, we tried to add the hairpin structure（Figure 5.A） found in the literature at the 5' end of triggers.
First, NUPACK and RNAfold were used to analyze the structure of the trigger with the 5' end hairpin, and then NUPACK was used to analyze the combination of the trigger with the 5' end hairpin and switches（Figure 5.B）.
Figure 5. The secondary structure of RNA that used in improving stability
(A) the secondary structure of 5' hairpin analyzed using NUPACK.(B) the secondary structure of 5' hairpin combining with trigger analyzed using NUPACK.
New Designed RNA Switch Logic Gate1—XOR
We wanted to develop various of logic gate using the RNA switches, but in the beginning, we don’t have an approach to design the XOR gate. So, we read a lot of literature and communicated with professors to find the solution. Finally, inspired by the NIMPLY gate [2], we designed the XOR gate.
We chose the trigger’s core sequence and added the nucleotide-binding domains at the triggers’ both ends to designed two triggers. When input one of these triggers, the switch can turn ON. And when input these two triggers at the same time, they can pair together and form a ring in the middle (Figure 6), so the switch will still be in OFF state.
The nucleotide-binding domains were designed using NUPACK. And the core sequences of trigger originate from previous literature.
Figure 6. The secondary structure of the pairing two triggers analyzed using NUPACK.
New Designed RNA Switch Logic Gate2—IMPLY
In order to enrich the types of RNA logic gates，we also designed the IMPLY gate. The IMPLY gate RNAs based on toehold switches and 3WJ repressors were generated by taking the core regulatory sequence of the toehold and the repressors running from the 5’ end of the binding domain through to the nucleotide immediately before the 21-nt linker sequence. Because 3WJ repressors’ core regulatory sequence had a length of 73 nt, spacers of 3n + 2, where n is a non-negative integer, were used to connect different toehold switches and 3WJ repressors together. Spacers of this length enabled successive repressor modules to remain in-frame through the full length of the gate RNA. In previous NAND gates, they chose 17-nt spacers and the results are good. So we also used 17-nt spacers inserted between toehold switch hairpin and 3WJ repressor hairpin. The spacers were designed using NUPACK to have single stranded secondary structures when flanked by the two hairpins (Figure 7). After created by NUPACK, we deleted those sequences that have the stop codons. The remaining sequences are tabulated and narrowed down further.
Then, we used RNA fold and RNA structure to test their structure and deleted some designed gates that have great structural changes in harpins. Finally, we designed 3 IMPLY gates.
Figure 7. The secondary structure of the IMPLY gates RNA analyzed using NUPACK, RNAfold and RNAstructure
The Subtractor
We used the NIMPLY gate to make up the subtractor. However, we only found one NIMPLY gate sequences in the literature [2]. In order to achieve the subtractor that can calculate different concentration range, we designed NIMPLY on the basis of that sequence.
The NIMPLY gate consists of a toehold switch and two input triggers that can pair together. The toehold and the trigger’s core sequences we used originated from the previous literature. Then we added the nucleotide-binding domains at both ends of the trigger’s core sequence. This sequence is named “trigger1”. The other one was designed for completely complementary pairing with trigger1 (Figure8) and not pairing with switch sequence.
These sequences were designed using NUPACK and were further screened using RNAfold and RNAstructure.
Figure 8. The secondary structure of the pairing two triggers analyzed using NUPACK.
Multiple-input
On the basis of single-input logic gates and two-input logic gates, we also designed the three-input logic gates and four-input logic gates. These enrich the types of our logic gates.
We combined the AND gate and 3WJ to realize the three-input logic gate. And we tested it with different input methods by using NUPACK. It showed as we expected.
Figure 9. The results of different input method to three-input gate analyzed by NUPACK
A)input trigger AND1 and trigger AND2 , both triggers are binding to the switch B）input trigger AND1, trigger AND2 and trigger 3wj, all triggers are binding to the switch C） input trigger 3wj, it is binding to the switch
We wanted to use two NIMPLY gates to compose a four-input gate, however, it didn’t work when we analyzed it in NUPACK. So we chose the OR gate and transformed the two toehold switch as a subtractor. And we tested different input method, it performed well in NUPACK.
Figure 10. The results of different input method to four-input gate analyzed by NUPACK
(A) input trigger A and trigger B, both triggers are binding to the switch. (B）input trigger A，trigger A is binding to the switch.(C)input trigger B, trigger B is binding to the switch. (D)input trigger A and trigger A*, A and A* are paired together. (E)input trigger B and trigger B*, B and B* are paired together.
References
[1] Meyer, S., Chappell, J., Sankar, S., Chew, R. and Lucks, J.B. (2016), Improving fold activation of small transcription activating RNAs (STARs) with rational RNA engineering strategies. Biotechnol. Bioeng., 113: 216-225. doi:10.1002/bit.25693
[2] Kim, J., Zhou, Y., Carlson, P. D., Teichmann, M., Chaudhary, S., Simmel, F. C., … Green, A. A. (2019). De novo-designed translation-repressing riboregulators for multi-input cellular logic. Nature Chemical Biology. doi:10.1038/s41589-019-0388-1
[3] Zadeh J N , Steenberg C D , Bois J S , et al. NUPACK: Analysis and design of nucleic acid systems[J]. Journal of Computational Chemistry, 2011, 32(1):170-173.