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].
(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 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 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")
|