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.

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 by parsing chain A of PDB structure 1gfl together with the header, which will be used later for visualizing the secoundary structure.

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

In [5]: gfp
Out[5]: <AtomGroup: 1gflA (1967 atoms)>

We want to use only Cα atoms for our calculations, so we select them into a new object calphas:

In [6]: calphas = gfp.ca

In [7]: calphas
Out[7]: <Selection: 'ca' from 1gflA (230 atoms)>

Build Hessian and calculate ANM modes

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α atoms) to ANM.buildHessian() method:

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

And calculate all anm modes as they will be needed for later:

In [10]: anm.calcModes(n_modes='all')

All modes are required to perform accurate mechanical stiffness calculations.

Stiffness Matrix Calculations

Mechanical stiffness calculations for the selected group of atoms can be performed using the calcMechStiff() function:

In [11]: stiffness = calcMechStiff(anm, calphas)

In [12]: stiffness
Out[12]: 
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.        ]])

To show the stiffness matrix as an image map use the following function:

In [13]: show = showMechStiff(stiffness, calphas, cmap='jet_r')

Note that ‘jet_r’ is the reverse of the jet colormap and is similar to the default coloring method of the VMD program.

The mean values of the mechanical stiffness matrix for each residue can be calculated using the showMeanMechStiff() function where the secoundary structure of the protein is drawn using header information.

In [14]: show = showMeanMechStiff(stiffness, calphas, header, 'A', cmap='jet_r')

Mechanical Stiffness in VMD

We can generate tcl files for visualizing mechanical stiffness with VMD using the writeVMDstiffness() function. Select one residue in indices ([3]) or series of residues ([3, 7] means from 3 aa to 7 aa inclusive) and a range of effective spring constant k_range ([0, 7.5]).

We provide gfp as well as calphas so VMD has information about the complete protein structure, which it can use for graphical representations.

In [15]: writeVMDstiffness(stiffness, gfp, [3,7], [0,7.5], filename='1gfl_3-7aa')

In [16]: writeVMDstiffness(stiffness, gfp, [3], [0,7], filename='1gfl_3aa')

A TCL file will be saved automatically and can be used later in VMD by running the following command line instruction. Results can be loaded automatically to VMD by setting keyword loadToVMD=True.

:: vmd -e 1gfl_3aa.tcl

or typing the following in the VMD TKConsole (VMD Main) for Linux, Windows and Mac users:

:: play 1gfl_3aa.tcl

The tcl file contains a method for drawing lines between selected pairs of residues, which are highlighted as spheres. The color of the line can be modified by changing the draw color red line in the output file. Only colors from VMD Coloring Method will work. Other changes can be done within VMD in the Graphical Representations menu.

The figure shows GFP results from writeVMDstiffness() function opened in VMD. Pairs of found residues LYS3-GLY116, LYS3-PRO211 and PRO211-ASN212 are shown as VDW spheres connected with red lines.

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 calcMechStiff().

LYS3    GLY116  6.91650667766
LYS3    PRO211  6.85989128668
LYS3    ASN212  6.69507284967
...

The range of spring constant for k_range can be checked as follows:

In [17]: calcStiffnessRange(stiffness)
Out[17]: (5.271418245781293, 17.38490944653242)

See also calcMechStiffStatistic() and calcStiffnessRangeSel() functions for more detailed analysis of the stiffness matrix.

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

In [18]: writeDeformProfile(stiffness, gfp, select='chain A and name CA', pdb_selstr='protein')

Calculate Distribution of Deformation

The distribution of the deformation in the distance contributed by each mode for a selected pair of residues has been described in [EB08], see Eq. (10) and plots are shown on Fig. (2).

These results can be plotted using plotting.showPairDeformationDist() or a list can be obtained using analysis.calcPairDeformationDist().

In [19]: d0 = calcPairDeformationDist(anm, calphas, 3, 132)

In [20]: show = showPairDeformationDist(anm, calphas, 3, 132)

Figure shows the plotted distribution for deformations between 3-132 residue in each mode k.

To obtain results without saving any file type:

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

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

In [23]: plot(d1[0], d1[1], 'k-', d2[0], d2[1], 'r-')
Out[23]: 
[<matplotlib.lines.Line2D at 0x7fd31191ca90>,
 <matplotlib.lines.Line2D at 0x7fd31191cb10>]