Team:Manchester/Model/PSO
Phosphate Starvation Operon
Achievements:
1. Successful compilation of parameters and ODEs for modelling response of native and recombinant phosphate starvation (PHO) regulon
2. Proved that the regulatory system response is sufficiently fast for fluctuating phosphate level
3. Predicted that the regulatory system is efficient enough to keep microcompartment synthesis to minimal during long period of starvation, as required to be implemented in real world scenario
4. Provided a framework for comparison of future experimental data
Introduction
The synthesis of microcompartments creates a massive metabolic burden that makes the bacteria “sick” and is very detrimental to cell growth, as observed during our experimental work. This would ultimately interfere with the real-life implementation of our bacteria in a continuous culture system.
In a continuous culture system it is important to balance biomass accumulation with PPK expression and phosphate accumulation. When phosphate levels are low, the cells should invest their resources into growth, rather than phosphate accumulation. To alleviate this problem, we designed and modelled a regulatory operon that aims to regulate the synthesis of microcompartment proteins under different external phosphate levels.
The successful implementation of our phosphate recovery system requires the external phosphate level to exceed a certain threshold. Once phosphate is above the threshold the bacteria will be begin to produce microcompartments to enhance phosphate accumulation. The trigger threshold would ideally be the concentration of phosphate deemed high in the water regulations. We inquired with regulatory bodies and found out that concentrations of phosphate above 1 mg/L (0.01 mM) would be unacceptable for treated effluent.
Thus, when the water is “clean” and the concentration of phosphate in the water dropped under the threshold, the bacteria would stop producing microcompartments and would be able to grow at wild type rates. Meanwhile, the response of the regulon has to be rapid enough for the synthesis of microcompartment to be turned on as the phosphate level in the water rises.
The envisaged design of the regulon involves the native Pho regulon in E. Coli, a recombinant Pho operon that encodes for TetR protein and a Tet operon that encodes for the microcompartment proteins. For simplicity we chose to have the microcompartment proteins controlled by one control circuit, rather than the three present in our Eut operon.
The regulation works as follows:
1. PhoR (R) binds with surface phosphate (Ps) to form repression complex (RR) in high phosphate level environment. TetR is not produced and EutS is produced as normal.
2. In low phosphate level, RR dissociates and releases some free PhoR (R).
3. Rfree rapidly autophosphorylates to form activated PhoR (RA).
4. RA phosphorylates PhoB to form Rfree and activated PhoB (BA).
5. BA binds to Pho operon to upregulate the transcription of PhoB, PhoR, tetR and AP.
6. TetR binds to Tet operon to repress the transcription of EutS mRNA, stopping the synthesis of microcompartments.
Note that BA also upregulates PhoB and PhoR, which makes the response to phosphate level intense and sharp. A detailed explanation of the Pho regulon was described in this paper(1).
Response of the native regulatory system where the product is alkaline phosphatase was simulated using Matlab and is depicted as below. This would be used to demonstrate the desired response pattern of TetR and EutS relative to native AP regulation.
Methods
The process of modelling was begun by researching on the native Pho regulon and the phosphate starvation response. The species involved were identified and the rate laws of each step in the system were developed. After that, efforts were put into designing and integrating operons that encode for tetR and eutS.
In order to simplify the modelling, a few assumptions were made. It was assumed that the transcription and translation rates were the same. The model does not include the pstSCAB-phoU operon in order to simplify the regulatory process. However, it was assumed that the formation of repression complex which would normally involve Ps, phoR, pstS and phoU to be dependent on Ps and 3 times the concentration of phoR, which would be reasonable because the concentrations of phoR, pstS and phoU are under the same regulation pathway.
| Species | Abbreviation |
|---|---|
| Surface phopshate | Ps |
| mRNA of phoR | mRNAR |
| mRNA of phoB | mRNAB |
| mRNA of TetR | mRNATetR |
| mRNA of EutS | mRNAEutS |
| mRNA of Alkaline Phosphatase | mRNAAP |
| phoR protein | R |
| phoB protein | B |
| TetR protein | TetR |
| EutS protein | EutS |
| Alkaline Phosphatase protein | AP |
| phoR and phosphate repression complex | RR |
| Phosphorylated phoR protein | RA |
| Phosphorylated phoB protein | BA |
| Parameters | Description |
|---|---|
| kTB | Basal transcription rate constant for Pho promoter |
| kTPB | Basal transcription rate constant for Pho promoter on plasmid |
| kTRB | Transcription rate constant for repressed Tet promoter |
| kTA | Activated transcription rate for Pho promoter |
| kH | Protein translation rate constant |
| kdmRNA | mRNA degradation rate constant |
| kdH | Protein degradation rate constant |
| KMpro | Dissociation constant for phoBA with Pho promoter |
| konRR | Association rate constant for repression complex |
| koffRR | Dissociation rate constant for repression complex |
| KI | Dissociation constant for repression complex |
| Ktet | Dissociation constant for TetR with Tet promoter |
| kra | Autophosphorylation rate of phoR |
| kba | Forward rate constant of phosphoryl transfer |
| krba | Reverse rate constant of phosphoryl transfer |
| kdba | Autodephosphorylation rate constant of phoBA |
| KMproP | Dissociation constant for phoBA with Pho promoter on plasmid |
BA available to each promoter is estimated as
\begin{equation*} [BA_{\textrm{each}}] = \frac{1}{4} \bigg\{[BA_{\textrm{total}}] - [BA_{\textrm{boundtoDNA}}]\bigg\} \end{equation*} \begin{equation*} [BA_{\textrm{boundtoDNA}}] = \frac{[DNA][BA_{\textrm{total}}]}{K_{MDNA}} = \frac{1}{3}[BA_{\textrm{total}}] \end{equation*}where
\begin{equation*} [DNA] = 6.64 \textrm{mM}, K_{MDNA} > 20 \textrm{mM} \end{equation*} \begin{equation*} [BA_{\textrm{each}}] = \frac{1}{4} \cdot \frac{2}{3}[BA_{\textrm{total}}]=0.167[BA_{\textrm{total}}] \end{equation*}Rate Laws:
Transcription of TetR
\begin{equation*} r_1 = k_{TA} \cdot \frac{0.167 [BA]}{K_{MproP} + 0.167 [BA]} + k_{TPB} \end{equation*}Degradation of TetR mRNA
\begin{equation*} r_2 = k_{dmRNA} [mRNA_{TetR}] \end{equation*}Degradation of R mRNA
\begin{equation*} r_3 = k_{dmRNA} [mRNA_{R}] \end{equation*}Translation of TetR
\begin{equation*} r_4 = k_{H} [mRNA_{TetR}] \end{equation*}Translation of R
\begin{equation*} r_5 = k_{H} [mRNA_{R}] \end{equation*}Degradation of TetR
\begin{equation*} r_6 = k_{dH} [TetR] \end{equation*}Degradation of R
\begin{equation*} r_7 = k_{dH} [R] \end{equation*}Degradation of RR
\begin{equation*} r_8 = k_{dH} [RR] \end{equation*}Transcription of EutS mRNA
\begin{equation*} r_9 = k_{TA} \cdot \frac{K_{Tet}}{K_{Tet} + [TetR]} + k_{TRB} \end{equation*}Degradation of EutS mRNA
\begin{equation*} r_{10} = k_{dmRNA} [mRNA_{EutS}] \end{equation*}Translation of EutS
\begin{equation*} r_{11} = k_{H} [mRNA_{EutS}] \end{equation*}Degradation of EutS
\begin{equation*} r_{12} = k_{dH} [EutS] \end{equation*}Dissociation of RR ⇌ R (+ pstS + phoU) + Ps
\begin{equation*} r_{13} = k_{offRR}[RR]-k_{onRR}[R]^3[Ps] \quad \textrm{where} \quad K_I = \frac{k_{offRR}}{k_{onRR}} \end{equation*}Transcription of AP, B and R mRNA
\begin{equation*} r_{14} = k_{TA} \cdot \frac{0.167 [BA]}{K_{Mpro} + 0.167 [BA]} + k_{TB} \end{equation*}Degradation of AP mRNA
\begin{equation*} r_{15} = k_{dmRNA} [mRNA_{AP}] \end{equation*}Translation of AP
\begin{equation*} r_{16} = k_{H} [mRNA_{AP}] \end{equation*}Degradation of AP
\begin{equation*} r_{17} = k_{dH} [AP] \end{equation*}Rapid Autophosphorylation of R
\begin{equation*} r_{18} = k_{ra} [R] \end{equation*}Degradation of RA
\begin{equation*} r_{19} = k_{dH} [RA] \end{equation*}Degradation of B mRNA
\begin{equation*} r_{20} = k_{dmRNA} [mRNA_{B}] \end{equation*}Translation of B
\begin{equation*} r_{21} = k_{H} [mRNA_{B}] \end{equation*}Degradation of B
\begin{equation*} r_{22} = k_{dH} [B] \end{equation*}Phosphorylation of B + RA ⇌ BA + R
\begin{equation*} r_{23} = k_{ba}[B][RA]-k_{rba}[BA][R] \end{equation*}Autodephosphorylation of BA
\begin{equation*} r_{24} = k_{dba}[BA] \end{equation*}Degradation of BA
\begin{equation*} r_{25} = k_{dH} [BA] \end{equation*}Ordinary Differential Equations (ODEs)
\begin{equation} \frac{d[R]}{dt} = r_{13}+r_5-r_7-r_{18}+r_{23} \end{equation} \begin{equation} \frac{d[RR]}{dt} = -r_{13}-r_8 \end{equation} \begin{equation} \frac{d[mRNA_{R}]}{dt} = r_{14}-r_3 \end{equation} \begin{equation} \frac{d[mRNA_{TetR}]}{dt} = r_1-r_2 \end{equation} \begin{equation} \frac{d[TetR]}{dt}= r_4-r_6 \end{equation} \begin{equation} \frac{d[mRNA_{EutS}]}{dt} = r_9-r_{10} \end{equation} \begin{equation} \frac{d[EutS]}{dt} = r_{11}-r_{12} \end{equation} \begin{equation} \frac{d[mRNA_{AP}]}{dt} = r_{14}-r_{15} \end{equation} \begin{equation} \frac{d[AP]}{dt} = r_{16}-r_{17} \end{equation} \begin{equation} \frac{d[RA]}{dt} = r_{18}-r_{19}-r_{23} \end{equation} \begin{equation} \frac{d[mRNA_{B}]}{dt} = r_{14}-r_{20} \end{equation} \begin{equation} \frac{d[B]}{dt} = r_{21}-r_{22}-r_{23}+r_{24} \end{equation} \begin{equation} \frac{d[BA]}{dt} = r_{23}-r_{24}-r_{25} \end{equation}The modelling of biological systems requires an accurate knowledge of the pathways and information on all kinetic parameters in order to make a reasonable prediction on the behaviour of the system.(2,3,4) However, the collection of parameters come with large uncertainties. Parameters reported from parameter fitting approaches suffer from the diverse experimental conditions due to the nature of extremely complex biological systems, e.g. measurements in vivo or in vitro, experimental noise, technical limitations, etc.
The issue of these uncertainties can be dealt through the application of ensemble modelling where reasonable sets of parameter values (“ensemble”) are sampled, based on probability distributions of likely values obtained from literature. This allows the predictions to accurately describe a population of all models, which are consistent with current knowledge about the system.
Ensemble modelling has been used with great success in many fields such as economics and engineering, but it has only been making its way into the field of synthetic biology with the help of the development of toolkits and methods. Although ensemble modelling is relatively novel method in SynBio, the consideration of uncertainties in well-defined probability distributions facilitates the modelling of extremely complicated biological systems where many parameters were obtained from various sources. Therefore, we are excited to incorporate ensemble modelling into prediction of our regulatory system. Rather than ignoring parameter variation, we are able to address the issue through modelling.
In order to maintain thermodynamic consistency during sampling of interdependent parameters, which in the case of our model are KI, konRR and koffRR, a multivariate distribution was used to generate an ensemble of parameters which are constrained by both uncertainties in the individual parameter and the relationship KI = koffRR / konRR. This extra effort was advantageous in that it generates the set of parameters interdependently, which fulfil the thermodynamic constraints imposed by the relationship(5) and this increases the ratio of successful and plausible iterations (plausibility filter explained below).
Ensemble modelling was done using parameters found from literature and weightings were assigned to the values based on the protocol described in the paper “Defining informative priors for ensemble modelling in systems biology”(5). Probability distributions of plausible values for each parameter were defined using log-normal distributions. Weightings of the parameters with their uncertainties and the summary of the distribution of parameters are documented in detail in appendix II. The aforementioned 13 ordinary differential equations (ODEs) are integrated in Matlab using ode15s solver which is suitable for stiff equations.
1000 iterations (simulations of the behaviour of the regulatory circuit using independent parameter samples) were run for each ensemble simulation. The results were filtered for:
1. no negative values for any of the species
2. TetR, EutS and AP concentration fluctuate
Results passing this quality / plausibility filter were plotted using median with 25th to 75th quantile -- the variability in results (confidence interval) is depicted in shades of colour.
Results
As Figure 2 suggests, the designed regulatory system works as expected in a constantly fluctuating phosphate concentration environment. The plot clearly shows that as the surface phosphate (Ps) concentration drops, the tetR concentration increases and consequently represses the transcription of EutS mRNA, leading to a drop in the concentration of EutS protein.
Bearing in mind that in a real-world scenario the phosphate level would not fluctuate rapidly, another simulation was run for a very long time with a low phosphate level to allow the system to reach equilibrium and then with a spike in phosphate for a considerable length of time. This would simulate how the Eut expression remains dormant and is activated following exposure of the system to high levels of phosphate.
The result of the simulation showed that given enough time at an elevated phosphate level, the increase in concentration of the microcompartment proteins would be more substantial.
Looking at this regulatory system from an alternative perspective, figure below shows how the regulatory system can be represented as a Boolean logic gate circuit. A highly regulated system would behave like a bistable switch (yes or no, 1 or 0) for the production of the end product, which in this case, the microcompartment proteins.
Conclusion
The phosphate starvation regulatory system was modelled by building from native Pho regulon and integrating recombinant operons which encode tetR and microcompartment proteins. A total of 13 ODEs were built with 17 kinetic parameters. Not only were uncertainties in the parameters considered by building an ensemble of parameters, but thermodynamic consistency between interdependent parameters was also maintained using a multivariate distribution.
The response pattern in the designed regulatory system was as expected. At low phosphate levels, the host organism (E. Coli) will remain dormant and will not produce microcompartments. Soon after the phosphate level rose past the threshold (1 mg/L), tetR level in the cell drops and microcompartments will be produced, enhancing the cell’s capability to accumulate polyphosphate. This regulatory system would be crucial for the implementation of our work on elevated polyphosphate accumulation by PPK encapsulation in microcompartment in real world scenario.
Further effort is needed to fit the parameters to get a minimal response at 4 uM and full response at 0.4 uM as well as to get [AP] = 25 uM when Ps = 0.4 uM as reported in experimental measurements from literature(6). Until then, it is not possible to relate the plots above to phosphate level data obtained from wastewater treatment plant, and to manipulate the maximal response of microcompartment synthesis at phosphate concentration above 0.01mM as mentioned in the introduction.
It is possible to fine-tune the regulatory system to be sharp and rapid. There are numerous ways to achieve this, such as a mutation of the pho operon on the plasmid for tighter regulation of tetR transcription or the choice of tet operon. Experimental validation of this system will certainly yield exciting results, which will further shape this model. Also, the pho operon could be interesting in mass production of desired proteins because the induction would be by phosphate starvation instead of other chemical species.
Acknowledgement
Grateful shout-out to Rainer’s group:
Aliah Hawari and especially Areti Tsigkinopoulou for helping extensively with the modelling from the start to the end!
References
- Wanner, B. (1993). Gene regulation by phosphate in enteric bacteria. Journal of Cellular Biochemistry, 51(1), pp.47-54.
- Achcar, F., Barrett, M. and Breitling, R. (2013). Explicit consideration of topological and parameter uncertainty gives new insights into a well-established model of glycolysis. FEBS Journal, 280(18), pp.4640-4651.
- Achcar, F., Kerkhoven, E., Bakker, B., Barrett, M. and Breitling, R. (2012). Dynamic Modelling under Uncertainty: The Case of Trypanosoma brucei Energy Metabolism. PLoS Computational Biology, 8(1), p.e1002352.
- Tsigkinopoulou, A., Baker, S. and Breitling, R. (2017). Respectful Modeling: Addressing Uncertainty in Dynamic System Models for Molecular Biology. Trends in Biotechnology, 35(6), pp.518-529.
- Tsigkinopoulou, A., Hawari, A., Uttley, M. and Breitling, R. (2017). Defining informative priors for ensemble modelling in systems biology. Unpublished manuscript.
- Van Dien, S. and Keasling, J. (1998). A Dynamic Model of the Escherichia coli Phosphate-Starvation Response. Journal of Theoretical Biology, 190(1), pp.37-49.
- Bernstein, J., Khodursky, A., Lin, P., Lin-Chao, S. and Cohen, S. (2002). 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, 99(15), pp.9697-9702.
- Maurizi, M. (1992). Proteases and protein degradation in Escherichia coli. Experientia, 48(2), pp.178-201.
- McCleary, W. (1996). The activation of PhoB by acetylphosphate. Molecular Microbiology, 20(6), pp.1155-1163.
- Kamionka, A. (2004). Two mutations in the tetracycline repressor change the inducer anhydrotetracycline to a corepressor. Nucleic Acids Research, 32(2), pp.842-847.
- Nelson, H. and Sauer, R. (1985). Lambda repressor mutations that increase the affinity and specificity of operator binding. Cell, 42(2), pp.549-558.
Appendix I
Appendix II
| Parameters | k values | error | weightings | Multiplicative(1)/Additive(0) | Source |
|---|---|---|---|---|---|
| TB | 6.41E-07 | NaN | 16 | 0 | (6) |
| TPB | 2.50E-07 | NaN | 16 | 0 | (6) |
| TRB | 8.20E-15 | NaN | 16 | 0 | assumed |
| TA *adjustable | 8.20E-02 | 4 | 16 | 1 | (6) |
| kH | 13.66 | NaN | 16 | 0 | (6) |
| kdRNA | 0.46 | NaN | 16 | 0 | (6) |
| 0.133 | NaN | 8 | 0 | (7) | |
| 0.122 | NaN | 8 | 0 | (7) | |
| kdH | 0.07 | NaN | 8 | 0 | (6) |
| 0.002 | NaN | 8 | 0 | (8) | |
| KMpro | 4.90E-07 | 8 | 32 | 1 | (9) |
| konRR | 1.00E+11 | 5 | 8 | 1 | assumed |
| koffRR | 4.30E-10 | 5 | 8 | 1 | assumed |
| KI *adjustable | 4.30E-23 | NaN | 32 | 0 | (6) |
| Ktet | 1.79E-07 | 6.39E-08 | 16 | 1 | (10) |
| 7.85E-07 | NaN | 4 | 0 | (11) | |
| 1.00E-03 | NaN | 16 | 0 | assumed for lower affinity mutation | |
| kra | 1.00E+05 | NaN | 32 | 0 | (6) |
| kba | 759 | NaN | 32 | 0 | (6) |
| krba | 759 | NaN | 32 | 0 | (6) |
| kdba | 18 | NaN | 32 | 0 | assumed |
| KMproP | 4.90E+00 | 8 | 32 | 1 | assumed |
| Parameter | Mode | CIfact | mu | sigma | Xmin | Xmax |
|---|---|---|---|---|---|---|
| TB | 6.38E-07 | 1.20999857 | -14.2557754 | 0.094883956 | 5.27E-07 | 7.72E-07 |
| TPB | 2.49E-07 | 1.20999857 | -15.197344 | 0.094883956 | 2.06E-07 | 3.01E-07 |
| TRB | 8.16E-15 | 1.20999857 | -32.4301813 | 0.094883956 | 6.75E-15 | 9.88E-15 |
| TA | 0.00168065 | 264.20671 | -3.62271232 | 1.663088749 | 6.36E-06 | 0.444038186 |
| kH | 13.5980968 | 1.20999857 | 2.6189328 | 0.094883956 | 11.2381098 | 16.4536776 |
| kdRNA | 0.27737535 | 3.66764267 | -0.95825833 | 0.569320055 | 0.075627691 | 1.017313653 |
| kdH | 0.03311486 | 35.1793226 | -1.91354971 | 1.222384352 | 0.000941316 | 1.164958275 |
| KMpro | 0.00030012 | 592.626217 | -4.80159551 | 1.819267643 | 5.06E-07 | 0.177858607 |
| konRR | 1.00E+11 | 24.9995002 | 26.6234751 | 1.137997851 | 4000079971 | 2.50E+12 |
| koffRR | 4.30E-10 | 24.9995002 | -20.2721968 | 1.137997851 | 1.72E-11 | 1.07E-08 |
| KI | 4.28E-23 | 1.20999857 | -51.4963812 | 0.094883956 | 3.54E-23 | 5.18E-23 |
| Ktet | 7.82E-07 | 4398.23255 | -9.34379245 | 2.172062414 | 1.78E-10 | 0.003438605 |
| kra | 99546.8284 | 1.20999857 | 11.5173864 | 0.094883956 | 82270.20353 | 120451.5198 |
| kba | 755.560428 | 1.20999857 | 6.63646273 | 0.094883956 | 624.4308448 | 914.2270351 |
| krba | 755.560428 | 1.20999857 | 6.63646273 | 0.094883956 | 624.4308448 | 914.2270351 |
| kdba | 17.9184291 | 1.20999857 | 2.89483271 | 0.094883956 | 14.80863664 | 21.68127356 |
| KMproP | 2.55932919 | 9.77122603 | 1.72192189 | 0.884407543 | 0.261925083 | 25.00778398 |
