Model Head Introduction

Each biological phenomenon is made up of many different parts and occurs in a complicated environment. Even while some parts can be researched in the lab, we are still unable to profoundly understand or describe all of the variables that affect the system. Mathematical models contribute to understanding the essence of events. "All models are incorrect, however some are helpful" George Box.

Euphoresis is a very complex yet unique project and its layer had to be simulated.
IPTG Inducible Promoters

Optimization and examination of IPTG induction
Safety Circuit Model

Optimization and parameter investigation of safety circuit and sensitivity

  LacI IPTG Promoter Regulation

Overview

Quorum sensing is the ability of microorganisms to communicate via signaling molecules that can diffuse in the extracellular environment. For biosafety concerns, the synthetic consortium was designed to be able to survive only in the presence of an extracellular factor in the medium and by codependency of its two microorganisms Nostoc and B. subtilis. Therefore, in our project we utilized IPTG inducible killswitch coupled with a quorum sensing system to prevent ‘cell survival’ outside the consortium.

For this purpose the team chose to utilize two Isopropyl β-d-thiogalactoside (IPTG) inducible promoters pGraC and pTrc-2 derivative promoter, whose activations lead to Nostoc’s and B. subtilis’ kill switches respectively. Both promoters are located in B. subtilis. However under the control of pGraC is the LuxI gene that produces the signaling molecule acyl homoserine lactone (AHL) and thus the extracellular factor IPTG also regulates the quorum sensing system that further induces Nostoc’s toxin/antitoxin kill switch. On the other hand, under control of pTrc-2 derivative is the TetR producing the repressor of B. subtilis kill switch. You can read more about our kill switch in our Design page.

Before proceeding to further analysis, it was crucial to break down the complex system of quorum sensing and kill switch mechanisms and make the appropriate assumptions to simplify model fluctuations and complexity starting with IPTG Inducible Promoters Regulation. This model provides the team a profound understanding of promoters activation regulation and the effect of change in concentration of its components. Furthermore we wanted to examine how fast the activation and repression of promoters occur so we would be able to analyze with accuracy the microorganism kill switches.

Model

Our system has two states: survival mode when there is still IPTG on the medium and both Nostoc and B. subtilis ensure their survival and starvation mode when there is no IPTG in the medium, BsrG toxin in B. subtilis is produced and VapB15 antitoxin in Nostoc exceed its toxin concentration VapC15. The model provides the team a profound understanding of promoters activation regulation and the effect of change in concentration of its components. Furthermore we wanted to examine how fast the activation and repression of promoters occur so we would be able to analyze with accuracy the microorganism kill switches. Last but not least, we performed a further analysis to investigate the effect of plasmid copy number and IPTG association constant kD in systems behavior, to provide the wet lab with proposed data. The system obey mass action kinetics and can be described by the following ODEs:

Ordinary Differential Equations

\begin{equation} (1)\quad \frac{d[\text{m}_L]}{dt} = k_{\text{m}_L} C_n - d_{\text{m}_L} [\text{m}_L] \end{equation}
\begin{equation} (2)\quad \frac{d[\text{L}]}{dt} = k_{\text{L}} [\text{m}_{\text{L}}] - d_{\text{L}} [\text{L}] \end{equation}
\begin{equation} (3)\quad \frac{d[\text{L}_2]}{dt} = k_{\text{d}} [\text{L}]^2 - k_{\text{d}}^{\text{r}} [\text{L}_2] - k_{\text{r}} [\text{L}_2\text{P}] - k_{\text{LI}} [\text{L}_2] [I]^2 + k_{\text{KL}}^{\text{r}} [\text{L}_2\text{I}] - d_{\text{L}_2} [\text{L}_2] \end{equation}
\begin{equation} (4)\quad \frac{d[\text{L}_2\text{I}]}{dt} = k_{\text{LI}} [\text{L}_2] [I]^2 - k_{\text{LI}}^{\text{r}} [\text{L}_2\text{I}] + k_{\text{a}} [\text{L}_2\text{P}] [I]^2 - k_{\text{a}}^{\text{r}} [\text{P}] [\text{L}_2\text{I}] - d_{\text{LI}} [\text{L}_2\text{I}] \end{equation}
\begin{equation} (5)\quad \frac{d[\text{I}]}{dt} = -2k_{\text{LI}} [\text{L}_2] [I]^2 + 2k_{\text{LI}}^{\text{r}} [\text{L}_2\text{I}] - k_{\text{a}} [\text{L}_2\text{P}] [I]^2 + 2k_{\text{a}}^{\text{r}} [\text{P}] [\text{L}_2\text{I}] + D_{\text{IPTG}} ([I_{\text{ex}}]-[I]) + 2d_{\text{LI}} [\text{L}_2\text{I}] \end{equation}
\begin{equation} (6)\quad \frac{d[\text{L}_2\text{P}]}{dt} = k_{\text{r}} [\text{P}][\text{L}_2] - k_{\text{r}}^{\text{r}} [\text{L}_2\text{P}] + k_{\text{a}}^{\text{r}} [\text{P}][\text{L}_2\text{I}] - k_{\text{a}}[\text{L}_2\text{P}][I]^2 \end{equation}
\begin{equation} (7)\quad \frac{d[\text{P}]}{dt} = k_{\text{r}}^{\text{r}} [\text{L}_2\text{P}] - k_{\text{r}} [\text{P}][\text{L}_2] + k_{\text{a}}[\text{L}_2\text{P}][I]^2 - k_{\text{a}}^{\text{r}} [\text{P}][\text{L}_2\text{I}] \end{equation}
\begin{equation} (8)\quad [P] + [\text{L}_2\text{P}] = C_n \end{equation}

