Introduction

Our team designed an integrated plan for rice sheath blight, including three steps: prevention of rice sheath blight by using Trichoderma atroviride to block the transmission of R. solani, early detection of rice sheath blight by electronic nose (ENose), precise detection by LAMP-LFD and biological control of Rhizoctonia solani by RNAi.

First, we would like to understand the process by which our product TACE releases spores of Trichoderma atroviride and propose suggestions for TACE improvement.

Second, a timer suicide switch model of Trichoderma atroviride was developed to simulate the dynamic change of the oscillator.

Thirdly, a rice sheath blight detection algorithm was established using machine learning methods. The response data of 6 gas sensors were collected by the ENose data acquisition equipment for processing. Then we used machine learning algorithms to classify the data, and achieve high accuracy.

Fourthly, a dynamic differential equation model was established to describe the dynamic process of shRNA production in E.coli and silencing of target genes in R.solani by using shRNA for biological control of R. solani.

Spore release model

The spores of Trichoderma atroviride can be released from the TACE in the water, and the concentration of Trichoderma atroviride spores in the water is related to the killing effect of Trichoderma atroviride on the R.solani, so we established a Trichoderma atroviride spore release model, and made the following assumptions for the model:

• Assumption 1: The effect of the interaction between spores can be ignored. Because spores have certain life activities, they may interfere with the release of nearby spores during the release process, but the interaction between spores is very weak, that is, spores can be considered as a single spore and do not form clump-like particles with neighboring spores.
• Assumption 2: The dynamic equilibrium can be achieved during the process of TACE absorption of spores. The TACE can absorb spores of Trichoderma atroviride in the solution, and can reach a dynamic equilibrium state after enough time.
• Assumption 3: Spores in the TACE cannot be fully released. When TACE releases spores in the solution, it will be affected by a variety of resistance, the most important of which is that in the process of releasing spores, the shape of TACE causes the void to change due to the release of spores, so that some of the original entry and exit gaps are blocked.
• Assumption 4: Trichoderma atroviride spores are evenly distributed in the solution. Trichoderma atroviride spores cannot be dissolved in the solution, and its own density is greater than the density of water, so the spores will precipitate at the bottom. In order to simplify the model, Trichoderma atroviride spores can be considered as evenly distributed in the solution.

The notations used in our release model are illustrated in Table 1.

Table 1. Notation

The absorbance `A` and the spore concentration `C` are consistent with Lambert-Beer law of light absorption:

$$A=abC \ (1.1)$$

Equation of spore concentration `C` and absorbance `A` can be obtained from Equation (1.1) :

$$C=\frac{A}{ab} \ (1.2)$$

After TACE reaches saturation, the spore release rate in the water satisfies the following rate equation:

$$\frac{dC}{dt}=\lambda e^{-kt} \ (1.3)$$

When `t_0=0 ` , the spores of Trichoderma atroviride were not released from TACE, so the concentration of Trichoderma atroviride spores in the solution was `C_0=0`. When `t\to\infty`, most of the spores in TACE have been released. There are still a small amount of spores that cannot be released due to changes in the internal structure of TACE, and the effective release rate is `\eta`. Then the concentration of Trichoderma atroviride spores in the solution at this time is `C_{\infty}=\eta C_{eq}`.

$$\begin{cases} \frac{dC}{dt}=\lambda e^{-kt}\\ C(0)=0\\ C(+\infty)=\eta C_{eq} \end{cases} \ (1.4)$$

Solve the differential equation (1.4) and substitute the above special solution into it, the following can be obtained:

$$\ln{(1-\frac{C}{\eta C_{eq}})=-kt} \ (1.5)$$

By sorting out Equation(1.5), the concentration `C` and time `t` of Trichoderma atroviride spore released by TACE swelling powder were obtained as follows:

$$C(t)=\eta C_{eq}(1-e^{-kt}) \ (1.6)$$

Equation (1.2) was substituted into equation (1.6) to obtain the relationship between the absorbance `A_{550}`of Trichoderma atroviride spore solution at 550 nm and time in the process of releasing Trichoderma atroviride spore from TACE:

$$A_{550}(t) = \eta A_{eq}(1-e^{-kt}) \ (1.7)$$

We measured the absorbance values of spore solutions of Trichoderma atroviride several times according to different release times, and fitted the release curve with experimental data, as shown in Figure 2.

Figure 2. Experimental data and fitting curves of TACE V1 (Left: 40 mesh, Right: 80 mesh) The fitting parameters and equations are shown in Table 2.

