Gaussian Network Model (GNM)

This example shows how to perform GNM calculations using an X-ray structure of ubiquitin. A GNM instance that stores the Kirchhoff matrix and normal mode data describing the intrinsic dynamics of the protein structure will be obtained. GNM instances and individual normal modes (Mode) can be used as input to functions in prody.dynamics module.

See [Bahar97] and [Haliloglu97] for more information on the theory of GNM.

[Bahar97]Bahar I, Atilgan AR, Erman B. Direct evaluation of thermal fluctuations in protein using a single parameter harmonic potential. Folding & Design 1997 2:173-181.
[Haliloglu97]Haliloglu T, Bahar I, Erman B. Gaussian dynamics of folded proteins. Phys. Rev. Lett. 1997 79:3090-3093.

Parse structure

We start by importing everything from the ProDy package:

In [1]: from prody import *

In [2]: from matplotlib.pylab import *

In [3]: ion() # turn interactive mode on

First we parse a PDB file by passing its identifier to parsePDB() function. Note that if file is not found in the current working directory, it will be downloaded.

In [4]: ubi = parsePDB('1aar')

In [5]: ubi
Out[5]: <AtomGroup: 1aar (1218 atoms)>

This file contains 2 chains, and a flexible C-terminal (residues 71-76). We only want to use Cα atoms of first 70 residues from chain A, so we select them:

In [6]: calphas = ubi.select('calpha and chain A and resnum < 71')

In [7]: calphas
Out[7]: <Selection: 'calpha and chai...and resnum < 71' from 1aar (70 atoms)>

See definition of “calpha”, “chain”, and other selection keywords in Atom Selections.

Note that, flexible design of classes allows users to select atoms other than alpha carbons to be used in GNM calculations.

Build Kirchhoff matrix

Instantiate a GNM instance:

In [8]: gnm = GNM('Ubiquitin')

We can build Kirchhoff matrix using selected atoms and GNM.buildKirchhoff() method:

In [9]: gnm.buildKirchhoff(calphas)

We can get a copy of the Kirchhoff matrix using GNM.getKirchhoff() method:

In [10]: gnm.getKirchhoff()
Out[10]: 
array([[ 11.,  -1.,  -1., ...,   0.,   0.,   0.],
       [ -1.,  15.,  -1., ...,   0.,   0.,   0.],
       [ -1.,  -1.,  20., ...,   0.,   0.,   0.],
       ..., 
       [  0.,   0.,   0., ...,  20.,  -1.,  -1.],
       [  0.,   0.,   0., ...,  -1.,  21.,  -1.],
       [  0.,   0.,   0., ...,  -1.,  -1.,  12.]])

Parameters

We didn’t pass any parameters, but GNM.buildKirchhoff() method accepts two of them, which by default are cutoff=10.0 and gamma=1.0, i.e. buildKirchhoff(calphas, cutoff=10., gamma=1.)

In [11]: gnm.getCutoff()
Out[11]: 10.0

In [12]: gnm.getGamma()
Out[12]: 1.0

Note that it is also possible to use an externally calculated Kirchhoff matrix. Just pass it to the GNM instance using GNM.setKirchhoff() method.

Calculate normal modes

We now calculate normal modes from the Kirchhoff matrix.

In [13]: gnm.calcModes()

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

Normal mode data

Get eigenvalues and eigenvectors:

In [14]: gnm.getEigvals().round(3)
Out[14]: 
array([  2.502,   2.812,   4.366,   5.05 ,   7.184,   7.65 ,   7.877,
         9.08 ,   9.713,  10.132,  10.502,  10.644,  10.888,  11.157,
        11.285,  11.632,  11.78 ,  11.936,  12.006,  12.218])

In [15]: gnm.getEigvecs().round(3)
Out[15]: 
array([[-0.064, -0.131, -0.245, ..., -0.256,  0.538, -0.   ],
       [-0.073, -0.085, -0.19 , ...,  0.006, -0.069,  0.032],
       [-0.076, -0.043, -0.135, ...,  0.017, -0.047,  0.018],
       ..., 
       [-0.092,  0.064,  0.105, ...,  0.032, -0.042,  0.006],
       [-0.07 ,  0.099,  0.054, ...,  0.031,  0.024, -0.014],
       [-0.081,  0.135,  0.124, ...,  0.013, -0.04 , -0.018]])

Get covariance matrix:

In [16]: gnm.getCovariance().round(2)
Out[16]: 
array([[ 0.08,  0.02,  0.01, ..., -0.01, -0.01, -0.01],
       [ 0.02,  0.02,  0.01, ..., -0.  , -0.  , -0.01],
       [ 0.01,  0.01,  0.01, ...,  0.  , -0.  , -0.  ],
       ..., 
       [-0.01, -0.  ,  0.  , ...,  0.01,  0.01,  0.01],
       [-0.01, -0.  , -0.  , ...,  0.01,  0.01,  0.02],
       [-0.01, -0.01, -0.  , ...,  0.01,  0.02,  0.05]])

Note that covariance matrices are calculated using the available modes in the model, which is the slowest 20 modes in this case. If the user calculates M modes, these M modes will be used in calculating the covariance matrix.

Individual modes

Normal mode indices start from 0, so slowest mode has index 0.

In [17]: slowest_mode = gnm[0]

In [18]: slowest_mode.getEigval().round(3)
Out[18]: 2.5019999999999998

In [19]: slowest_mode.getEigvec().round(3)
Out[19]: 
array([-0.064, -0.073, -0.076, -0.112, -0.092, -0.143, -0.164, -0.205,
       -0.24 , -0.313, -0.192, -0.152, -0.066, -0.07 , -0.025, -0.031,
        0.001, -0.006, -0.015,  0.027,  0.042,  0.055,  0.063,  0.09 ,
        0.09 ,  0.069,  0.132,  0.175,  0.145,  0.121,  0.195,  0.218,
        0.158,  0.217,  0.245,  0.214,  0.225,  0.171,  0.2  ,  0.151,
        0.102,  0.043, -0.029, -0.064, -0.072, -0.086, -0.09 , -0.078,
       -0.057, -0.011,  0.016,  0.061,  0.058,  0.043,  0.029,  0.013,
        0.004,  0.011, -0.013, -0.037, -0.05 , -0.059, -0.07 , -0.094,
       -0.094, -0.099, -0.097, -0.092, -0.07 , -0.081])

By default, modes with 0 eigenvalue are excluded. If they were retained, slowest non-trivial mode would have index 6.

Plot results

ProDy plotting functions are prefixed with show. Let’s use some of them to plot data:

Contact Map

In [20]: showContactMap(gnm);
../../_images/enm_analysis_gnm_contact_map.png

Cross-correlations

In [21]: showCrossCorr(gnm);
../../_images/enm_analysis_gnm_cross_corr.png

Slow mode shape

In [22]: showMode(gnm[0]);

In [23]: plt.grid();
../../_images/enm_analysis_gnm_mode.png

Square fluctuations

In [24]: showSqFlucts(gnm[0]);
../../_images/enm_analysis_gnm_sqflucts.png