Survival mode

promoter regulation with IPTG
Figure 1: LacI in the presence of IPTG

LacI is constitutively expressed in B. subtilis and produces the LacI protein, which undergoes further dimerization. Then the LacI dimer binds to pGraC and pTrc-2 derivative promoter and represses transcription of LuxI and TetR respectively. However when IPTG is still in the medium, it diffuses into the cell and inhibits the repression by binding to the LacI dimer and thus allowing the expression of genes of interest.


Results

promoter regulation with IPTG- LacI without IPTG
Figure 2: LacI without the presence of IPTG

Due to its high diffusion rate, IPTG diffuses instantly into the cell. LacI IPTG complex (orange) excess over LacI (blue) and LacI dimer (green), indicates that within seconds LacI dimer (the repressor) is occupied due to invasion of IPTG in the cell resulting in the formation of the LacI IPTG complex, that increases (second plot). There are two derepression mechanisms according to literature. First IPTG depression mechanism occurs from binding to free LacI dimer, while the second one from binding directly to the repressed promoter. Since the data show no existence of a repressed promoter (third plot) we can assume that the free LacI dimer concentration is not sufficient for repressing. In contrast, the activated promoter almost reaches equilibrium (third plot) while the protein of interest is increasing (first plot). We can now hypothesize for the kill switch and quorum sensing models that activation of IPTG inducible promoters happen instantly. Because of the time scale difference between promoter activation (seconds) and protein production (minutes) we can use the Quassy Steady State Approximation (QSSA) to construct formulas that describe the expression of genes of interest with accuracy.


Starvation mode

When IPTG is depleted, the LacI dimer is free to repress the genes of interest.

promoter regulation with IPTG
Figure 3: LacI without the presence of IPTG

Results

Promoter regulation without IPTG
Figure 4: Promoter regulation without IPTG

From the simulation in IPTG starvation conditions, we can observe that no LacI IPTG complex was formed (second plot). Once LacI reaches equilibrium LacI dimer starts increasing (second plot). Hence the promoter is fully repressed (third plot) and no protein production occurs(first plot). From the above data we can now hypothesize that once the microorganisms have no IPTG available in the medium, the repression will allow the production of genes that regulate quorum sensing and kill switch mechanisms.

Transcription LacI: geneLacI → mRNALacI

Degradation mRNA LacI: mRNALacI → ∅

Translation LacI: mRNALacI → LacI

Degradation LacI: LacI → ∅

Dimerization LacI: LacI + LacI ↔ LacI2

Degradation LacI dimer: LacI2 → ∅

LacI IPTG complex formation: LacI2 + IPTG ↔ LacI2IPTG

Degradation LacI IPTG complex: LacI2IPTG → ∅

Promoter Regulation: PromoterUnoccupied + LacI2IPTG ↔ PromoterOccupied

IPTG Diffusion: IPTGExtracellular → IPTGIntracellular

Table 1: Compartment symbols
Compartment Symbol
mRNALacI mL
LacI L
LacI2 L2
LacI2IPTG L2I
PromoterOccupied PO
PromoterUnoccupied Pu
IPTGExtracellular Iex
IPTGIntracellular I

Plasmid CN and Promoter association constant over variation of IPTG inducer concentration

Results

Inspired by Alejandro Vignoni presentation on iGEM 2020 Measurement Webinars, we aim to understand in depth the effect of modifications possible to perform and measure by wet lab.We therefore want to further investigate the effects of different plasmid copy numbers and IPTG concentrations. The goal of this analysis was to find a balance between those two variations to achieve maximal possible and coordinated production of kill switch activators.

IPTG plasmid CN variation
Figure 5: IPTG plasmid CN variation

For the analysis we use three values of plasmid copy numbers (8 blue, 15 black and 30 pink) that correspond to low, medium and high plasmid vectors respectively. The concentrations used for IPTG variation were taken from the range of values provided by the wet lab that ensures no IPTG negative effects on cell growth. As IPTG concentration decreases, there is significant difference in LuxI production while for the highest possible IPTG concentration low, medium and high plasmid vectors appear no big difference. The distinction in TetR production is less significant when IPTG decreases and TetR protein rages in certain values of Y axis. In order to achieve maximal production that has the minimum possible difference between the two products LuxI and TetR when IPTG decreases, it is proposed to use a medium plasmid vector according to the data.

IPTG kDp variation t=30s
Figure 6: IPTG kDp variation t=30s

Accordingly, the promoter affinity with IPTG constant kD was investigated in different IPTG concentrations. Promoter affinity with the IPTG inducer is presented in blue, pink and black color for high, medium and weak affinity respectively, while the IPTG variation is presented in different line types according to the diagram's footnote. From the results above, we conclude that gene transcription and translation rates have a major impact in TetR and LuxI formation, since the variation in affinity seems to have considerable difference in LuxI in contrast to TetR.

