Molecular Dynamics Analysis BCKDHA

From Bioinformatikpedia
Revision as of 17:34, 5 September 2011 by Reisinger (talk | contribs) (Protein)

A brief check of results

To verified that the simulations finished properly we first use the command

  • gmxcheck -f wt.xtc

How many frames are in the trajectory file and what is the time resolution?

  • frames: 2001
  • time resolution: 5ps

How long did the simulation run in real time (hours), what was the simulation speed (ns/day) and how many years would the simulation take to reach a second?

  • real time: 9h27:35
  • simulation speed: 25.370 ns/day
  • simulation speed: 107991 years/second

Which contribution to the potential energy accounts for most of the calculations?

  • potential energy: -1.24431e+06

Visualization of results

To get a pdb file to be able to visualize the model with pymol we used the Swiss army knife gromacs tool trjconv:

  • trjconv -s wt.tpr -f wt.xtc -o protein.pdb -pbc nojump -dt 10

Quality assurance

Energy calculations

To calculate the different energies we used the command:
g_energy -f wtMD.edr -o energy.xvg
After submitting this command we had to choose the energy which should calculated.

  • Pressure: 13
  • Temperature: 12
  • Potential: 9
  • Total Energy: 11

Pressure

Figure2: Plot of the pressure of the MD simulation
Energy Average Err.Est. RMSD Tot-Drift (bar)
Pressure 1.01601 0.015 71.2152 -0.0706383

In Figure 2 the pressure of the molecular dynamic system is shown. The average value is 1.0161 bar which is shown in the table above. But as we can see the pressure ranges from about -250 bar to 250 bar. Since there is such a big range of 500 bar we are not sure if this value of 1.0161 bar is the equilibrium which should be reached or only the arithmetic average of all the values.

Temperature

Figure3: Plot of the temperature of the MD simulation
Energy Average Err.Est. RMSD Tot-Drift (K)
Temperature 297.941 0.0047 0.954498 0.00557078

In Figure 3 the temperature of the MD simulation is shown. As we can see it ranges between 294 K and 302 K so it has a very small deviation of the average value of 297.9 K. Since there is only such a small fluctuation we can see that the temperature in the system is quite stable which means that it reached an equilibrium.

Potential

Figure4: Plot of the potential of the MD simulation
Energy Average Err.Est. RMSD Tot-Drift (kJ/mol)
Potential -1.24431e+06 66 1041.57 -463.992

Figure 4 shows the potential of the md system. As we can see in the picture the potential ranges from -1.24e+06 kJ/mol to -1.25e+06 kJ/mol. Although this is a very huge range of 10000 we can see that all in all the potential is very low. This low potential indicates that the protein is quite stable. Since the structure of a protein is responsible for the functino of a protein this stability is important for the function of the protein.

Total Energy

Figure5: Plot of the total energy ofthe MD simulation
Energy Average Err.Est. RMSD Tot-Drift (kJ/mol)
Total Energy -1.02119e+06 65 1279.76 -459.819

The low potential energy already indicated that the total energy of the system has to be quite low. By looking at Figure 5 we can see that the energy is a bit higher than the potential energy but it is still very low. Additionally there is less variation in the energy since it ranges between -1.017e+06 kJ/mol and -1.025e+06 kJ/mol. Again we can say that such a low energy stands for a stable protein which indicates that the simulation was correct.

Minimum distance between periodic boundary cells

To calculate the minimum distance we used the command
g_mindist -f wtMD.xtc -s wtMD.tpr -od minimal-periodic-distance.xvg -pi
After submitting this command we chose group 1 to calculate the minimum distance for the whole protein.

Figure6: Minimum distance between periodic boundary cells

The shortest periodic distance is 1.40945 (nm) at time 6090 (ps) between atoms 25 and 6490.

Root mean square fluctuations

For calculating the RMSF of a protein each atom of this protein is compared with the calculated average stucture of the protein. By comparing them it is possible to find out how much it varies from its average position and so the flexibility of this region can be calculated. Regions with a high fluctuation are more flexible than regions with a low one.

To calculate the minimum distance we used the command
g_rmsf -f wtMD.xtc -s wtMD.tpr -o rmsf-per-residue.xvg -ox average.pdb -oq bfactors.pdb -res
After submitting this command we had to choose the group we want the RMSF to be calculated for:

  • Protein: 1
  • C-alpha: 3
RMSF for protein RMSF for C-alpha
Figure7: RMSF for protein
Figure8: RMSF for C-alpha

Figure 7 where the RMSF for the whole protein is plotted points out that mainly the beginning of the protein features the most fluctuation. This indicates that this region of the protein is the one with the most flexibility. In the middle part of the protein there are nearly no peaks which means that this part is fixed. In the end there is a peak which is a bit higher than the rest of the protein. This could suggest that this region is also a bit flexible.
In Figure 8 we did the RMSF calculations not for the whole protein but only for the c-alpha atoms. Those atoms are the central carbon atoms of the protein which means that in this case calculation only the backbone is considered. But by comparing the two plots above with each other we can see that in both cases the beginning of the protein is the part which has definitly the huges fluctuation of the average structure. So not only the side chains differ from the average structure but also the backbone which indicates for a strong flexibility.

