Elastic network analysis

Here we use a Gaussian network model to characterise conformational states of a trajectory.

Last executed: Dec 27, 2020 with MDAnalysis 1.0.0

Last updated: January 2020

Minimum version of MDAnalysis: 1.0.0

Packages required:

• MDAnalysis ([MADWB11], [GLB+16])

• MDAnalysisTests

Optional packages for visualisation:

• matplotlib

Note

The elastic network analysis follows the approach of ([HKP+07]). Please cite them when using the MDAnalysis.analysis.gnm module in published work.

import MDAnalysis as mda
from MDAnalysis.tests.datafiles import PSF, DCD, DCD2
from MDAnalysis.analysis import gnm
import matplotlib.pyplot as plt
%matplotlib inline

Loading files

The test files we will be working with here feature adenylate kinase (AdK), a phosophotransferase enzyme. ([BDPW09])

u1 = mda.Universe(PSF, DCD)
u2 = mda.Universe(PSF, DCD2)

Using a Gaussian network model

Using a Gaussian network model to represent a molecule as an elastic network, we can characterise the concerted motions of a protein, and the dominance of these motions, over a trajectory. The analysis is applied to the atoms in the selection. If two atoms are within the cutoff distance (default: 7 ångström), they are considered to be bound by a spring. This analysis is reasonably robust to the choice of cutoff (between 5-9 Å), but the singular value decomposition may not converge with a lower cutoff.

We use the GNMAnalysis class (API docs) for the analysis.

nma1 = gnm.GNMAnalysis(u1,
                      select='name CA',
                      cutoff=7.0)
nma1.run()

The output is saved in nma1.results: the time in picoseconds, the first eigenvalue, and the first eigenvector, associated with each frame.

len(nma1.results)

98

nma2 = gnm.GNMAnalysis(u2,
                      select='name CA',
                      cutoff=7.0)
nma2.run()

Unlike normal mode analysis, Gaussian network model analysis uses only a single eigenvalue to represent the rotation and translation of each frame. The motion with the lowest positive eigenvalue represents the dominant motion of a structure. The frequency of this motion is the square root of the eigenvalue.

Plotting the probability distribution of the frequency for the first eigenvector can highlight variation in the probability distribution, which can indicate trajectories in different states.

Below, we plot the distribution of eigenvalues. The dominant conformation state is represented by the peak at 0.06.

eigenvalues1 = [res[1] for res in nma1.results]
eigenvalues2 = [res[1] for res in nma2.results]

histfig, histax = plt.subplots(nrows=2, sharex=True, sharey=True)
histax[0].hist(eigenvalues1)
histax[1].hist(eigenvalues2)

histax[1].set_xlabel('Eigenvalue')
histax[0].set_ylabel('Frequency')
histax[1].set_ylabel('Frequency');

When we plot how the eigenvalue varies with time, we can see that the simulation transitions into the dominant conformation and stays there in both trajectories.

time1 = [res[0] for res in nma1.results]
time2 = [res[0] for res in nma2.results]
linefig, lineax = plt.subplots()
plt.plot(time1, eigenvalues1, label='DCD')
plt.plot(time2, eigenvalues2, label='DCD2')
lineax.set_xlabel('Time (ps)')
lineax.set_ylabel('Eigenvalue')
plt.legend();

DCD and DCD2 appear to be in similar conformation states.

Using a Gaussian network model with only close contacts

The MDAnalysis.analysis.gnm.closeContactGNMAnalysis class provides a version of the analysis where the Kirchhoff contact matrix is generated from close contacts between individual atoms in different residues, whereas the GNMAnalysis class generates it directly from all the atoms. In this close contacts class, you can weight the contact matrix by the number of atoms in the residues.

nma_close = gnm.closeContactGNMAnalysis(u1,
                                        select='name CA',
                                        cutoff=7.0,
                                        weights='size')
nma_close.run()

eigenvalues_close = [res[1] for res in nma_close.results]

plt.hist(eigenvalues_close)
plt.xlabel('Eigenvalue')
plt.ylabel('Frequency');

time_close = [res[0] for res in nma_close.results]
ax = plt.plot(time_close, eigenvalues_close)
plt.xlabel('Time (ps)')
plt.ylabel('Eigenvalue');

References

[1] Oliver Beckstein, Elizabeth J. Denning, Juan R. Perilla, and Thomas B. Woolf. Zipping and Unzipping of Adenylate Kinase: Atomistic Insights into the Ensemble of Open↔Closed Transitions. Journal of Molecular Biology, 394(1):160–176, November 2009. 00107. URL: https://linkinghub.elsevier.com/retrieve/pii/S0022283609011164, doi:10.1016/j.jmb.2009.09.009.

[2] Richard J. Gowers, Max Linke, Jonathan Barnoud, Tyler J. E. Reddy, Manuel N. Melo, Sean L. Seyler, Jan Domański, David L. Dotson, Sébastien Buchoux, Ian M. Kenney, and Oliver Beckstein. MDAnalysis: A Python Package for the Rapid Analysis of Molecular Dynamics Simulations. Proceedings of the 15th Python in Science Conference, pages 98–105, 2016. 00152. URL: https://conference.scipy.org/proceedings/scipy2016/oliver_beckstein.html, doi:10.25080/Majora-629e541a-00e.

[3] Benjamin A. Hall, Samantha L. Kaye, Andy Pang, Rafael Perera, and Philip C. Biggin. Characterization of Protein Conformational States by Normal-Mode Frequencies. Journal of the American Chemical Society, 129(37):11394–11401, September 2007. 00020. URL: https://doi.org/10.1021/ja071797y, doi:10.1021/ja071797y.

[4] Naveen Michaud-Agrawal, Elizabeth J. Denning, Thomas B. Woolf, and Oliver Beckstein. MDAnalysis: A toolkit for the analysis of molecular dynamics simulations. Journal of Computational Chemistry, 32(10):2319–2327, July 2011. 00778. URL: http://doi.wiley.com/10.1002/jcc.21787, doi:10.1002/jcc.21787.