Table 1: Table of Parameters
Biochemical Reaction Symbol Rate Value Ref.
LacI IPTG Promoters Regulation
Transcription rate mRNA LacI kmL 0.47 min -1 Registry
Translation rate LacI kL 15 min-1 [3]
Degradation rate mRNALacI dmL 0,462 min-1 [3]
Degradation rateLacI dL 0,2 min-1 [3]
Dimerization association rateLacI kd 50 nM/min [3]
Dimerization association rateLacI krr 12 min-1 [3]
Repression LacI2 dissociation rate krr 2.4 min-1 [3]
Association rate LacI2 IPTG kLI 30x10-7nM-2min [3]
Association constant LacI2IPTG kLI 3x10-7nM-2min [3]
Dissociation rate LacI2IPTG krLI 12 min-1 [3]
Derepressionassociation rate LacI2by IPTG ka 3x10-7nM-2min [3]
Derepression dissociation rate LacI2by IPTG ka 4.8x103nM-2min-1 [3]
Degradation rate LacI2 dLI2 0.2 min-1 [3]
Degradation rate LacI2IPTG complex kind 3x10-7nM-2min [3]
Induction dissociation constant LacI2IPTG krind 4,8x103nM-2min-1 [3]
Diffusion IPTG DIPTG 0,92 min-1 [3]
Kill Switch B.subtilis
Transcription rate mRNA TetR kmT 3.61 min-1 est.
Translation rate TetR kT 10 min-1 [7]
Degradation rate mRNA TetR dmT 2 min-1 [7]
Degradation rate TetR dT 8 min-1 [7]
Transcription rate Toxin BsrG kmTox 4.26 min-1 est.
Translation rate Toxin BsrG kTox 1 min-1 ass.
Degradation rate Toxin BsrG dTox 0,02 min-1 ass.
Degradation rate mRNA Toxin BsrG dmTox 4 min-1 BioNumbers
Induction constant pTet01 KD 0.006 [7]
Hill Coefficient n 2 [7]
Basal expression TetR βTetR 0.92 wet
Basal expression Toxin BsrG βBsrG 0.7 wet
Quorum Sensing Kill switch Nostoc
AHL Transport rate constant of AHL (Gram Negative) DC 3 min-1 BioNumbers
AHL Diffusion rate of AHL through cell membrane DB 2 min-1 [1]
Dissociation rate LuxR AHL monomer krRA 10 min-1 [1]
Dissociation constant LuxR AHL monomer KDRA 100 molecules [1]
Association constant LuxR AHL monomer kRA 0.1 min-1molecules-1 krRA/KDRA
Dissociation rate LuxR AHL dimer kr(RA)2 1 min-1 [1]
Dissociation constant LuxR AHL dimer Dr(RA)2 20 molecules [1]
Association rate LuxR AHL dimer kr(RA)2 0.05 min-1molecules-1 kr(RA)2/Dr(RA)2
Synthesis rate of AHL S 0.08 min-1 [1]
Degradation rate AHL intracellular dAin 0.057 min-1 [1]
Degradation rate AHL extracellular dAex 0.04 min-1 [1]
Degradation rate LuxR AHL monomer dRA 0.156 min-1 [1]
Degradation rate LuxR AHL dimer d(RA)2 0.017 min-1 [1]
Degradation rate mRNA LuxR dmR 0.247 min-1 wet
Degradation rate mRNA LuxI dml 0.047 min-1 [1]
Degradation rate LuxR dR 0.2 min-1 [1]
Degradation rate LuxI dI 0.02 min-1 [1]
Transcription rate LuxR kmR 2.52 min-1 est.
Translation rate LuxR kR 2 min-1 [1]
Transcription rate LuxI kml 3.96 min-1 est.
Translation rate LuxI kI 2 min-1 [1]
Translation rate LuxI kmB 0.37 min-1 est.
Translation rate Antitoxin VapB15 kB 1 min-1 ass.
Transcription rate Toxin VapC15 kmC 0.84 min-1 est.
Translation rate Toxin VapC15 kC 1 min-1 ass.
Degradation rate mRNA Antitoxin VapB15 dmB 4 min-1 BioNumbers
Degradation rate Antitoxin VapB15 dB 0.02 min-1 ass.
Degradation rate mRNA Toxin VapC15 dmC 4 min-1 BioNumbers
Degradation rate Toxin VapC15 dC 0.02 min-1 ass.
Basal expression Antitoxin VapB15 βB 0.52-1 wet
Basal expression LuxI βI 0.8-1 wet
Other Parameters
Plasmid Copy Number Nostoc (medium) CnN 13 molecules/min wet
Plasmid Copy Number B.subtilis (medium) CnB 18 molecules/min-1 wet
Typical transcription elongation rate tr 90 nts/sec BioNumbers
RNA polymerase II footprint f 40 nts BioNumbers
Bacillus subtilis cell volume Vc,B 0.9 μm3 BioNumbers
Cyanobacterium cell volume Vc,N 1.8 μm3 BioNumbers

Transcription rate estimation

First we collected from the literature data about the typical transcription rate (tr) and RNA polymerase II footprint (f). Then, we utilized the formula below to estimate the transcription rate of the genes of interest using the maximal transcription rate of 100000 a.u. provided by Salislab.

