7.1.3. Distances¶

Distances can be a very useful observable to track conformational changes. They can often be directly related to real experimental observables such as NOEs from NMR experiments or distances from cross-linking or FRET experiments.

Here we calculate a simple distance between two C_α atoms as an approximation to the distance between two chromophores attached to the correspoding residues isoleucine 52 (I52) and lysine 145 (K145) used in a FRET experiment [Henzler-Wildman2007].

First we need to create an index file containing the two groups:

mkdir -p analysis/dist/I52_K145 && cd analysis/dist/I52_K145
gmx make_ndx -f ../../../MD/md.tpr -o I52_K145.ndx

Use interactive commands like the following [1]

keep 0
del 0
r 52 & a CA
name 0 I52
r 145 & a CA
name 1 K145
q

to generate the index file I52_K145.ndx.

Then run gmx distance and compute the distance between the two atoms:

printf "I52\nK145\n" | gmx distance -s ../../../MD/md.tpr -f ../../../MD/md.xtc -n I52_K145.ndx -o dist.xvg

The dist.xvg file contains the distance in nm for each time step in ps, which can be plotted [2].

Timeseries of the distance between |Calpha| atoms of I52 and K145.

Timeseries of the distance between the C_α atoms of I52 (NMP domain) and K145 (LID domain).

(You can also use the centered and fitted trajectory md_fit.xtc as an input instead of md.xtc to make sure that the distance calculation does not contain any jumps due to periodic boundary effects, or use gmx mindist.)

See also

[Beckstein2009] for a discussion of FRET distances in AdK.

Footnotes

[1]

Note that one has to be careful when selecting residue ids in make_ndx. It is often the case that a PDB file does not contain all residues, e.g. residues 1–8 might be unresolved in the experiment and thus are missing from the PDB file. The file then simply starts with residue number 9. Gromacs, however, typically renumbers residues so that they start at 1. Thus, in this hypothetical case, a residue that might be referred to in the literature as “residue 100” might actually be residue 92 in the simulation (\(N^\mathrm{sim}_\mathrm{res} = N^\mathrm{PDB}_\mathrm{res} - (\mathrm{min} N^\mathrm{PDB}_\mathrm{res} - 1)\)). Thus, if you wanted to select the C_α atom of residue 100 you would need to select r 92 & a CA in make_ndx.

[2]

To plot in Python, make sure to not write xvg legend information to the output file (using the -xvg none flag)

printf "I52\nK145\n" | \
     gmx distance -s ../../../MD/md.tpr -f ../../../MD/md.xtc \
                  -n I52_K145.ndx -o dist.xvg -xvg none

so that you can easily read the data with numpy.loadtxt():

import matplotlib.pyplot as plt
import numpy

t,d,x,y,z = numpy.loadtxt("dist.xvg", unpack=True)

fig = plt.figure(figsize=(5,2.5))
ax = fig.add_subplot(111)
fig.subplots_adjust(bottom=0.2)

ax.fill_between(t, d, color="orange", linestyle="-", alpha=0.1)
ax.plot(t, d, color="orange", linestyle="-", label="I52-K145")

ax.set_xlabel("time $t$ (ps)")
ax.set_ylabel(r"C$_\alpha$ distance (nm)")
ax.legend(loc="best")

fig.savefig("d_I52_K145_ca.png", dpi=300)
fig.savefig("d_I52_K145_ca.svg")
fig.savefig("d_I52_K145_ca.pdf")