Table 2. Fitting parameters and fitting equations.

As can be seen from Figure 2, at 100 minutes, the absorbance value of the solution after spore release by absorbing 40 mesh TACE is `A_{550}(100)=1.309`, and that of the solution after spore release by absorbing 80 mesh TACE is `A_{550}(100)=1.452`. Most of the spores were released in 100 minutes by the two different sizes of TACE. At the same time, the amount of spores released by the absorption 80 mesh TACE was more than 40 mesh TACE, indicating that the spore release performance of absorption 80 mesh TACE was better than 40 mesh TACE. From the `\eta A_{eq}` column in Table 2, maximum spore release of 80 mesh TACE was higher than that for 40 mesh TACE when `t\to\infty`.

From the above analysis, it can be seen that the first generation of TACE products released Trichoderma atroviride spore for a short time, which is not conducive to the release of Trichoderma atroviride spore in the rice field for a long period of time. On the other hand, its total release is in a low amount , which is not conducive to the formation of local high concentration areas at the water surface. Therefore, we iterated on TACE product, replaced the new wrapping material, and made the product into pill shape.

We obtained spore release concentration data for TACE V2 with a diameter of 3 mm and 6 mm from wetlab, and then we fitted the experimental data using the spore release model (Figure 3). According to the fitting curve, although the spore release of the 3mm product is less than the 6mm product, it can achieve the optimal release in a short time. We analyzed the results with the wet lab members and decided to change the ratio of TACE V2.

Figure 3. Experimental data and fitting curves of TACE V2(Left: 3mm, Right: 6mm).

Table 3. Fitting parameters and fitting equations.

ShB detection model

The ShB detection model processes the gas sensor data collected by ENose, and uses pre-processing data as the input of the neural network to classify whether rice suffer from rice sheath blight(ShB). The main modules include data preprocessing module, data conversion module, classification algorithm module and result output module. The flowchart is shown in Figure 4.

Figure 4. Flowchart of ENose ShB detection model.

The data processing module is the basis of all the following modules. First, using linear interpolation to replace the missing values. Then, the Kalman filter is used to filter the sensors data, which is helpful to reduce the error of sensors. Set `i(i=1,2,\cdots ,6)` for the gas sensor id, and their observations are given by `z_{k_{i}}`, then the data is filtered using Kalman filter in equation(2.1).

After filtering the raw data using the Kalman filter, the time series data is shown in Figure 5.It can be seen in Figure 5, the noise in the time series data is reduced.

Figure 5. Comparison of time series data before and after Kalman filtering.

To make full use of the gas sensors data, we use sliding window algorithm to cut off the data. Set `X` for the sensors data after Kalman filtering, sliding window size is `\alpha`, and sliding step size is `\beta`，then the number of sliping window is `m=\frac{n-\alpha}{\beta}`, the sensors time series length is`n`, then after the sliding window is cut off, the data is `X`, which is shown in equation(2.2).

The data conversion module uses Gramian Angular Field(GAF) algorithm to encode time series as images. GAF encodes time series into images by polar coordinates based matrix and it can preserve absolute temporal correlation. The time series `x'_{j}` after filtering is first normalized to between 0 and 1, which is defined in equation (2.3).

Then, angular cosine and the time stamp are used to encode the rescaled data into polar coordinates. We use the Gramian Angular Summation Field (GASF) to generate images, and the GASF is defined in Equation(2.4).

The flow chart of GASF algorithm is shown in Figure 6.

Figure 6. The process of encoding time series as images using GAF algorithm.

After converting the data, 8000 gas datasets (Clean air: 1600; Healthy rice: 1600; Rice with rice sheath blight: 1600; Rice with blast: 1600) is randomly shuffled and then an 80:20 split is performed to construct the training set and test set. We perform data augmentation (such as random rotation, random horizontal flip) on the training set, and then normalize the data set to train ResNet18. After training 30 epochs, the training accuracy of our model is 99.38%, and the test accuracy is 99.87%(Figure 7).

Figure 7. Confusion matrix of rice disease classification. (0=Clean air, 1=Healthy rice, 2=Rice with rice sheath blight, 3=Rice with blast).

According to the test results, our algorithm can not only distinguish whether rice has ShB, but also distinguish ShB and rice blast. Nevertheless, we believe that our detection model has some shortcomings. For example, environmental data such as temperature and humidity are not put into the network for training, but rice diseases are often associated with high temperature and humidity.

RNA Silence model