\begin{equation} k_{\text{transcription}} = \frac{tr}{f} \cdot \frac{TIR_{\text{salis}}}{100000} \end{equation}
  1. Y. Boada, A. Vignoni and J. Picó, "Multiobjective Identification of a Feedback Synthetic Gene Circuit," in IEEE Transactions on Control Systems Technology, vol. 28, no. 1, pp. 208-223, Jan. 2020, doi: 10.1109/TCST.2018.2885694.
  2. Maass S, Sievers S, Zühlke D, Kuzinski J, Sappa PK, Muntel J, Hessling B, Bernhardt J, Sietmann R, Völker U, Hecker M, Becher D. Efficient, global-scale quantification of absolute protein amounts by integration of targeted mass spectrometry and two-dimensional gel-based proteomics. Anal Chem. 2011 Apr 1;83(7):2677-84. doi: 10.1021/ac1031836. Epub 2011 Mar 11. PMID: 21395229.
  3. Stamatakis M, Mantzaris NV. Comparison of deterministic and stochastic models of the lac operon genetic network. Biophys J. 2009 Feb;96(3):887-906. doi: 10.1016/j.bpj.2008.10.028. Erratum in: Biophys J. 2009 Nov 4 ;97(9):2651. PMID: 19186128; PMCID: PMC2996131.
  4. Stekel, D.J., Jenkins, D.J. Strong negative self regulation of Prokaryotic transcription factors increases the intrinsic noise of protein expression. BMC Syst Biol 2, 6 (2008). https://doi.org/10.1186/1752-0509-2-6
  5. Roee Amit, Hernan G. Garcia, Rob Phillips, Scott E. Fraser, Building Enhancers from the Ground Up: A Synthetic Biology Approach, Cell, Volume 146, Issue 1, 2011, Pages 105-118, ISSN 0092-8674, https://doi.org/10.1016/j.cell.2011.06.024.
  6. Tran DTM, Phan TTP, Doan TTN, Tran TL, Schumann W, Nguyen HD. Integrative expression vectors with Pgrac promoters for inducer-free overproduction of recombinant proteins in Bacillus subtilis. Biotechnol Rep (Amst). 2020 Oct 9;28:e00540. doi: 10.1016/j.btre.2020.e00540. PMID: 33163371; PMCID: PMC7599426.

Assumptions of the model

  1. The operation of the expression system is the same for both pGraC and pTrc-2 derivative promoters and is studied under the same regulation.
  2. The amount of species transformed by the reactions only depends on the current amount of species, the rates at which these reactions proceed, and their stoichiometry.
  3. It models what happens in a single cell.


  Safety Circuit Model

Overview

For biosafety concerns, the team decided to design two kill switch mechanisms responsible for each microorganism cell death. Additionally, in order to link the cell density to cell survival, a quorum sensing circuit needed to be coupled. In B. subtilis the BsrG gene is put under control of the pTrc-2 promoter which is repressed by the TetR protein when IPTG is present. On the other hand, in Nostoc the VapB15 gene is under control of a quorum sensing inducible promoter, while the VapC15 gene expresses constitutively VapB15 antitoxin’s toxin. The quorum sensing signaling molecule N-Acyl homoserine lactone (AHL) is produced by the LuxI gene located in B. subtilis and controlled by its IPTG inducible promoter. Then, AHL diffuses through B. subtilis cell membrane into the extracellular environment and when reaches high concentrations into Nostoc’s cells. AHL forms a dimer complex with the LuxR gene that is constitutively expressed in Nostoc and induces antitoxin expresion. When the cells escape from the microspheres, IPTG is absent resulting in the expression of B. subtilis toxin leading to cell death. Additionally, Nostoc’s cell death is achieved since there is no more AHL induction due to LuxI repression in IPTG inducer absence and toxin levels exceed the antitoxin ones. This adds an extra layer of safety, as Nostoc is now dependent on B. subtilis for AHL, while B. subtilis from IPTG making them unable to escape alone.

We aim to provide the team with an in-depth understanding of that complex system as well as data that can improve systems behavior and particular points that need to be carefully examined in the lab.

Ordinary Differential Equations

The KS QS system involves much interaction of compounds between B. subtilis and Nostoc, thus, we use differential equations to describe the rate of change of each compound. The kinetics of the system can be described by the law of mass action. However, according to the results from the IPTG promoter regulation, Hill functions are derived for LuxI and TetR genes (1), (2) in respect to the Quasi Steady State Approximation (QSSA). Additionally, following the same reasoning, we can also utilize QSSA for BsrG and VapB15 production and form Hill function equations (3), (8) due to time scale difference in protein and mRNA production. Hence we are able to describe systems behavior in a more accurate way while reducing its complexity.

B.subtilis:

\begin{equation} (1)\quad \frac{d[\text{LuxI}]}{dt} = k_{\text{m}_I} k_I \frac{C_{nB}}{d_{\text{m}_I} d_I} \left( \beta_I + \frac{(1-\beta_I ) [\text{IPTG}]^n}{k_{\text{pGrac}} \frac{C_{nB}}{[\text{IPTG}]^n} + [\text{IPTG}]^n} \right) \end{equation}
\begin{equation} (2)\quad \frac{d[\text{TetR}]}{dt} = \frac{k_{\text{m}_T} k_T C_{nB}}{d_{\text{m}_T} d_T} \left( \beta_T + \frac{(1-\beta_T) [\text{IPTG}]^n}{k_{\text{pTrc-2}} \frac{C_{nB}}{[\text{IPTG}]^n} + [\text{IPTG}]^n} \right) \end{equation}
\begin{equation} (3)\quad \frac{d[\text{BsrG}]}{dt} = \frac{k_{\text{m}_B} k_B C_{nB}}{d_{\text{m}_G} d_G} \left( \beta_G + \frac{(1-\beta_G) (K_{D} C_{nB})^n}{(K_{D} C_{nB})^n + [\text{TetR}]^2} \right) \end{equation}

Nostoc:

