ANM calculations are carried out in the same way as GNM calculations except that the ANM uses the Hessian matrix rather than the Kirchhoff matrix. There are also additional analysis options arising from the modes being 3D, such as visualization thereof in VMD using the normal mode wizard (NMWiz) and the comparison of ANM modes with the transition vectors between experimental structures.
from prody import *
from numpy import *
from matplotlib.pyplot import *
%matplotlib inline
confProDy(auto_show=False)
confProDy(auto_secondary=True)
@> ProDy is configured: auto_show=False @> ProDy is configured: auto_secondary=True
First we parse a structure from the PDB. We will use the open structure of adenylate kinase and also parse the closed structure for comparison afterwards.
open_aa = parsePDB('4ake', compressed=False)
@> PDB file is found in working directory (4ake.pdb). @> 3459 atoms and 1 coordinate set(s) were parsed in 0.08s. @> Secondary structures were assigned to 279 residues.
open_ca = open_aa.select('calpha and chain A')
open_ca
<Selection: 'calpha and chain A' from 4ake (214 atoms)>
closed_aa = parsePDB('1ake', compressed=False)
@> PDB file is found in working directory (1ake.pdb). @> 3804 atoms and 1 coordinate set(s) were parsed in 0.11s. @> Secondary structures were assigned to 306 residues.
closed_ca = closed_aa.select('calpha and chain A')
closed_ca
<Selection: 'calpha and chain A' from 1ake (214 atoms)>
Next, we instantiate an object of ANM class for the open structure.
anm_open_ca = ANM('open_AKE')
We then build the Hessian matrix and calculate normal modes.
anm_open_ca.buildHessian(open_ca)
@> Hessian was built in 0.16s.
anm_open_ca.calcModes()
@> 20 modes were calculated in 0.09s.
We can show mode shapes, square fluctuations, and cross-correlations as we did with GNM. Let's start with the first nonzero mode:
showMode(anm_open_ca[0]);
legend();
showSqFlucts(anm_open_ca[0]);
If we want to calculate mean square fluctuations using computed modes, then we do not need to select any specific mode.
showSqFlucts(anm_open_ca); # MSFs based on "computed modes"
showCrossCorr(anm_open_ca[0]);
# cross-correlation based on "computed modes"
showCrossCorr(anm_open_ca);
ANM modes can be visualized with the VMD plugin NMWizard. This plugin can be found by navigating from VMD Main menu --> Analysis --> Normal Mode Wizard.
In order to visualize the ANM modes in NMWizard, we need to write them into a .nmd file using the method writeNMD
.
However, before saving them in .nmd format, let's extend the ANM model based on $C^\alpha$-atoms to all heavy atoms by the method extendModel
.
anm_open_aa, atoms_all = extendModel(anm_open_ca, open_ca, open_aa)
writeNMD('4ake_anm_aa', anm_open_aa, atoms_all)
'4ake_anm_aa.nmd'
# !vmd
First we will align the open and closed structures before visualizing them together in VMD. Since the chains of these two structures match, we do not need to find the mathcing atoms in these chains. Let's first look at the RMSD between the resolved structures.
calcRMSD(closed_ca, open_ca)
75.04657622215835
Alignment of the closed structure onto the open one can be easily done using the function superpose
:
aligned_closed_ca, T = superpose(closed_ca, open_ca)
RMSD between these structures is minimized after this alignment as we will see below:
calcRMSD(aligned_closed_ca, open_ca)
7.130700152775837
writePDB('1akeA_ca_alg.pdb', aligned_closed_ca)
'1akeA_ca_alg.pdb'
writePDB('4akeA_ca.pdb', open_ca)
'4akeA_ca.pdb'
# !vmd
We can easily calculate the deformation vector describing the change between two structures of the protein AKE, and then systematically compare it with ANM modes. Because the chains of these structures match and we have already aligned them, it is straightforward to calculate the deformation vector using calcDeformVector
.
defvec = calcDeformVector(open_ca, aligned_closed_ca)
We then show the overlap or correlation cosine between the ANM modes and the deformation vector. The cumulative overlap is the square root of the sum of squared overlaps.
showOverlap(defvec.getNormed(), anm_open_ca);
showCumulOverlap(defvec.getNormed(), anm_open_ca, 'r');
We observe that the first nonzero mode in the above figure, whose index is 0, overlaps best with the transition described by the deformation vector. Therefore, we can use this mode to generate a trajectory by the method traverseMode
. This takes steps in both directions starting from the provided structure to generate conformers along the mode chosen.
Let's first generate a trajectory along this mode based on $C^\alpha$ atoms, and then calculate RMSD between the last conformer in the trajectory and the closed structure to see how close we can get to the latter. Remember RMSD between aligned structures is 7.13 $\mathring{A}$.
traj_ca = traverseMode(anm_open_ca[0], open_ca, rmsd=2.0)
traj_ca.setAtoms(open_ca)
@> Parameter: rmsd = 2.00 A @> Parameter: n_steps = 10 @> Step size is 0.20 A RMSD @> Mode is scaled by 0.5118772591655905.
calcRMSD(aligned_closed_ca, traj_ca[-1])
5.662906247164603
Let's generate a similar trajectory for all atom this time by using the previously generated all-atom ANM model, anm_open_aa. The next step is to use traverseMode
again to generate conformers along the extended mode.
traj_aa = traverseMode(anm_open_aa[0], atoms_all)
traj_aa.setAtoms(atoms_all)
@> Parameter: rmsd = 1.50 A @> Parameter: n_steps = 10 @> Step size is 0.15 A RMSD @> Mode is scaled by 1.0679493809502514.
writePDB('4ake_model_aa.pdb', traj_aa)
'4ake_model_aa.pdb'
We can now visualize it in VMD.
# !vmd 4ake_mode1_aa.pdb