MBG wiki | RecentChanges | Blog | 2024-03-19 | 2024-03-18

More on eigenvectors and eigenvalues

Calculation of eigenvalues and eigenvectors is the first step for performing a large number of other analyses. Some of these are discussed below.

NOTE WELL : For these analyses to be valid, you have to be able to show that your simulation sufficiently sampled the molecule's configurational space. See Eigenvectors and eigenvalues, sufficient sampling issues.

Having a physical analogue in your mind for what PCA analysis is doing is not easy. A nice simplified but revelant example can be found in wikipedia [1]

Determination of collective motions

If you do a plot of the eigenvalues (xmgr carma.eigenvalues.dat) you will see that the top ~10% of eigenvectors accounts for the majority of the observed motion. You can visualise the motion due to each selected eigenvector with carma and vmd. Here it goes :

carma -v -cov -eigen -proj 1 10 1 CA.dcd CA.psf
echo  "Amplitude limits are " `sort -r -n -k 2 carma.fluctuations.dat | tail -1 | awk '{print $2}'` `sort -n -k 2 carma.fluctuations.dat | tail -1 | awk '{print $2}'`

The file carma.fluctuations.dat is organised as follows : successive lines correspond to successive frames from the DCD file. Successive columns correspond to successive eigenvectors (in decreasing order of corresponding eigenvalues). The numbers are the amplitudes with which the given eigenvectors contribute to the given frames (ie structures). Large numbers (in absolute value) imply that the corresponding eigenvector contributes significantly to the corresponding structure. The unix shell line above determines the range of amplitudes observed in the given run.

carma -verb -write -cov -eigen -play 1 <amp1> <amp2> CA.dcd CA.psf
vmd CA.psf carma.play.dcd

While you are inside VMD, you can save a pdb file containing snapshots from the carma.play.dcd (every, say, 200 frames) and use the pdb to prepare a figure with rasmol (or vmd, pymol, …) :

 Motion of HrcQb due to principal component

Essential dynamics

Once you know the principal components (and assuming that sufficient sampling has been demonstrated), you can -at least in principle- use these components to obtain all structures that are consistent with them. Carma is rather primitive in its approach : it can produce a pseudo-trajectory file containing a number of structures (but any number, eg 5 or 10 million) that are consistent with the observed (from the simulation) components. Note that the succession of structures in this file has no physical meaning (there is no such thing as 'time'), and that the amplitudes for the various eigenvectors are chosen randomly and uniformly from the observed fluctuations (contained in the carma.fluctuations.dat file). Performing the calculation is easy :

carma -v -cov -eigen -proj 1 50 1 CA.dcd CA.psf
carma -v -cov -eigen -arti 50 20000 CA.dcd CA.psf

will produce a DCD file with the name carma.arti.dcd containing 20000 structures consistent with the 50 largest components.

You should view the resulting structures with caution (and you should never try to obtain frequencies from it) : because carma choses amplitudes randomly and uniformly the sampling of conformations is unphysical. You can convince yourself for this by performing (on the carma.arti.dcd file) an analysis similar to the one discussed in Eigenvectors and eigenvalues, sufficient sampling issues. For example, a 1-2 plot would look like this

 projection of fluctuations on the 1-2 eigenvector plane for ARTI which is convincingly unnatural. Note also that the structures corresponding to the corners of the distribution shown above can be highly distorted.

Determination of correlated motions' fluctuations

When you calculate the average structure (with xplor) you also get the atomic fluctuations. These fluctuations do not differentiate between being correlated or not. After principal component analysis, it is possible to prepare graphs that show the amount of correlated fluctuations (per residue per eigenvector). Doing it :

#include <stdio.h>
#include <math.h>

main()
{
        float   x, y, z;
        int     c;

        while( ( c = scanf("%f %f %f", &x, &y, &z)) == 3 )
                printf("%10.8f\n", sqrt( x*x + y*y + z*z));
        if ( c > 0 )
                fprintf(stderr, "ERROR. Not triplets ???\n");
}

mkfluct < 1.dat > 1.fluct
xmgr 1.fluct

 Correlated fluctuations per residue per eigenvector

Entropy estimation

It is possible to use the eigenvalues of the mass weighted covariance matrix to calculate an upper bound to solute's configurational entropy (but including all atoms, even hydrogens). carma can do the trick, but because the results have not been validated, we will skip this for the time being.