Our team has built detailed and specific models of each module around our project design, but the modeling is not only a quantitative description of biological processes. We also hope to establish a close relationship between the modeling and the experiment as the experiment provides the modeling data to optimize parameters while the modeling gives some guidance to the experiment or solve some problems that the experiment cannot deal with.
In the smell-sensing module, after the bacteria cells sense the smell signal, it will produce another signal LuxI that initiate the memory module. The initial amount of LuxI was obtained from the experimental data and was used in the next module. In the memory module, in order to quantify the amount of taRNA, we constructed a series of kinetic models and an image. In the reproduction module, we performed kinetic analysis of IPTG-induced crRNA production and got the amount of crRNA over time and IPTG concentration. The crRNA along with the taRNA from the memory module form an AND gate, which is possible to be quantitatively characterized by the dynamics models of transcription and translation. In the erasure module, we conducted the kinetic analysis of the AiiO induced by ATc and the degradation of AHL by AiiO, aiming to obtain the amount of AHL over time and ATc concentration. In addition to studying the kinetic models of these four modules, we also constructed a genome-scale metabolic model for the entire organism to simulate its true metabolism and growth.
Each model is closely connected with the experiment.
The models we built above are expected to solve the four problems we encountered in the experiment:
· How to control the AND gate that produces benzyl alcohol according to people's wishes? (AND gate model)
· How to make our assay produce more benzyl alcohol? (Metabolic model)
· Is it possible to achieve the memory erasure? (Erasure model)
· Which appropriate methods should we choose to control leakage of the memory module? (Control leakage model)
And our model description will revolve around these four questions.
. AND Gate Model
The AND gate depends on benzyl alcohol and IPTG. The bacteria cells sense benzyl alcohol and convert it into the initial AHL, which will trigger memory module. The AHL accumulated by the memory module combines with the LuxR expressed constitutively to form a complex which activates taRNA. The activated crRNA is induced by IPTG. Therefore, to obtain the quantitative description of AND gate, we divided the problem into lots of small circuits as follows:
o The circuit of initial LuxI production ▶
o The circuit of AHL production ▶
o The circuit of LuxR constitutive expression ▶
o The circuit of complex formation ▶
o The circuit of LuxI production controlled by promoter luxPR ▶
o The circuit of gene taRNA activation controlled by promoter luxPR ▶
o The circuit of gene crRNA activation induced by IPTG ▶
o The circuit of taRNA-crRNA complex formation ▶
* The circuit of initial LuxI production back to AND gate model→
The LuxI produced by the sensing module will trigger the initiation of memory module. But how did we obtain the production of LuxI from the sensing module?
Since kinetic parameters of the sensing module are difficult to find, It is so hard for us to obtain the production of LuxI from sensing different benzyl alcohol concentrations by modeling. However, the experiment measured the Fluorescence/OD600 in sensing different concentrations of benzyl alcohol and we converted this data to the concentrations of LuxI, which will be passed to the next module.
|Concentrations of benzyl alcohol(g/L)||Fluorescence/OD600 (a.u.)||Concentrations of LuxI (nM)|
The quantitative conversion between the concentration of LuxI and Fluorescence/OD600 is obtained from the literature and the concentrations of LuxI in sensing different benzyl alcohol concentrations will be used in the AHL production .
* The circuit of AHL production back to AND gate model→
The reactions of AHL production catalyzed by LuxI are described as follows:
|KAHL||Synthesis rate of AHL by LuxI||min-1||0.04 |
|D||Diﬀusion rate of AHL through the cell membrane||min-1||2 |
|KLuxI||Degradation rate of LuxI||min-1||0.027 |
|d1||Degradation rate of AHL intracellular||min-1||0.057 |
|d2||Degradation rate of AHL in medium||min-1||0.04 |
|Vcell||Typical volume of E. coli.||µL/cell||1.1*10-9 |
|Vex||Typical volume of microﬂuidic device||µL||1*10-3 |
|N||The numbers of the bacteria||2*108|
We modeled reactions (r.1-r.5) with the following set of ODEs. The initial LuxI from the last circuit was used to quantitatively characterize the initial AHL.
According to our characterization on luxPR, the concentration of AHL must be higher than 10 nM for triggering the memory module. However, we can know from Figure 5 that the initial AHL concentration has been more than 10 nM without sensing benzyl alcohol. So the modeling predicted the experiment that it is necessary to control the leakage of sensing module.
* The circuit of LuxR constitutive expression back to AND gate model→
The detailed reactions are described as follows:
|KLuxR||Translation rate of mRNA LuxR||min-1||10 |
|kmLuxR||Degradation rate of mRNA LuxR||min-1||0.247 |
|KmLuxR||Constitutive transcription rate of mRNA LuxR||min-1||3.0050 [1,5,6]|
|kLuxR||Degradation rate of LuxR||min-1||0.2 |
The model predictions use parameter values from paper and the steady-state value of LuxR constitutive expression will be used in the circuit of complex formation.
* The circuit of complex formation back to AND gate model→
The reactions are described as follows:
|k-1||Dissociation rate of [LuxR-AHLin]||min-1||10 |
|kd1||Dissociation constant of [LuxR-AHLin]||nM||150.9581 |
|k-2||Dissociation rate of dimer [LuxR-AHLin]2||min-1||1 |
|kd2||Dissociation constant of [LuxR-AHLin]2||nM||30.1916 |
|kc||Degradation rate of [LuxR-AHLin]||min-1||0.156 |
|kc2||Degradation rate of [LuxR-AHLin]2||min-1||0.017 |
We achieved the quantification of the dimer by using a series of ODEs (f.4-f.7). The LuxR concentration we took is the steady-state value of the constitutive expression. The total amount of AHL is composed of the initial AHL from sensing module and the accumulated AHL from memory module. The complex forms a dimer, which will bind the regulatory regions of two downstream promoters luxPR. The activated luxPR can lead to the activation of gene taRNA and the production of LuxI.
* The circuit of LuxI production controlled by promoter luxPR back to AND gate model→
The reactions are described as follows:
|kmLuxI||Degradation rate of mRNALuxI||min-1||0.247 |
|KLuxI||Translation rate of mRNALuxI||min-1||3.09 |
|KmLuxI||Plasmid copy number times LuxI transcription rate||nM*min-1||23.3230 |
|k-3||Dissociation constant of [LuxR-AHLin]2 to the promoter luxPR||min-1||10 |
|kd3||Dissociation rate of [LuxR-AHLin]2 to the promoter luxPR||nM||200 |
|α||Basal expression of LuxI||0.01 |
We combined the quantification of the dimer from previous circuit, the circuit of LuxI catalyzing AHL and the reactions of this circuit to model the following set of ODEs:
A. The production of LuxI from the positive feedback circuit varies with time.
B. The production of AHL from the positive feedback circuit varies with time.
The accumulation of LuxI and AHL was caused by positive feedback and reached stability over time (Figure 12). This part of the model could not only predict the feasibility of the memory module, but also realize the quantification of various substances, which laid a foundation for the quantitative description of AND gate. In addition, the steady-state values of LuxI and AHL will be used in the erasure module.
* The circuit of gene taRNA activation controlled by promoter luxPR back to AND gate model→
The model of this circuit and the model of the previous circuit have the same process, except that the same lux box.
The previous circuits have realized the quantification of LuxI and AHL in the momory module. We can get the number of lux box that combined with the dimer and the amount of activated gene taRNA. The steady-state value of activated taRNA will be passed to the next module.
* The circuit of gene crRNA activation induced by IPTG back to AND gate model→
The reactions are described as follows:
|KmLacI||lacI transcription rate constant||min-1||0.247 |
|kmLacI||lacI mRNA degradation rate constant||min-1||0.462 |
|KLacI||LacI monomer translation rate constant||min-1||15 |
|kLacI||LacI monomer degradation rate constant||min-1||0.2 |
|K2I||LacI dimerization rate||nM-1*min-1||50 |
|k-2I||LacI dimer dissociation rate constant||min-1||10-3 |
|kI2||LacI dimer degradation rate constant||min-1||0.2 |
|KIR||Association rate constant for 1st derepression mechanism||nM-2*min-1||3*10-7 |
|k-IR||Dissociation rate constant for 1st derepression mechanism||min-1||12 |
|kI2R2||Repressor-inducer degradation constant||min-1||0.2 |
|KIO||Association rate constant for repression||nM-1*min-1||960 |
|k-IO||Dissociation rate constant for repression||min-1||2.4 |
|KIRO||Association rate constant for 2nd derepression mechanism||nM-2*min-1||3*10-7 |
|k-IRO||Dissociation rate constant for 2nd derepression mechanism||nM-2*min-1||4.8*103 |
|D1||IPTG passive diffusion constant||min-1||0.92 |
A. The concentration of activated gene crRNA over time in different IPTG concentrations and different promoters. Due to the little difference in steady-state value, we performed partial amplification to obtain Figure B, C. (The promoter strength of J23110, J23114 and J23113 is becoming stronger in turn)
B. The concentration of activated gene crRNA over time in different IPTG concentrations of promoter J23113.
C. The concentration of activated gene crRNA over time in different strength of promoters of 20nM IPTG.
We obtained that the variation of activated crRNA with IPTG concentrations is very little. The results show that we can't get the quantitative effect but the qualitative effect of activated crRNA induced by IPTG. To achieve the different activated crRNA concentrations, we tried different promoter strengths to express LacI. Then, we found that LacI has a certain effect on the amount of activated crRNA even if it is not appreciable. But it can obviously affect the AND gate because a few activated crRNA can result in significant influence, which we can summarize from the next circuit.
* The circuit of hairpin structure formation back to AND gate model→
DNA sequences crRNA (Dcr) and taRNA (Dta) are transcribed into crRNA (Rcr) and taRNA (Rta). The hairpin structure of crRNA was not only opened slowly, but also forms a complex (Ract) by hybridizing with taRNA, which was opened faster. The open of the hairpin structure can make the protein expressed. We defined the ρthon/off = protein expression in the presence of Dta/ protein expression in the absence of Dta to quantificationally describe the degree of AND gate .
The reactions are described as follows:
|rcrtx||the transcription rate of crRNA||min-1|
|rtatx||the transcription rate of taRNA||min-1|
|racttl||the translation rate of crRNA_taRNA||min-1|
|rcrtl||the translation rate of crRNA||min-1|
|kx||the rate of reaction x||min-1|
|Kx||the Michaelis−Menten constants of reaction x||nM|
We modeled reactions (r.27-r.31) with the following set of ODEs to quantitatively characterize AND gate.
|Ktx||the Michaelis− Menten constants of transcription||nM||4.2 |
We use MATLAB to describe the system of differential equations (f.20-f.26) and the simulation results are shown as the following.
The results above indicate that the increase of crRNA preferentially promotes the opening of the AND gate, which results in producing more proteins. We got the steady-state value of activated taRNA (5.4 nM) from the circuit of gene taRNA activation. This value was used in the formula (f.26) to obtain the quantification of AND gate varying with crRNA (Figure 18).
In our project, the amount of taRNA could reach the steady-state over time after the bacteria cells sensing the smell, so the ρthon/off of AND gate depends on the amount of activated crRNA induced by IPTG. However, the model of the previous circuit shows that the concentrations of IPTG have little effects on crRNA. So we can try different promoters to express LacI to achieve different effects of AND gate according to people's wishes.
By combining the models described above, we've solved the problem that how to control the AND gate to produce benzyl alcohol quantitatively.
1. Accumulating taRNA rather than crRNA in the positive feedback circuit will result in better quantitative effect, which provides tutorial advice to the design of our project. From Figure 14, we can infer that the positive feedback will reach a steady state as a constant which can't be varied by human's will. From Figure 18, we can figure out that the activation degree of the AND gate is more dependent on the amount of crRNA. On account of this, we should regard crRNA as the controllable variable while regard taRNA as the accumulated substance.
2. Through the model, we proposed that controlling the production of crRNA by IPTG induction is irrational (Figure 16). We suggested to vary the intensity of the promoter that expresses LacI or take other more quantitative inducing methods, in order to control the production of benzyl alcohol more quantitatively.
. Metabolic Model
Genome-scale metabolic models (GEMs) are becoming increasingly popular recently to predict metabolic fluxes for various systems-level metabolic studies. GEMs are the reconstruction of many organisms' entire metabolic networks and describe gene-protein-reaction (GPR) relationships. Applications of GEMs reconstructed for various organisms have received considerable success . In our project, we use the tool COnstraints-Based Reconstruction and Analysis for Python (COBRApy) and Escher to reconstruct genome-scale model of metabolic networks in E. coli to achieve our goals [10,11].
The followings are our goals of this model:
1. Visualize the related metabolic pathways of benzyl alcohol production, so that we can easily find possible methods we can take to produce more benzyl alcohol.
2. Simulate the metabolic fluxes of E. coli in different growth medium and guide the selection of additional nutrients in the medium during the experiment.
3. Simulate the amounts of benzyl alcohol production under different gene-knockout strains by building genome-scale models of metabolic networks of E. coli. This can give some guidance and verification to our experiments of reproducation module.
1. The visualization of benzyl alcohol-related metabolic pathways
We put the map in jpg format here and you can move the mouse to zoom in and view all benzyl alcohol-related metabolic pathways.
mouse over the image：
Figure 20. The red box labeled genes are what we want to knock out, and the purple box labeled metabolites are what we want to produce more.
From the map, we can easily know that these genes tyrB, aspC tyrA, trpE and tnaA knockout can make organism produce more phenylpyruvate, a substrate of benzyl alcohol production . We can also observe that these genes control the synthesis of some essential amino acids, from which we suspect that the growth of the bacteria will be affected if these genes are knocked out. We will verify this conclusion in the next step.
2. Simulating growth media
· What will happen to the growth of bacteria if tyrB, aspC, tyrA, trpE and tnaA are knocked out?
|Types of medium||LB||D||Y||W||WY||DY||WD||WDY|
The types of medium as follows:
LB: LB medium without adding other amino acids; D: LB medium adding Aspartic acid; Y: LB medium adding Tyrosine; W: LB medium adding Tryptophan; WY: LB medium adding Tyrosine and Tryptophan; DY: LB medium adding Tyrosine and Aspartic acid; WD: LB medium adding Tryptophan and Aspartic acid; WDY: LB medium adding Tyrosine, Aspartic acid Tryptophan.
As a result of the modeling (Table 8), at least all L-tyrosine, L-tryptophan, L-aspartic acid have to be added in the medium, even without one of them, the growth rate of the bacteria will be 0. This result guides the experimental group that at least three supplements, L-tyrosine, L-tryptophan, and L-aspartic acid are added to the medium.
3. Simulating Deletions
· Could the bacteria produce benzyl alcohol if there is no gene knockout?
· How to make our bacteria produce more benzyl alcohol?
Our model firstly predicted that the bacteria could not produce benzyl alcohol without gene knocked out, which was consistent with the results of the experiment. In addition, the model also predicted how to produce more benzyl alcohol. Details are as follows:
Objective value is the total flux of all reactions calculated by FBA and indicates the growth rate of the bacteria. Benzyl alcohol value solved by FVA is the maximum flux of benzyl alcohol transported to the extracellular region, which can reflect the amount of benzyl alcohol production. C: aspC knocked out; B: tyrB knocked out ; E: trpE knocked out; A: tyrA knocked out; BC: tyrB and aspC knocked out; BCA: tyrB , aspC and tyrA knocked out; BCE: tyrB , aspC and trpE knocked out; BCEA: tyrB , aspC, tyrA and trpE knocked out; BCEAT: tyrB , aspC tyrA, trpE and tnaA knocked out.
From left to right: line 1 is marker; mdlB+, mdlC+ and vcm17+ are the results of gene mdlB, mdlC and vcm17 after the bacteria were induced with 0.01mM IPTG for 3 and half hours; mdlB-, mdlC- and vcm17- are the results of gene mdlB, mdlC and vcm17 without the induction of IPTG under the same condition; mdlB’, mdlC’ and vcm17’ are blank controls (without templates in PCR ). The results show that whether induced or not, all three genes have been transcribed.
From the results of model (Figure 21), we can conclude that E. coli doesn’t produce benzyl alcohol if there is no gene knockout, while produces sufficient benzyl alcohol if tyrB, aspC and trpE are knocked out. Meanwhile, knockout of genes also affects the growth of the bacteria.
In our experiments, firstly the organism that have three genes vcm17, mdlB and mdlC but without genes knocked out. It had been verified by HPLC that no benzyl alcohol was reproduced. However, the result of reverse transcription proved that three genes were all transcribed (Figure 22). The results indicate that benzyl alcohol cannot be produced due to insufficient substrate without knocking the genes, which is consistent with the results of our model.
1. The visualization of benzyl alcohol-related metabolic pathways
We obtain the metabolic model of Escherichia coli str. K-12 substr. MG1655 (iJO1366.json) from BiGG database. According to the previous design introduction, in order to make E. coli produce benzyl alcohol, it is necessary to use phenylpyruvate as a substrate and add three enzyme genes to E. coli. In order to simplify the model, we added a simple reaction in the metabolic model, directly producing benzyl alcohol from phenylpyruvate, without showing the three enzymes in detail.
We used Escher, COBRA and COBRApy to visualize the metabolic model. In our project benzyl alcohol is produced with phenylpyruvate as a substrate, so how to produce more phenylpyruvate is an important issue. We first used COBRApy to add a simple reaction to produce benzyl alcohol from phenylpyruvate as a substrate in the downloaded metabolic model iJO1366. Then we visualize benzyl alcohol-related metabolic pathways using Escher's Edit function and export it as PNG. In addition, we can also save the file as html format, which contains various functions of Escher. Through Escher you can get the details of the metabolites, reactions and genes in the BiGG database by moving the mouse and you can also edit the maps. The html format files and the codes can be obtained and viewed ( click here).
2. Simulating growth media
Due to the knockout of certain genes, the necessary amino acid cannot be synthesized and the growth of the bacteria is 0. Therefore, we simulate the growth environment of the bacteria in the model, and simulate the change of nutrients in the medium by changing the metabolic flow of nutrients to predict the growth of bacteria.
In order to produce more benzyl alcohol, more phenylpyruvate, a substrate for the synthesis of benzyl alcohol should be firstly produced, so we need to knock out the genes involving in the metabolic pathway that consumes phenylpyruvate. From the results of the previous visualization, we obtained the genes that need to be knocked out, which is also consistent with the paper.
In our model, we use COBRApy simulate deletions to get the difference in the amount of benzyl alcohol production. Simulations using flux balance analysis (FBA) can be solved using Model.optimize(), by which we can gather values for all reactions and metabolites. While FBA will not always give unique solution, because multiple flux states can achieve the same optimum, lux variability analysis (FVA) finds the ranges of each metabolic flux at the optimum. So in our model，we use FVA to obtain the flux of a single reaction and the objective. The simulated environment of the bacteria in this part of our model we used temporarily were consistent with the literature. The fermentation medium contained the following components: D-glucose (20 g/L), L-tyrosine (0.2 g/L), L-tryptophan (0.1 g/L), L-aspartic acid (3 g/L).In the next section we will verify the effect of different media on bacterial growth.
By constructing the metabolic model, we’ve worked out the problem that how to produce more benzyl alcohol.
1. In order to produce enough benzyl alcohol, at least three genes, tyrB, aspC and trpE, should be knocked out. However, this approach will have negative impact on bacterial growth.
2. At least three kinds of amino acids, L-tyrosine, L-tryptophan and L-aspartic, should be added into culture media for normal bacterial growth.
. Erasure Model
To achieve the erasure of the memory of bacteria, we designed the erasure module induced by ATc in our project. But we didn’t know whether the degradation of AHL by AiiO could eliminate the positive feedback of quorum sensing. Therefore, the erasure model was built by us to predict the effect of this erasure design.
* The circuit of AiiO expression induced by ATc
The reactions described as follows:
|KmRNA||transcription rate constant of the promoter||s-1|
We modeled reactions (r.31-r.41) with the following set of ODEs to quantitatively characterize the AiiO expression induced by ATc.
|Ptet_copy||plasmid copy number||s-1||4|
|KX||disassociation rate constant of TetR and ATc||s-1||0.36 |
|Kd||disassociation rate constant of TetR and DNA||s-1||0.1 |
|β||Original (unrepressed) transcription rate constant of the promoter||s-1||0.023 |
|kmRNATetR||degradation rate constant of mRNATetR||s-1||0.009 |
|kTetR||degradation rate constant of TetR||s-1||0.631 |
|KTetR||degradation rate constant of TetR||s-1||235.5 |
|kAiiO||degradation rate constant of AiiO||s||0.00017 |
|KAiiO||translation rate constant of AiiO||s-1||125 |
|kATc||degradation rate constant of ATc||s-1||0.00077 |
|n||Hill coefficient||2 |
From the results above, we can get that the different concentrations of ATc will have a significant impact on the concentration of AiiO, which will degrade the AHL.
* The circuit of AHL degradation
We used the concentration of AHL that reaches the steady-state from the memory module as the initial value of the erasure module. The parameters of AiiO degradation we took according to the paper .
|UAiiO||Enzyme activity of AiiO protein||nmol/(mg*min)||27.8 |
From the results above, we can know that the AHL concentration can be lowered whatever the concentration of ATc is, which show that the memory can be erased. However, when the ATc concentration is 1nM, the degradation of AHL needs more time.
By constructing the erasure module, we concluded that the expression of AiiO driven by ATc can reach the erasure effect.
. Control Leakage Model
To handle the leakage of the positive feedback circuit well, we make three alternative plans, and the detailed circuit information is described in the previous design. But we don't know which plan is more appropriate, so we build models of the two plans in order to make a better choice.
A. The degradation of AHL by luxI.
B. We tag an LVA flag on the C terminal of LuxI,hoping to decompose LuxI more quickly.
The concentration of LuxI was considered as the leakage of luxPR. The change of AHL concentrations under LuxI and LuxI-LVA were compared by modeling.
|KLuxI-LVA||degradation of LuxI-LVA||0.189 ||min-1|
The results above indicate that when LVA flag is tagged on the C terminal of LuxI, the decrease of AHL concentration can be realized significantly. However, such a concentration of AHL can still trigger the positive feedback. According to our characterization on luxPR, the concentration of AHL must be lower than 10-1 nM for inhibiting the positive feedback.
* Plan Ⅱ:
In this circuit, we need to solve two following problem:
1. When no smell molecules are detected, which RBS should be chosen so that LuxI leaked by luxPR promoter cannot trigger the positive feedback;
2. After sensing smell molecules, which RBS should be chosen so that LuxI generated by luxPR can trigger the positive feedback.
The circuits related in this part have been introduced in AND gate model. The difference is that the expression parameter of the promoter is set to 0.01 to describe its leakage, in order to meet our needs.
1. Due to the strong degradation ability of AiiO, the model inferred that no matter how to change RBS, the concentration of AHL would decrease from 0, which meant the positive feedback couldn’t be triggered.
2. Under the condition of sensing benzyl alcohol, on behalf of the former model named “The circuit of initial AHL production”, we took 150nM as the initial value of AHL concentration. According to the model built on different RBSs, we made a heatmap as follows:
From the map we can tell that the concentration of AHL is more dependent on the RBS on the upstream of AiiO. Since only if the concentration of AHL has reached 102 nM can the positive feedback be triggered, we can take B0030, B0029 or B0035 as the RBS on the upstream of AiiO.
By constructing the control leakage model, we can figure out that Plan I cannot realize the leakage control effect. Plan II can realize the leakage control effect under some circumstances. However, if the RBSs are not chosen properly in the Plan II, excessive inhibition will be formed, resulting in failures in the initiation of the positive feedback when sensing smells.
 Verdaasdonk J S , Lawrimore J , Bloom K . Determining absolute protein numbers by quantitative fluorescence microscopy[J]. Methods in cell biology, 2014, 123C:347-365.
 Weber M , Buceta J . Dynamics of the quorum sensing switch: stochastic and non-stationary effects[J]. BMC Systems Biology, 2013, 7(1).
 Boada Y , Vignoni A , Pico J . Engineered control of genetic variability reveals interplay between quorum sensing, feedback regulation and biochemical noise[J]. ACS Synthetic Biology, 2017, 6(10).
 Bernhard O̸. Palsson, Goldenfeld N . An Introduction to Systems Biology: Design Principles of Biological Circuits and Systems Biology: Properties of Reconstructed Networks[J]. Physics Today, 2007, 60(6):63-64.
 Wu C H , Lee H C , Chen B S . Robust synthetic gene network design via library-based search method[J]. Bioinformatics, 2011, 27(19):2700-2706.
 Iverson S , Haddock T L , Beal J , et al. CIDAR MoClo: Improved MoClo Assembly Standard and New E. coli Part Library Enables Rapid Combinatorial Design for Synthetic and Traditional Biology[J]. ACS Synthetic Biology, 2015:151019092657002.
 Stamatakis M , Mantzaris N V . Comparison of Deterministic and Stochastic Models of the lac Operon Genetic Network[J]. Biophysical Journal, 2009, 96(3):887-906.
 Senoussi A, Lee Tin Wah J, Shimizu Y, et al. Quantitative Characterization of Translational Riboregulators Using an in Vitro Transcription–Translation System[J]. ACS synthetic biology, 2018, 7(5): 1269-1278.
 Zhang C, Hua Q. Applications of genome-scale metabolic models in biotechnology and systems medicine[J]. Frontiers in physiology, 2016, 6: 413.
 Ebrahim A, Lerman J A, Palsson B O, et al. COBRApy: constraints-based reconstruction and analysis for python[J]. BMC systems biology, 2013, 7(1): 74.
 King Z A, Dräger A, Ebrahim A, et al. Escher: a web application for building, sharing, and embedding data-rich visualizations of biological pathways[J]. PLoS computational biology, 2015, 11(8): e1004321.
 Pugh S, McKenna R, Halloum I, et al. Engineering Escherichia coli for renewable benzyl alcohol production[J]. Metabolic Engineering Communications, 2015, 2: 39-45.
 Sun Z, Ning Y, Liu L, et al. Metabolic engineering of the L-phenylalanine pathway in Escherichia coli for the production of S-or R-mandelic acid[J]. Microbial cell factories, 2011, 10(1): 71.
 Mei-Chao Z , Ya-Nan C , Bin Y , et al. Characteristics of quenching enzyme AiiO-AIO6 and its effect on Aeromonas hydrophila virulence factors expression[J]. Journal of Fisheries of China, 2011.