\begin{equation} (4)\quad \frac{d[\text{mRNA LuxR}]}{dt} = k_{\text{m}_R} Cn_{N} - d_{\text{m}_R} [\text{mRNA LuxR}] \end{equation}
\begin{equation} (5)\quad \frac{d[\text{LuxR}]}{dt} = k_{R} [\text{mRNA LuxR}] - d_{R} [\text{LuxR}] \end{equation}
\begin{equation} (6)\quad \frac{d[\text{AR}]}{dt} = k_{\text{AR}} [\text{LuxR}][\text{AHL}_{\text{in}}] - k_{\text{AR}}^r[\text{AR}] - d_{\text{AR}}[\text{AR}] \end{equation}
\begin{equation} (7)\quad \frac{d[\text{AR}_2]}{dt} = k_{\text{AR}_2} [\text{AR}]^2 - k_{\text{AR}_2}^r [\text{(AR)}_2] - d_{\text{AR}_2} [\text{(AR)}_2] \end{equation}
\begin{equation} (8)\quad \frac{d[\text{VapB}15]}{dt} = \frac{k_{\text{m}_B} k_B Cn_{N}}{d_{\text{m}_B} d_B} \left( \beta_B + \frac{(1-\beta_B) [\text{AHL}_{\text{in}}]^n}{\frac{k_{\text{AR}} Cn_{N}}{[\text{AHL}_{\text{in}}]^n} + [\text{AHL}_{\text{in}}]^n} \right) \end{equation}
\begin{equation} (9)\quad \frac{d[\text{mRNA VapC}15]}{dt} = k_{\text{m}_C} Cn_{N} - d_{\text{m}_C} [\text{mRNA VapC}15] \end{equation}
\begin{equation} (10)\quad \frac{d[\text{VapC}15]}{dt} = k_C [\text{mRNA VapC}15] - d_C [\text{VapC}15] \end{equation}

Quorum Sensing:

\begin{equation} (11)\quad \frac{d[\text{AHL}]}{dt} = S [\text{LuxI}] - D \left( [\text{AHl}] - \left[ \text{AHL}_{\text{ex}} \right] - \left[ \text{AHL}_{\text{in}} \right] \right) - d_{\text{AHL}} [\text{AHL}] \end{equation}
\begin{equation} (12)\quad \frac{d[\text{AHL}_{\text{ex}}] }{dt} = D \left( [\text{AHL}] - \left[ \text{AHL}_{\text{ex}} \right] - \left[ \text{AHL}_{\text{in}} \right] \right) \frac{V_{\text{c}}}{V_{\text{e}}} - T \left( [\text{AHL}_{\text{in}}] - [\text{AHL}_{\text{ex}}] \right) \frac{V_{\text{c}}}{V_{\text{e}}} - d_{\text{AHL}_{\text{ex}}} [\text{AHL}_{\text{ex}}] \end{equation}
\begin{equation} (13)\quad \frac{d[\text{AHL}_{\text{in}}]}{dt} = T([\text{AHL}_{\text{in}}] - [\text{AHL}_{\text{ex}}]) - d_{\text{AHL}}[\text{AHL}_{\text{in}}] \end{equation}

Understanding Systems Kinetics

Visualisation of safety circuit
Figure 7: Visualisation of safety circuit

Once IPTG is present in B. subtilis, LuxI and TetR proteins production occurs (1), (2). BsrG toxin that is responsible for B. subtilis cell death, is repressed and we expect low levels of BsrG production due to the gene's basal expression (3). Additionally, LuxI protein is formed and the AHL signaling molecule starts synthesizing and diffuses into the extracellular environment (11). Once it is transported to Nostoc (13), AHL LuxR dimer is formed (6) and dimerized (7) because of LuxR transcription (4) and expression (5). Then the induction of antitoxin VapB15 begins (8) and ensures cell survival overcoming its toxins VapC15 production (9),(10).

Quorum Sensing Results

IPTG plasmid CN variation
Figure 8: IPTG Plasmid CN variation

From the simulation results we observe that AHL (blue) diffuses relatively fast, and reaches high concentration in the extracellular environment (black) while diffusing in Nostoc also happens rapidly (pink), possibly due to high synthesis rate of AHL in B. subtilis and high extracellular levels. Between 1-3 minutes of diffusion AHL levels in Nostoc reach extremely high amounts and continue their increasing trend.

AHL diffusion stimulation t=1min
Figure 9:
 AHL diffusion stimulation t=1min
AHL diffusion stimulation t=5min
Figure 10:
 AHL diffusion stimulation t=5min
AHL diffusion stimulation t=20min
Figure 11:
AHL diffusion stimulation t=20min

Kill Switch Mechanisms

TetR BsrG simulation
Figure 12: TetR BsrG simulation

In the first plot, B. subtilis TetR (purple) production is increasing rapidly leading to BsrG (green) repression and as expected there is only a small, negligible amount of the toxin available due to its gene basal expression. In the second plot we can observe that just in a few seconds, Nostoc is expressing LuxR protein (red) and LuxR AHL complex is formed (yellow) at an accelerated rate because of AHL’s invention into the cell. Antitoxin VapB15 (purple) is produced due to its gene basal expression, a fact that is further confirmed by the Toxin Antitoxin plot. There we can observe the low VapC15 levels in red (again due VapC15 gene basal expression) and the dramatic rise of VapB15 at the same second as in the second plot.

TetR BsrG simulation
Figure 13:
 TetR BsrG simulation
Nostoc quorum sensing simulation t=12sec
Figure 14:
 Nostoc quorum sensing simulation t=12sec