RNA interference is a process that can specifically and selectively disrupt the expression of target genes. We developed a dynamic differential equation to describe this process. The production and silencing of shRNA were occurred in E.coli and R.solani, respectively. For the production process, the plasmid with shRNA sequences was inserted into the cytoplasm of Escherichia coli and transcribed to obtain Primary shRNA (pri-shRNA). Pri-shRNA was treated with Drosha/DGCR8 complex to form precursor shRNA(pre-shRNA), which was extracted to obtain pre-shRNA. For the silencing process, R.solani extracellular pre-shRNA is transported into the cytoplasm, loaded onto the Dicer/TRBP/PACT complex, and subsequently further processed into mature shRNA. Mature shRNA is associated with the Argonaute protein containing RISC, which mediates the cleavage and degradation of the target gene mRNA and functions as RNA interference.

Figure 8. Flowchart of shRNA production and silencing process.

To establish the RNA silencing model, we make the following assumptions for the model:

• Assumption 1: The gene copy number is constant at 40. In general, the transcription rate of shRNA plasmids in Escherichia coli is affected by a variety of elements, and it is difficult for us to quantitatively analyze due to many factors such as experimental environment. To simplify the model, we assume that the gene copy number is constant.
• Assumption 2: The enzymatic catalytic processes in the model all satisfy Michaelis-Menten equation. Michaelis-Menten equation are developed based on the kinetic theory of enzymatic reactions, and making assumptions about the enzymatic reactions involved here helps us to focus on the kinetic properties of the enzyme.
• Assumption 3: Not all pre-shRNA in E. coli cytoplasm can enter R.solani. On the one hand, there is a loss in extracting pre-shRNA from Escherichia coli. On the other hand, the pre-shRNA could not all be absorbed by R.solani.
• Assumption 4: The total available amount of Ago&RISC is assumed to be a constant. During silence, Ago&RISC exists in three forms: Ago&RISC, Ago&RISC combining shRNA, Ago&RISC combining shRNA and target mRNA. Ago&RISC combining shRNA and target mRNA can be decomposed in order to downregulate target mRNA. We assume that the total available amount of Ago&RISC is 100.

In the production process, the biochemical reactions involved are as follows:

First, the copy number of shRNA on the plasmid vector was `c`indicates that pri-shRNA is obtained by transcription at a constant rate `\alpha`, with `d_{pri\_shRNA}`denotes the degradation rate of pri-shRNA. The kinetic constants of Michaelis-Menten equation of Drosha/DGCR8 are `V_{Max,DD}`and `K_{m,DD}`, respectively, then the differential equation for the production of pri-shRNA is as follows:

pri-shRNA was loaded onto Drosha/DGCR8 compound and processed to form pre-shRNA. Using `d_{pre\_shRNA}` to represent the degradation rate of pre-shRNA, the differential equation for the production of pre-shRNA is as follows:

For the silencing process, the biochemical reactions involved are as follows:

After entering R.solani, the pre-shRNA binds to the Dicer/TRBP/PACT compound and is further processed into ds-siRNA. The percentage of pre-shRNA produced by E.coli that entered R.solani was `\eta`, and the Michaelis-Menten equation kinetic constants of Dicer/TRBP/PACT compound were `V_{Max,Dicer}` and `K_{m,Dicer}` , respectively. The degradation rate of pre-shRNA into R.solani was `d_{Pre\_shRNA_R}`. The differential equation is as follows:

pre-shRNA is processed to get ds-siRNA, and ds-siRNA binds to RISC is then further processed to give siRNA, and the RISC·siRNA compound mediates mRNA cleavage and degradation. We denote the degradation rate of ds-siRNA by `v_{b,RISC}`and `v_{da,RISC}` respectively. We can get the following differential equation:

We assume that the total available amount of RISC is a constant `RISC+RISC \cdot siRNA+RISC\cdot siRNA \cdot mRNA=K_{total}`, where, Only `RISC\cdot siRNA \cdot mRNA` compound can be degraded, and the rate of mRNA degradation is denoted by `d_{mRNA}`, and the rate of compound degradation is denoted by `d_{Rsm}`. We can get:

The parameters used in the model are as follows:

We used Python to simulate the model, and it can be seen from the simulation results (Figure 9) that the mRNA concentration in the simulation curve using RNAi reached the maximum at about 1 day, after which the mRNA concentration decreased and finally stabilized at 47.06 molec. The simulation curve without RNAi showed that the mRNA concentration would eventually stabilize at 67.51 molec. According to the inhibition ratio, as the siRNA complex degrades mRNA, the inhibition rate of RNAi also increases to 30.29%. In conclusion, from the simulation results, the RNA silencing model well describes the production and action process of shRNA, and most importantly, it can guide the design of RNAi wet experiments.

