Team:Lambert GA/Model
OBJECTIVE
While researching potential parts for a phosphate biosensor, Lambert iGEM discovered the Pho Regulon, a mechanism that naturally exists in E. coli to manage inorganic phosphate (Pi) levels within cells. It is regulated by a phosphate-sensitive signaling pathway that responds to levels of extracellular Pi to transcribe its regulatory genes [1]. Based on The Pho regulon and the pathogenesis of Escherichia coli [2], the Pho Regulon pathway has a threshold of 4uM Pi. Concentrations greater than 4uM inhibit expression of regulatory genes, while concentrations less than 4uM do not.
Part BBa_K2447000, created by the 2017 NUS Singapore iGEM team, utilizes the Pho Regulon signaling pathway to control GFP expression of the BioBrick under varying phosphate concentrations. NUS Singapore improved the part to increase its sensitivity, causing the part to be able to detect ranges of phosphate concentrations significantly greater than 4uM.
Before cloning and characterizing the part as a phosphate biosensor, Lambert iGEM developed a deterministic Ordinary Differential Equation model to simulate the biosensor’s activity over time and explore whether the 4uM threshold still limits the biosensor. Rather than a stochastic or exploratory model, a deterministic ODE model can predict specific concentrations of species within the model based on parameter values and initial conditions.
Lambert iGEM’s model creates ordinary differential equations of all the reactions of the signaling pathway in a single E. coli cell and evaluates them with an initial input of extracellular phosphate concentration to predict GFP expression.
MODEL DEVELOPMENT
During a virtual meeting with Styczynski Research Group from Georgia Institute of Technology, Dr. Mark Styczynski advised the team on initial steps for creating the model. He recommended searching for rate constants for the reactions of the signaling pathway in existing literature as the best method for modeling.
Over the summer, Lambert iGEM had virtual biosensor collaboration meetings with the Queens and University of Rochester iGEM teams to troubleshoot each others’ biosensor project ideas. The University of Rochester team gave input on the model development and suggested splitting the pathway into modules - modeling them separately - and then connecting the pieces, since the Pho Regulon signaling pathway is so complex.
DIAGRAM OF SIGNALING PATHWAY
For building the model, the team followed the process outlined in A Tutorial on Mathematical Modeling of Biological Signaling Pathways [3]. First, using MATLAB Simbiology software, Lambert iGEM diagrammed each reaction of the signaling pathway, from the entering of Pi into the cell to final GFP expression.
Figure 1. Diagram of the reactions of the Pho Regulon signaling pathway in a single E. coli cell, created in MATLAB Simbiology software.
BIOCHEMICAL REACTIONS
The 22 biochemical reactions in the model are shown as follows:
| Number | Reaction | Parameter(s) | Description |
|---|---|---|---|
| 1 | [outer membrane].PstS + [inner membrane].Pi + [outer membrane].Pi -> [outer membrane].PstS_Pi | kf*[outer membrane].PstS*[inner membrane].Pi*[outer membrane].Pi | Pi (inorganic phosphate) binds to the transport protein PstS, forming the bound Pst_Pi complex. |
| 2 | [outer membrane].PstS_Pi -> [outer membrane].PstS + [outer membrane].Pi | kf*[outer membrane].PstS_Pi | After transporting Pi across the periplasmic space, PstS dissociates from Pi. |
| 3 | [outer membrane].[PstA/Cp] + [outer membrane].Pi -> [inner membrane].Pi + [outer membrane].[PstA/C] | kf*[outer membrane].[PstA/Cp]*[outer membrane].Pi | The phosphorylated transmembrane channel PstA/C actively transports Pi across the inner membrane. |
| 4 | [inner membrane].Pi + [inner membrane].PhoU -> [inner membrane].PhoU_Pi | kf*[inner membrane].Pi*[inner membrane].PhoU | Pi binds to PhoU in the cytoplasm, forming the bound PhoU_Pi complex. |
| 5 | [inner membrane].Pi + [inner membrane].PhoU -> [inner membrane].PhoU_Pi | kf*[inner membrane].Pi*[inner membrane].PhoU | PhoU and Pi dissociate. |
| 6 | [inner membrane].PhoU_Pi + [inner membrane].PstB <-> [inner membrane].[inactive PstB] | kf*[inner membrane].PhoU_Pi*[inner membrane].PstB - kr*[inner membrane].[inactive PstB] | The bound PhoU_Pi complex inhibits and therefore inactivates the permease PstB in a reversible reaction. |
| 7 | [outer membrane].[PstA/C] + [inner membrane].PstB -> [outer membrane].[PstA/Cp] | kf*[outer membrane].[PstA/C]*[inner membrane].PstB | Active PstB phosphorylates the transmembrane channel PstA/C. |
| 8 | [inner membrane].PhoU_Pi + [inner membrane].DiPhoR -> [inner membrane].[stabilized PhoR] | kf*[inner membrane].PhoU_Pi*[inner membrane].DiPhoR | The bound PhoU_Pi complex stabilizes the histidine kinase PhoR. |
| 9 | [inner membrane].DiPhoR <-> [inner membrane].DiPhoRp | kf*[inner membrane].DiPhoR - kr*[inner membrane].DiPhoRp | The PhoR dimer autophosphorylates in a reversible reaction. |
| 10 | [inner membrane].DiPhoRp <-> [inner membrane].DiPhoRpp | kf*[inner membrane].DiPhoRp - kr*[inner membrane].DiPhoRpp | PhoR autophosphorylates again in a reversible reaction. |
| 11 | [inner membrane].DiPhoRpp + [inner membrane].PhoB <-> [inner membrane].DiPhoRpp_PhoB | kf*[inner membrane].DiPhoRpp*[inner membrane].PhoB - kr*[inner membrane].DiPhoRpp_PhoB | PhoR binds to the transcription factor PhoB in a reversible reaction. |
| 12 | [inner membrane].DiPhoRpp_PhoB -> [inner membrane].PhoBp + [inner membrane].DiPhoRp | kf*[inner membrane].DiPhoRpp_PhoB | PhoR phosphorylates PhoB, and they dissociate. |
| 13 | [inner membrane].PhoB + [inner membrane].DiPhoRp <-> [inner membrane].DiPhoRp_PhoB | kf*[inner membrane].PhoB*[inner membrane].DiPhoRp - kr*[inner membrane].DiPhoRp_PhoB | PhoR binds to a different PhoB protein in a reversible reaction. |
| 14 | [inner membrane].DiPhoRp_PhoB -> [inner membrane].DiPhoR + [inner membrane].PhoBp | kf*[inner membrane].DiPhoRp_PhoB | PhoR phosphorylates PhoB, and they dissociate. PhoR is now back to its original unphosphorylated dimer. |
| 15 | [inner membrane].DiPhoR + [inner membrane].PhoBp <-> [inner membrane].DiPhoR_PhoBp | kf*[inner membrane].DiPhoR*[inner membrane].PhoBp - kr*[inner membrane].DiPhoR_PhoBp | PhoR binds to phosphorylated PhoB in a reversible reaction. |
| 16 | [inner membrane].DiPhoR_PhoBp -> [inner membrane].DiPhoR + [inner membrane].PhoB | kf*[inner membrane].DiPhoR_PhoBp | PhoR, acting as a bifunctional enzyme, dephosphorylates PhoB, and they dissociate. |
| 17 | [inner membrane].PhoBp <-> [inner membrane].DiPhoBpp | kf*[inner membrane].PhoBp - kr*[inner membrane].DiPhoBpp | Phosphorylated PhoB dimerizes in a reversible reaction. |
| 18 | [inner membrane].DiPhoBpp + [inner membrane].pPhoB <-> [inner membrane].[active pPhoB] | kf*[inner membrane].DiPhoBpp*[inner membrane].pPhoB - kr*[inner membrane].[active pPhoB] | PhoB binds to the promoter of the Pho Regulon genes, activating the promoter. |
| 19 | [inner membrane].[active pPhoB] -> [inner membrane].mRNA + [inner membrane].[active pPhoB] | kf*[inner membrane].[active pPhoB] | Transcription of the downstream genes occurs. |
| 20 | [inner membrane].mRNA -> null | kf*[inner membrane].mRNA | mRNA degradation occurs. |
| 21 | [inner membrane].mRNA -> [inner membrane].GFP + [inner membrane].mRNA | [translation_GFP.kf]*[inner membrane].mRNA | Translation of GFP occurs. |
| 22 | [inner membrane].GFP -> null | kf*[inner membrane].GFP | GFP degradation occurs. |
Table 1. Model reactions of the Pho Regulon signaling pathway. Reactions 9-17 were obtained from the article "Quantifying dynamic mechanisms of auto-regulation in Escherichia coli with synthetic promoter in response to varying external phosphate levels" [4], which modeled the signaling pathway from PhoR activity through gene expression.
INITIAL VALUES OF SPECIES
The initial value of each species is shown below. Protein abundance data was gathered from Datanator [5], as well as from Quantifying dynamic mechanisms of auto-regulation in Escherichia coli with synthetic promoter in response to varying external phosphate levels [4].
| Species | Initial Value |
|---|---|
| PstS | 481 |
| Pi | 0 micromoles |
| PstS_Pi | 0 |
| PstA/C | 197 |
| PstA/Cp | 0 |
| DiPhoR | 24 |
| DiPhoRp | 0 |
| DiPhoRpp | 4.0000e-8 |
| PhoB | 35 |
| DiPhoRpp_PhoB | 0 |
| pPhoB | 1 |
| Active pPhoB | 0 |
| mRNA | 0 |
| PstB | 250 |
| GFP | 0 |
| PhoU | 126 |
| PhoU_Pi | 0 |
| Inactive PstB | 0 |
| Stabilized PhoR | 0 |
Table 2. Initial values of model species.
ODE EQUATIONS
Lambert iGEM determined that the Mass Action kinetic law was the most appropriate for describing the reactions of the model [3], and generated 24 ODE equations of the reactions, as shown below:
| ODE Equations |
|---|
| d([outer membrane].Pi)/dt = -([PstS_Pi binding].kf*PstS*[inner membrane].Pi*[outer membrane].Pi) + ([PstS_Pi dissociation].kf*PstS_Pi) - ([Pi transport].kf*[PstA/Cp]*[outer membrane].Pi) |
| d(PstS)/dt = 1/[outer membrane]*(-([PstS_Pi binding].kf*PstS*[inner membrane].Pi*[outer membrane].Pi) + ([PstS_Pi dissociation].kf*PstS_Pi)) |
| d(PstS_Pi)/dt = 1/[outer membrane]*(([PstS_Pi binding].kf*PstS*[inner membrane].Pi*[outer membrane].Pi) - ([PstS_Pi dissociation].kf*PstS_Pi)) |
| d([PstA/C])/dt = 1/[outer membrane]*(-([PstA/C phosphorylation].kf*[PstA/C]*PstB) + ([Pi transport].kf*[PstA/Cp]*[outer membrane].Pi)) |
| d(PstB)/dt = 1/[inner membrane]*(-([PstA/C phosphorylation].kf*[PstA/C]*PstB) - ([PstB inhibition].kf*PhoU_Pi*PstB-[PstB inhibition].kr*[inactive PstB])) |
| d([PstA/Cp])/dt = 1/[outer membrane]*(([PstA/C phosphorylation].kf*[PstA/C]*PstB) - ([Pi transport].kf*[PstA/Cp]*[outer membrane].Pi)) |
| d([inner membrane].Pi)/dt = -([PstS_Pi binding].kf*PstS*[inner membrane].Pi*[outer membrane].Pi) + ([Pi transport].kf*[PstA/Cp]*[outer membrane].Pi) - ([PhoU_Pi binding].kf*[inner membrane].Pi*PhoU) + ([PhoU_Pi dissociation].kf*PhoU_Pi) |
| d(PhoU)/dt = 1/[inner membrane]*(-([PhoU_Pi binding].kf*[inner membrane].Pi*PhoU) + ([PhoU_Pi dissociation].kf*PhoU_Pi)) |
| d(PhoU_Pi)/dt = 1/[inner membrane]*(([PhoU_Pi binding].kf*[inner membrane].Pi*PhoU) - ([PhoR stabilization].kf*PhoU_Pi*DiPhoR) - ([PstB inhibition].kf*PhoU_Pi*PstB-[PstB inhibition].kr*[inactive PstB]) - ([PhoU_Pi dissociation].kf*PhoU_Pi)) |
| d([inactive PstB])/dt = 1/[inner membrane]*(([PstB inhibition].kf*PhoU_Pi*PstB-[PstB inhibition].kr*[inactive PstB])) |
| d([stabilized PhoR])/dt = 1/[inner membrane]*(([PhoR stabilization].kf*PhoU_Pi*DiPhoR)) |
| d(DiPhoR)/dt = 1/[inner membrane]*(-(([autophosphorylation(1)].kf*DiPhoR)*[inner membrane]-([autophosphorylation(1)].kr*DiPhoRp)*[inner membrane]) + (([PhoB_phosphorylation(2)].kf*DiPhoRp_PhoB)*[inner membrane]) - ([DiPhoR_PhoBp binding].kf*DiPhoR*PhoBp-[DiPhoR_PhoBp binding].kr*DiPhoR_PhoBp) + ((PhoB_dephosphorylation.kf*DiPhoR_PhoBp)*[inner membrane]) - ([PhoR stabilization].kf*PhoU_Pi*DiPhoR)) |
| d(DiPhoRp)/dt = 1/[inner membrane]*((([autophosphorylation(1)].kf*DiPhoR)*[inner membrane]-([autophosphorylation(1)].kr*DiPhoRp)*[inner membrane]) - (([autophosphorylation(2)].kf*DiPhoRp)*[inner membrane]-([autophosphorylation(2)].kr*DiPhoRpp)*[inner membrane]) + (([PhoB_phosphorylation(1)].kf*DiPhoRpp_PhoB)*[inner membrane]) - ([DiPhoRp_PhoB binding].kf*PhoB*DiPhoRp-[DiPhoRp_PhoB binding].kr*DiPhoRp_PhoB)) |
| d(DiPhoRpp)/dt = 1/[inner membrane]*((([autophosphorylation(2)].kf*DiPhoRp)*[inner membrane]-([autophosphorylation(2)].kr*DiPhoRpp)*[inner membrane]) - ([DiPhoRpp_PhoB binding].kf*DiPhoRpp*PhoB-[DiPhoRpp_PhoB binding].kr*DiPhoRpp_PhoB)) |
| d(PhoB)/dt = 1/[inner membrane]*(-([DiPhoRpp_PhoB binding].kf*DiPhoRpp*PhoB-[DiPhoRpp_PhoB binding].kr*DiPhoRpp_PhoB) - ([DiPhoRp_PhoB binding].kf*PhoB*DiPhoRp-[DiPhoRp_PhoB binding].kr*DiPhoRp_PhoB) + ((PhoB_dephosphorylation.kf*DiPhoR_PhoBp)*[inner membrane])) |
| d(DiPhoRpp_PhoB)/dt = 1/[inner membrane]*(([DiPhoRpp_PhoB binding].kf*DiPhoRpp*PhoB-[DiPhoRpp_PhoB binding].kr*DiPhoRpp_PhoB) - (([PhoB_phosphorylation(1)].kf*DiPhoRpp_PhoB)*[inner membrane])) |
| d(PhoBp)/dt = 1/[inner membrane]*((([PhoB_phosphorylation(1)].kf*DiPhoRpp_PhoB)*[inner membrane]) + (([PhoB_phosphorylation(2)].kf*DiPhoRp_PhoB)*[inner membrane]) - ([PhoB dimerisation].kf*PhoBp-[PhoB dimerisation].kr*DiPhoBpp) - ([DiPhoR_PhoBp binding].kf*DiPhoR*PhoBp-[DiPhoR_PhoBp binding].kr*DiPhoR_PhoBp)) |
| d(DiPhoRp_PhoB)/dt = 1/[inner membrane]*(([DiPhoRp_PhoB binding].kf*PhoB*DiPhoRp-[DiPhoRp_PhoB binding].kr*DiPhoRp_PhoB) - (([PhoB_phosphorylation(2)].kf*DiPhoRp_PhoB)*[inner membrane])) |
| d(DiPhoBpp)/dt = 1/[inner membrane]*(([PhoB dimerisation].kf*PhoBp-[PhoB dimerisation].kr*DiPhoBpp) - ([DiPhoBpp_pPhoB binding].kf*DiPhoBpp*pPhoB-[DiPhoBpp_pPhoB binding].kr*[active pPhoB])) |
| d(DiPhoR_PhoBp)/dt = 1/[inner membrane]*(([DiPhoR_PhoBp binding].kf*DiPhoR*PhoBp-[DiPhoR_PhoBp binding].kr*DiPhoR_PhoBp) - ((PhoB_dephosphorylation.kf*DiPhoR_PhoBp)*[inner membrane])) |
| d(pPhoB)/dt = 1/[inner membrane]*(-([DiPhoBpp_pPhoB binding].kf*DiPhoBpp*pPhoB-[DiPhoBpp_pPhoB binding].kr*[active pPhoB])) |
| d([active pPhoB])/dt = 1/[inner membrane]*(([DiPhoBpp_pPhoB binding].kf*DiPhoBpp*pPhoB-[DiPhoBpp_pPhoB binding].kr*[active pPhoB])) |
| d(mRNA)/dt = 1/[inner membrane]*(((transcription.kf*[active pPhoB])*[inner membrane]) - ((mRNA_degradation.kf*mRNA)*[inner membrane])) |
| d(GFP)/dt = 1/[inner membrane]*(([translation_GFP.kf]*mRNA) - (GFP_degradation.kf*GFP)) |
Table 3. Model ODE equations. These ODE equations describe the concentration of each species in the model over time.
RATE CONSTANTS
Next, for reaction rate constants, the same article [4] provided some of the values, but Lambert iGEM still lacked many of the constants for the other reactions. However, research to find the missing constants was inconclusive, as the metabolomics of the Pho Regulon are incomplete. The team turned to parameter estimation, using Simbiology Model Analyzer to fit NUS Singapore’s characterization data and estimate the missing parameters, with guidance from Ms. Yue Han of Styczynski Research Group. Since the characterization data simply reflects phosphate concentration to final GFP expression with no data for GFP expression over time, it was difficult to estimate a global set of parameters that would simulate GFP expression dependent on phosphate concentration. Instead, Lambert iGEM selected one characterization data point (one phosphate concentration and one GFP expression value) at a time to estimate parameters with a local solver, Isqnonlin, and then analyzed the relationship between each set of parameters. To confirm that the estimation was accurate, Lambert iGEM also estimated the same parameters described in the article Quantifying dynamic mechanisms of auto-regulation in Escherichia coli with synthetic promoter in response to varying external phosphate levels [4], and these estimations not only fell within the stated physiological ranges, but also matched the given values closely. Values for parameters across each set were uniform, except for GFP translation.
Phosphate concentration to each corresponding GFP translation rate constant was plotted, and an equation for the GFP translation rate constant dependent on phosphate concentration was derived.
[translation_GFP.kf] = -121517.16788599*log([outer membrane].Pi) + 778114.76297209
The rest of the estimated rate constants are shown below:
| Reaction Number | Name(s) | Estimated Value |
|---|---|---|
| 1 | [PstS_Pi binding].kf | 1 |
| 2 | [PstS_Pi dissociation].kf | 1 |
| 3 | [Pi transport].kf | 1 |
| 4 | [PhoU_Pi binding].kf | 1 |
| 5 | [PhoU_Pi dissociation].kf | 1 |
| 6 | [PstB inhibition].kf | 1 |
| [PstB inhibition].kr | 1 | |
| 7 | [PstA/C phosphorylation].kf | 1.000147861 |
| 8 | [PhoR stabilization].kf | 0.999999972 |
| 9 | [autophosphorylation(1)].kf | 25.33326162 |
| [autophosphorylation(1)].kr | 8.118317422 | |
| 10 | [autophosphorylation(2)].kf | 25.36529935 |
| [autophosphorylation(2)].kr | 8.116903976 | |
| 11 | [DiPhoRpp_PhoB binding].kf | 100.0309931 |
| [DiPhoRpp_PhoB binding].kr | 44.9040607 | |
| 12 | [PhoB_phosphorylation(1)].kf | 21.37675216 |
| 13 | [DiPhoRp_PhoB binding].kf | 100.0030047 |
| [DiPhoRp_PhoB binding].kr | 94.94854646 | |
| 14 | [PhoB_phosphorylation(2)].kf | 21.43778755 |
| 15 | [DiPhoR_PhoBp binding].kf | 100.0360146 |
| [DiPhoR_PhoBp binding].kr | 34.99807536 | |
| 16 | [PhoB_dephosphorylation].kf | 12.97103369 |
| 17 | [PhoB dimerisation].kf | 100.1002524 |
| [PhoB dimerisation].kr | 24.89215776 | |
| 18 | [DiPhoBpp_pPhoB binding].kf | 10000.00234 |
| [DiPhoBpp_pPhoB binding].kr | 1000.000968 | |
| 19 | [transcription].kf | 0.13 |
| 20 | [mRNA_degradation].kf | 0.004067274 |
| 21 | [translation_GFP.kf] | See equation above |
| 22 | [GFP_degradation].kf | 0.00968 |
Table 4. Estimated rate constants for each biochemical reaction.
ASSUMPTIONS
The model assumes that all rate constants are constant and are not affected by environmental factors except for GFP translation. It also assumes that the input phosphate concentration remains constant throughout simulation of GFP expression. For simplicity, the only gene expressed at the end of the pathway is GFP, while in E. coli, regulatory genes of the Pho Regulon signaling pathway are also expressed.
RESULTS
With complete ODE equations, Lambert iGEM used Simbiology Model Analyzer to simulate GFP expression over time with different inputs of phosphate concentrations - inputs used by NUS Singapore in their team’s characterization, as well as inputs that would be used by Lambert iGEM for new characterization.
Figure 2. Graph of the model’s simulated GFP expression to phosphate concentration.
The model confirms that phosphate concentrations above the 4uM threshold described by The Pho regulon and the pathogenesis of Escherichia coli [2] will still result in levels of GFP expression. Therefore, cloning part BBa_K2447000 for use as a phosphate biosensor would be sensitive enough to detect the range of phosphate concentrations commonly found in aquaponics systems: 0uM to 100uM.
The model also follows the same log relationship as the original characterization data from 2017 NUS Singapore, shown below:
Figure 3. Graph of the original characterization data from 2017 NUS Singapore.
Compared to NUS Singapore, Lambert iGEM’s model estimates smaller values of GFP expression because the model represents the GFP expression of one single E. coli cell, whereas NUS Singapore’s data is from a culture of cells.
COMPARISON TO LAMBERT IGEM’S CHARACTERIZATION DATA
Having created a model of the Pho Regulon signaling pathway and simulated its activity, Lambert iGEM successfully cloned part BBa_K2447000 and added characterization with a range of phosphate concentrations found in aquaponics systems, as shown below:
Figure 4. Lambert iGEM’s characterization data of Part BBa_K2447000.
Lambert iGEM’s characterization data matches the model even more closely than NUS Singapore’s characterization data.
REFERENCES
[1] Santos-Beneit, F. (2015). The Pho regulon: a huge regulatory network in bacteria. Frontiers in Microbiology, 6. https://doi:10.3389/fmicb.2015.00402.
[2] Crépin, S., Chekabab, S., Bihan, G. L., Bertrand, N., Dozois, C. M., & Harel, J. (2011). The Pho regulon and the pathogenesis of Escherichia coli. Veterinary Microbiology, 153(1-2), 82-88. https://doi:10.1016/j.vetmic.2011.05.043.
[3] Zi, Z. (2012). A Tutorial on Mathematical Modeling of Biological Signaling Pathways. Methods in Molecular Biology Computational Modeling of Signaling Networks, 41-51. https://doi:10.1007/978-1-61779-833-7_3.
[4] Uluşeker, C., Torres-Bacete, J., García, J. L., Hanczyc, M. M., Nogales, J., & Kahramanoğulları, O. (2019). Quantifying dynamic mechanisms of auto-regulation in Escherichia coli with synthetic promoter in response to varying external phosphate levels. Scientific Reports, 9(1). https://doi:10.1038/s41598-018-38223-w.
[5] Karr Lab, I. (n.d.). Data for modeling and simulating cells. Retrieved October 25, 2020, from https://datanator.info/.