Toxin Antitoxin simulation
Figure 15:
Toxin Antitoxin simulation

The data below can be summarized and combined for examination in the plot below. In IPTG presence both kill switch’s toxins VapC15 and and BsrG are at low levels while the antitoxin already overcomes the basal expression of the gene. That results from successful activation of LuxI and TetR, the genes responsible for the safety circuit. Finally, AHL is produced from B. subtilis diffuses outside the cell and some amount of AHL is already present in Nostoc, verifying the successful VapG15 induction.

KS overview simulation
Figure 16: KS overview simulation

Copy Number RBS variation

We additionally inquired into the impacts of various plasmid copy numbers and RBS variation in BsrG toxin of B. subtilis and VapB15 antitoxin in Nostoc, in an effort to inform the wet lab team about adjustments for systems improvement and to accurately interpret the results.


B. subtilis BsrG toxin:


BsrG plasmid CN RBS variation
Figure 17: BsrG plasmid CN RBS variation

For the analysis we use three values to describe the relative RBS strength following iGEM registry’s characterizations (0.2 blue, 0.6 pink and 1 black) that correspond to low, medium and high RBS strength respectively. The copy numbers used for high, medium and low plasmid vectors are embedded in the plot's memo. BsrG is repressed so it is not reaching high levels. However the increasing rate of its production is valuable to predict its behavior when IPTG is depleted. According to the data, in order to achieve an immediate production rate to ensure most rapid cell death, it is proposed to use either a high plasmid vector with high or medium RBS strength or a medium plasmid vector with strong RBS.


Nostoc VapB15 antitoxin:


VapB15 RBS plasmid CN variation
Figure 18: VapB15 RBS plasmid CN variation

Following the same reasoning of data interpretation, the RBS relative strength and plasmid copy number were examined. In the case of antitoxin, corresponding to the intended properties of the system, a medium antitoxin production rate is acquired so that the VapB15 always overcome VapC15 levels, but also when AHL is depleted, the VapB15 levels can easily weaken. Therefore it is proposed to utilize either a medium plasmid copy number with medium RBS strength, low plasmid copy number with high or medium RBS strength or high copy number with low RBS strength.

Sensitivity analysis

According to model predictions, our system's behavior in IPTG presence appears to have the desired outcome and can be further optimized with modifications in some of its parameters. Sensitivity analysis was a vital next step for our investigation since we want to design a robust kill switch circuit for Euphoresis. Our aim was to investigate the sensitivity of our system to low levels of the IPTG inducer and the signaling molecule AHL, especially when these fundamental elements are diminished. We performed sensitivity analysis to determine which system traits had the most effects on the sensitivity of our system. The outcomes demonstrated that our kill switch circuit has outstanding sensitivity, which is precisely aligned with the objectives of Euphoresis. Due to the great sensitivity, our circuit can be strongly activated when AHL and IPTG are scarce, which ultimately indicates effective cell death. Additionally, sensitivity analysis offered insightful data regarding the time frame in which this sensitivity is most evident. Our capacity to convey to the wet lab what to expect when monitoring these parameters contributes to applying this information to guide our research in the lab. This cooperative strategy makes sure that our wet lab team can interpret the data more precisely, which will ultimately result in more effective experimental design and a better comprehension of the behavior of our system. Due to its ability to shed light on the crucial variables affecting a system's performance, sensitivity analysis has thus been crucial in helping us reach our project's objectives.

B. subtilis BsrG toxin

The sensitivity analysis was accomplished using Sobol indices, a reliable method for assessing the impact of input variables on system outcomes. We created a substantial set of 5000 samples to guarantee the validity of our investigation. We were able to completely examine IPTG variation in depletion circumstances and gain reliable sensitivity insights for BsrG behavior because of the extensive number of samples. The behavior of BsrG is essential to ensure efficient cell death. Our simulation conducted for 60 minutes, which allowed us to monitor the system's dynamic response to the IPTG inducer when its stocks become scarce. From the results, we observe low fluctuation of BsrG behavior in the first few minutes, while the sensitivity increases and reaches constant maximal sensitivity (0.8 - 1) for the next 50 minutes which starts to decrease in the next 10 minutes indicating the complete depletion of IPTG and therefore the continued production of BsrG which is no longer dependent on IPTG stocks. The data extracted from this analysis are perfectly in line with the objectives of Euphoresis and point the wet lab to carefully monitor the variation of BsrG over 60 minutes.

First-order Sobol index maintains all other parameters constant while measuring the specific influence of IPTG depletion on the BsrG behavior. The overall influence of IPTG variation input on BsrG output, including both its direct effect and any interactions with other parameters, is measured by overall-order Sobol index. The fact that the first-order and total-order Sobol indexes exhibit comparable results over time indicates that the sensitivity analysis should include other parameters such as the association constant of IPTG to TetR and TetR association constant to BsrG as well as parameters affecting those two genes promoter and RBS strengths. Due to time limitation and computational cost of such analysis to be performed using the extensive sample range of 5000 samples to achieve maximal accuracy, the analysis has not been conducted and remains to be carried out in the future. Sobol first order index indicates that IPTG inducer influence could become insignificant after 50 minutes. In other words, during the first 50 minutes, the direct effects of IPTG concentration on the generation of BsrG toxin remain relatively steady. After this time, the system possibly obtained a steady state with minimal sensitivity to variations in IPTG levels.

