Linalool Production by MEV Pathway
Overview
Our Linalool production model highlights the production of linalool by the mevalonate pathway in our bacteria, in the presence of the quorum sensing circuit in it.
The mevalonate pathway uses Acetyl CoA and Aceto-acetyl CoA as the starting materials for production of linalool with the help of a set of enzymes catalysing the conversion of one intermediate to another.
The MEV pathway is given below:-
The mevalonate pathway is incorporated inside the bacteria with the help of the MEV plasmid expressing the enzymes involved in the pathway. The genetic circuit of the MEV plasmid is given below :-
Objective
We determine the concentration of linalool produced by the E.coli in the bacterial medium in presence of the quorum sensing circuit. The simulation is performed through a system of several ODEs to determine the linalool concentration above the shoe insole with time.
Model
We use simple enzyme kinetics to determine the production of Linalool starting from Acetyl-CoA and enumerate the concentrations of the intermediates. We thus generate a system of 12 ODEs governing the rate of production of each of the molecules and simulate them using Python to determine the concentration of the produced Linalool with time, both beneath the shoe insole and above the shoe insole.
Our ODE system is given below :-
NOTE: The parameters used for the quorum sensing model are the most suitable parameters initially set for obtaining the best desired result from the QS model, and not the parameters obtained after fitting our QS model with the Wetlab data.
Explanation of the equations :-
The bacterial population, “N’, used here is the population curve obtained in presence of the quorum sensing circuit inside the bacteria, to effectively quantify the amount of Linalool produced in the presence of quorum sensing.
Here, “E” is the concentration of each of the enzymes catalysing the reactions involved in the mevalonate pathway and produced by the MEV plasmid.
IMPORTANT ASSUMPTIONS :-
We have assumed here that all the rate of gene expression in the MEV plasmid for the production of the enzymes is the same for all the enzymes as all the enzymes are downstream of the same promoter.
We have assumed that most of the enzymes have the similar half lives, and thus, same rates of degradation. This assumption has been made to reduce the number of ODEs for the model and make it more compact.
We assume that linalool diffuses from inside the shoe insole out to the external environment above the insole. The diffusion is governed by just the concentration gradient of linalool above and below the insole for the ODE model.
“ACA” refers to acetyl-coA, which is already present inside the cell at a certain concentration, and thus, it has its own intrinsic rate of production.
“Li” is the linalool concentration inside the shoe insole, that is, concentration of linalool directly above the bacterial medium. It gradually diffuses out of the shoe insole into the external environment through the semipermeable membrane of the shoe insole at a diffusion constant of “D”. The rate of diffusion is directly proportional to the concentration gradient between the linalool inside the shoe insole “Li” and the linalool outside the shoe insole “Lo”.
“Lo” refers to the linalool that has diffused out of the shoe insole and into the external environment. Apart from diffusing out of the shoe insole, it also has a rate of evaporation, “v”, as linalool is a volatile compound and evaporates at room temperature.
Plots :-
We solve the ODE system using Python and simulate it to determine the concentration of Linalool produced with time, in the presence of the quorum sensing circuit :-
Explanation in Oscillatory Dynamics
Fig: Population Dynamics with time
Fig: Linalool production with time
We observe that the oscillation observed in the cell population is consistent with the concentration of linalool produced over every time slot of simulation, that is, when the oscillating bacterial population reaches its peak population, the oscillatory linalool concentration also achieves its peak due to more bacteria producing linalool at that point in time, and vice versa when the bacterial oscillation reaches its crest.
Stochastic Model
We built a stochastic model for the concentration of linalool produced by the bacteria having the quorum sensing circuit in it, by simulating the model using the Gillespie Algorithm, where we treat each reaction/occurrence as an event and determine which event takes place in each timestamp. We then add/subtract an entity from the current number of molecules of that entity for that timestamp, based upon how the number of molecules of that entity changes upon one occurrence of that particular event.
The reason for building the stochastic model is to study the effect of noise in our Linalool production circuit, and observe whether noise can enhance or diminish the amount of linalool produced by the bacteria, because any biological system in contact with the outside environment will be influenced by several external factors which can subject the circuit a variety of noise levels, and hence, studying any biological system in the presence of external noise levels become extremely important.
The linalool production model In the presence of quorum sensing consists of 16 more events as stated below, in addition to the 19 events considered while simulating the quorum sensing model :-
We simulate the events using the Gillespie Algorithm to obtain the stochastic linalool concentration curve produced by the bacteria with time. The plot is given below :-
Chemical Master Equation :-
We thus determine the Master Equation of our MEV circuit governing the probabilities of each state in the circuit.
P(n21,t) = P(Lo,t) External Linalool Master Equation,
which gives us the probability of having [Lo] concentration of linalool on the shoe insole at a certain time t.
To determine P(Lo,t), we run numerous realizations of our stochastic model using the Gillespie Algorithm to plot the Stationary Distribution of Linalool concentration over the time interval, 50 < t (in hours) < 550, for the bacterial population with quorum sensing and find out the probabilities of the amount of linalool produced being at a particular concentration during the time interval.
The Stationary Distribution of linalool concentration, within 50 < t (in hours) < 550, is given below :-
Observation :-
We observe that external noise doesn’t create much of an effect over the amount of linalool produced by the bacteria. It however makes the oscillatory dynamics in linalool production a bit more vigorous, making the oscillation more prominent that what is observed in the deterministic model. Except that, the amount of linalool produced is observed to be almost the same for both the stochastic and deterministic models, indicating that noise doesn’t have any effect on linalool production by our bacteria.
Parameter Space Analysis
We perform parameter space analysis for the parameters ‘ke’ and ‘de’ to determine the how the linalool production varies upon changing the parameters ‘ke’ and ‘de’ for examination of the optimum parameter space for our model.
ke – Promoter Strength
de – Degradation Rate of each Enzyme
We therefore plot a 10 X 10 array in Python, with the 2 axes representing the 2 parameters:-
X axis – enzyme degradation rate (de) ,
Y axis – promoter strength (ke)
With each grid indicating a certain parameter space of the form: (de , ke).
NOTE: The ke and de values used here are taken from the most suitable parameter space for obtaining the best desired result from the quorum sensing model
An Interesting Observation :-
We observe that the amount of linalool produced remains around the desired value of 1750 nM, which we have chosen as the optimum amount for linalool production, over a wide range of parameter space, thereby indicating that even large fluctuations in the degradation rate of the enzymes, as well as fluctuations in the promoter activity, shouldn’t have too much off an effect on the amount of linalool produced by the circuit.
Linalool Diffusion Model
We animate the diffusion of linalool produced by the bacteria living on the bacterial medium inside the shoe insole, through the semipermeable membrane of the insole into the outside environment of the insole, using the Gray-Scott diffusion model to demonstrate the entire diffusion process of linalool through the insole material in a better way other than simplistic ODEs.
The linalool diffusion model through the insole membrane follows the throupled Reaction-Diffusion equation system as given below :-
We took a 3 X 26 hexagonal grid and used the finite difference method to solve the coupled reaction-diffusion PDE system. The reason for choosing such a grid system is stated below :-
1) The bottom row in the grid system represents the space just above the bacterial medium, where the linalool produced directly by the bacteria, “Li”, gets accumulated.
2) The middle row represents the insole material/membrane, through where internal linalool produced by the bacteria, “Li”, travels to reach the outer environment above it.
3) The top row represents the space just above the shoe insole, that is, the outside environment, where the internal Linalool, “Li”, diffuses out to from inside the insole. This outer Linalool, “Lo”, also has its own rate of evaporation, “v”, by which is slowly gets lost into the environment.
This is therefore a cross sectional view of the shoe insole!
The model is simulated over a timespan of 7 < t (in hours) < 44 to generate the animation of the diffusion of linalool through the insole material with time. The animation is given below.
An Interesting Observation :-
Note that the linalool diffusion occurs from inside the insole out to the external environment above the insole. We observe that after a period of around 36 hours, the linalool concentration both inside and outside the insole reaches equilibrium, and proceeds to remain constant until the bacterial population onsets its death phase.
Conclusion
We observe, through in silico, that the bacteria produces linalool in an efficient manner. Due to the presence of the quorum sensing circuit in the bacteria, linalool gets produced in a sustained manner over a long period of time, due to enhanced population life of the bacteria, approximately over more than 25 days, rather than only getting produced for 10 days in case of a normal bacterial population.
We therefore verify our goal of producing linalool through the mevalonate pathway inside our bacteria, and produce it in a sustained manner with the help of quorum sensing.
References
https://2019.igem.org/Team:XJTU-CHINA
https://2021.igem.org/Team:FAFU-CHINA
https://chem.libretexts.org/Bookshelves/General_Chemistry/Book%3A_Structure_and_Reactivity_in_Organic_Biological_and_Inorganic_Chemistry_(Schaller)/IV%3A__Reactivity_in_Organic_Biological_and_Inorganic_Chemistry_2/02%
3A_Mathematical_Tools_in_Reaction_Kinetics/2.06%3A_Characterization-__The_Mathematics_Behind_Enzyme_Kinetics">
https://2023.igem.wiki/iiser-kolkata/
https://2021.igem.org/Team:HZAU-China/Model
Quorum Sensing for Population Regulation
Overview
Our quorum sensing model presents a novel approach towards controlling the cell population at a level way below the carrying capacity by blocking the production of a protein named DNAA, a key initiator molecule for cell division. This is in contrast with the traditional method of population regulation through quorum sensing, where a cell toxin gene is introduced to promote cell death once toxic metabolites accumulated due to cell death.
Our Quorum Sensing Circuit is given below:-
Objective
We determine the population curve for normal E.coli growth in a bacterial medium and then compare it to E.coli growth in presence of the Quorum Sensing circuit inside the it. The simulation is performed through a system of several ODEs to determine the population curve of the bacteria with time.
Model
Reference model :-
Our initial target is to determine a general bacterial growth curve without any transformation. We thus generate a system of 2 ODEs governing the rate of change of substrate (S) with time and population (N) with time.
For population dynamics determination, we have modified the Logistic Growth Equation to enhance the analysis of our model, by incorporating in a new death term which is shown and explained below eventually.
For substrate consumption, we have employed the Monod Substrate Model.
We initially set the most suitable parameters targeted towards obtaining the best desired result from the reference model. The value of carrying capacity K=700 is chosen as a general feed bacterial concentration, which is basically the CFU value scaled down by a particular factor for efficient model simulation.(These parameters are not fitted to the wetlab data. The fitted parameters to the wetlab data are given in the Results section.)
Growth Equation Modification :-
We know that any bacterial population passes through 4 phases, the lag phase, the log phase, the stationary phase, and the death phase. However, the shortcoming of the logistic growth equation is that it can never determine the death phase of the bacteria. However, for efficient prediction of the longevity of our shoe insole, we need to determine how long can the medium sustain the bacterial population, and hence, determination of the death phase in the bacterial population becomes a crucial factor for us.
We know that the death phase of a bacterial population in a medium arises mainly due to reduction in the substrate/food concentration in the medium with time, as well as due to the accumulation of toxic secondary metabolites in the medium released during and after a bacteria dies. Essentially, the concentration of toxic secondary metabolites in the medium increases with time. We thus substract the following term,
to the growth equation to demonstrate the bacterial death caused upon by reduction in substrate concentration and the accumulation of toxic secondary metabolites with time, to eventually bring in the death phase in the population curve.
🡪 This is the secondary metabolites accumulation term, which is demonstrated by the Hill function, to highlight the fact that the accumulation of secondary metabolites in the medium peaks up after a certain point in time, thereby bringing about the death phase in the population. This term is proportional to the population N because greater the population of the bacteria, greater will be the amount of secondary metabolites released into the medium.
🡪 This is the death term due to reduction in substrate concentration. It denotes that as the substrate concentration decreases, rate of bacterial death become greater, basic Hill function, and once the substrate concentration goes below a particular threshold, death becomes really rapid in the population, bringing about the death phase.
Since death induction due to the secondary metabolite accumulation event and due to the substrate reduction event are independent of each other, therefore, we multiply both the terms to generate the overall death rate of the bacteria, and determine when the death phase comes in the population.
Thus, our overall death term for death phase incorporation becomes :-
Hence, the normal logistic equation is :-
And we introduce the death term in the population to bring about the death phase by modifying it as such :-
Plots :-
We solve the ODE system using Python and simulate it to determine bacterial growth curve, as well as the substrate concentration curve, with time :-
Quorum Sensing Model :-
We generate a system of 9 ODEs to model the gene expression of our QS circuit by employing simple gene kinetics to determine the population (N) and substrate (S) curves with time in presence of the quorum sensing circuit.
We modify our newly developed growth equation accordingly to incorporate the change in the bacterial growth rate due the binding of the sense mRNA of DNAA gene with its antisense part, to block the production of the DNAA protein, to determine the effect of the antisense mRNA upon bacterial growth.
Our ODE system is :-
We initially set the most suitable parameters targeted towards obtaining the best desired result from the quorum sensing model. (These parameters are not fitted to the wetlab data. The fitted parameters to the wetlab data are given in the Results section.)
Modifications in the growth equation :-
In our initial logistic growth model,
The antisense mRNA of DNAA gene binds with its sense counterpart, thereby inhibiting the production of the DNAA protein. As the production of DNAA gets hindered, cell division also gets hindered, and as a consequence, the rate of population growth reduces. Thus, we need to modify the Growth term to incorporate the reduction in population growth due to the presence of the antisense mRNA.
We henceforth modify the growth term as such :-
The idea is that, once the antisense mRNA reaches above a threshold concentration, it binds to the sense mRNA of the DNAA gene, thereby inhibiting the production of the DNAA protein. Once this occurs, population growth also reduces, which is incorporated by the Hill function in the population growth term.
Plots :-
We solve the ODE system using Python and simulate it to determine bacterial growth curve, as well as the substrate concentration curve, with time in the presence of the quorum sensing circuit :-
Observation :-
We observe that our bacterial population oscillates (oscillation is due to the negative feedback loop in the QS circuit) at a cell concentration level of around 250, quite below the carrying capacity of the original population, which is 700. As a result, due to a lower population concentration, the substrate consumption, as well as the accumulation of toxic secondary metabolites in the medium, slows down in case of the bacterial population with quorum sensing and thus, the population life of the bacteria with the quorum sensing circuit gets extended upto a period of around 25 days, quite longer than population life of the original bacterial population, which is around 10 days.
We are thus able to forecast the increase of bacterial population life in the medium by incorporating the quorum sensing circuit inside the bacteria.
Explanation of Oscillatory Dynamics :-
We simulate the bacterial population up to 4 different time slots to demonstrate each event occurring in that time slot.
Fig 1: t = 13.5 sec
Antisense concentration below threshold, so population increasing. Cell cycle not arrested. TetR does not inhibit production of LuxI to produce AHL.
Fig 2: t = 14.5 sec
Antisense concentration crosses threshold, so cell cycle gets arrested. Population growth stops, and starts to decrease due to only cell death occurring. TetR starts to inhibit production of AHL by hindering LuxI production.
Fig 3: t = 17.5 sec
Antisense concentration goes below threshold due to tetR inhibiting the production of AHL by hindering LuxI production, as AHL is the inducing molecule for the pLux promoter synthesizing the antisense mRNA. So cell cycle arrest gets removed. Population growth restarts. Now, the inhibition of AHL production gets removed as tetR concentration also reduces with reduction in the concentration of the antisense as both of them are expressed by the pLux promoter. AHL thus again starts inducing the pLux promoter to produce the antisense mRNA.
Fig 4: t = 20 sec
Antisense concentration again increases and reaches beyond the threshold, thereby again arresting the cell cycle and stopping population growth. Thus, population again comes down.
Fig 5: t = t` sec
This entire cycle gets repeated over time, thereby producing the oscillatory dynamics in the cell population.
An Interesting Observation :-
If observed carefully, it is seen that the time interval between a maxima in the population and its adjacent minima is approximately 2 hours. Thus, we can conclude that the cell cycle remains arrested over a span of 2 hours due to the sense-antisense binding, after which the binding gets removed, and so the arrest, thereby explaining the fact that only cell death occurs within that 2 hours, making the population decline during that time interval.
GFP Characterization of the QS Circuit :-
We first characterize the backbone circuit for quorum sensing, given below, using GFP.
For this, we only use equations (1) to (5) from the QS model and introduce a new ODE to determine the concentration of GFP as follows :-
d[GFP]/dt = lC + kc*N[A2R2]^10/(α^10+A2R2^10) – dGFP[GFP]
The plots indicate that GFP production only starts after AHL has reached a certain threshold concentration, thereby indicating AHL binds to the pLux promoter for gene expression only after it has reached a certain threshold concentration.
We now characterize our main quorum sensing circuit with the antisense using GFP as given below to determine the oscillatory dynamics of the colour intensity expressed by GFP, and hence verify the oscillation in gene expression due to the negative feedback loop.
For this, we use the Equations (1) to (7) from the QS Model and introduce the ODE for GFP, to the system as follows :-
Parameters remain the same.
The oscillation in GFP expression hence verifies the negative feedback loop in the quorum sensing circuit.
Stochastic Model
We also built 2 stochastic models for the bacterial population, one with quorum sensing and another without quorum sensing, by simulating the models using the Gillespie Algorithm, where we treat each reaction/occurrence as an event and determine which event takes place in each timestamp. We then add/subtract an entity from the current number of molecules of that entity for that timestamp, based upon how the number of molecules of that entity changes upon one occurrence of that particular event.
The reason for building the stochastic model is to study the effect of noise in our Quorum Sensing Circuit, and to observe whether noise can destroy or enhance the oscillatory behaviour of the circuit, because any biological system in contact with the outside environment will be influenced by several external factors which can subject the circuit a variety of noise levels, and hence, studying any biological system in the presence of external noise levels become extremely important.
The population model without quorum sensing consists of 4 events as follows :-
We simulate the events using the Gillespie Algorithm to obtain the stochastic population curve for the bacteria with time, for a general bacterial population without quorum sensing. The plot is given below :-
The population model with quorum sensing consists of 19 events as follows :-
We simulate the events using the Gillespie Algorithm to obtain the stochastic population curve for the bacteria with time, having the quorum sensing circuit in it. The plot is given below :-
Chemical Master Equation :-
We thus determine the Master Equation of our Quorum Sensing circuit governing the probabilities of each state in the circuit.
The Master Equation is thus defined as :-
To determine P(N,t), we run numerous realizations of our stochastic model using the Gillespie Algorithm to plot the Stationary Distribution of the number of cells in the population over the time interval, 24 < t (in hours) < 220, for the general bacterial model without quorum sensing and find out the probabilities of the bacterial population being at a particular concentration during the time interval.
The Stationary Distribution for the model with Quorum sensing, within 24 < t (in hours) < 220, is given below :-
The Population Master Equation peaks at around N = 690, thereby being consistent with the bacterial population curve without the quorum sensing circuit, where the stationary phase population density is around 700.
To determine P(N,t), we thus run numerous realizations of our stochastic model using the Gillespie Algorithm to plot the Stationary Distribution of the number of cells in the population over the time interval, 50 < t (in hours) < 600, for our bacterial model with quorum sensing and find out the probabilities of the bacterial population being at a particular concentration during the time interval.
The Stationary Distribution for the model with Quorum sensing, within 50 < t (in hours) < 600, is given below :-
The Population Master Equation peaks at around N = 260, thereby being consistent with the bacterial population curve having the quorum sensing circuit, where the average stationaryL phase population density is around 260.
Observation :-
We observe that external noise doesn’t create much of an effect over our bacterial population with the quorum sensing circuit. It however makes the oscillatory dynamics a bit more vigorous, making the oscillation more prominent that what is observed in the deterministic model. Except that, our bacterial population is observed to oscillate around the same mean of around 260 for both the stochastic and deterministic models, indicating that noise doesn’t have any effect on controlling the population threshold level, but only has an effect upon increasing the oscillatory behaviour of the circuit.
Parameter Sensitivity Analysis
Parameter sensitivity analysis is conducted to identify whether and how the parameters have influence on the bacterial population density during the various phases of its growth.
We define the Global Sensitivity Coefficient as :-
We calculate S over the time interval, 0 < t (in hours) < 700, using a resolution of 10,000 timestamps in between. The parameters for which we perform sensitivity analysis are :-
The Sensitivity Coefficients (S) of the different parameters are plotted below :-
Parameter Space Analysis
We perform parameter space analysis for the parameter spaces (ki, kr) and (ki, kc) to determine the optimal working range for the parameters to achieve the desired output from our QS circuit.
ki – Promoter strength of pTet promoter
kr – Promoter strength of pTac promoter
kc – Promoter strength of pLux promoter
CHOOSING OPTIMAL POPULATION LEVEL :-
We observe that as the carrying capacity of the bacterial population decreases, the substrate consumption and accumulation of secondary metabolites also slows down, and hence, the death phase of the population also comes after a longer time, increasing the population life. However, the mean carrying capacity cannot be lowered beyond a certain level, as then, the amount of linalool production will also reduce due to very less number of cells.
Thus, we need to choose an optimum mean value of cell population at the stationary phase, around which the population will oscillate, to ensure optimum linalool production, as well as prolonged population life.
We choose this optimum mean population level at the stationary phase to be 1/3rd of the original carrying capacity, K, of the population, so that the stationary phase, and hence population life, gets increased 3 times of the initial duration.
Therefore, mean optimal population level, M = K/3 = 233.
We then calculate the Mean Deviation of N from M for different values of the parameter ordered paris: (ki, kr) and (ki, kc), between the time interval 50 < t (in hours) < 160.
We define the Mean Deviation from M, μ, as :-
The colour coding represents the mean deviation from M, μ. The bluer region represents a less value of μ, and thus, a more optimal and desirable working range for the parameter spaces to achieve the desired working of the circuit.
An Interesting Observation :-
We observe that at a higher strength of the pTet promoter (ki), the strength of the pTac promoter (kr) should be lowered to achieve the desired output and vice versa. It therefore reflects that simply increasing the strengths of both the pTet and pTac promoters is not recommendable for achieving our desired output from the circuit.
Similarly, at a higher strength of the pTet promoter (ki), the strength of the pLux promoter (kc) should also be lowered to achieve the desired output and vice versa. Therefore, we can again conclude that simply increasing the strengths of both the pTet and pLux promoters is not recommendable for achieving our desired output from the circuit. However, at very high and very low strengths of the pTet promoter, the strength of the pLux promoter should be maintained at a constant value to obtain the desired output.
Bacterial Diffusion Model
We animate the growth of bacteria, having the quorum sensing circuit in it, in the medium with time using the Gray-Scott Diffusion Model.
The bacterial growth follows the coupled Reaction-Diffusion equation system as given below :-
We took a 50 X 50 hexagonal grid and used the finite difference method to solve the coupled reaction-diffusion PDE system. The model is simulated over a timespan of 7 < t (in hours) < 44 to generate the animation of bacterial diffusion in the population with time. The animation is given below.
Observation :-
Note that the bacterial diffusion occurs during a time interval, and then stops occurring for sometime, and again restarts after a time interval. This can be explained by the constant oscillation in the cell population.
When the population is increasing with time, bacterial diffusion occurs into the surrounding medium. After that, when the population starts to decrease with time, bacterial diffusion stops occurring within that time interval when the population decrease, and the animation becomes static. After that, the population again starts increasing, and so does bacterial diffusion, making the animation dynamic again. And the cycle eventually continues with time.
Results
We now try to fit our math model with the wetlab results. For this, we used the wetlab data obtained from [0.1mM IPTG + AHL] induction. We took this data because 0.1nM IPTG provides the least stress condition and is most suitable for proper bacterial growth.
For comparing the results, we took the CFU data for the bacterial population from the wetlab result, because CFU gives the most reliable cell number for a bacterial culture.
From the wetlab CFU results, it was observed that the carrying capacity for the control bacterial population (without the QS Circuit) was around 1,04,80,00,000. The drylab model has the carrying capacity for the control bacterial cell concentration to be 700.
Hence, to convert the bacterial cell concentration (N) to the actual bacterial population in the culture (to relate it with CFU data), we use the following conversion factor for the bacterial cell concentration :-
CFU (Population) = (1,04,80,00,000/700) * N (Cell Concentration)
We accordingly fit the wetlab CFU data with the math model and plot the CFU vs t for the fitted math model for both the Control and the QS Circuit (Hifi) bacterial populations. The fitted curves are given below :-
The fitting parameters for the Control population, fitted using the Reference model, are given below :-
The fitting parameters for the Hifi/QS Circuit population, fitted using the Quorum Sensing model, are given below :-
NOTE: We have kept the degradation rates of all the molecules (proteins, mRNA) unchanged.
Efficiency of Fit :-
CONTROL :-
The deviation from the wetlab data for the control population for each timestamp is given below :-
Standard Deviation of the fitted Drylab model from the Wetlab data = 6.318*107.
HIFI/QS CIRCUIT :-
The deviation from the wetlab data for the population with the QS circuit (Hifi) for each timestamp is given below :-
Standard Deviation of the fitted Drylab model from the Wetlab data = 2.949*107.
Order of Population (CFU data) = 108.
Order of Standard Deviation = 107.
Order of relative deviation from wetlab data = 107/108 = 0.1
We thereby observe that our drylab model has a relative deviation of just around 0.1 from the wetlab data. We hence conclude that our model fit to the experimental data is quite a good fit when seen from the perspective of biological systems.
Long Term Prediction :-
The prediction made by our math model for the above fitted parameters, upto a time of 10,000 minutes, for both the control population and the QS circuit (Hifi) population, is given below :-
Observations :-
From the Wetlab data and its corresponding model fit, we observe that our bacterial population with Quorum Sensing is indeed being maintained at a population level below the original carrying capacity, thereby increasing the bacterial population life in the medium for the population with quorum sensing upto 9000 minutes, longer than the population life for the control population which is 8000 minutes.
We also observe a slight oscillation in the population with quorum sensing from the fitted math model curve, at around 200 to 1500 mins, thereby verifying the negative feedback loop of the QS circuit.
Conclusion
We observe that the bacterial population oscillates in the presence of the quorum sensing circuit, thus verifying the negative feedback loop of the circuit. But the most important observation is that the bacterial population gets controlled at a level quite lower than the carrying capacity of the control population. Thereby, the consumption of the bacterial medium also slows down due to the presence of less bacteria in the medium, as well as slower accumulation of toxic secondary metabolites, and thus, the media life gets extended up to a longer period the presence of the quorum sensing circuit, quite longer than a normal media life.
We are thus able to achieve our goal of making the bacterial population last longer in the medium by maintaining it at a certain threshold level quite below its original carrying capacity, thereby ensuring that our target molecule, Linalool, gets produced in a sustainable manner over a prolonged period of time.
References
https://link.springer.com/article/10.1007/s11538-016-0160-6
https://www.nature.com/articles/s41540-021-00196-4
https://pubmed.ncbi.nlm.nih.gov/31612240/
https://2018.igem.org/Team:Tsinghua
https://2017.igem.org/Team:Wageningen_UR/Results/Quorum_Sensing
https://www.nature.com/articles/nature08753
https://2020.igem.org/Team:Fudan/Model
https://2021.igem.org/Team:NCTU_Formosa/Efficiency_Optimization_Model
Molecular Docking and Deep Learning
Overview
We synthesize antisense mRNAs, of the DnaA protein, of different lengths to calculate the length dependent binding efficiency of the antisense mRNAs with the sense mRNA and determine the best length of antisense for sense-antisense binding to achieve the target of controlling cell population at our desired level.
We therefore use IntaRNA for predicting RNA-RNA interaction, Alphafold 3 for determining the structure of the sense-antisense complex, and take help of DegScore, a deep learning algorithm, to predict the half lives of the antisense mRNAs of different lengths.
Objective
We derive a relation of the binding efficiency (ΔG) of the sense mRNA of DnaA with the antisense mRNAs based upon their different lengths and also a relation between the half life and length of the antisense mRNAs to determine the length dependent degradation rate of the antisense mRNAs. We then use the relations to develop a model of the DnaA concentration and the bacterial population as a function of the antisense of different lengths and find out which length gives us the most desired outcome.
RNA-RNA Docking
We synthesize the antisense mRNAs of 5 different lengths :-
1) 50 bp
2) 150 bp
3) 275 bp
4) 400 bp
5) 500 bp
We then determine the ΔG of binding of the antisense mRNAs with the sense mRNA of DnaA based upon different lengths and come up with a relation between the ΔG of binding and the length of the antisense mRNAs.
The software used for predicting the ΔG of interaction is IntaRNA and the software used to predict structure of the sense-antisense complex is Alphafold 3.
The ΔG vs L curve for the sense-antisense binding and the corresponding fitting parameters are given below :-
Fig: Here, x-axis (x1) represents length of the antisense (in bp) and y-axis (y1) represents ΔG (in Kcal/mol) of interaction
The fitted equation is thus :-
ΔG = a1L + b1 , with a1 = -1.50269 kcal/mol/bp , b1 = -8.64668 kcal/mol
where, L = Length of the antisense in base pairs (bp)
ΔG = Free energy of interaction (in kcal/mol)
The correlation constant, r = -0.9998 indicates a very efficient fitting of the equation with the values. The higher negative values of ΔG of interaction in case of larger L indicates a very efficient and stable sense-antisense binding for the antisense mRNAs of greater lengths.
Deep Learning
We predict the half lives of the antisense mRNAs of the 5 different lengths (50 bp, 150 bp, 275 bp, 400 bp, 500 bp) using DegScore, a deep learning model which uses ridge regression trained on Eterna Roll-Your-Own structure data to predict RNA degradation.
DegScore Functionality :-
Degscore models degradation at a given nucleotide i as a linear function of nucleotides surrounding i :-
Obtaining length dependent degradation rates :-
The DegScore() is a callable function in Python intrinsic in the DegScore library of the deep learning model. This function takes input the mRNA sequence and its secondary structure (in dots and brackets form) and returns the degradation scores of each nucleotide in the sequence, as well as the half life of the entire mRNA sequence.
We determined the secondary structures of our antisense mRNAs of 5 different lengths using RNAfold, and gave the sequence and structure of each mRNA as input to the Degscore() function called in Python to determine the half life of that mRNA sequence and eventually, obtain the relation between half life and length for the mRNA sequences.
The code for giving the 50 bp antisense as input to the Degscore() function and its corresponding output is given below :-
The Degradation Score obtained at each nucleotide for the 50 bp antisense is given below :-
The t1/2 vs L curve for antisense mRNAs and the corresponding fitting parameters are given below :-
Fig: Here, x-axis (x1) represents length of the antisense (in bp) and y-axis (y1) represents half life (in hrs) of the corresponding antisense mRNA
The fitted equation is thus :-
t1/2 = a2Lb2 , with a2 = 26.9932 hr-1/bpb2 , b2 = -0.639596
where, L = Length of the antisense in base pairs (bp)
t1/2 = Half life of the antisense mRNA in hours (hrs)
The correlation constant, r = -0.9998 indicates a very efficient fitting of the equation with the values. The higher values of t1/2 in case of smaller L indicates that antisense mRNAs of smaller lengths, though having lower ΔG of interaction, are more stable and thus, can persist for a longer time inside the cell.
Model
We use the ODE system from our quorum sensing model and introduce 2 new ODEs to the system, given below :-
The parameters ‘m’ and ‘r3’ have their previous meanings and values from the QS model.
We write m as a function of L (Length of the antisense mRNA) as :-
m = -ΔG/C
m = -(a1*L+b1)/C
We therefore rewrite Equation (11) as :-
We, hence, also write Equation (8) as :-
Now, our QS model, by default, was prepared taking into consideration the antisense mRNA was of length 400 bp. Thus, to determine the value of C, we choose the original value of m from our QS model and put L=400. We determine C = 407.3333 nM-1.
We also modify Equation (7) from our QS model to incorporate the length dependent (L) degradation rate of the antisense mRNA :-
NOTE: All old parameters have their previous meanings from the QS model, and their values are most suitable parameter values initially chosen for obtaining the best desired result from our QS circuit.
Model explanation :-
Equation (10) represents transcription of the DnaA mRNA and Equation (11) represents translation of the mRNA into the DnaA protein.
In Equation (11), we have used the Hill term,
, to demonstrate the repression in the translation of the M into DnaA caused by the binding of the antisense mRNA (CCA) to the sense mRNA (M). We have written :-
m = -(a1*L+b1)/C
because greater the length of the antisense, more negative is the ΔG of interaction, and hence, more efficient will be the sense-antisense binding. A higher value of m makes the Hill function reach 0 at lower concentration of the repressor molecule, indicating a more efficient repression. Thus, a larger antisense length, and hence a more negative ΔG, will ensure that the repression starts occurring at an even lower concentration of [CCA], thereby demonstrating more efficient and quicker repression of translation by antisense mRNAs of longer lengths.
We have used the same ‘m = -(a1*L+b1)/C ’ in both equations (8) and (11) because both of them indicate repression in growth/expression due to the antisense mRNA, CCA, concentration. And hence, the factor by which the DnaA production gets reduced should ideally be the same factor by which the growth of the bacterial population also gets reduced.
We have simply modified Equation (7) by writing dCCA as a function of the antisense mRNA length. Previously, it was only the degradation rate of the 400 bp antisense molecule.
Plots :-
We solve the ODE system using Python and plot the DnaA concentration vs time obtained from repression by each the 5 different lengths of the antisense mRNAs :-
Fig: DnaA concentration for the 5 different antisense mRNA lengths
We also plot the bacterial concentration (N) vs time for the 5 cases :-
Observation :-
We observe that as the length of the antisense mRNA becomes larger, the bacterial population achieves a higher carrying capacity, but as a higher carrying capacity means that the substrate consumption as well as secondary metabolites accumulation will be at higher rates, the death phase of the population also comes faster as the length of the antisense mRNA increases.
Parameter Space Analysis
We perform parameter space analysis for the parameter space (L, kc) to determine the best length of the antisense mRNA and its corresponding pLux promoter strength for most efficient working of our QS circuit.
L – Length of the antisense mRNA for DnaA
kc – pLux promoter strength
We use M=233 as the optimal mean population level and define the mean deviation from M, μ, for this case as :-
The colour coding represents the mean deviation from M, μ. The bluer region represents a less value of μ, and thus, a more optimal and desirable working range for kc at different values of L to achieve maximum efficiency of the circuit.
Analysis for the Fitted QS model :-
We now use the fitted parameters for the quorum sensing model and set M=500 to determine the best possible antisense mRNA length for our actual QS circuit used in experiments (wetlab data, fitted to math model) to achieve the best desired output of controlling the population at a mean level of 500, lower than the carrying capacity of 700.
The corresponding heat map for our working QS circuit, using the fitted parameters, and just varying kc and L, is given below :-
(All things have the same meanings as the previous heat map)
The parameters used in this case are parameters obtained after fitting our QS model with the [0.1nM IPTG + AHL] induction Wetlab experimental data.
Fig: Heat map of mean deviation from M=500, μ, for our working QS circuit used in the wetlab using the fitted parameters. Only kc and L are varied.
Interesting Observations :-
BEST SUITABLE PARAMETER MODEL :-
We observe that the pLux promoter strength only varies by a factor of 102 as antisense mRNA length varies from 50 bp to 600 bp for maximum efficiency of the circuit. It becomes almost constant for antisenses of lengths higher than 300 bp.
We find out that our circuit works most efficiently for 2 parameter spaces :-
kc = 12 nM/hr, L = 300, 350, 400, 450, 500 bp (kc = 6 nM/hr and L = 400 bp are the default circuit parameters)
kc = 120 nM/hr, L = 50, 100 bp.
One important observation we can arrive at from the above data is that as the length of the antisense mRNA decreases, the pLux promoter should be increased in strength to achieve the desired output from the circuit, as the smaller antisense mRNA has a less ΔG of binding with the sense mRNA of the DnaA protein.
FITTED QS MODEL :-
We observe that as the best efficiency of the circuit is observed at around a pLux promoter strength (kc) of 0.01 nM/min. If observed carefully along the row with kc = 0.01 nM/min, the darkest blue grid is obtained at an antisense length of 300 bp.
We hence claim that for our working QS circuit in the wetlab, the best efficiency can be observed for an antisense mRNA of length 300 bp at a pLux promoter strength of 0.01 mM/min.
For our fitted pLux promoter strength of kc = 0.0001 nM.min, the best efficiency is found at an antisense length of around 400-500 bps.
Another interesting observation is that at a particular value of kc, the efficiency of the circuit does not vary significantly for different lengths of the antisense mRNA.
‘AntiRNA’ Algorithm
We saw that there were weren’t any such algorithms available which could predict the best length of the antisense mRNA for sense-antisense binding to achieve the best population control output from our designed Quorum Sensing Circuit.
We, therefore, proceeded upon to develop ‘AntiRNA’, an algorithm developed by incorporating certain deep learning, structure prediction and molecular docking softwares like DegScore, Seqfold and IntaRNA, as well as our modified Quorum Sensing mathematical model for bacterial population determination, to predict the best antisense mRNA length for sense-antisense binding to achieve the desired and most efficient population regulation with our designed QS circuit.
The AntiRNA functionality is given below :-
For more information, visit AntiRNA GitLab
Prediction by AntiRNA
We ran the algorithm by giving the required inputs to obtain the best antisense mRNA length for the DnaA protein for most efficient target achievement.
All ODE parameters were set to default. (The initially set most suitable parameters to obtain the best desired result from our Quorum Sensing model are the default parameters here!)
The output of the run is given below :-
Observation :-
We observe that the best antisense length predicted by AntiRNA is 300 bp, and the Heat Map also shows the same for the given parameters (default).
The 1D Heat map obtained from the algorithm is observed to be consistent with the 2D Heat Map given above for kc = 2 nM/hr, which is the default value of kc in the model.
Conclusion
We therefore observe that our model is successfully able to predict the length dependent binding of the antisense mRNAs of different lengths with the sense mRNA of DnaA to achieve different population threshold levels with our QS circuit.
We assert that for our working QS circuit with fitted wetlab data, the best antisense mRNA length is 300 bps, obtained at a pLux promoter strength of 0.01 mM/min.
Our designed algorithm, AntiRNA is also quite robust and efficiently predicts the best antisense mRNA length for achieving the desired population threshold with the QS circuit.
Such a method of model designing and integration of the math model into our designed algorithm promises to be an extremely robust and efficient method for circuit designing and predicting the different outcomes the circuit can give upon experimentation. This method can thus be characterised as a novel form of approach for in genetic circuit and model designing, combing computational software outputs with mathematical equations to achieve an integrated model to predict appropriate behaviours at the molecular, organismal, as well as population levels.
Additional Work
Linalool is known to effectively bind with several proteins on the surface of the fungal and bacterial membranes, as well as inside it, thereby distorting the structure as well as permeability of the fungal and bacterial membranes and in the process, disrupting their metabolism, thereby halting their growth.
We chose 6 such bacterial and fungal proteins and fatty acids with which linalool is known to bind effectively and used their PDB structures to determine the efficiency of binding of linalool with these proteins and how stable is this linalool-protein complex using the HDOCK Docking WebServer (PDB structure of linalool was generated through software) :-
1) Exo-B-(1,3)-Glucanase from Candida Albicans in complex with Castanospermine at 1.85 A (PDB ID – 1EQC)
2) Crystal structure of BmrR bound to puromycin (PDB ID – 3Q3D)
3) Crystal structure of CviR bound to antagonist chlorolactone (CL) (PDB – 3QP5)
4) Crystal structure of Bacillus subtilis YmdB (PDB – 4B2O)
5) Structure of an integral membrane delta(14)-sterol reductase (PDB – 4QUV)
6) Ergosterol (PDB generated through software)
The docked structures are given below :-
Fig: Docking Score of the various complexes formed with linalool
Inference :-
We observe that the Docking Scores/Energies of the complexes are well negative in value (hence, negative ΔG for the entire formation process), thereby indicating that the complex formation of these molecules with linalool is feasible and the complexes are stable enough to be formed. Thus, we verify that linalool indeed binds with fungal and bacterial membrane proteins/fatty acids and in the process, halting their growth, thereby justifying the antimicrobial properties of linalool.
Linalool is particularly efficient in binding with 1QEC, and membrane Ergosterol, thereby indicating that it is a particularly effective antifungal molecule against Candida albacins, the fungus responsible for most foot fungal infections.
References
https://2021.igem.org/Team:IISER_Kolkata
https://2023.igem.wiki/iiser-kolkata/
https://www.nature.com/articles/s41540-021-00196-4
https://2023.igem.wiki/tsinghua/
https://www.nature.com/articles/s42256-022-00571-8
https://www.nature.com/articles/s41467-022-28776-w
https://academic.oup.com/bioinformatics/article/24/24/2849/196716?login=false
https://www.nature.com/articles/s41467-023-42528-4
https://www.nature.com/articles/s41467-021-21194-4
https://academic.oup.com/nar/article/45/W1/W365/3829194?login=false
https://bmccomplementmedtherapies.biomedcentral.com/articles/10.1186/s12906-024-04457-7
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9377531/
https://academic.oup.com/nar/article-lookup/doi/10.1093/nar/gkx279
Wayment-Steele, H., Kim., D.S. DegScore. (2022) Zenodo. https://doi.org/10.5281/zenodo.7130659
https://www.biorxiv.org/content/10.1101/2021.03.29.437587v1
http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi
https://www.nature.com/articles/s41586-024-07487-w
https://alphafoldserver.com/about
https://www.biorender.com/
Media Dependent Kill Switch
Overview
To prevent live GMO release into the environment, we propose a media-dependent bacterial kill switch, such that the kill switch remains dormant as long as the bacteria lives on the medium. However, as soon as it escapes out of the medium, the kill switch gets activated, immediately killing the bacteria.
Thus, through this kill switch, we seek to ensure that no GMO remains alive once it has left the medium, thereby presenting GMO contamination of the environment.
The genetic circuit for the media-dependent kill switch is given below :-
Working of the Circuit
We chose IPTG as the repressing molecule governing the action of the kill switch as most of the promoters in our genetic circuits will be IPTG induced, and thus it would prove to be an effective media-incorporated molecule for the functioning of our media-dependent kill switch.
If the bacteria leaves the medium, there will be no IPTG present to repress the action of the kill switch, and hence, the kill switch will get activated, leading to the death of the bacteria.
Objective
We determine the robustness of our kill switch by calculating the time required to kill the bacteria once it has escaped out of the medium, as well as realize the probability of bacterial death once the kill switch has got activated.
Model
We generate the ODE system for the model :-
Explanation of the equations :-
The pCl/Bla promoter producing LacI is getting partially repressed by the CIλ protein. So, we have modified the Hill Term for repression of the promoter as,
to demonstrate that in the presence of optimal concentration of CIλ, the
expression of the pCl/Bla promoter does not completely stop, but
instead, gets reduced to 1/3rd of the original expression.
LacI repressed the pLac promoter producing CIλ and ccdA. However, IPTG suppresses this action of LacI upon the pLac promoter and thus, we use the following Hill Term,
to incorporate the repression of pLac by LacI and repression of the effect of LacI by IPTG. Thus, IPTG indirectly induces the expression of the pLac promoter.
The probability that the bacteria dies depends upon the product of the concentrations of ccdB and mf-lon. Greater the concentration of each molecule, greater is the probability of cell death. Therefore, we have used the following function to relate the probability of cell death upon the ccdB and mf-lon concentrations :-
Plots :-
We varied the concentration of IPTG as :-
to inculcate the fact that the bacteria escapes out of the medium at time . So, before time , the IPTG concentration available to the bacteria is I0, but after time τ, it abruptly drops to 0 as IPTG is no longer available to the bacteria after τ.
We thus solve the ODE system using Python and plot Pkill vs t to determine the probability of bacterial death once it has escaped out of the medium at time τ, using 3 different values of τ.
τ = 15 mins :-
τ = 30 mins :-
τ = 50 mins :-
An Interesting Observation:-
We observe that the probability that the bacteria gets killed after escaping from the medium surges to the value of 1 within almost 10 mins of the bacteria leaving the medium. Thus, we can conclude that bacterial death almost surely occurs within the first 10 mins after the bacteria has left the medium.
Now, this 10 mins is quite less than the division time for E.coli, which is around 20 mins, so we can say once the bacteria has left the medium, it dies even before it can divide into 2 daughter cells. Hence, our kill switch not only kills the bacteria, it also ensures that the bacteria gets killed even before it can start multiplying outside the medium, which accordingly verifies the robustness of our kill switch.
Integration with Quorum Sensing Model :-
We define the probability of a bacteria escaping the medium as :-
Here, N and S are the bacterial population and substrate concentrations obtained through simulation of our quorum sensing model. Less the substrate and greater the population, greater is the probability for the bacteria to leave the medium and escape outside.
We thus plot Pesc vs t using Python and obtain :-
Therefore, probability of the occurrence of bacterial death at any point in time can be defined as :-
Pdeath = Pkill. Pesc
We assume that the bacteria escapes the medium at time τ. Hence, we plot Pdeath vs t for 4 different values of τ :-
Bacterial Escape Model
We animate the probability of bacterial escape and eventual death of the bacteria due to the kill switch getting activated using the Gray-Scott reaction-diffusion model.
For that, we use the 2 X 26 hexagonal grid (lower row represents the bacterial medium and upper row represents the environment outside the insole) and use the finite difference method to solve the coupled PDEs given below :-
ASSUMPTIONS:-
1) We assume that the bacteria escapes only when Pesc > 0.74.
2) We assume that the bacteria escapes in a direction perpendicular to the medium (along the vertical axis of the grid)
From the Kill Switch ODE model, we observe that probability of bacterial death becomes almost 1 within 10 mins of the bacterial escape. Thus, we model the animation to demonstrate that the bacteria dies within 10 mins after it has escaped out of the medium.
The model is simulated over a timespan of 7 < t (in hours) < 44 to generate the animation of the bacterial escape from the medium and eventual death due to kill switch with time. The animation is given below.
Bottom row: Yellow represents high bacterial concentration
Top row: Dark green represents an escaped bacteria
Observation :-
We observe that the bacteria escapes out of the population once the cell concentration becomes greater than 650, that is, almost reaches the carrying capacity, thus verifying the Probability of escaping the medium.
The bacteria dies within 10 minutes after it has escaped out of the medium, thus verifying the kill switch ODE model.
Conclusion
We thus observe that our Kill Switch, being governed by the availability of IPTG in the surrounding environment, provides an exhaustively robust and efficient method of cell death outside the bacterial medium, thereby ensuring that no GMO contamination in the environment takes place as even if the GMO has left the medium, it will no longer be able to live upto the time to undergo division, as the kill switch would have already done its job by then, as observed by the model where the probability of bacterial death surges to 1 within just 10 minutes of the bacterial escape, thereby not even providing the bacteria the time to divide.
Our model for the media dependent bacterial kill switch presents a novel approach towards kill witch modelling as it takes a coupled deterministic-probabilistic approach to determine the occurrence of bacterial death, as well as determine the chances of the bacteria escaping out of the medium and insole.
References
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4718764/
https://2020.igem.org/Team:Fudan/Design/kill_switch
https://2023.igem.wiki/iiser-kolkata/
https://2018.igem.org/Team:Oxford/Model/Killswitch
