Recentering and Rewrapping Coordinates

As in any simulation conducted with periodic boundary conditions, molecules may appear "broken" or may "jump" back and forth across the box. To recenter the protein and rewrap the molecules within the unit cell to recover the desired rhombic dodecahedral shape, invoke trjconv:

gmx trjconv -s md_0_10.tpr -f md_0_10.xtc -o md_0_10_center.xtc -center -pbc mol -ur compact

Choose "Protein" for centering and "System" for output. Note that centering complexes (protein-ligand, protein-protein) may be difficult for longer simulations involving many jumps across periodic boundaries. In those instances (particularly in protein-protein complexes), it may be necessary to create a custom index group to use for centering, corresponding to the active site of one protein or the interfacial residues of one monomer in a complex.

To extract the first frame (t = 0 ns) of the trajectory, use trjconv -dump with the recentered trajectory:

gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o start.pdb -dump 0

For even smoother visualization, it may be beneficial to perform rotational and translational fitting. Execute trjconv as follows:

gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o md_0_10_fit.xtc -fit rot+trans

Choose "Backbone" to perform least-squares fitting to the protein backbone, and "System" for output. Note that simultaneous PBC rewrapping and fitting of the coordinates is mathematically incompatible. Should you wish to perform fitting (which is useful for visualization, but not necessary for most analysis routines), carry out these coordinate manipulations separately, as indicated here.

Analyzing Protein-Ligand Interactions and Ligand Dynamics

This tutorial cannot possible cover all analysis methods that you may wish to perform. A few basic operations will be illustrated here.

The 2-propylphenol ligand can engage in a hydrogen bond with the Gln102 side chain. The GROMACS hbond module can easily be employed to calculate the number of hydrogen bonds between any groups of atoms, but in our case, the only values will be 1 or 0. For a more detailed look at how the ligand is interacting with Gln102, we will compute the distance between the hydroxyl group of JZ4 and the carbonyl O atom of Gln102. For a hydrogen bond to be formed, the typical criterion is that the donor and acceptor atoms will be separated by a distance ≤ 3.5 Å (0.35 nm). Use the distance module to calculate the distance over the course of the trajectory, using command-line selection syntax (see gmx help selections for examples and more syntax).

gmx distance -s md_0_10.tpr -f md_0_10_center.xtc -select 'resname "JZ4" and name OAB plus resid 102 and name OE1' -oall

The average distance is 0.31 ± 0.05 nm, consistent with the formation of a hydrogen bond.

The second criterion usually applied in determining the presence of a hydrogen bond is the angle formed among the donor, hydrogen, and acceptor atoms. There are different conventions for calculating the angle. In the GROMACS hbond module, the angle is defined as hydrogen-donor-acceptor, and this angle should be ≤ 30°. To perform this analysis, first create index groups for the donor atoms (which must include both the donor heavy atom and the bonded hydrogen) and the acceptor atom:

gmx make_ndx -f em.gro -o index.ndx
...
 > 13 & a OAB | a H12
 (creates group 23)
 > 1 & r 102 & a OE1
 (creates group 24)
 > 23 | 24
 > q

Analyze the angle formed by these three atoms using the angle module:

gmx angle -f md_0_10_center.xtc -n index.ndx -ov angle.xvg

Note that angle takes no -s argument, unlike most other GROMACS analysis modules. Angle calculations do not require mass or periodicity information, atom names, etc. so the trajectory is processed without a .tpr or coordinate file. The value returned by the command should be roughly 23 ± 9°. This outcome may be somewhat unexpected, as the index group we constructed was in the order of OAB, H12, OE1, which would correspond to the donor-hydrogen-acceptor distance, which we expect to be close to linear (~180°). Inspect the contents of the index file and you will find:

[ JZ4_&_OAB_H12_Protein_&_r_102_&_OE1 ]
1616 2624 2636

make_ndx has sorted the atom numbers automatically from low to high, thus the outcome of the calculation is the acceptor-donor-hydrogen angle, the same angle that the hbond module would have calculated. So the result is consistent with formation of a hydrogen bond, since it is ≤ 30°. To get the desired angle of donor-hydrogen-acceptor, we would have to manually edit the index group in a text file to reorder the atom numbers (2624 2636 1616). Re-running the angle calculation with this index group yields an average value of 147 ± 11°.

Finally, we may be interested in quantifying how much the ligand binding pose has changed over the course of the simulation. To compute a heavy-atom RMSD of just JZ4, create a new index group for it:

gmx make_ndx -f em.gro -n index.ndx
...
 > 13 & ! a H*
 > name 26 JZ4_Heavy
 > q

Execute the rms module, choosing "Backbone" for least-squares fitting and "JZ4_Heavy" for the RMSD calculation. By doing so, the overall rotation and translation of the protein is removed via fitting and the RMSD reported is how much the JZ4 position has varied relative to the protein, which is a good indicator of how well the binding pose was preserved during the simulation.

gmx rms -s em.tpr -f md_0_10_center.xtc -n index.ndx -tu ns -o rmsd_jz4.xvg

The calculated RMSD should be about 0.1 nm (1 Å), indicating only a very small change in the ligand's position.

Protein-Ligand Interaction Energy

To quantify the strength of the interaction between JZ4 and T4 lysozyme, it may be useful to compute the nonbonded interaction energy between these two species. GROMACS has the ability to decompose the short-range nonbonded energies between any number of defined groups. It is important to note that this quantity is NOT a free energy or a binding energy. In fact, most force fields are not parametrized in such a way that this quantity is actually physically meaningful. CHARMM is parametrized to specifically target quantum mechanical interaction energies with water, so it is intrinsically balanced against meaningful quantities, and as such the interaction energy can be useful.

Calculation of an interaction energy is carried out via the energygrps keyword in the .mdp file. Despite being an .mdp keyword, interaction energy calculations should not be considered part of a normal simulation. Decomposition of the short-range energies is incompatible with running on a GPU and also slows the calculation down unnecessarily. The mdrun module does not need to do this extra work to perform a valid simulation. As such, only compute interaction energies as a part of your analysis, not your dynamics. Create a new .tpr file from an .mdp file that has energygrps = Protein JZ4 defined, like this one:

gmx grompp -f ie.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o ie.tpr

Next, invoke mdrun with the -rerun option to recalculate energies from the existing simulation trajectory:

gmx mdrun -deffnm ie -rerun md_0_10.xtc -nb cpu

Note the use of -deffnm to read ie.tpr and write all output files to ie.* as their file names. The -rerun option takes the name of the trajectory for which you want to recompute energies, and -nb cpu tells mdrun to only attempt to run on CPU hardware and ignore any GPU that might be available. As stated above, this type of calculation cannot be performed on a GPU. The rerun should be very fast, completing in just a few minutes.

Extract the energy terms of interest via the energy module. The terms we are interested in are Coul-SR:Protein-JZ4 and LJ-SR:Protein-JZ4.

gmx energy -f ie.edr -o interaction_energy.xvg

The average short-range Coulombic interaction energy is -20.5 ± 7.4 kJ mol^-1 and the short-range Lennard-Jones energy is -99.1 ± 7.2 kJ mol^-1. It may be tempting to draw conclusions from the relative magnitudes of these quantities, but even though CHARMM was parametrized against explicit water interaction energies, decomposition of interaction energies further into these components is not necessarily real. There is no way to experimentally verify these quantities, so it is impossible to know whether they are meaningful. The total interaction energy, however, is useful in this case. That value (after propagating the error according to the standard formula for addition of two quantities) is -119.6 ± 10.3 kJ mol^-1.