Sensitivity analysis about the effect of IPTG decleration using first and second sobol index on BsrG over time
Sensitivity analysis about the effect of IPTG decleration using first and second sobol index on BsrG over time
Figure 19:Sensitivity analysis about the effect of IPTG decleration using first and second sobol index on BsrG over time

Nostoc VapB15 antitoxin

Accordingly we calculated VapB15 sensitivity over IPTG inducer and AHL signaling molecule depletion utilizing Sobol indices technique over 5000 samples to achieve maximum accuracy. Therefore we investigated with accuracy IPTG and AHL direct impact to the VapB15 behavior through an extensive sample set. High sensitivity of the antitoxin to the above-mentioned factors is crucial for the system to efficiently and immediately deactivate the safety circuit leading to cell death. According to the results, IPTG reduction exhibits high sensitivity over VapB15 production in a time frame of just 10 minutes before it starts to decrease, which perfectly aligns with Euphoresis objectives. Therefore we may conclude that in that 10 minutes wet lab should very carefully examine both VapB15 gene expression as well as its inducer LuxI for outcome’s verification. Furthermore, analyzing VapB15 sensitivity output over AHL reduction, we observe a high sensitivity (over 1) for the first few minutes while up to 40 minutes VapB15 behavior reaches a steady state still considerably dependent on AHL produced by B. subtilis (0.8 - 1). This indicates that as long as AHL is still in the extracellular environment due to AHL stock produced from B. subtilis VapB15 behavior remains steady. We therefore proposed to further investigate LuxI and VapB15 dynamic responses by the wet lab team.

Again the first order index examines the direct parameters effect on VapB15 while keeping the other parameters steady and total order index investigates both direct and indirect effects considering the rest of the system parameters. In the particular analysis no other kinetic parameter was taken into account due to time limitations and computational burden of performing this analysis with accuracy and extensive data set. Therefore, as the results indicate, total and first order indexes exhibit similar sensitivity outcomes and a further analysis is needed to achieve optimal understanding of systems parameters influence.

Sensitivity analysis about the effect of IPTG and AHL decleration using first and second sobol index on VapB15 over time
pp
Figure 20:Sensitivity analysis about the effect of IPTG and AHL decleration using first and second sobol index on VapB15 over time

Overall Model Sensitivity

The outcome of VapB15 and BsrG behavior over IPTG reduction is interpreted in the same plot under 60 minutes in order to help the team genuinely and efficiently the differences when tuning the same parameter, IPTG concentration in this case. Sobol's first order index is high (0.8-1) for VapB15 in the first few minutes while up to 60 minutes it is considered insignificant (0-0.2). BsrG exhibits stable, high sensitivity according to first order index (0.8-1) until 50 minutes that starts diminishing linearly. Hence, the kill switch mode is hyper sensitive in low levels of IPTG as we intended to achieve.

VapB15 RBS plasmid CN variation
Figure 21: VapB15 RBS plasmid CN variation

Assumptions of the model

  1. Mean cell volume is a constant. We assume that the cell volume is constant,thus, the cell surface area is also constant.
  2. Το IPTG diffuses rapidly into the cell and its concentration is stable.
  3. We regard the transmembrane diffusion of AHL as a simple free diffusion that can be described by Fick's law. Thus the flux of transmembrane molecules is directly proportional to the gradient concentration. Fick's First Law: \(J = -D \frac{d[S]}{dx}\)
  4. All degradations occurring in the system can be described with first order reactions.

