Team:Queens Canada/Static-Dynamic

In Silico Protein Stability


An immediate concern of ours was whether our biosensor constructs would be structurally stable in a realistic environment. Given the relatively small size of our sensing-proteins, the addition of two relatively large fluorophores to achieve FRET was thought to lead to thermodynamic, and steric instability that would compromise the stability of the sensing-protein, resulting in impaired biosensor function. To understand the potential fluorophore-induced stress placed on our constructs: Fluorescent biosensor protein for phosphate (FBP-Pi), Fluorescent biosensor protein for potassium (FBP-K), Fluorescent biosensor protein for glucose (FBP-Glc), Fluorescent biosensor protein for PTH (FBP-PTH), and Fluorescent biosensor protein for FGF23 (FBP-FGF23), we performed nanoscale molecular dynamic simulations. These in silico simulations provided us with valuable insights into the nanoscale interactions of our biosensor constructs in different solvent systems - allowing us to assess each biosensor’s stability.


How did we perform nanoscale molecular dynamic simulations?

All the simulations were performed on the software GROMACS 2020.2 with the NVIDA CUDA toolkit. The methodology employed to generate our models was adapted from Justin Lemkul’s GROMACS Tutorials.

Step 1. Create Fluorescent Biosensor Protein Constructs in PyMOL and Chimera

We used Chimera and PyMOL, both molecular modelling software’s, to create our fluorescent biosensor protein constructs before performing any advanced dynamic modelling. This approach of static modelling first gave us a chance to look for surface residues to perform mutagenesis and assess whether our protein is potentially compatible with FRET – dependent on the distance between fluorophores on the binding proteins tails.

Phosphate binding protein fluorescent construct. Protein structures were obtained from the RCSB Protein Data Bank and a construct was made using Chimera software. Torsion angles between fluorescent and binding proteins are adjusted for display purposes.. mNeonGreen (green), binding protein (blue), and mCherry (red) are all displayed using PyMOL. A. Phosphate binding protein with a cartoon preset coupled with two fluorophores (mNeonGreen and mCherry). B. Construct in a stick preset to provide insights regarding the construct volume in 3D space.

Step 2. Generate the Protein Topology.

After importing the FBP construct files into GROMACS, the protein topology was generated. The topology contains the defining characteristics of a protein within a simulation, including nonbonded (atom types, charges) and bonded (bonds, angles, dihedrals) parameters. It also contains a force field, which is a collection of equations and associated constants that attempt to recapitulate the physical characteristics of a protein. The PLS-AA/L all-atom force field was used for our models.

Step 3. Create and Solvate a Box (Environment Generation).

For our modelling purposes a simple aqueous system was generated. This involved defining a box that would contain our FBP constructs, and then systematically filling the box with water molecules.

Step 4. Neutralizing System Net Charge.

The solvated system we had created contains our charged FBP construct. This net charge can lead to several artifacts and inaccuracies in the simulation; therefore, it was neutralized by adding ions. We systematically added Na+ and Cl- ions to the solvation to neutralize any charge that was present in the protein.

Step 5. Energy Minimization.

We have now created a solvated, and electroneutral system. To ensure there were no steric clashes or misfolding in our constructs we performed structural relaxation through energy minimization.

Step 6. Isothermal-Isochoric Equilibration

To avoid system collapse and begin performing the molecular dynamic simulations, the solvent and ions were equilibrated around the protein. This was done by equilibrating the system on the basis of temperature to ensure proper orientation of the protein.

Step 7. Isothermal-Isobaric Equilibration

This step is almost identical to Step 6 except the system was equilibrated and stabilized on the basis of pressure and density instead of temperature.

Step 8. Run Nanoscale Molecular Dynamic Simulation.

This is the final step. After our system has been energy minimized and equilibrated, we can start collecting data on the nanoscale dynamics of each atom in our system. This provides us information on the interactions of our protein constructs in the simulation system, allowing the interpretation of thermodynamic and structural stability.

Results & Interpretation

How will generating these models help us better understand our models?

Several nanoscale molecular dynamic simulations were performed on all our constructs: FBP-Pi, FBP-K, FBP-Glc, FBP-PTH, and FBP-FGF23. We simulated all our constructs on the software GROMACS 2020.02 and created nanoscale models of these proteins in a simple water and 0.1M NaCl solution. By doing so, we were able to visualize the trajectory of these proteins over one nanosecond.

