Analysing a molecular dynamics trajectory: solvated chignolin

This tutorial runs a short molecular dynamics simulation of chignolin (PDB 1UAO, a 10-residue β-hairpin) in a box of water at constant pressure, and analyzes the trajectory with the analysis tools of EasyHybrid: box density, alignment, RMSF, Ramachandran plot, secondary structure (DSSP), principal components and the radial distribution function of water. The tools are described in the User Guide pages of the same names.

1. Load the solvated system

File > Open examples/pkl/1UAO_OPLS_folded_solvated.pkl (also used in the Molecular dynamics basics tutorial): chignolin with OPLS-AA, 665 TIP3P waters and two Na+ ions in a periodic cubic box of 27.5 Å, 2135 atoms.

2. Run a short NPT simulation

Open Simulate > Molecular Dynamics and choose the Leap Frog Dynamics integrator, which allows pressure coupling. Set 10 000 steps of 0.001 ps (10 ps), a temperature of 300 K, tick Pressure Coupling with a pressure of 1 atm, and save the trajectory every 100 steps (101 frames). Click Run; the run took about 10 minutes on a desktop computer. When it finishes, import the trajectory (the Import... button of the Finished! dialog).

A 10 ps run is far too short for real conclusions; it is enough to learn the tools. Use nanoseconds for a study.

3. Box density: is the system equilibrated?

Open Analysis > Box Density Analysis, choose the trajectory and click Run (Figure 1). The density starts at about 1.0 g/cm3, rises during the first frames while the box adjusts to the pressure, and then fluctuates around 1.04 g/cm3 (standard deviation 0.012 g/cm3). The log file lists the volume and density of every frame and their statistics.

4. Align the trajectory and compute the RMSF

The fluctuations of the protein must be measured after removing its overall rotation and translation. Open Analysis > Align Trajectory, choose the trajectory, align By the backbone (N, CA, C, O) on frame 0, and click Align.

Then open Analysis > RMSF Analysis, choose the aligned trajectory, tick Backbone atoms only, group by Residue and click Run (Figure 2). The backbone RMSF stays between 0.24 and 0.32 Å along the hairpin and rises at the C-terminal Gly10 (0.41 Å): the termini are the most flexible part, as expected.

5. Ramachandran plot

Open Analysis > Ramachandran Plot, choose the trajectory, tick Use all trajectory frames and click Run (Figure 3). The 808 points (8 residues × 101 frames) fall in the β region (the two strands), in the αR region (the turn) and, for glycine, in the left-handed region. Hover over a point to identify the residue and frame.

6. Secondary structure and principal components

Open Analysis > Trajectory Analysis..., choose the trajectory, keep Ignore water and ions, select Secondary structure (DSSP) and click Run (Figure 4). Tyr2–Asp3 and Thr8–Trp9 are strands and Pro4–Gly7 a turn in every frame: the β-hairpin is stable during the run.

Select Principal components (PCA) with CA atoms and click Run (Figure 5). The first component describes 31.7% of the fluctuations of the CA atoms and the first two 50.4%. Animate PC1 creates an object that moves along the first component, a simple way to see the main collective motion of the hairpin.

7. The structure of the water: O–O radial distribution function

Use the trajectory as imported, not the aligned one: the alignment rotates the whole box, and the periodic images no longer match. Select the water oxygens and add them to the selection list (see the Radial distribution functions tutorial), then open Analysis > RDF Analysis (g(r)) with that group as Group A and Group B and click Run (Figure 6). The first peak, at about 2.8 Å, is the first hydration shell of liquid water; beyond about 8 Å g(r) approaches 1. The decrease after 14 Å appears because r max (15 Å) is larger than half the box length: keep r max below half of the box.

Going further

• Run longer simulations and compare the RMSF of the two halves, or use Clustering and Hydrogen bonds in Trajectory Analysis.

• Compute the carboxylate–water g(r) of Asp3 and Glu5 (see the Radial distribution functions tutorial).