Figure 9. Simulation results of RNA silencing model.

Oscillator model

In order to prevent engineering bacteria from escaping into the environment, we designed a timed suicide switch (Figure 10) reference to Goodwin Oscillator, but the relevant properties need to be analyzed by model. Here, we hope to verify the concept of timed suicide switch by establishing a suicide switch model of E.coli, and provide theoretical guidance for practical experiments of this type of suicide switch.

Figure 10. Oscillator design(Left: Oscillator; Right: Effector).

In order to establish the oscillator model, the following assumptions are necessary for the model:

• Assumption 1: The amount of DNA in the same cell is constant and the physiological behavior is same. The DNA in the same cell can be expressed at the corresponding gene fragment.
• Assumption 2: Transcription and translation processes are carried out under saturation conditions. In the process of transcription and translation, polymerases, ribosomes, amino acids and nucleotides are present in large amounts.
• Assumption 3: Degradation of protein and mRNA as well as reactive degradation. That is, proteins and mRNA can be degraded directly without intermediate products.
• Assumption 4: Transcription rate can be modelled by Hill equation. That is because the transcription rate of each gene is determined by the concentration of its protein product, and the repressor protein binds to the regulatory region of the gene faster than transcription and translation.
• Assumption 5: The Hill coefficient (`n`) approximates the number of cooperative ligand binding sites on the receptor, and ligand molecules bind to a receptor simultaneously.
• Assumption 6: The translation rate of each gene is equal, which is `K_2`.

The expression processes of the 3 genes lacI,tetR and cI are described as follows:

Where `i=TetR,LacI,cI`, concentration of `DNA_i` is `[DNA_i]`, concentration of `mRNA_i`is `[mRNA_i]`, concentration of protein expression product is `[p_i]`, `k_1`is transcription rate of `mRNA_i`, `k_2` is translation rate of `mRNA_i`, `d_1`is the rate of `mRNA_i` decay, and `d_2` is the rate of `p_i` decay.

For the transcription process, the ODE for transcription rate is described as follows:

The transcription rate `K_1` in equation (4.2) can be modeled using Hill Equation because the transcription rate of `mRNA_i`depends on the protein `p_i `concentration. There are 3 different states of `p_i`, which are the unbound state `(P)`, dimer state when unbound `(P_2)` and the binding state to the operon receptor `(P_2 * O)`. The quantitative relationship between the 3 states is `P\gt\gt P_2\gt\gt P_2*O `.

Therefore, set `\gamma` is maximal repressor binding rate, `\gamma_0` is leakiness of repressor, `K_b` is the number of repressor molecules required to reduce expression by 50%, and `n` is Hill coefficient, `k_1` is described by the Hill equation：

To carry out dimension reduction and simplify the analysis of this system, we normalize the ODEs:

• 1. We write `p_i` in units of `K_b^{-1/n}`, write `m_i` in units of `\frac{d_2}{TK_b^{\frac{1}{n}}}`, where `T` is transcription rate, and write time in units of `\frac{1}{d_1}` , then the rewritten ODEs is: Where `\alpha=\frac{\gamma K_b^{\frac{1}{n}}T}{d_1 d_2}` , `\alpha_0=\frac{\gamma_0 K_b^{\frac{1}{n}}T}{d_1 d_2}` .
• 2. We take `\lambda` as a free parameter, and `p_i`can be rewritten as `p_i^{'}=\frac{p_i}{\lambda}` , define `\beta=d_2,\ d_2=\frac{k_2}{\lambda}, \lambda=\frac{k_2}{d_2}` , redefine the ODE:

The ODEs of the oscillator are defined as:

We set the oscillator model parameters as:

First, we set the initial number of mRNA molecules and the number of protein molecules to be both 0. The oscillator model was simulated for 1000 minutes (Figure 11). Note that only the curve representing the number of cI protein molecules is visible, as all plotted concentrations are identical and overlap. From the simulation results, it can be seen that the peak value is 288 molecules at 8 minutes, and the steady state is reached at 28 minutes, where the number of protein molecules is maintained at 210 molecules.

Figure 11. Oscillator simulation result.