Pymol analysis of average and bfactors

These average.pdb file was produced automatically during the calculation of the RMSF because it is needed for comparisons. This file contains the average structure of the protein. Because of the option -oq the bfactor.pdb file was produced additionally. In this file the temperature factors (bfactors) are calculated and added to a reference structure. Normally the parts of the protein which are most flexible have also the highest temperature. To find out if this is the fact in our case we used pymol to analyze the average structure and the bfactors structure. Additionally we compared the predicted average structure with the original structure of our protein.


Protein

1u5b/average 1u5b/bfactors bfactors/average
Figure9: Alignment of the original structure with the average structure
Figure10: Alignment of the original structure with the structure containing the bfactors
Figure11: Alignment of the average structure with the structure containing the bfactors
RMSD: 1.169 RMSD: 0.377 RMSD: 1.422


To find out how accurate the calculated average structure is we aligned it with the original structure of out protein (1u5b). As we can see in Figure 9 and additionally because of the RMSD value of 1.169 the superposition of the two structures is not covering perfectly. The middle part of the protein is aligned quite good. The most deviating parts are the two ends of the structures on the left and on the right side of the picture.

C-alpha

1u5b/average 1u5b/bfactors bfactors/average
Figure11: Alignment of the original structure with the C-alpha atoms of the average structure
Figure12: Alignment of the original structure with the C-alpha atoms of the structure which contains the bfactors
Figure13: Alignment of the C-alpha atoms of the average structure with the C-alpha atoms of the structure containing the bfactors
RMSD: 0.955 RMSD: 0.300 RMSD: 0.993

Radius of gyration

To calculate the radius of gyration we used the command
g_gyrate -f wtMD.xtc -s wtMD.tpr -o radius-of-gyration.xvg
After submitting this command we chose group 1 to calculate the radius of gyration for the whole protein

Radius off gyration MD BCKDHA.png

Structural analysis

First we had to use the command
trjconv -f wtMD.xtc -o wtMD_nojump.xtc -pbc nojump
This is important because the protein possibly jumpes out of the box so the trajectory has to be rebuild. This has the effect that the particles are back in the center.

Solvent accessible surface area

To calculate the solvent accessible surface area we used the command
g_sas -f wtMD_nojump.xtc -s wtMD.tpr -o solvent-accessible-surface.xvg -oa atomic-sas.xvg -or residue-sas.xvg
After submitting this command we had to choose two groups. Both times we chose protein.

Solvent accessible surface SAS over time per residue SAS over time per atom
Solvent-accessible-surface MD BCKDHA.png
Residue-sas MD BCKDHA.png
Atomic-sas MD BCKDHA.png

Hydrogen bonds

To calculate the hydrogen bonds between protein and protein and between protein and water we used the commands

echo 1 1 | g_hbond -f wtMD_nojump.xtc -s wtMD.tpr -num hydrogen-bonds-intra-protein.xvg
echo 1 12 | g_hbond -f wtMD_nojump.xtc -s wtMD.tpr -num hydrogen-bonds-protein-water.xvg

protein and protein protein and water
Hydrogen bonds intra protein MD BCKDHA.png
Hydrogen-bonds-protein-water MD BCKDHA.png

salt bridges

Ramachandran plot

To calculate the ramachandran plot we used the command
g_rama -f wtMD_nojump.xtc -s wtMD.tpr -o ramachandran.xvg

Ramachandran plot MD BCKDHA.png

Analysis of dynamics and time-averaged properties

RMSD matrix

To calculate the RMSD matrix we used the command
g_rms -s wtMD.tpr -f wtMD_nojump.xtc -f2 wtMD_nojump.xtc -m rmsd-matrix.xpm -dt 10
After submitting this command we had to choose two groups. Both times we chose protein.

Rmsd-matrix MD BCKDHA.jpg
  • Minimum: 0.000
  • Maximum: 0.579
  • Average: 0.329

cluster analysis

To calculate the cluster we used the command
echo 6 6 | g_cluster -s wtMD.tpr -f wtMD_nojump.xtc -dm rmsd-matrix.xpm -dist rmsd-distribution.xvg -o clusters.xpm -sz cluster-sizes.xvg -tr cluster-transitions.xpm -ntr cluster-transitions.xvg -clid cluster-id-over-time.xvg -cl clusters.pdb -cutoff 0.1 -method gromos -dt 10

For the analysis of the first two cluster we used pymol.

Cluster MD BCKDHA.png
Rmsd-distribution MD BCKDHA.png


542 cluster
RMSD cluster1+2: 0.709

internal RMSD

To calculate the internal RMSD we used the command
g_rmsdist -s wtMD.tpr -f wtMD_nojump.xtc -o distance-rmsd.xvg
After submitting this command we chose group 1 to calculate the internal RMSD for the protein

Distance-rmsd MD BCKDHA.png
  • rmsmax: -1000
  • rmscmax: inf