Anisotropic Network Model (ANM)

This example shows how to perform ANM calculations, and retrieve normal mode data. An ANM instance that stores Hessian matrix (and also Kirchhoff matrix) and normal mode data describing the intrinsic dynamics of the protein structure will be obtained. ANM instances and individual normal modes (Mode) can be used as input to functions in dynamics module.

See [Doruker00] and [Atilgan01] for more information on the theory of ANM.

[Doruker00]Doruker P, Atilgan AR, Bahar I. Dynamics of proteins predicted by molecular dynamics simulations and analytical approaches: Application to a-amylase inhibitor. Proteins 2000 40:512-524.
[Atilgan01]Atilgan AR, Durrell SR, Jernigan RL, Demirel MC, Keskin O, Bahar I. Anisotropy of fluctuation dynamics of proteins with an elastic network model. Biophys. J. 2001 80:505-515.

Parse structure

We start by importing everything from the ProDy package:

In [1]: from prody import *

In [2]: from pylab import *

In [3]: ion()

We start with parsing a PDB file by passing an identifier. Note that if a file is not found in the current working directory, it will be downloaded.

In [4]: p38 = parsePDB('5uoj')

In [5]: p38
Out[5]: <AtomGroup: 5uoj (3138 atoms)>

We want to use only Cα atoms, so we select them:

In [6]: calphas = p38.select('protein and name CA')

In [7]: calphas
Out[7]: <Selection: 'protein and name CA' from 5uoj (343 atoms)>

We can also make the same selection like this:

In [8]: calphas2 = p38.select('calpha')

In [9]: calphas2
Out[9]: <Selection: 'calpha' from 5uoj (343 atoms)>

To check whether the selections are the same, we can try:

In [10]: calphas == calphas2
Out[10]: True

Note that, ProDy atom selector gives the flexibility to select any set of atoms to be used in ANM calculations.

Build Hessian

We instantiate an ANM instance:

In [11]: anm = ANM('p38 ANM analysis')

Then, build the Hessian matrix by passing selected atoms (351 Cα’s) to ANM.buildHessian() method:

In [12]: anm.buildHessian(calphas)

We can get a copy of the Hessian matrix using ANM.getHessian() method:

In [13]: anm.getHessian().round(3)
Out[13]: 
array([[11.348, -1.575,  0.585, ...,  0.   ,  0.   ,  0.   ],
       [-1.575, 10.709,  2.121, ...,  0.   ,  0.   ,  0.   ],
       [ 0.585,  2.121,  7.943, ...,  0.   ,  0.   ,  0.   ],
       ...,
       [ 0.   ,  0.   ,  0.   , ...,  4.221, -1.615,  2.557],
       [ 0.   ,  0.   ,  0.   , ..., -1.615, 10.354, -1.451],
       [ 0.   ,  0.   ,  0.   , ...,  2.557, -1.451,  8.425]])

Parameters

We didn’t pass any parameters to ANM.buildHessian() method, but it accepts cutoff and gamma parameters, for which default values are cutoff=15.0 and gamma=1.0.

In [14]: anm.getCutoff()
Out[14]: 15.0

In [15]: anm.getGamma()
Out[15]: 1.0

Note that it is also possible to use an externally calculated Hessian matrix. Just pass it to the ANM instance using ANM.setHessian() method.

Calculate normal modes

Calculate modes using ANM.calcModes() method:

In [16]: anm.calcModes()

Note that by default 20 non-zero (or non-trivial) and 6 trivial modes are calculated. Trivial modes are not retained. To calculate a different number of non-zero modes or to keep zero modes, try anm.calcModes(50, zeros=True).

Normal modes data

In [17]: anm.getEigvals().round(3)
Out[17]: 
array([0.156, 0.249, 0.36 , 0.863, 0.995, 1.215, 1.294, 1.47 , 1.568,
       1.808, 1.944, 2.083, 2.239, 2.326, 2.494, 2.627, 2.777, 2.821,
       2.882, 3.095])

In [18]: anm.getEigvecs().round(3)
Out[18]: 
array([[ 0.037, -0.025,  0.038, ..., -0.013, -0.006, -0.004],
       [-0.017, -0.08 ,  0.024, ..., -0.003, -0.008,  0.005],
       [ 0.055,  0.02 ,  0.053, ..., -0.083, -0.013,  0.01 ],
       ...,
       [ 0.023, -0.087, -0.017, ...,  0.152, -0.044,  0.035],
       [ 0.016, -0.01 ,  0.003, ...,  0.038, -0.035,  0.06 ],
       [ 0.056, -0.012, -0.011, ..., -0.071, -0.024,  0.068]])

You can get the covariance matrix as follows:

In [19]: anm.getCovariance().round(2)
Out[19]: 
array([[ 0.02,  0.01,  0.01, ...,  0.01,  0.  ,  0.02],
       [ 0.01,  0.04, -0.02, ...,  0.02,  0.  ,  0.  ],
       [ 0.01, -0.02,  0.05, ..., -0.01,  0.01,  0.02],
       ...,
       [ 0.01,  0.02, -0.01, ...,  0.29,  0.05, -0.07],
       [ 0.  ,  0.  ,  0.01, ...,  0.05,  0.02, -0.  ],
       [ 0.02,  0.  ,  0.02, ..., -0.07, -0.  ,  0.07]])

Covariance matrices are calculated using the available modes (slowest 20 modes in this case). If the user calculates M slowest modes, only they will be used in the calculation of covariances.

Individual modes

Normal mode indices in Python start from 0, so the slowest mode has index 0. By default, modes with zero eigenvalues are excluded. If they were retained, the slowest non-trivial mode would have index 6.

Get the slowest mode by indexing ANM instance as follows:

In [20]: slowest_mode = anm[0]

In [21]: slowest_mode.getEigval().round(3)
Out[21]: 0.156

In [22]: slowest_mode.getEigvec().round(3)
Out[22]: array([ 0.037, -0.017,  0.055, ...,  0.023,  0.016,  0.056])

Write NMD file

ANM results in NMD format can be visualized using Normal Mode Wizard VMD plugin. The following statement writes the slowest 3 ANM modes into an NMD file:

In [23]: writeNMD('p38_anm_modes.nmd', anm[:3], calphas)
Out[23]: 'p38_anm_modes.nmd'

Note that slicing an ANM objects returns a list of modes. In this case, slowest 3 ANM modes were written into NMD file.

View modes in VMD

First make sure that the VMD path is correct

In [24]: pathVMD()
Out[24]: '/home/exx/miniconda3/envs/py27/bin/vmd'

# if this is incorrect use setVMDpath to correct it
In [25]: viewNMDinVMD('p38_anm_modes.nmd')

This will show the slowest 3 modes in VMD using NMWiz. This concludes the ANM example. Many of the methods demonstrated here apply to other NMA models, such as GNM and EDA.

Advanced visualization in jupyter notebooks

You can visualize structures and modes determined from ANM or GNM calculations in jupyter notebooks using another python module, py3Dmol. It is a java-script library that can visualize structural elements with light weight customization.

You can find an example notebook.