We set the initial mRNA molecules and protein molecules to `m_{lacI}(t = 0)=m_{tetR}(t = 0) = m_{cI}(t = 0)=0`, `p_{tetR} = 10, p_{lacI} = p_{cI} = 0 `, The oscillator model was then simulated for 1000 minutes, and the simulation results are shown in Figure 12. It is clear that after a period of time in the oscillator system, the system reaches a steady state with the peak value is 3000 molecules of the 3 proteins and the peak-to-peak period is 175 min.

Figure 12. Results of 1000 min simulation of oscillator model

In order to change the peak amplitude of the oscillator model and the peak when the steady state is reached, we also explored the effect of different initial conditions on the oscillator model. In addition, we used 3D representation that plotted the simulation of the oscillator system for 1000 minutes (Figure 13), to fully understand the behavior of the oscillator system. Figure 13(a) shows that the system will rapidly reach the steady state under different initial protein numbers for 1000 minutes, unless all concentrations are the same, the system will begin to approach the limit cycle, and the steady state simulation is shown in red. Figure 13(b) shows the changes in TetR and LacI protein concentrations during the simulation. Figure 13(c). shows the curve of the region near the steady state near the limit cycle of Figure 13(b). Figure 13(d) shows that when the oscillator system is greater than 800 minutes, except for the steady state (indicated in red), it is related to the limit cycle. Figure 13(e) shows that for more than 200 minutes, the closer the initial conditions are to the same, the longer it takes for the system to reach the limit cycle.

Figure 13. 3D representation of oscillator model

After exploring the behavior of the oscillator with different parameters, we explore the dynamics of the mazEF system. mazEF is a toxin-antitoxin module located on the Escherichia coli chromosome and that of some other bacteria, including pathogens. MazF is a stable toxin that is capable of exerting toxic effects to kill cells. MazE is an unstable antitoxin that binds to mazF to form mazEF to keep cells alive, and its cellular concentration drops more rapidly than that of MazF, leaving MazF to exert its toxic effect, leading to cell death. The oscillator is coupled to the mazEF system, and based on equation (4.6), equation (4.8) is added to describe the dynamic change process of mazEF.

Where all parameters are identical as outlined above(4.7) except for the parameters of mazEF. We set
`k_1=0.03 [min^{-1}],k_2=10.0 [min^{-1}],K=40,d_1 =0.25[min^{-1}],d_2=1.155\times10^{-2}[min^{-1}]` and `Kas=2.0\times10^{-3}[min^{-1}]`.

We set different initial parameters of mRNA molecules and protein molecules for 1000 minutes of simulation, as shown in Figure 14. The oscillator model simulates oscillations in the number of repressor, mazE and mazF protein molecules. Notably, the mazF minimum amplitude corresponds to the cI protein minimum amplitude, which is due to the fact that both are controlled by the TetR operon. In addition, it can be seen that the number of mazF protein molecules does not decrease to 0 like other genes, because mazF protein has a long half-life and is more stable in cells. In the simulation, it can be seen that after a period of time, the content of mazF exceeds the content of mazE, and the toxin in the cells gradually plays a role, thus causing cell death. For example, in Figure 14(b), the cells will accumulate toxin proteins in the cells at 437 minutes, 626 minutes, 735 minutes and 806 minutes, thus causing cell death.

Figure 14. Simulation results of oscillator model under different initial conditions

References

