Normal Mode Algebra

This part shows how to use some handy features of Mode objects.

ANM Calculations

We will compare modes from two ANMs for the same protein, but everything applies to comparison of ANMs and PCAs (as long as they contain same number of atoms).

Let’s get started by getting ANM models for two related protein structures:

In [1]: from prody import *

In [2]: str1 = parsePDB('1p38')

In [3]: str2 = parsePDB('1r39')

Find and align matching chains

In [4]: matches = matchChains(str1, str2)

In [5]: match = matches[0]

In [6]: ch1 = match[0]

In [7]: ch2 = match[1]

Minimize RMSD by superposing ch2 onto ch1:

In [8]: ch2, t = superpose(ch2, ch1)  # t is transformation, already applied to ch2

In [9]: calcRMSD(ch1, ch2)
Out[9]: 0.8984016339868081

Get ANM models for each chain

In [10]: anm1, ch1 = calcANM(ch1)

In [11]: anm2, ch2 = calcANM(ch2)

In [12]: anm1[0]
Out[12]: <Mode: 1 from ANM 1p38>

Let’s rename these ANM instances, so that they print short:

In [13]: anm1.setTitle('1p38_anm')

In [14]: anm2.setTitle('1r39_anm')

This is how they print now:

In [15]: anm1[0]
Out[15]: <Mode: 1 from ANM 1p38_anm>

In [16]: anm2[0]
Out[16]: <Mode: 1 from ANM 1r39_anm>

Calculate overlap

We need NumPy in this part:

In [17]: from numpy import *

Multiplication of two Mode instances returns dot product of their eigenvectors. This dot product is the overlap or cosine correlation between modes.

Let’s calculate overlap for slowest modes:

In [18]: overlap = anm1[0] * anm2[0]

In [19]: overlap
Out[19]: -0.9840211954504443

This shows that the overlap between these two modes is 0.98, which is not surprising since ANM modes come from structures of the same protein.

To compare multiple modes, convert a list of modes to a numpy.array():

In [20]: array(list(anm1[:3])) * array(list(anm2[:3]))
Out[20]: 
array([-0.9840211954504443, -0.9815834854497119, -0.9913578118318815],
      dtype=object)

This shows that slowest three modes are almost identical.

We could also generate a matrix of overlaps using numpy.outer():

In [21]: outer_product = outer(array(list(anm1[:3])), array(list(anm2[:3])))

In [22]: outer_product
Out[22]: 
array([[-0.9840211954504443, -0.14494461667593267,
        -0.0021711558324741245],
       [0.1483667882791774, -0.9815834854497119, 0.0807736109528223],
       [0.01043287216285443, -0.08407811447297674, -0.9913578118318815]],
      dtype=object)

This could also be printed in a pretty table format using printOverlapTable():

In [23]: printOverlapTable(anm1[:3], anm2[:3])
Overlap Table
                      ANM 1r39_anm
                    #1     #2     #3
ANM 1p38_anm #1   -0.98  -0.14   0.00
ANM 1p38_anm #2   +0.15  -0.98  +0.08
ANM 1p38_anm #3   +0.01  -0.08  -0.99

Scaling

Mode instances can be scaled, but after this operation they will become Vector instances:

In [24]: anm1[0] * 10
Out[24]: <Vector: (Mode 1 from ANM 1p38_anm)*10>

Linear combination

It is also possible to linearly combine normal modes:

In [25]: anm1[0] * 3 + anm1[1] + anm1[2] * 2
Out[25]: <Vector: (((Mode 1 from ANM 1p38_anm)*3) + (Mode 2 from ANM 1p38_anm)) + ((Mode 3 from ANM 1p38_anm)*2)>

Or, we could use eigenvalues for linear combination:

In [26]: lincomb = anm1[0] * anm1[0].getEigval() + anm1[1] * anm1[1].getEigval()

It is the name of the Vector instance that keeps track of operations.

In [27]: lincomb.getTitle()
Out[27]: '((Mode 1 from ANM 1p38_anm)*0.148971269751) + ((Mode 2 from ANM 1p38_anm)*0.24904210757)'

Approximate a deformation vector

Let’s get the deformation vector between ch1 and ch2:

In [28]: defvec = calcDeformVector(ch1, ch2)

In [29]: abs(defvec)
Out[29]: 16.687069727870362

Let’s see how deformation projects onto ANM modes:

In [30]: array(list(anm1[:3])) * defvec
Out[30]: 
array([-5.608605947838815, 2.1539336595932643, -3.1370160919865846],
      dtype=object)

We can use these numbers to combine ANM modes:

In [31]: approximate_defvec = sum((array(list(anm1[:3])) * defvec) *
   ....:                          array(list(anm1[:3])))
   ....: 

In [32]: approximate_defvec
Out[32]: <Vector: ((-5.60860594784*(Mode 1 from ANM 1p38_anm)) + (2.15393365959*(Mode 2 from ANM 1p38_anm))) + (-3.13701609199*(Mode 3 from ANM 1p38_anm))>

Let’s deform 1r39 chain along this approximate deformation vector and see how RMSD changes:

In [33]: ch2.setCoords(ch2.getCoords() - approximate_defvec.getArrayNx3())

In [34]: calcRMSD(ch1, ch2)
Out[34]: 0.8209600870337735

RMSD decreases from 0.89 A to 0.82 A.