# Elastic network analysis¶

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

Last executed: Feb 06, 2020 with MDAnalysis 0.20.2-dev0

Last updated: January 2020

Minimum version of MDAnalysis: 0.17.0

Packages required:

Optional packages for visualisation:

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


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.

:

nma1 = gnm.GNMAnalysis(u1,
selection='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,
selection='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 for res in nma1.results]
eigenvalues2 = [res for res in nma2.results]

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

histax.set_xlabel('Eigenvalue')
histax.set_ylabel('Frequency')
histax.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 for res in nma1.results]
time2 = [res 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,
selection='name CA',
cutoff=7.0,
weights='size')
nma_close.run()

:

eigenvalues_close = [res for res in nma_close.results]

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

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