Figure 1.1. The nanoscale trajectory of fluorescent biosensor protein for phosphate (FBP-Pi).

Figure 1.2. The nanoscale trajectory of fluorescent biosensor protein for potassium (FBP-K).

Figure 1.3. The nanoscale trajectory of fluorescent biosensor protein for glucose (FBP-Glc).

Figure 1.4. The nanoscale trajectory of fluorescent biosensor protein for parathyroid hormone (FBP-PTH).

Figure 1.5. The nanoscale trajectory of fluorescent biosensor protein for fibroblast growth factor 23 (FBP-FGF23).

Assessment of the atomic trajectory of our proteins over one nanosecond allowed us to understand and visualize the stability of our biosensor constructs. Despite this simulation being fairly rudimentary in the landscape of possibilities provided by molecular dynamics, it provided our team with valuable insight towards structural properties our proteins, including (1) the process of energy minimization, (2) root-mean-square-deviation, and (3) the radius of gyration. This data is represented through Figures 2-4, respectively.

Figure 2. Structural relaxation of biosensor constructs through energy minimization.

To ensure our fluorescent biosensor protein (FBP) constructs were devoid of steric clashes and inappropriate geometries, we performed structural relaxation through energy minimization via the steepest decent algorithm (1). Although this does not provide specific insights towards the stability of the FBP construct, it ensures the protein system is stable and ready for MD simulations.

Figure 3. Assessment of protein structure compactness through radius of gyration measurement.

The radius of gyration of a protein serves as an indicator of protein structure compactness – and by extension, its stability. It has an essential role in structural comparison of the constructs as it provides insight into how the proteins would fluctuate if expressed. The radius of gyration is concerned how regular secondary structures are compactly packed into the 3D structure of a protein, with the protein domains of four major classes; α, β, α/β, and α + β, having characteristic radius of gyration. This determines the protein structure compactness, allowing us to infer how stable our proteins are folded. The results presented in Figure 3. suggest that all of our constructs, potentially with the exception of FBP-FGF23, are stable.

Figure 4. Assessing conformational variation through root-mean-square deviation (RMSD).

Root-mean-square-deviation (RMSD) is indicative of the variation in a protein’s conformation over time via measuring the distance between coordinates (i.e. the average distance between backbone atoms of a protein). We conducted this analysis to assess if there were any major, unfavorable conformational changes in our constructs that may confer a loss of function. Given by the low RMSD values (<0.5 nm) we can infer that our constructs are in a relatively stable conformation.

Future Directions

While our nanoscale molecular dynamics simulations were informing, they only provided us a glance into how our constructs actually behave given the duration of simulation was only one nanosecond. Therefore, the main direction we will take when looking to model our constructs in the future will be to increase the simulation time. To do this, we will need access to large CPU clusters to provide the computational power required to simulate these systems for 10, or 100 nanoseconds.

Further, due to the lack of computational power we were unable to simulate scenarios such as ligand-binding, or protein-protein docking with the coiled-coil interacting with the cysteine residues of the main construct. These simulations would provide insight towards how the binding of a ligand would affect the protein conformation, along with how interaction of the coiled-coil with the construct would affect thermodynamic and stability properties of our proteins.


  1. Lobanov, M. Yu., Bogatyreva, N. S., and Galzitskaya, O. V. (2008) Radius of gyration as an indicator of protein structure compactness. Mol Biol. 42, 623–628.
  2. Sargsyan, K., Grauffel, C., and Lim, C. (2017) How Molecular Size Impacts RMSD Applications in Molecular Dynamics Simulations. J. Chem. Theory Comput. 13, 1518–1524.
  3. J.A. Lemkul (2018) "From Proteins to Perturbed Hamiltonians: A Suite of Tutorials for the GROMACS-2018 Molecular Simulation Package, v1.0" Living J. Comp. Mol. Sci. In Press. GitHub.
  4. The PyMOL Molecular Graphics System, Version 2.0 Schrödinger, LLC, 2019,
  5. Humphrey, W., Dalke, A. and Schulten, K., "VMD - Visual Molecular Dynamics", J. Molec. Graphics, 1996, vol. 14, pp. 33-38.