[1] He J, Zhong C, Mi J. Modeling of drug release from bioerodible polymer matrices[J]. Drug Delivery, 2005, 12(5): 251-259.
[2] Goutelle S, Maurin M, Rougier F, et al. The Hill equation: a review of its capabilities in pharmacological modelling[J]. Fundamental & clinical pharmacology, 2008, 22(6): 633-648.
[3] Ainurofiq A, Choiri S. Drug release model and kinetics of natural polymers-based sustained release tablet[J]. Latin Am J Pharm, 2015, 34: 1328-1337.
[4] Zhang W, Liu T, Brown A, et al. The Use of Electronic Nose for the Classification of Blended and Single Malt Scotch Whisky[J]. IEEE Sensors Journal, 2022, 22(7): 7015-7021.
[5] Qi P F, Meng Q H, Zeng M. A CNN-based simplified data processing method for electronic noses[C]//2017 ISOCS/IEEE International Symposium on Olfaction and Electronic Nose (ISOEN). IEEE, 2017: 1-3.
[6] Liu T, Zhang W, Li J, et al. A Multiscale Wavelet Kernel Regularization-Based Feature Extraction Method for Electronic Nose[J]. IEEE Transactions on Systems, Man, and Cybernetics: Systems, 2022.
[7] Zhang W, Liu T, Ye L, et al. A novel data pre-processing method for odour detection and identification system[J]. Sensors and Actuators A: Physical, 2019, 287: 113-120.
[8] Cui H, Dong X, Shang K. An Improved Method for Long-term Drift Compensation of Electronic Nose with Batch Control[C]//2022 IEEE 12th International Conference on Electronics Information and Emergency Communication (ICEIEC). IEEE, 2022: 145-148.
[9] Borowik P, Adamowicz L, Tarakowski R, et al. Development of a Low-Cost Electronic Nose for Detection of Pathogenic Fungi and Applying It to Fusarium oxysporum and Rhizoctonia solani[J]. Sensors, 2021, 21(17): 5868.
[10] Eytan Zlotorynski. Insights into the kinetics of microrna biogenesis and turnover.NatureReviews Molecular Cell Biology, 2019.
[11] Rafał Krzysztoń, Daniel Woschée, Anita Reiser, Gerlinde Schwake, Helmut H. Strey, andJoachim O. Rädler. Single-cell kinetics of sirna-mediated mrna degradation.Nanomedicine:Nanotechnology, Biology and Medicine, 21:102077, 2019.
[12] Jiajia Guo, Jianjun Wang, and Genxing Xu. Pharmacokinetic perspective and the deliveryof sirna in vivo.Pharmaceutical and Clinical Research, 18(4):363–365, 2010.
[13] Shlomi Reuveni, Isaac Meilijson, Martin Kupiec, Eytan Ruppin, and Tamir Tuller. Genome-scale analysis of translation elongation with a ribosome flow model.PLOS ComputationalBiology, 7(9):1–18, 09 2011.
[14] Heather A Ferguson, Jennifer F Kugel, and James A Goodrich. Kinetic and mechanisticanalysis of the rna polymerase ii transcription reaction at the human interleukin-2 pro-moter11edited by k. yamamoto.Journal of Molecular Biology, 314(5):993–1006, 2001.
[15] Nancy J. Pokrywka and David S. Goldfarb. Nuclear export pathways of trna and 40sribosomes include both common and specific intermediates.Journal of Biological Chemistry,2021/10/02.
[16] D.N. Fouillen A. et al. Bouvette, J. Korkut. High-yield production of human dicer bytransfection of human hek293-ebna1 cells grown in suspension.BMC Biotechnol 18, 76(2018)., 2018.
[17] Antoine et al Baudrimont. Multiplexed gene control reveals rapid mrna turnover.Scienceadvances, 3, 2017.
[18] Elowitz M.B. and Leibier S. (2000) A synthetic oscillatory network of transcriptional regulators. Nature, volume 403(6767); pp. 335–338. ISSN 00280836.
[19] Gonze D. and Ruoff P. (2020) The Goodwin Oscillator and its Legacy. Acta Biotheoretica; pp. 1–18. ISSN 15728358.
[20] Danino T., Mondrag´on-Palomino O., Tsimring L., et al. (2010) A synchronized quorum of genetic clocks. Nature, volume 463(7279); pp. 326–330. ISSN 00280836.
[21] M¨uller S., Hofbauer J., Endler L., et al. (2006) A generalized model of the repressilator. Journal of Mathematical Biology, volume 53(6); pp. 905–937. ISSN 03036812.
[22] Gonze D. and Abou-Jaoud´e W. (2013) The Goodwin Model: Behind the Hill Function. PLoS ONE, volume 8(8); p. 69 573. ISSN 19326203.
[23] Andersen J.B., Sternberg C., Poulsen L.K., et al. (1998) New unstable variants of green fluorescent protein for studies of transient gene expression in bacteria. Applied and Environmental Microbiology, volume 64(6); pp. 2240–2246. ISSN 00992240.
[24] Purcell O., Savery N.J., Grierson C.S., et al. (2010) A comparative analysis of synthetic genetic oscillators. Journal of The Royal Society Interface, volume 7(52); pp. 1503–1524. ISSN 1742-5689.
[25] Strelkowa N. and Barahona M. (2010) Switchable genetic oscillator operating in quasi-stable mode. Journal of the Royal Society Interface, volume 7(48); pp. 1071–1082. ISSN 17425662.
[26] Potvin-Trottier L., Lord N.D., Vinnicombe G., et al. (2016) Synchronous long-term oscillations in a synthetic gene circuit. Nature, volume 538(7626); pp. 514–517. ISSN 14764687.