Table 2: Table of Parameters
Biochemical Reaction Symbol Rate Value Ref.
LacI IPTG Promoters Regulation
Transcription rate mRNA LacI kmL 0.47 min -1 Registry
Translation rate LacI kL 15 min-1 [3]
Degradation rate mRNALacI dmL 0,462 min-1 [3]
Degradation rateLacI dL 0,2 min-1 [3]
Dimerization association rateLacI kd 50 nM/min [3]
Dimerization association rateLacI krr 12 min-1 [3]
Repression LacI2 dissociation rate krr 2.4 min-1 [3]
Association rate LacI2 IPTG kLI 30x10-7nM-2min [3]
Association constant LacI2IPTG kLI 3x10-7nM-2min [3]
Dissociation rate LacI2IPTG krLI 12 min-1 [3]
Derepressionassociation rate LacI2by IPTG ka 3x10-7nM-2min [3]
Derepression dissociation rate LacI2by IPTG ka 4.8x103nM-2min-1 [3]
Degradation rate LacI2 dLI2 0.2 min-1 [3]
Degradation rate LacI2IPTG complex kind 3x10-7nM-2min [3]
Induction dissociation constant LacI2IPTG krind 4,8x103nM-2min-1 [3]
Diffusion IPTG DIPTG 0,92 min-1 [3]
Kill Switch B.subtilis
Transcription rate mRNA TetR kmT 3.61 min-1 est.
Translation rate TetR kT 10 min-1 [7]
Degradation rate mRNA TetR dmT 2 min-1 [7]
Degradation rate TetR dT 8 min-1 [7]
Transcription rate Toxin BsrG kmTox 4.26 min-1 est.
Translation rate Toxin BsrG kTox 1 min-1 ass.
Degradation rate Toxin BsrG dTox 0,02 min-1 ass.
Degradation rate mRNA Toxin BsrG dmTox 4 min-1 BioNumbers
Induction constant pTet01 KD 0.006 [7]
Hill Coefficient n 2 [7]
Basal expression TetR βTetR 0.92 wet
Basal expression Toxin BsrG βBsrG 0.7 wet
Quorum Sensing Kill switch Nostoc
AHL Transport rate constant of AHL (Gram Negative) DC 3 min-1 BioNumbers
AHL Diffusion rate of AHL through cell membrane DB 2 min-1 [1]
Dissociation rate LuxR AHL monomer krRA 10 min-1 [1]
Dissociation constant LuxR AHL monomer KDRA 100 molecules [1]
Association constant LuxR AHL monomer kRA 0.1 min-1molecules-1 krRA/KDRA
Dissociation rate LuxR AHL dimer kr(RA)2 1 min-1 [1]
Dissociation constant LuxR AHL dimer Dr(RA)2 20 molecules [1]
Association rate LuxR AHL dimer kr(RA)2 0.05 min-1molecules-1 kr(RA)2/Dr(RA)2
Synthesis rate of AHL S 0.08 min-1 [1]
Degradation rate AHL intracellular dAin 0.057 min-1 [1]
Degradation rate AHL extracellular dAex 0.04 min-1 [1]
Degradation rate LuxR AHL monomer dRA 0.156 min-1 [1]
Degradation rate LuxR AHL dimer d(RA)2 0.017 min-1 [1]
Degradation rate mRNA LuxR dmR 0.247 min-1 wet
Degradation rate mRNA LuxI dml 0.047 min-1 [1]
Degradation rate LuxR dR 0.2 min-1 [1]
Degradation rate LuxI dI 0.02 min-1 [1]
Transcription rate LuxR kmR 2.52 min-1 est.
Translation rate LuxR kR 2 min-1 [1]
Transcription rate LuxI kml 3.96 min-1 est.
Translation rate LuxI kI 2 min-1 [1]
Translation rate LuxI kmB 0.37 min-1 est.
Translation rate Antitoxin VapB15 kB 1 min-1 ass.
Transcription rate Toxin VapC15 kmC 0.84 min-1 est.
Translation rate Toxin VapC15 kC 1 min-1 ass.
Degradation rate mRNA Antitoxin VapB15 dmB 4 min-1 BioNumbers
Degradation rate Antitoxin VapB15 dB 0.02 min-1 ass.
Degradation rate mRNA Toxin VapC15 dmC 4 min-1 BioNumbers
Degradation rate Toxin VapC15 dC 0.02 min-1 ass.
Basal expression Antitoxin VapB15 βB 0.52-1 wet
Basal expression LuxI βI 0.8-1 wet
Other Parameters
Plasmid Copy Number Nostoc (medium) CnN 13 molecules/min wet
Plasmid Copy Number B. subtilis (medium) CnB 18 molecules/min-1 wet
Typical transcription elongation rate tr 90 nts/sec BioNumbers
RNA polymerase II footprint f 40 nts BioNumbers
Bacillus subtilis cell volume Vc,B 0.9 μm3 BioNumbers
Cyanobacterium cell volume Vc,N 1.8 μm3 BioNumbers

Transcription rate estimation

First we collected from the literature data about the typical transcription rate (tr) and RNA polymerase II footprint (f). Then, we utilized the formula below to estimate the transcription rate of the genes of interest using the maximal transcription rate of 100000 a.u. provided by Salislab.

\begin{equation} k_{\text{transcription}} = \frac{tr}{f} \cdot \frac{TIR_{\text{salis}}}{100000} \end{equation}
  1. Y. Boada, A. Vignoni and J. Picó, "Multiobjective Identification of a Feedback Synthetic Gene Circuit," in IEEE Transactions on Control Systems Technology, vol. 28, no. 1, pp. 208-223, Jan. 2020, doi: 10.1109/TCST.2018.2885694.
  2. Maass S, Sievers S, Zühlke D, Kuzinski J, Sappa PK, Muntel J, Hessling B, Bernhardt J, Sietmann R, Völker U, Hecker M, Becher D. Efficient, global-scale quantification of absolute protein amounts by integration of targeted mass spectrometry and two-dimensional gel-based proteomics. Anal Chem. 2011 Apr 1;83(7):2677-84. doi: 10.1021/ac1031836. Epub 2011 Mar 11. PMID: 21395229.
  3. Stamatakis M, Mantzaris NV. Comparison of deterministic and stochastic models of the lac operon genetic network. Biophys J. 2009 Feb;96(3):887-906. doi: 10.1016/j.bpj.2008.10.028. Erratum in: Biophys J. 2009 Nov 4 ;97(9):2651. PMID: 19186128; PMCID: PMC2996131.
  4. Stekel, D.J., Jenkins, D.J. Strong negative self regulation of Prokaryotic transcription factors increases the intrinsic noise of protein expression. BMC Syst Biol 2, 6 (2008). https://doi.org/10.1186/1752-0509-2-6
  5. Roee Amit, Hernan G. Garcia, Rob Phillips, Scott E. Fraser, Building Enhancers from the Ground Up: A Synthetic Biology Approach, Cell, Volume 146, Issue 1, 2011, Pages 105-118, ISSN 0092-8674, https://doi.org/10.1016/j.cell.2011.06.024.
  6. Tran DTM, Phan TTP, Doan TTN, Tran TL, Schumann W, Nguyen HD. Integrative expression vectors with Pgrac promoters for inducer-free overproduction of recombinant proteins in Bacillus subtilis. Biotechnol Rep (Amst). 2020 Oct 9;28:e00540. doi: 10.1016/j.btre.2020.e00540. PMID: 33163371; PMCID: PMC7599426.
forest