MechStiff tutorial

This is an example how to use Mechanical Stiffness Calculations implemented in ProDy package.

The theory of MechStiff has been describe in:

Eyal E, Bahar I (2008) Toward a molecular understanding of the anisotropic response of proteins to external forces: insights from elastic network models Biophys J 94: 3424-35. PMID: 18223005; PMC2292382

http://www.ncbi.nlm.nih.gov/pubmed/18223005?ordinalpos=5&itool=EntrezSystem2.PEntrez.Pubmed.Pubmed_ResultsPanel.Pubmed_RVDocSum



In [1]:

    
%matplotlib inline
from prody import *
import matplotlib.pylab as plt



In [2]:

    
gfp, header = parsePDB('1gfl', header=True)









    



@> Connecting wwPDB FTP server RCSB PDB (USA).
@> 1gfl downloaded (1gfl.pdb.gz)
@> PDB download via FTP completed (1 downloaded, 0 failed).
@> 3950 atoms and 1 coordinate set(s) were parsed in 0.04s.



In [3]:

    
gfp









    Out[3]:





<AtomGroup: 1gfl (3950 atoms)>



In [4]:

    
calphas = gfp.select('protein and chain A and name CA')



In [5]:

    
calphas









    Out[5]:





<Selection: 'protein and chain A and name CA' from 1gfl (230 atoms)>

Building Hessian matrix using ANM model

We instantiate an ANM instance and we are building Hessian matrix by passing selected atoms.



In [6]:

    
anm = ANM('GFP ANM analysis')



In [7]:

    
anm.buildHessian(calphas, cutoff=13.0)









    



@> Hessian was built in 0.09s.



In [8]:

    
anm.getHessian().round(3)









    Out[8]:





array([[ 9.017,  2.741, -0.979, ...,  0.   ,  0.   ,  0.   ],
       [ 2.741,  5.32 ,  1.181, ...,  0.   ,  0.   ,  0.   ],
       [-0.979,  1.181,  5.663, ...,  0.   ,  0.   ,  0.   ],
       ..., 
       [ 0.   ,  0.   ,  0.   , ...,  4.794, -0.052, -3.557],
       [ 0.   ,  0.   ,  0.   , ..., -0.052,  3.885,  1.827],
       [ 0.   ,  0.   ,  0.   , ..., -3.557,  1.827,  7.321]])

Mechanical Stiffness Matrix Calculations

We can perfom MechStiff calculations for selected group of 230 Ca's atoms.



In [9]:

    
anm.buildMechStiff?



In [10]:

    
anm.buildMechStiff(calphas)









    



@> 690 modes were calculated in 0.17s.
@> Calculating stiffness matrix.
@> Stiffness matrix calculated in 0.23s.
@> The range of effective force constant is: 5.27141824578 to 17.3849094465.



In [11]:

    
anm.getStiffness()









    Out[11]:





array([[ 0.        ,  8.54744466,  8.55133411, ...,  8.73270017,
         8.36529536,  7.42692504],
       [ 8.54744466,  0.        ,  8.89990284, ...,  8.98097054,
         8.61849698,  7.77257692],
       [ 8.55133411,  8.89990284,  0.        , ...,  8.87863489,
         8.55675667,  7.84521673],
       ..., 
       [ 8.73270017,  8.98097054,  8.87863489, ...,  0.        ,
         8.87731778,  8.36590599],
       [ 8.36529536,  8.61849698,  8.55675667, ...,  8.87731778,
         0.        ,  9.52997339],
       [ 7.42692504,  7.77257692,  7.84521673, ...,  8.36590599,
         9.52997339,  0.        ]])

How to take a part of MechStiff matrix:



In [12]:

    
cut = anm.getStiffness()[5:15, 20:40]



In [13]:

    
plt.imshow(cut)
plt.colorbar()









    Out[13]:





<matplotlib.colorbar.Colorbar at 0x7f31b35de990>

To obtain matrix with effective spring constant values use:



In [14]:

    
showMechStiff(anm, calphas, 'jet_r')









    



@> 690 modes were calculated in 0.16s.
@> Calculating stiffness matrix.
@> Stiffness matrix calculated in 0.22s.
@> The range of effective force constant is: 5.27141824578 to 17.3849094465.






    Out[14]:





(<matplotlib.image.AxesImage at 0x7f31af918690>,
 <matplotlib.colorbar.Colorbar at 0x7f31af8a4c50>)

Mean value of spring constant is also significant and it can be obtain by using showMeanMechStiff() function. Arrows and trianglars represents beta strands and helixes. This information can be also seen on the 3D structure of protein using writeDeformProfile().



In [15]:

    
showMeanMechStiff(anm, calphas, header, 'A', 'jet_r')









    



@> The range of mean effective force constant is: 8.48582065351 to 12.5870266927.






    Out[15]:





<function matplotlib.pyplot.show>

Mechanical Stiffness in VMD

Mechanical Stiffness results can be seen in VMD program using:

writeVMDstiffness() - shows pair of residues with selected range of spring constant
writeDeformProfile() - deformability profile will be loaded to VMD program



In [16]:

    
pdb = gfp.select('chain A')



In [17]:

    
writeVMDstiffness(anm, pdb, [3,7], [0,7.5], filename='1gfl_3-7aa', loadToVMD=False)









    



@> Creating VMD file.






    Out[17]:





<closed file '1gfl_3-7aa.tcl', mode 'w' at 0x7f31b3545e40>



In [18]:

    
writeVMDstiffness?



In [19]:

    
writeVMDstiffness(anm, pdb, [3], [0,7], filename='1gfl_3', loadToVMD=False)









    



@> Creating VMD file.






    Out[19]:





<closed file '1gfl_3.tcl', mode 'w' at 0x7f31b3545ed0>



In [20]:

    
anm.getStiffnessRange()









    Out[20]:





(5.2714182457812999, 17.384909446532447)



In [21]:

    
writeDeformProfile(anm, pdb, selstr='chain A and name CA', pdb_selstr='protein', loadToVMD=False)









    



@> PDB file with deformability profile has been saved.
@> Creating TCL file.

Calculate Distribution of Deformation

Distribution of the deformations in the distance d contributed by each mode k in the presence of extensional forces applied to residues i and j. In this example it will be between residue nr 3 and 132.



In [22]:

    
showPairDeformationDist(anm, calphas, 3, 132)









    



@> Deformation was calculated in 0.63s.






    Out[22]:





<function matplotlib.pyplot.show>

How to put multiple results on one plot

The results can be used directly from the function. Below the example how to show to results on one plot.



In [23]:

    
import matplotlib
import matplotlib.pylab as plt



In [24]:

    
D1 = calcPairDeformationDist(anm, calphas, 3, 212)









    



@> Deformation was calculated in 0.65s.



In [25]:

    
D2 = calcPairDeformationDist(anm, calphas, 132, 212)









    



@> Deformation was calculated in 0.75s.



In [26]:

    
matplotlib.rcParams['font.size'] = '16'
fig = plt.figure(num=None, figsize=(12,8), facecolor='w')
plt.plot(D1[0], D1[1], 'k-', D2[0], D2[1], 'r-')
plt.xlabel('mode (k)', fontsize = '18')
plt.ylabel('d$^k$' '($\AA$)', fontsize = '18')
plt.grid()
plt.show()



In [ ]: