Model
Explain your model's assumptions, data, parameters, and results in a way that anyone could understand.
Explain your model's assumptions, data, parameters, and results in a way that anyone could understand.
After conducting laboratory experiments with hlFAB-GFP and receiving a negative result, specifically no GFP expression in the presence of PFOA, we sought to understand the molecular interactions at play between FAB and PFOA. Our goal was to determine if modifications to hlFAB could enhance its binding affinity for PFOA, thereby improving its lower detection limit. To achieve this, we employed computational approaches to simulate the binding process, evaluate its effectiveness, and explore potential optimizations.
We initiated the process by using AutoDock Vina, a high-throughput molecular docking tool, to model the interaction between PFOA and the binding domain of hlFAB. After docking, we selected the best pose—based on binding energy and spatial orientation—for further investigation. This docked complex was then imported into Amber for molecular dynamics (MD) simulations. The purpose of these simulations was to validate that PFOA remained stable in the FAB binding domain, confirming a successful binding interaction. After validating the complex, we employed the Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) method to quantitatively assess the strength of the binding. MMPBSA works by separating the ligand (PFOA) from the receptor (FAB) and calculating the binding free energy (ΔG). A more negative ΔG indicates a stronger binding affinity. The MMPBSA method is grounded in the following free energy equation:
ΔG_bind = ΔE_vdw + ΔE_elec + ΔG_solv - TΔS
Where each term represents the van der Waals energy, electrostatic energy, solvation free energy, and the entropy contribution, respectively. By leveraging a custom Python script, we were able to further analyze charge contributions to the binding interactions, helping us to pinpoint key residues driving the binding process. Additionally, determining the ΔG allowed us to accurately calculate the dissociation constant (Kd) and predict the rate constants for binding kinetics, enhancing our understanding of FAB's performance.
In our next step, we sought to optimize the FAB binding by introducing mutations. Using ChimeraX's rotamer tool, we performed rational mutations of key residues to explore whether these changes could improve the binding affinity for PFOA. We also created a computational pipeline to automate the MMPBSA workflow, streamlining the process of evaluating binding affinities for various FAB mutations and PFOA. This pipeline significantly enhanced our ability to iterate on the design of FAB mutants and evaluate their performance efficiently.
WTo create the computational model of the hlFAB-PFOA complex, we utilized AlphaFold to predict the 3D structure of the FAB-GFP fusion protein, based on the amino acid sequence published in Dr. Berger’s research. The structure of PFOA was retrieved from its PubChem entry as an SDF file, and the atomic charges for PFOA were calculated using Antechamber, a widely accepted tool for small molecule parameterization. This step ensured that the ligand was accurately represented in our simulations.
For the docking stage, we selected the entire FAB binding domain to explore all possible binding orientations using AutoDock Vina. After the docking results were generated, the top-ranked pose was extracted and split into three files: the receptor (FAB), the ligand (PFOA), and the complex (FAB-PFOA). These files were then fed into our MMPBSA pipeline.
The pipeline we developed includes multiple stages to ensure accurate free energy calculations. These stages include:
Topology Creation: LEaP was used to parameterize the receptor-ligand complex.
Energy Minimization: A minimization step was conducted to relax any unfavorable steric clashes in the system.
Heating:The complex was heated to 300K, mimicking physiological temperature.
Density Equilibration:The system’s density was equilibrated to ensure accurate simulation dynamics.
Equilibration:Additional steps ensured the system was properly equilibrated.
Production Run: A 10 ns MD production run was performed to gather trajectory data.
MMPBSA Analysis: The MMPBSA method was applied to calculate the binding free energy.
Data Curation: A custom Python script was used to organize the raw data into a format that could be easily analyzed.
We used standard parameters for most stages to ensure that the protocol is reproducible and could be applied to future studies involving other ligands or receptors. For more in-depth explanation of the pipeline, check the methods in the Dry-Lab Experiments page.
Mann MM, Berger BW. A genetically-encoded biosensor for direct detection of perfluorooctanoic acid. Sci Rep. 2023 Sep 13;13(1):15186. doi: 10.1038/s41598-023-41953-1. PMID: 37704644; PMCID: PMC10499884.
The initial docking of hlFAB with PFOA was successful, providing a solid foundation for further analysis. Using Autodock Vina, the docking procedure predicted an initial binding free energy (ΔG) of -11.3016 kcal/mol for the wild-type hlFAB bound to PFOA. This value reflects the strength of the interaction between the protein and the ligand, where a more negative ΔG indicates a stronger, more favorable binding.
To better understand the significance of this result, we calculated the dissociation constant (Kd) from the ΔG using the following equation:
ΔG = RT * ln(Kd)
Where:
R, ideal gas constant, (1.987 * 10^-3 cal/mol·K)
Temperature (300K)
ΔG is the binding free energy (-11.3016 Kcal/mol)
By rearranging the equation to solve for ΔG:
e(ΔG / RT) = Kd
Substituting in the values, we get:
5.836 nano MolarThis calculation gives a dissociation constant (Kd) of approximately 5.836 nM, which indicates strong binding. Since Kd means the concentration required of the ligand to bind to half of all the receptors, a low Kd, especially in the nanomolar range, signifies that the ligand (PFOA) binds tightly to the protein (hlFAB), making it an effective candidate for detection in our experimental work.
After running the first MMPBSA, we generated a graph from the python script, in order to see the charge decomposition:
We used this graph to determine likely residues to mutate.
To track the impact of various mutations on the binding affinity, we created a table comparing the ΔG values and calculated Kd for each mutation. These results show how each mutation influences the interaction between hlFAB and PFOA, and whether the changes enhance or weaken binding strength.
| Mutation | ΔG (kcal/mol) | Kd (M) | Effect on Binding Strength | Why we chose this mutation |
|---|---|---|---|---|
| Wild type + PFOA | -11.3016 | 5.8359 x 10^-9 | N/A | Baseline |
| Wild type + Palmitic Acid | -2.4944 | 0.0152 | N/A | To compare PFOA to it’s natural ligand. |
| ILE 52 to ARG + PFOA | -10.0965 | 4.4065 x 10^-8 | Slightly worse | Due to this residue having a very small contribution (thin slice of the pie graph) to the bind, changing this residue to a positively charged resiude (ARG) would increase the binding strength to PFOA, since PFOA has a net charge of -1. |
| ILE 308 to ARG + PFOA | -13.1878 | 2.4654 x 10^-10 | Great improvement. | Due to this residue having a very small contribution (thin slice of the pie graph) to the bind, changing this residue to a positively charged resiude (ARG) would increase the binding strength to PFOA, since PFOA has a net charge of -1. |
| SER 349 to THR + PFOA | -3.1394 | 0.0051 | Severly worse | In the binding pocket, There is a Serine and a Threonine that interact with the
oxygens on PFOA, and they form hydrogen bonds between FAB and PFOA. So the idea was to make the two
residues that interact with PFOA the same.
|
| PHE 50 to ARG + PFOA | -6.9105 | 9.2318 x 10^-6 | Severly worse | Due to this residue having a very small contribution (thin slice of the pie graph) to the bind, changing this residue to a positively charged resiude (ARG) would increase the binding strength to PFOA, since PFOA has a net charge of -1. |
| THR 351 to SER + PFOA | -12.2177 | 1.2550 x 10^-9 | Slightly worse | In the binding pocket, There is a Serine and a Threonine that interact with the
oxygens on PFOA, and they form hydrogen bonds between FAB and PFOA. So the idea was to make the two
residues that interact with PFOA the same.
|
This table organizes the results of each mutation, highlighting how the mutations affect the protein's ability to bind with PFOA. For example, the ILE 308 mutation shows an improvement in binding, while the SER 349 mutation significantly weakens the interaction. These findings provide clear direction on which mutations could improve the lower detection limit of hlFAB, and helps identify key residues such as SER 349 and THR 351.
To get an idea of how our system would play out and determine the effect size of different perturbations on the system, we modeled our genetic constructs in the Virtual Cell software using the Gibson-Bruck stochastic solver. The stochastic solver is more accurate for nanoscale reactions than deterministic differential equations, and Virtual Cell provides a view-friendly user interface for designing reaction networks.
Whenever possible, we extracted rate constants from existing literature. Due to inconclusive experimental results, all our rate constants were extracted from literature or estimated. All our models are publicly available on VCell and are also saved to our team’s software gitlab. You will need to download Virtual Cell to see our models. See our Engineering page for more details on how we constructed and refined our models, and see our Contributions page for how to access our models.
| Reaction Name | Kf Value | Kr Value | Annotations |
|---|---|---|---|
| PFAS_DiffusionAndMessenger | 0.001 s-1 | 0.001 s-1 | https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8913905/ it is not well known if PFAS interact with proteins for cell transport. in absence of better evidence, we assume it acts by simple diffusion |
| PFAS_Bind_Syntran | 1.0 s-1.μM-1 | 7.195E-5 s-1 | https://www.sciencedirect.com/science/article/pii/S016041201932759X Kd is 71.95 uM kr/kf=71.95*10^-6 if kf=1 kr=71.95*10^-6 https://pubs.acs.org/doi/abs/10.1021/es304030x this says kd=107 ??M about the same |
| Syntran_Degr | 9.62704417444E-6 s-1 | 0 s-1 | https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6660218/ protein half lives range from 2-20 hours in e. coli. since our e coli are likely in the stationary phase, their protein turnover rate is much lower and therefore is closer to 20 hours. k=ln(2)/(20*60*60)=9.62704417444E-6 |
| SynTRan_Production | 0.00155709342561 s-1 | 0 s-1.μM-1 | https://bionumbers.hms.harvard.edu/bionumber.aspx?s=y&id=100059&ver=26&hlid=56413 12 for slower rates synthetic transcription factor is 1734/3 aa (12.0 / (1734.0 / 3.0)) https://www.nature.com/articles/s41467-023-41176-y at slow growth rates elongation can fall to 9 aa so 9 /(1734/3)=0.0155709342561 also says translation initiation time is about 10 s so .1*0.0155709342561=0.00155709342561 |
| Production | 0.00763650248 s-1 | 0 s-1.μM-1 | https://book.bionumbers.org/what-is-faster-transcription-or-translation/ transcription 40-80 nt/sec. syntrans is 1734 nt rbs is 12 40/1746=0.0229095074456 https://bionumbers.hms.harvard.edu/bionumber.aspx?id=111997&ver=5 initiation about 20 genes per minute so 1/3 gene per second 0.0229095074456/3=0.00763650248187 |
| degrade_synt_pfas_cplx | 9.62704417444E-6 s-1 | 0 s-1.μM-1 | same as protein degradation |
| mRNA_Degr | 3.85E-4 s-1 | 0 s-1 | https://genomebiology.biomedcentral.com/articles/10.1186/s13059-022-02811-x https://www.frontiersin.org/journals/microbiology/articles/10.3389/fmicb.2020.02111/full https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6660218/ Half life of mRNA is about 0.5 hour, so k=ln(2)/(0.5*60*60) |
| Syntran_bind_pLex | 4 s-1.μM-1 | 1 s-1 | https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5873372/ In figure 4 our promoter has strength of about 80%. assuming that full strength means all DNA is taken, you could roughly approximate K as 80/20 = 4/1=kf/kr |
| Leak | 4.553734E-4 s-1 | 0 s-1.μM-1 | https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5873372/ inducibility is 40 times leak, so rate of leak transcription is roughly 1/40 of induced 0.01821493624/40=0.0004553734 |
| make GFP_mRNA | 0.01821493624 s-1 | 0 s-1.μM-1 | https://book.bionumbers.org/what-is-faster-transcription-or-translation/ transcription 40-80 nt/sec. gfp + rbs is 732 nt 40/732=0.0546448087432 https://bionumbers.hms.harvard.edu/bionumber.aspx?id=111997&ver=5 initiation about 20 genes per minute so 1/3 gene per second 0.0546448087432/3=0.01821493624 |
| Degr_GFP_mRNA | 0.01821493624 s-1 | 0 s-1.μM-1 | https://genomebiology.biomedcentral.com/articles/10.1186/s13059-022-02811-x https://www.frontiersin.org/journals/microbiology/articles/10.3389/fmicb.2020.02111/full https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6660218/ Half life of mRNA is about 0.5 hour, so k=ln(2)/(0.5*60*60) |
| SF_GFP_degr | 1.83E-5 s-1 | 0 s-1 | Was halflife of 2 hrs but substituted for 4.15/hr or 0.00115/s From https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4162443/ |
| Reaction Name | Kf Value | Kr Value | Annotations |
|---|---|---|---|
| PFAS_DiffusionAndMessenger | 0.001 s-1 | 0.001 s-1 | https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8913905/ it is not well known if PFAS interact with proteins for cell transport. in absence of better evidence, we assume it acts by simple diffusion (probably slowly) |
| Fab_GFP_Production | 0.0116652085156 s-1 | 0 s-1.μM-1 | https://book.bionumbers.org/what-is-faster-transcription-or-translation/ transcription 40-80 nt/sec. fabgfp is 1143 nt with rbs 40/1143=0.0349956255468 https://bionumbers.hms.harvard.edu/bionumber.aspx?id=111997&ver=5 initiation about 20 genes per minute so 1/3 gene per second 0.0349956255468/3=0.0116652085156 |
| FAB_GFP_mRNA_degredation | 3.85E-4 s-1 | 0 s-1 | https://genomebiology.biomedcentral.com/articles/10.1186/s13059-022-02811-x https://www.frontiersin.org/journals/microbiology/articles/10.3389/fmicb.2020.02111/full https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6660218/ Half life of mRNA is about 0.5 hour, so k=ln(2)/(0.5*60*60) |
| FAB_GFP_mRNA translation | 0.00238726790451 s-1 | 0 s-1.μM-1 | https://bionumbers.hms.harvard.edu/bionumber.aspx?s=y&id=100059&ver=26&hlid=56413 12nt/s for slower rates https://www.nature.com/articles/s41467-023-41176-y at slow growth rates elongation can fall to 9 aa fab gfp is 377 aa so 9/377=0.0238726790451 above paper also says translation initiation time is about 10 s so =0.00238726790451 |
| FAB_GFP_Protein_degr | 9.62704417444E-6 s-1 | 0 s-1 | https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6660218/ protein half lives range from 2-20 hours in e. coli. since our e coli are likely in the stationary phase, their protein turnover rate is much lower and therefore is closer to 20 hours. k=ln(2)/(20*60*60)=9.62704417444E-6 |
| PFAS_and_FAB_GFP_ bind | 192.4 s-1.μM-1 | 9.99E-4 s-1 | calculated using graph and table from paper: https://www.nature.com/articles/s41598-023-41953-1 |
| PFAS_FAB_GFP_degr | 9.62704417444E-6 s-1 | 0 s-1.μM-1 | https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6660218/ protein half lives range from 2-20 hours in e. coli. since our e coli are likely in the stationary phase, their protein turnover rate is much lower and therefore is closer to 20 hours. k=ln(2)/(20*60*60)=9.62704417444E-6 |
| Reaction Name | Kf Value | Kr Value | Annotations |
|---|---|---|---|
| PFAS_DiffusionAndMessenger | 0.1 s-1 | 0.1 s-1 | This is a rough first approximation of the unknown mechanisms for activating prmA |
| PFAS_bind_pRMA | 1 s-1.μM-1 | 0.01 s-1 | Taken from PFAS_Detectorv4_pillow, original had no source |
| Make_GFP_mRNA | 50⁄1200 s-1 | 0 s-1.μM-1 | estimated 643 bp, transcription rate of 10bp/sec |
| pRMA_Leak | 50⁄3600 s-1 | 0 s-1.μM-1 | Leakiness rate constant taken from PFAS Detectorv4_pillow, original had no source |
| Degr_GFP_mRNA | 3.85E-4 s-1 | 0 s-1 | https://genomebiology.biomedcentral.com/articles/10.1186/s13059-022-02811-x https://www.frontiersin.org/journals/microbiology/articles/10.3389/fmicb.2020.02111/full https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6660218// Half life of mRNA is about 0.5 hour, so k=ln(2)/(0.5*60*60) |
| Make_SF_GFP | 0.04201680672 s-1 | 0 s-1.μM-1 | From the paper https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4162443/ GFP protein production from gene to immature (?) GFP with active promoter was 150/hr and another 3.6 /hr to get to the mature version I will use 150/hr = 150/hr *1hr/3600 sec= 0.04166/sec |
| SF_GFP_degr | 1.83E-5 s-1 | 0 s-1 | https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6660218/ https://www.mdpi.com/2073-4409/12/14/1846 |
| SF_GFP_dimerization | 0.83333 s-1 | 0.01666 s-1 | https://pubmed.ncbi.nlm.nih.gov/18601531/ |
| SF_GFP_dimer_degr | 3.33333333333E-5 s-1 | 0 s-1 | Rate constant used from PFAS Detectorv4_pillow, from similar dimerization degradation reaction, original had no source |
| Construct Name | Biomodel Name | Creator of Biomodel | Short Description |
|---|---|---|---|
| Prma_GFP Construct | pRMA_GFP construct 1 v1.3 | Pillow123 | This construct explored the use of the pRMA promoter in response to PFAS to produce SF_GFP |
| FAB_GFP Construct | Construct 3_Cycle3_dgl | dglVcell | This construct explored the FAB protein mechanism in expressing FAB_GFP in response to PFAS |
| Synthetic_Transcription_Factor Construct | FAB_GFP_Construct_v1.1_Douglas | dglVcell | In this construct we modeled the dynamics of GFP under the control of a potentially PFAS inducible synthetic transcription factor |
This simulates GFP production when under the control of the prmA promoter, which has been shown to upregulate transcription in the presence of PFAS. This is a model of our part BBa_K5114823. Because the mechanism of the prmA promoter is poorly understood, we approximated its response to PFAS as the DNA binding to PFAS molecules which increases its rate of transcription. This model also takes into account the leakiness of the promoter. Experimental evidence suggests PFAS presence only increases expression of prmA controlled proteins by about 3 fold.
This graph demonstrates the dynamics of the simple genetic circuit. The presence of GFP in the absence of PFAS is due to the basal production rate of the promoter. As shown, the amount of GFP by the end of the simulation increases as the amount of PFAS added increases. There appears to be a significant increase in the amount of GFP produced compared to basal production in the presence of as little as 0.01 uM (micromolar) of PFAS. As the amount of initial PFAS increases, the increase in GFP production slows down, likely due to saturation of DNA with PFAS. Thus, the highest curves represent maximal GFP production rate.
This model served as our first experience with modeling genetic circuits to detect PFAS, however we were skeptical that prmA would work in our chassis of E. coli because the transcriptional machinery that regulate the prmA promoter appears to be unique to Rhodococcus, as we found in a BLAST search. Additionally, chemicals very rarely bind directly to DNA to modulate transcriptional activity, so this model should not be taken as entirely true.
hlFAB-GFP is a protein made of circularly permuted GFP and human liver fatty acid binding protein, created by Mann and Berger in 2023. It acts as a PFAS sensor by fluorescence upon binding with PFAS. When no PFAS is bound, water can infiltrate the hydrophobic barrel of GFP that shields its chromophore and quench the fluorescence. When PFAS is bound to the hlFAB domain, less water can infiltrate and fluorescence is stronger. This models our part BBa_K5114228.
Simulations were carried out for 100 minutes (6000) seconds at various concentrations of initial PFAS. The environment was made to be 1000 um^3 in volume because that is the expected amount of water “available” to each individual cell (max density of E. coli is around 10^9 cells/ml). Cell volume was set to 1 um^3 based on established cell sizes of E. coli at stationary phase. Simulations were carried out at non-steady-state and steady-state (concentration of FAB_GFP and FAB_GFP_mRNA are stable) conditions, found by deterministic ODE modeling. Regulations for PFAS are on the parts per trillion (ppt) level, so we used 1 ppt as our minimum concentration. Assuming 1 ppt=1 ng/L and the molar weight of PFOA, a common model PFAS, is 414 g/mol, that means 1 ppt is approximately 2E-6 uM of PFOA. Therefore, we simulated from 1E-6 uM initial PFAS in the environment up to 1E-2 uM, which would correlate to roughly 5 ppb.
Interestingly, steady state modeling does not yield significant differences in maximal fluorescence. This is likely because the fluorescence is actually limited by the amount of PFAS available to bind to hlFAB-GFP. This can be seen in the graphs below. Since the environment is set to be 1000 times larger in volume than the cell, the final concentration of bound hlFAB-GFP is 1000 times the initial concentration of PFAS, indicating nearly every molecule of PFAS is bound to a hlFAB-GFP.
From 1E-6 to 1E-4, the stochastic nature of the simulation is evident: each step up on the graph represents
one PFAS molecule diffusing into the cell from the environment and binding to one molecule of hlFAB-GFP.
Thus, our model indicates that biosensors that directly fluoresce upon binding to PFAS are limited by the
amount of PFAS available to them, and not necessarily by the binding affinity of the protein. Future work
should focus on increasing the effective fluorescence of the molecule.
A synthetic estradiol transcription factor (synTrans or STF for short) was originally created by fusing the hormone binding domain of human estrogen receptor alpha with the DNA binding domain of the Lex transcription factor and the Herpes-Simplex Virus protein VP16 to recruit ribosomes. The STF has been successfully used in yeast as a gratuitous transcription factor. Since PFAS has been shown to be an agonist for human estrogen receptor, we attempted to express it in E. coli to see if we could use the STF as a transcription factor inducible by PFAS. This is a model for our parts TBD and TBD.
Cell and environment volumes were kept at 1 and 1000 um^3, respectively. Initial PFAS concentrations were also kept the same. Simulations were carried out in ideal settings that were not at steady state and had no promoter leakage (expression of GFP without the STF binding to the DNA), as well as in more realistic conditions where there was promoter leakage and the concentrations of the STF and its mRNA were constant.
Compared to hlFAB-GFP, this construct appears to be much more sensitive to small quantities of PFAS. Each initial PFAS concentration resulted in much more GFP produced compared to the hlFAB-GFP. The increase in GFP production as initial PFAS concentration increases does not appear to increase logarithmically as did the hlFAB-GFP construct. However, the graphs imply the amount of GFP produced has yet to reach a steady state, indicating that longer exposure times to PFAS may produce more differentiated levels of GFP. When considering GFP production in the presence of leaking promoters, it appears that it is more difficult to distinguish concentrations below 1E-4, as the shapes of the graph are very similar. It is possible that longer incubation time will make for more differentiated graphs, although more testing is required. The differentiability of the graphs will also depend on the sensitivity of the fluorimeter used. Steady state conditions appear to slightly increase the rate of GFP at any given time, likely because there is more free STF available to bind to PFAS as soon as PFAS diffuses into the cell. In summary, a synthetic transcription factor inducible by PFAS will likely produce much more fluorescence than something similar to hlFAB-GFP, with the fluorescence more limited by time than by the amount of PFAS available. However, promoter leakage must be kept to a minimum to better separate very low concentrations of PFAS, or be incubated for longer periods of time. Our modeling demonstrates that steady state conditions do not significantly affect the fluorescence time or maximal fluorescence for hlFAB-GFP, and affects the STF by simply increasing the rate at which PFAS can be bound and GFP can be produced.