Team:TU-Eindhoven/Model
Model
“Simplicity is the soul of efficiency” – Austin Freeman
Welcome to the modeling page! To get more insight into the reactions of the two-component system and the production of the gas vesicles we built a kinetic model. This model will be used to assess the protein cascades, transcription, translation, and interactions, as well as optimize our concept by tweaking input parameters. On this page, we will give a quick overview of the model, a more in-depth explanation of how the model was constructed, and discuss the results.
The Overview
The model that was built is a deterministic, mass-action kinetics-based model, using ordinary differential equations (ODEs) to simulate the activity of the two-component TtrS/R system and subsequent activation of the acoustic reporter genes (ARG). This model aims to explore the effect of different input parameters on the expression of gas vesicle proteins. By varying the amount of tetrathionate (Ttr) that is added to the model, we can estimate and predict the behavior of our system. This provides us with insight into what the dose-response behavior is expected to be like, what our theoretical limit of detection might be, and how long our signal lags behind initial detection of the biomarker. In addition, the number of gas vesicles (GVs) that are produced can be assessed and related to the generated ultrasound signal.
We worked with a computational modeling program called SimBiology. While working with SimBiology we discovered that the online documentation for this app is quite limited, and only consists of basic tutorials for getting started with SimBiology. Therefore, we wanted to make a guide that includes tips and tricks that can’t be found online easily, but are very useful if you use this program for iGEM. This guide can be found on Contributions.
The Method
Introduction
The kinetic model was built using the MATLAB add-on SimBiology and simulated using the add-on SimBiology Analyzer. These apps provide tools for modeling, simulating, and analyzing dynamic systems based on ODEs which are calculated automatically. The kinetics of our system can be described through the ODEs corresponding to the chemical reactions listed below. We are dealing with a stiff model because the reactions in our system contain time constants varying by several orders of magnitude. Therefore, we chose to simulate our model using the Sundials solver, which is a solver suited for stiff systems of ODEs.
Figure 1: Simplified Schematic of the TtrS/R system in combination with the expression of ARG
The system that we used in our bacteria includes several components, which we all tried to include in our final model. A schematic representation of the components that we modeled can be seen in Figure 1. In addition to this full model, an intermediate model was made to validate the accuracy of the TtrS/R system, containing GFP expression instead of GV expression.
The first component of the model comprises the activation of the lac promoter ($P_{lac}$) by IPTG (Figure 1 step 1). Activation of this promoter is modeled to lead to the transcription and translation of TtrS, which is the sensing protein, and TtrR, the regulator protein. In reality, TtrS and TtrR are present on separate plasmids. Transcription of the TtrS mRNA is in fact regulated by IPTG, however, the transcription of the TtrR mRNA is in reality regulated by doxycycline instead. This mechanism is explained in more detail at Proof of Concept. In the lab, no doxycycline was used to induce transcription of TtrR, as the leakiness of the promoter itself resulted in sufficient TtrR production. For more information on this see Engineering Success segment Doxycycline. In SimBiology, we chose to simplify the transcription of TtrR and TtrS by modeling them under the same promoter, meaning we model a single transcription reaction for both regulators. Because no doxycycline is needed for induction, this simplification will have a limited impact on the accuracy of the model.
Lac promoter activation is modeled similarly to the model by Stamatakis and Mantzaris [1]. In this component of the model three reactions are used to describe the activation of the lac promoter. The lac inhibitor ($LacI$) is modeled as a constant tetramer species that binds to the active promoter ($P_{lac}^*$), thereby inactivating it. When IPTG is added, two reactions can occur. In the first, IPTG binds to free lac inhibitor, making it unable to bind to the active promoter, resulting in effective promoter activation. IPTG is also able to bind to the lac inhibitor in its bound state, causing it to release from the promoter, resulting in promoter activation. We assumed that the inhibitor concentration is consistent, meaning production and degradation are not included.
Transcription is modeled as a two-step process: binding of RNA polymerase ($RNAP$) to the promoter, and transcription of the complex towards mRNA transcripts (Figure 1 step 2). The translation reactions of TtrR and TtrS are modeled separately as they have separate ribosome binding sites and mRNA lengths. This process is, again, modeled as a two-step process: complexation of the ribosome ($R$) and the mRNA transcript ($T_{TtrR}$ and $T_{TtrS}$), and production of the corresponding proteins ($TtrR$ and $TtrS$) (Figure 1 step 3). After translation, the mRNA transcript and ribosome are released and can be reused.
Further, the signaling cascade induced by tetrathionate ($Ttr$) is modeled (Figure 1 step 4-7). The cascade is modeled as binding of tetrathionate to the sensing protein TtrS, inducing autophosphorylation of TtrS (Figure 1 step 4). TtrR can then bind to either phosphorylated ($TtrS_p$), or dephosphorylated TtrS. The binding of unphosphorylated TtrS and TtrR is not depicted in the schematic but was included in the model. The binding of phosphorylated TtrS to TtrR results in a complex which in turn results in phosphorylation of TtrR ($TtrR_p$) (Figure 1 step 5). Dephosphorylation of phosphorylated species is modeled independently of phosphatases, as it is not known whether these play a role in the dephosphorylation of TtrR and TtrS. Phosphorylated TtrR dimerizes into a complex ($TtrR_p:TtrR_p$) (Figure 1 step 6), which then binds to the promoter of ARG ($P_{ttrB}$). Unphosphorylated TtrR could possibly also form this dimer, however, whether this occurs, if this unphosphorylated dimer can act as a promoter, and to what degree that would be, remains to be elucidated [2]. Thus we chose not to include this behavior in our model. Once the phosphorylated TtrR dimer binds to the promoter region, the promoter becomes activated (Figure 1 step 7), allowing RNA polymerase to bind. RNA polymerase binding is followed by the transcription and translation of ARG, modeled likewise to the transcription and translation of the two-component system (Figure 1 step 8-9).
Finally, the expression and assembly of the gas vesicle proteins are modeled (Figure 1 step 10). The ARG gene codes for 12 different gas vesicle proteins, which may all play a role in GV formation. However, the process of GV formation is not fully elucidated yet. Therefore, we chose to model GV formation as a self-assembling process, where a sufficient amount of gas vesicle proteins will lead to an irreversible assembly of a gas vesicle. It is known that gas vesicle protein A (GvpA) is the major constituent of the gas vesicle wall and that gas vesicle protein C (GvpC) is often located at the exterior surface of the GV, thereby strengthening the entire structure [3,4]. However, only GvpA is essential for the formation [4], therefore, we chose to model the GV formation solely through the assembly of GvpA. In each ARG gene, two repeats of the GvpA protein are present, and it is estimated that 46 GvpA proteins are needed to form a GV with an average size of about 200 nm [3]. Thus, we model that a GV is formed when 46 GvpA proteins are present.
As rate parameters were not available in the literature for all reactions, some assumptions were made. Most of the rates are retrieved from literature, corresponding to that exact, or an equivalent chemical reaction. Rates that could not be retrieved from literature were estimated. These rates were determined by varying them and thereby fitting the simulation to the data retrieved from literature and the lab. In the model, all proteins and mRNA transcripts are assumed to degrade with a certain degradation rate. This degradation rate is identical for all proteins and identical for all mRNA transcripts. All mass-action chemical reactions, rate constants, and initial concentrations that were used in the model can be found in the tables below.
During the development of the model, we went through several iterations of the design cycle. You can read all about this at Engineering Success, segment Computational Modeling. The complete model can be found on our Github repository.
$ LacI + P_{lac}^* \overset{k_1^b}{\underset{k_1^u}{\rightleftharpoons}} LacI:P_{lac} $
$LacI + 2 IPTG \overset{k_2^b}{\underset{k_2^u}{\rightleftharpoons}} LacI:IPTG $
$LacI:P_{lac} + 2 IPTG \overset{k_3^b}{\underset{k_3^u}{\rightleftharpoons}} LacI:IPTG + P_{lac}^*$
$P_{lac}^* + RNAP \overset{k_4^b}{\underset{k_4^u}{\rightleftharpoons}} P_{lac}^*:RNAP $
$P_{lac}^*:RNAP \xrightarrow{k_{tx}} T_{TtrR} + T_{TtrS} + P_{lac}^* + RNAP $
$T_{TtrR} + R \overset{k_5^b}{\underset{k_5^u}{\rightleftharpoons}} T_{TtrR}:R $
$T_{TtrR}:R \xrightarrow{k_{tl}} T_{TtrR} + TtrR + R $
$T_{TtrS} + R \overset{k_6^b}{\underset{k_6^u}{\rightleftharpoons}} T_{TtrS}:R $
$T_{TtrS}:R \xrightarrow{k_{tl}} T_{TtrS} + TtrS + R $
$TtrS + Ttr \overset{k_7^b}{\underset{k_7^u}{\rightleftharpoons}} TtrS_p + Ttr $
$TtrS_p \xrightarrow{k_{dephos}} TtrS $
$TtrR + TtrS \overset{k_8^b}{\underset{k_8^u}{\rightleftharpoons}} TtrR:TtrS $
$TtrR + TtrS_p \overset{k_9^b}{\underset{k_9^u}{\rightleftharpoons}} TtrR:TtrS_p $
$TtrR:TtrS_p \overset{k_{10}^b}{\underset{k_{10}^u}{\rightleftharpoons}} TtrR_p + TtrS $
$TtrR_p \xrightarrow{k_{dephos}} TtrR $
$2 TtrR_p \overset{k_{11}^b}{\underset{k_{11}^u}{\rightleftharpoons}} TtrR_p:TtrR_p$
$TtrR_p:TtrR_p + P_{ttrB} \overset{k_{12}^b}{\underset{k_{12}^u}{\rightleftharpoons}} P_{ttrB}^* $
$P_{ttrB}^* + RNAP \overset{k_{13}^b}{\underset{k_{13}^u}{\rightleftharpoons}} P_{ttrB}^*:RNAP $
$P_{ttrB}^*:RNAP \xrightarrow{k_{tx}} P_{ttrB}* + RNAP + T_{ARG} $
$T_{ARG} + R \overset{k_{14}^b}{\underset{k_{14}^u}{\rightleftharpoons}} T_{ARG}:R$
$T_{ARG}:R \xrightarrow{k_{tl}} T_{ARG} + R + Gvp $
Degradation Reactions
$T_{ARG} \xrightarrow{k_{deg1}} \emptyset$
$T_{TtrS} \xrightarrow{k_{deg1}} \emptyset$
$T_{TtrR} \xrightarrow{k_{deg1}} \emptyset$
$TtrR \xrightarrow{k_{deg2}} \emptyset$
$TtrS \xrightarrow{k_{deg2}} \emptyset$
$Gvp \xrightarrow{k_{deg2}} \emptyset$
$TtrS_p \xrightarrow{k_{deg2}} \emptyset$
$TtrR_p \xrightarrow{k_{deg2}} \emptyset$
$TtrR_p:TtrR_p \xrightarrow{k_{deg2}} \emptyset$
$TtrS:TtrR \xrightarrow{k_{deg2}} \emptyset$
$TtrS_p:TtrR \xrightarrow{k_{deg2}} \emptyset$
-
Stamatakis M, Mantzaris N. Comparison of Deterministic and Stochastic Models of the lac Operon Genetic Network. Biophysical Journal. 2009;96(3):887-906.
-
Merk L, Shur A, Pandey A, Murray R, Green L. Engineering Logical Inflammation Sensing Circuit for Gut Modulation. 2020;.
-
Offner S, Ziese U, Wanner G, Typke D, Pfeifer F. Structural characteristics of halobacterial gas vesicles. Microbiology. 1998;144(5):1331-1342.
-
Hill A, Salmond G. Microbial gas vesicles as nanotechnology tools: exploiting intracellular organelles for translational utility in biotechnology, medicine and the environment. Microbiology. 2020;166(6):501-509.
-
Arkin A, Ross J, McAdams H. Stochastic Kinetic Analysis of Developmental Pathway Bifurcation in Phage λ-Infected Escherichia coli Cells. Genetics. 1998;149(4):1633-1648.
-
Stamatakis M, Mantzaris N. Comparison of Deterministic and Stochastic Models of the lac Operon Genetic Network. Biophysical Journal. 2009;96(3):887-906.
-
Újvári A, Martin C. Thermodynamic and Kinetic Measurements of Promoter Binding by T7 RNA Polymerase. Biochemistry. 1996;35(46):14574-14582.
-
Kierzek A, Zaim J, Zielenkiewicz P. The Effect of Transcription and Translation Initiation Frequencies on the Stochastic Fluctuations in Prokaryotic Gene Expression. Journal of Biological Chemistry. 2001;276(11):8165-8171.
-
Gao R, Stock A. Quantitative Kinetic Analyses of Shutting Off a Two-Component System. mBio. 2017;8(3).
-
Elf J, Li G, Xie X. Probing Transcription Factor Dynamics at the Single-Molecule Level in a Living Cell. Science. 2007;316(5828):1191-1194.
-
Wong O, Guthold M, Erie D, Gelles J. Interconvertible Lac Repressor–DNA Loops Revealed by Single-Molecule Experiments. PLoS Biology. 2008;6(9):e232.
-
Bernstein J, Khodursky A, Lin P, Lin-Chao S, Cohen S. Global analysis of mRNA decay and abundance in Escherichia coli at single-gene resolution using two-color fluorescent DNA microarrays. Proceedings of the National Academy of Sciences. 2002;99(15):9697-9702.
-
Porter S, Armitage J. Phosphotransfer in Rhodobacter sphaeroides Chemotaxis. Journal of Molecular Biology. 2002;324(1):35-45.
The Results
Introduction
By varying the dose of tetrathionate (Ttr) that is added to the model, we were able to create dose-response curves that describe the relationship between the exposure to Ttr and the magnitude of protein expression after a set exposure time. The exposure time was established by determining the time point at which the steady-state in protein expression was reached, which was at t=2500 min. We created two dose-response curves: one for the expression of GFP (Figure 1a), and one for the expression of gas vesicles (Figure 1b).
Figure 1: a) Dose-response curve of GFP expression in response to varying doses of Ttr. b) Dose-response curve of gas vesicle expression in response to varying doses of Ttr
When comparing the dose-response curve of GFP upon tetrathionate induction, which was created using the first model version, to the dose-response curve that was obtained in the lab, as well as to that from literature, we can observe great similarity between the curves. The experimentally defined curve shows a dose-response between 0.01 and 0.1 mM, which is similar to the results of Daeffler et al. [1], and a half-maximal effective concentration (EC50), the concentration that induces a response halfway between the baseline and maximum [2], of 0.0491 ± 0.0023 mM. The model-derived curve shows slightly less sensitivity, with a dose-response between 0.01 and 1 mM, and an EC50 of 0.2098 ± 0.004 mM. Further, we calculated and compared the Hill coefficient, which is a measure for the cooperativity in a binding process. The Hill coefficient calculated with the model nearly matches the experimentally obtained Hill coefficient by Daeffler et al. [1], yielding 1.199 ± 0.025 and 1.5 ± 0.1 respectively.
The difference in dose-response is likely due to inaccuracies in our rate parameters, or due to simplifications that were made in the model. Because of the correspondence of our model results with the lab results, it can be validated that our model is accurate, meaning the extended model (including the addition of gas vesicle assembly) can also be assumed to be accurate. This extended model is, however, expected to show a similar deviation in gas vesicles expression compared to emperically established results in vitro.
The model including gas vesicle assembly yields a dose-response curve with the number of gas vesicles for different inducer concentrations. The EC50 was calculated to be 0.2099 ± 0.004 mM. Unfortunately, little is known on the exact tetrathionate concentration in the gut during active inflammations, thus further research needs to be done before we can relate important metrics such as Limit-of-Detection and confirm the efficacy of our system in vivo. Further, the Hill coefficient of this curve is equal to 1.200 ± 0.025. The Hill coefficient being greater than 1.0 indicates positive cooperativity between tetrathionate and TtrS.
Figure 2: Simulation of the gas vesicle expression over time at a Ttr dose of 0.2099 mM
Using the EC50 as a dose, because pathological concentrations of tetrathionate were not published in literature, we simulated the gas vesicle expression over time to determine the optimal time for ultrasound measurement after ingestion of the pill. This simulation is shown in Figure 2. The curve illustrates that gas vesicle production keeps increasing for about 25 hours after induction with Ttr. A steady-state production is reached after 25 hours, at a count of 116 gas vesicles per E. coli cell. Bourdeau et al. [1] have shown that an average of 100 gas vesicles per cell can result in an adequate ultrasound signal in vivo, although this is also dependent on the number of cells that are measured. From this information, we can conclude that after 25 hours, sufficient gas vesicles are present for a proper ultrasound measurement in vivo at a Ttr concentration of 0.2099 mM. This suggests that our system produces enough signal in humans for proper ultrasound pictures.
-
Bourdeau R, Lee-Gosselin A, Lakshmanan A, Farhadi A, Kumar S, Nety S et al. Acoustic reporter genes for noninvasive imaging of microorganisms in mammalian hosts. Nature. 2018;553(7686):86-90.
-
Neubig R, Spedding M, Kenakin T, Christopoulos A. International Union of Pharmacology Committee on Receptor Nomenclature and Drug Classification. XXXVIII. Update on Terms and Symbols in Quantitative Pharmacology. Pharmacological Reviews. 2003;55(4):597-606.