Mechanical Stiffness Calculations

This example shows how to perform mechanical resistance calculation for GFP protein (1gfl) and visualize the results using Matplotlib library and VMD program.

An ANM instance that stores Hessian matrix and normal mode data describing the intrinsic dynamics of the protein structure will be used as an input (model) as well as cooridinates of protein structure (coords, pdb).

See [EB08] for more information about the theory of mechanical resistance calculations and more examples.

[EB08](1, 2) Eyal E., Bahar I. Toward a Molecular Understanding of the Anisotropic Response of Proteins to External Forces: Insights from Elastic Network Models. Biophys. J. 2008 94:3424-34355.

Parse structure

We start by importing everything from the ProDy package:

In [1]: from prody import *

In [2]: from pylab import *

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

We start with parsing a PDB file by passing an identifier (1gfl). Note that if a file is not found in the current working directory, it will be downloaded directly from Protein Data Bank database. We will include a header of PDB file which will be used in the next function.

In [4]: gfp, header = parsePDB('1gfl', header=True)

In [5]: gfp
Out[5]: <AtomGroup: 1gfl (3950 atoms)>

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

In [6]: calphas ='protein and chain A and name CA')

In [7]: calphas
Out[7]: <Selection: 'protein and chain A and name CA' from 1gfl (230 atoms)>

Note that, ProDy atom selector gives the flexibility to select any set of atoms to be used in calculations. To obtain more information about selection see: AtomGroup.

Build Hessian

In the next step we instantiate an ANM instance:

In [8]: anm = ANM('GFP ANM analysis')

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

In [9]: anm.buildHessian(calphas, cutoff=13.0)

Those two actions are required to perform mechanical stiffness calculations. Optained ANM model will be used to mechanical striffness calculations. We can get a copy of the Hessian matrix using ANM.getHessian() method:

In [10]: anm.getHessian().round(4)
array([[ 9.0167,  2.7414, -0.9794, ...,  0.    ,  0.    ,  0.    ],
       [ 2.7414,  5.3204,  1.181 , ...,  0.    ,  0.    ,  0.    ],
       [-0.9794,  1.181 ,  5.6629, ...,  0.    ,  0.    ,  0.    ],
       [ 0.    ,  0.    ,  0.    , ...,  4.7941, -0.0517, -3.5568],
       [ 0.    ,  0.    ,  0.    , ..., -0.0517,  3.885 ,  1.8274],
       [ 0.    ,  0.    ,  0.    , ..., -3.5568,  1.8274,  7.321 ]])

Stiffness Matrix Calculations

Mechanical stiffness calculations for selected group of atoms can be performed using ANM.buildSM() method:

In [11]: anm.buildMechStiff(calphas)
ImportError                               Traceback (most recent call last)
<ipython-input-11-1412d885312d> in <module>()
----> 1 anm.buildMechStiff(calphas)

/home/jkrieger/anaconda2/lib/python2.7/site-packages/prody/dynamics/anm.pyc in buildMechStiff(self, coords, n_modes, kbt)
    421         sm = np.zeros((n_atoms, n_atoms), np.double)
--> 422         from .smtools import calcSM
    423'Calculating stiffness matrix.')

ImportError: /home/jkrieger/anaconda2/lib/python2.7/site-packages/prody/dynamics/ file too short

In [12]: anm.getStiffness()

Mechanical stiffness matrix is avaliable using ANM.getStiffness() method. To save stiffness matrix as an image map use following function:

In [13]: showMechStiff(anm, calphas, 'jet_r')

Note that ‘jet_r’ will reverse the colormap of image map which will be similar to coloring method of VMD program.

Mean value of mechanical stiffness matrix can be calculated using showMeanStiff() function where the secoundary structure of protein is drawing using header information.

In [14]: showMeanMechStiff(anm, calphas, header, 'A', 'jet_r')

Mechanical Stiffness in VMD

To generate tcl file for VMD program with mechanical striffness calculations use writeVMDstiffness() method. Select one residue in indices ([3]) or series of residues ([3, 7], means from 3 aa to 7 aa including) and a range of effective spring constant k_range ([0, 7.5]). This faunction required also pdb with complete protein structure which will be used in VMD representation. If calphas instead of full protein structure will be used in this function, the representation of protein in VMD program will not be accurate. In this example we considered chain A therefore suitable selection will be used:

In [15]: pdb ='chain A')
In [16]: writeVMDstiffness(anm, pdb, [3,7], [0,7.5], filename='1gfl_3-7aa',
   ....:                                                      loadToVMD=False)

In [17]: writeVMDstiffness(anm, pdb, [3], [0,7], filename='1gfl_3')

Results will be loaded automatically to VMD. Use loadToVMD=False to change it. TCL file will be saved automatically and can be used later by using linux command line:

:: vmd -e 1gfl_3aa.tcl

or in VMD TKConsole (VMD Main) for Linux, Windows and Mac users: :: play 1gfl_3aa.tcl

Tcl file contains drawing line method between selected pairs of residues which are highlighted as a VDW spheres. Color of the line can be modified by changing draw color red line in output file. Only colors from VMD Coloring Method will worked. Other changes can be done in VMD program in Graphical Representations section.


GFP results from vmdfile.writeVMDstiffness() method opened VMD. Pair of found residues: LYS3-GLY116, LYS3-PRO211 and PRO211-ASN212 are shown as VDW sphesres connected with red line.

Additionally, 1gfl_3aa.txt file is created. It contains a list of residue pairs with the value of effective spring constant (in a.u. because kbT=1) obtained from ANM.buildSM() method.

LYS3    GLY116  6.91650667766
LYS3    PRO211  6.85989128668
LYS3    ASN212  6.69507284967

The range of spring constant for k_range can be check:

In [18]: anm.getStiffnessRange()

See also ANM.getMechStiffStatistic() and ANM.getStiffnessRangeSel() function for detailed analysis of stiffness matrix.

The results of mean value of mechanical stiffness calculation can be seen in VMD program using:

In [19]: writeDeformProfile(anm, pdb, selstr='chain A and name CA',\
   ....:                                pdb_selstr='protein')

Calculate Distribution of Deformation

Distribution of the deformation in the distance contributed by each mode for selected pair of residues has been described in [EB08], see Eq. (10) and plots are shown on Fig. (2). The results can be received using plotting.showPairDeformationDist() to obtain a plot or analysis.calcPairDeformationDist() to receive a list with data that can be modified.

In [20]: calcPairDeformationDist(anm, calphas, 3, 132)

In [21]: showPairDeformationDist(anm, calphas, 3, 132)

Distribution of the deformation plot between 3-132 residue in each mode k.

To obtain results without saving any file typed:

In [22]: d1 = calcPairDeformationDist(anm, calphas, 3, 212)

In [23]: d2 = calcPairDeformationDist(anm, calphas, 132, 212)

In [24]: print d1[0], d1[1]

In [25]: plot(d1[0], d1[1], 'k-', d2[0], d2[1], 'r-')