How to compute the RMSF using GROMACS
The gmx rmsf module
GROMACS is one of the most widely used software for conducting molecular simulations. On this website, you will find different GROMACS tutorials that will help you to run protein simulations.
However, once your simulation is complete, you should also be interested to continue with the analysis of the resulting trajectory. Among different analyses at your disposal, one that often holds great importance is the Root Mean Square Fluctuation, fondly known as RMSF. In this article, you’ll learn about the concept of RMSF and gain practical insights into computing it using GROMACS’ dedicated gmx rmsf
module.
What is the RMSF?
The Root Mean Square Fluctuation (RMSF) is an important structural feature that one can study in Molecular Dynamics (MD). It measures the fluctuation of an atom or a group of atoms (e.g., a protein residue) along the course of a simulation, and is useful to identify the most mobile regions of our system providing valuable insights on the structural dynamics of the system.
The RMSF ($\rho$) is computed as the square root of the variance of the fluctuation around the average position. In other words, we compute the root mean square of the deviations of the atomic positions from their mean positions over a series of time steps of a molecular dynamics simulation.
For a generic atom $i$ we have that:
$$ \rho_i = \sqrt {\big \langle ( r_i - \langle r_i \rangle)^2 \big \rangle } $$
where $r_i$ represent the coordinates of atom $i$, and $\langle r_i \rangle$ its ensemble average position.
The interpretation of the RMSF value is quite straightforward:
-
A region characterized by high RMSF values indicates high divergence from the average position and therefore high structural mobility.
-
The opposite is true for regions displaying low RMSF values which will not deviate much from the average position, meaning that they are more rigid during the simulation.
RMSD vs RMSF
You may notice that the basic concept of this feature is somewhat similar to the concept of Root Mean Square Deviation (RMSD). However, the information that we retrieve from the two analyses is different.
- The RMSD shows how much the molecule of interest deviates from its initial structure over the course of the simulation. Therefore, it is used as an indicator of convergence and is generally plotted as a function of time. When the RMSD curve becomes flat the system is not changing significantly anymore and it can be considered equilibrated.
- The RMSF determines the fluctuation of an atom (or a group of atoms) during the simulation. This means that it is generally used to check the flexibility of residues within a protein during the course of the simulation. For this reason, it needs to be plotted against the residue number to measure the contributions of different amino acids to the overall molecular motion.
Typically, one would perform an RMSD analysis to verify that the value reached a plateau and that the structure has equilibrated.
Subsequently, you may want to perform the analysis of the RMSF during the equilibrated portion of the trajectory to observe which residues fluctuate the most even when the structure is stable.
The gmx rmsf command
You can use this tool to compute the RMSF of different residues during an MD simulation.
|
|
- Call the program with
gmx
- Select the
rmsf
module - Select the
-f
flag and provide the trajectory file from your simulation (simulation.xtc
) - Choose the reference structure (
simulation.tpr
) with the-s
flag. This is required to communicate GROMACS the masses and connectivity of each atom. - Call the
-o
flag and decide how you want to name the output xvg file (rmsf.xvg
) - The
-res
flag is used to compute the RMSF for each residue.
You can use the -b
and -e
flags to compute the RMSF for a specific time frame.
- Pick the starting time (1000 ps) for the RMSD calculation using the
-b
flag - Select the ending time (10000 ps) using the
-e
option
In protein simulation, we generally compute the RMSF of the $C_{\alpha}$ which gives us the fluctuation of the backbone for each residue. At the prompt, select the group C-alpha
.
The RMSF will be computed from the average structure in the analyzed time frame and you will find the output in the rmsf.xvg
file.
As usual, you can plot the resulting xvg file with your favorite tool. Here I showed you how you can plot it with Grace or Python.
Grace
|
|
![RMSF as a function of time in GROMACS](/gromacs_rmsf/rmsf.webp)
Python
|
|
![](/gromacs_rmsf/rmsf_p.webp)
The spikes in RMSF correspond to the residues having wider fluctuation along the trajectory.
B-factor analysis with GROMACS
You can also use the gmx rmsf
command to perform the B-factor analysis.
In this analysis, the RMSF values are converted into B-factors, offering a valuable measure of atom fluctuation within your system. These B-factor values are then efficiently writtein to a new PDB file, where they are conveniently placed in the last column (for an in-depth look at PDB formatting).
The usefulness of this approach becomes apparent when you utilize tools like PyMOL to import the PDB file and apply color schemes based on B-factor values to your molecular structure. This dynamic visualization offers a direct means of evaluating atom fluctuations within your molecular system.
To perform this analysis you will need to add the -oq
flag to the previous command and decide how to name the output file (bfac.pdb
).
|
|
The outcome is a PDB file enriched with B-factor data, ready for import into visualization tools like PyMOL or VMD. This file allows you to visually identify the most dynamic regions within your system, offering valuable insights into molecular mobility.