Tutorial:Overlap Script

From Harlem

In this tutorial we are going to use the module of intermolecular interactions and a coarse grained model of alpha-helical membrane proteins. Alpha-helices are represented by rigid bodies, the interactions between aminoacids of different helices are given by DFIRE core. We use volume constraints between helices with the purpose to avoid overlapped conformations. The molecule is Glycophorin A (1AF0.pdb) but we only take into account the alpha-helix transmembrane domains from the crystal structure.

You can download the initial files here Media:crystal_2.zip

1. We load the Crystal_2.hlm file in Harlem and run the script RunMC_Empirical.py using Harlem. This a Monte Carlo (MC) simulation script, in this case each rigid body will be translated and rotated during 10000 MC steps and each random conformation will be accepted according to the Metropolis Monte Carlo method.
2. From the previous step, a trajectory file all_mc_pts.dat is generated. It contains the coordinates of the center of mass and orientation (expressed in quaternions) of each helix during the MC simulation.
3. In this example we are interested in calculate the root medium square deviation (rmsd of CA atoms) of the different conformations generated by Metropolis MC. We can calculate the rmsd, relative to the crystal structure of glycophorin A, and pairwise energy running overlap_trajectory.py.
4. Plot rmsd versus energy.


The arrow indicates the minimum energy conformation.
5. From the graph we can identify the minimum energy conformation for the rigid body model of glycophorin A.
6. Homework, create a python script that saves the minimum energy conformation from the MC trajectory. You can use as a basis the overlap_trajectory.py
7. The minimum energy conformation can be saved as minimum.pdb and the crystal structure as initial.pdb. We load both conformations on Harlem, the image will be like this:



Crystal and minimum energy conformation of glycophorin A (grey and cyan ribbons respectively)

8. We can apply the next command on Harlem command view:

overlapmol($MINIMUM$*.ca)($initial$*.ca)

9. The value of rmsd is 4.126 A and the final figure is:



10. It is possible to calculate the energy value of the minimum conformation by applying the next script:

pmset = GetCurMolSet()
empmod = pmset.GetEmpiricalMod(1)
empmod.Initialize()
total = empmod.ScoreEnergy()