Difference between revisions of "Molecular Dynamcis analysis"
(→Interaction Energy: Coulomb) |
(→Interaction Energy: Van der Waals) |
||
| Line 605: | Line 605: | ||
We also take a look at the non covalent interaction energies between the protein and the solvent |
We also take a look at the non covalent interaction energies between the protein and the solvent |
||
| − | [[File: |
+ | [[File:1a6z_282b_vanderwaals.png|300px|thumb|center|Figure 3.10: Van der Waals energy of the simulation over time.]] |
Energy Average Err.Est. RMSD Tot-Drift |
Energy Average Err.Est. RMSD Tot-Drift |
||
------------------------------------------------------------------------------- |
------------------------------------------------------------------------------- |
||
| − | LJ-SR:Protein-non-Protein - |
+ | LJ-SR:Protein-non-Protein -2086.53 9.1 137.417 9.01909 (kJ/mol) |
LJ-14:Protein-non-Protein 0 0 0 0 (kJ/mol) |
LJ-14:Protein-non-Protein 0 0 0 0 (kJ/mol) |
||
| + | The Van der Waals energy describes the close range interaction of non covalent bound atoms. Figure 3.10 shows the interaction energy of the protein with the solvent. Here we see a smaller fluctuation compared to the coulomb energy, but also here, a convergence is not obervable. The energy level is a bit lower as the one of the wild-type. As LJ-14:Protein-non-Protein is just 0, we did not include this term to the figure. |
||
| − | |||
| − | The Van der Waals energy describes the close range interaction of non covalent bound atoms. Figure 2.10 shows the interaction energy of the protein with the solvent. Here we see a smaller fluctuation compared to the coulomb energy, but also here, a convergence is not obervable. The energy level is a bit lower as the one of the wild-type. As LJ-14:Protein-non-Protein is just 0, we did not include this term to the figure. |
||
====Minimum distances between periodic images==== |
====Minimum distances between periodic images==== |
||
Revision as of 15:11, 17 September 2011
by Robert Greil and Cedric Landerer
Contents
- 1 Wildtype
- 2 Mutation 8a [C282Y]
- 3 Mutation 8b [C282S]
- 4 References
Wildtype
First of all, we checked the resulting file with gmxcheck.
- gmxcheck -f ref_md.tpr.xtc
Result
Reading frame 0 time 0.000 # Atoms 68601 Precision 0.001 (nm) Last frame 2000 time 10000.000 Item #frames Timestep (ps) Step 2001 5 Time 2001 5 Lambda 0 Coords 2001 5 Velocities 0 Forces 0 Box 2001 5
The Simulation toked 6h33:50 ans the simulation speed was 36.564 ns/day. So, to reach 1 second of simulation, we had to wait around 75061 years. But this is a bit to long for this Project, so we just used the results we got. The potential energy was fluctuating about -9.185e+05 kJ/mol with a range of about 0.15e+04 kJ/mol. These information are given in the different log-files provided by the simulation. To create the pdb file for the visualization, we used the command trjconv -s ref_md.tpr -f ref_md.tpr.xtc -o protein.pdb -pbc nojump -dt 10.
To create the images, we saved each frame in PyMol<ref>The PyMOL Molecular Graphics System, Version 1.2r3pre, Schrödinger, LLC</ref> as an image, which we converted into gif format by for file in *.png; do convert "$file" "$(basename $file .png).gif"; done. Than we were able to use gifsicle to create an animated gif with the command gifsicle *.gif -loop.
As one can see in Figure 1.1 and Figure 1.2, we have mostly a motion in space. The motion within the protein, like it is shown by the normal mode analysis is not identifiable. There is no movement of the beta-sheet or helical region against each other.
Quality assurance
For the quality assurance, we checked the convergence of different relevant parameters like the temperature, pressure or density. If they reach a thermal equilibrium, we can be sure our simulation has not to be extended any further.
Convergence of energy terms
To create the necessary files to check the convergence of energy terms we used the command:
g_energy -f ref_md.tpr.edr -o xvg/xxxx.xvg.
Temperature
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Temperature 297.942 0.0055 1.11471 0.0143589 (K)
Figure 1.3 shows the temperature over the simulation time. The energy fluctuates between 294K (20.85°C) and 302K (28.85°C) around a mean of about 298K (24.85°C). So, the protein is simulated in a temperature range at which a hypothermia would be the consequence. As the temperature stays the whole simulation time in the same range, we can rate this as a stable region, but the fluctuation is about 8K. So the region is quite large and as there is no trend observable, so we can not see any convergence.
Pressure
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Pressure 0.995939 0.014 83.256 0.0113633 (bar)
Figure 1.4 shows the pressure during the simulation. The pressure is not cnovereged as it ranges from -300 bar to 300 bar, which corresponds to about 600 atmospheres. The normal pressure in a human cell is about 0.37 bar<ref>D. A. T. Dick and Leah M. Lowenstein Osmotic Equilibria in Human Erythrocytes Studied by Immersion Refractometry</ref>, so the pressure during the simulation is not in the physiological range.
Energy
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Potential -919346 21 889.917 -54.8929 (kJ/mol) Kinetic En. 164147 3 614.133 7.91086 (kJ/mol) Total Energy -755199 19 1096.12 -46.9806 (kJ/mol)
Figure 1.5 shows the energy terms of the simulation. Included are the potential, the kinetic and the total energy. The total energy is calculated by substracting the potential energy from the kinetic energy. All energy terms reach the end of the simulation with about the same value as at the begining of the simulation. So, the fluctuation is just in a small range.
Volume
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Volume 694.068 0.012 0.512013 0.00522608 (nm^3)
The simulation volume shown in Figure 1.6 moves around 694 nm^3. A convergence is not obersvable as the fluctuation is about 4 nm^3 during the whole simulation. But as the range of fluctuation is very small, we would classify this as a stable volume.
Density
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Density 999.461 0.018 0.737307 -0.00750088 (kg/m^3)
The density of the simulation stays all the time in the same range (Figure 1.7), with just a small divergence. Therefor we can assume a stable density for the simulation.
Box
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Box-X 9.93815 5.8e-05 0.00244379 2.49479e-05 (nm) Box-Y 9.93815 5.8e-05 0.00244379 2.49479e-05 (nm) Box-Z 7.02734 4.1e-05 0.00172803 1.76268e-05 (nm)
The size of the box containing the simulation volume is very constant (Figure 2.8).
Interaction Energy: Coulomb
Here we take a look at the interaction energies between the protein and the solvent.
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Coul-SR:Protein-non-Protein -21213.5 39 373.942 197.347 (kJ/mol) Coul-14:Protein-non-Protein 0 0 0 0 (kJ/mol)
The coulomb energy is not converged at all like it is shown in Figure 1.9. So the high fluctuation means that we have a high fluctuation interaction with the surrounding solvet. That could mean, that the physiological environment is not matched at all. This could also be the reason for the strong movement in space. So the electrostatic interferences pushes the protein in one direction. As Coul-14:Protein-non-Protein is just 0, we did not include this term to the figure.
Interaction Energy: Van der Waals
We also take a look at the non covalent interaction energies between the protein and the solvent
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- LJ-SR:Protein-non-Protein -2237.6 18 140.778 93.3922 (kJ/mol) LJ-14:Protein-non-Protein 0 0 0 0 (kJ/mol)
The Van der Waals energy describes the close range interaction of non covalent bound atoms. Figure 1.10 shows the interaction energy of the protein with the solvent. Here we see a smaller fluctuation compared to the coulomb energy, but also here, a convergence is not obervable. As LJ-14:Protein-non-Protein is just 0, we did not include this term to the figure.
Minimum distances between periodic images
To create the corresponding .xvg files, we used the command g_mindist -f ref_md.tpr.xtc -s ref_md.tpr -od xvg/min_dist_XXXX.xvg -pi.
Figure 1.11 shows the minimal periodic distance on the protein calculated by the simulation. Figure 1.12 shows the minimal periodic distance on the C-alpha groups. In both cases, the boxes 1,2 and 3 have the same values over time, so just box 3 is visible in the plots.
For the C-alpha groups, the shortest periodic distance is 2.6075 (nm) at time 6015 (ps), between atoms 2513 and 4306. For the whole protein, the shortest periodic distance is 1.74506 (nm) at time 6120 (ps),between atoms 2532 and 3799. As the shortest periodic distance is whithin the whole protein smaller than just in the C-alpha groups, the smallest distance must be placed in a side chain. Also the fluctuation interval becomes smaller if we look at the whole protein. If the distance becomes shorter than the cut-off distance for the electrostatic interactions, the energy level would be increase dramaticly.
Root mean square fluctuations
To get the root mean square<ref>http://en.wikipedia.org/wiki/Root_mean_square_fluctuation</ref> (RMSF) fluctuations of the protein, we used the command:
g_rmsf -f ref_md.tpr.xtc -s ref_md.tpr -o xvg/rmsf_per_res.xvg -ox average.pdb -oq bfactors.pdb -res
Figure 1.13 indicates that the HFE protein is very stable. Just 2 smaller regions about residue 89 and about residue 196 show a larger flexibility. This is not surprising as, according to UniProt, both regions are within a loop.
Figure 1.14 shows the average structure superimposed with the b-factor structur. Both structures are collored according to there b-factors. Blue reagions are less flexible than red region. We can observe a higher difference in structure in regions with higher b-factors. The possible functional region around the three helices is shown as very stable.
Convergence of RMSD
To create all necessary files for the root mean square deviation (RMSD) calculation , we used the commands:
trjconv -f ref_md.tpr.xtc -o traj_nojump.xtc -pbc nojumpg_rms -f traj_nojump.xtc -s ref_md.tpr xvg/rmsd_XXX_vs_start.xvgused option 1,1 and option 1,4.echo 1 | trjconv -f traj_nojump.xtc -s ref_md.tpr -o protein.xtcg_rms -f protein.xtc -s average.pdb -o xvg/rmsd_all_vs_average.xvgused option 1,1 and option 1,4.
In all cases, no plateau is reached. The difference between the comparison of the starting structure and the average structure with the all atom and backbone model is clearly visible. Whereas the RMSD against the starting structure is continuously increasing (Figure 1.15, 1.17), the RMSD against the average structure shows an oscillating behavior(Figure 1.16, 1.18). But in both cases (starting and average structure) the comparisson with the all atom and backbone show the same behavior (Figure 1.15, 1.17 and Figure 1.16, 1.18).
The starting structure should be the better measure for convergence because, if the RMSD against the starting structure remains nearly constant, an equilibrium is reached, while an constant RMSD against the average structure could also mean an equal change in both directions over time, like becoming a structure closer to the average and a becoming a structure with a higher RMSD.
Convergence of radius of gyration
Gyration<ref>http://en.wikipedia.org/wiki/Radius_of_gyration</ref> is a measure of the size of an object, a surface, or an ensemble of points.
g_gyrate -f ref_md.tpr.xtc -s ref_md.tpr -o xvg/radius_gyration.xvg
There is no convergens of the gyration observable (Figure 1.19). The gyration moves in a range from 2.25 to 2.5. The X and Y directed gyrations show an comparable behavior, while the gyration in Z-direction behave totaly different.
Structural analysis: properties derived from configurations
Solvent accessible surface area
To get an idea which residues are accessible by the solvent, we calculated the Solvent accessible surface area (SASA)<ref>http://en.wikipedia.org/wiki/Accessible_surface_area</ref> by using the command:
g_sas -f traj_nojump.xtc -s ref_md.tpr -o xvg/solvent_acc_surf.xvg -oa xvg/atomic_sas.xvg -or xvg/residue_sas.xvg
Figure 1.20 shows the total solvent accessable surface area over time. The area is not changing over time, so we have no strong movement within the protein. There is no element adjustment observable. Interestingly, an equal amount of hydrophobic and hydrophilig residues are accessable. Figure 1.21 show the average accessable area of all atoms. Figure 1.22 show the same for each residue. In both cases, a core region wich is lies protected within the protein is not observable.
Hydrogen bonds
We also looked at the hydrogen bonds within the protein and the one which connect the protein with the surrounding solvent. Therefor we used the commands:
echo 1 1 | g_hbond -f traj_nojump.xtc -s ref_md.tpr -num xvg/hydro_bonds_intra_prot.xvgecho 1 12 | g_hbond -f traj_nojump.xtc -s ref_md.tpr -num xvg/hydro_bonds_prot_water.xvg
Figure 1.23 show that the number of hydrogen binds is very stable, also the number of pairs which can form a hydrogen bond is fluctuating less that is most likely because there is no strong movement of the different secondary structure elements against each other. So, there are just less hydrogen bonds which brake up or form new. Figure 1.24 shows the number of hydrogen bonds between the protein and the surrounding solvent. The number of hydrogen bonds here is about 3 times larger than within the protein. This is not surprising as the most residues are accessable by the solvent. In this case, the fluctuation is much stronger than within the protein. Most likely because of the movement of the protein in space within the solvent.
Ramachandran (phi/psi) plot
The Ramachandran plot<ref>http://en.wikipedia.org/wiki/Ramachandran_plot</ref> visulizes the phi and psi angles of the amino acid back bone. To get a plot of the hfe protein, we used the command:
g_rama -f traj_nojump.xtc -s ref_md.tpr -o xvg/ramachandran.xvg
 Figure 1.26: Standard Ramachandran plot<ref>http://www.cryst.bbk.ac.uk/PPS95/course/3_geometry/rama.html</ref>. |
We can observe a lot of angle combinations (Figure 1.25) which are not part of the standard ramachandran plot (Figure 1.26). Specialy in the positive phi range. Most residues are part of beta-sheets and alpha helices. Also, the HFE protein contains, according to Figure 1.25, many left handed alpha helices.
Analysis of dynamics and time-averaged properties
Root mean square deviations again
g_rms -s ref_md.tpr -f traj_nojump.xtc -f2 traj_nojump.xtc -m rmsd-matrix.xpm -dt 10xpm2ps -f rmsd-matrix.xpm -o rmsd-matrix.eps -rainbow blue
Figure 1.27 shows the rmsd over time, we see that the rmsd decreases over time in both dimensions. That means, that the movment after 1200 ps is less strong then before. The protein converges after 1200 ps.
Cluster analysis
Here we analysed cluster distribution and the size of the clusters. We used the command:
echo 6 6 | g_cluster -s ref_md.tpr -f traj_nojump.xtc -dm rmsd-matrix.xpm -dist xvg/rmsd_disstribution.xvg -o clusters.xpm -sz xvg/cluster_sizes.xvg -tr cluster_transitions.xpm -ntr xvg/cluster_transitions.xvg -clid xvg/cluster_id_over_time.xvg -cl clusters.pdb -cutoff 0.1 -method gromos -dt 10
470 clusters were found by g_cluster. The most clusters consists of just one structure. Just 2 clusters constits of 7 structures (Figure 1.28). The RMSD distribution shown in Figure 1.29 shows that the most structures have an RMSD larger than 0.2. This indicates a flexible structure.
Distance RMSD
To get the dRMSD, we used the command:
g_rmsdist -s ref_md.tpr -f traj_nojump.xtc -o xvg/distance_rmsd.xvg
Figure 1.30 shows the same behavior as Figure 1.15. Also no convergence is observable.
Mutation 8a [C282Y]
First of all, we checked the resulting file with gmxcheck.
- gmxcheck -f 282a_md.tpr.xtc
Result
Reading frame 0 time 0.000 # Atoms 68604 Precision 0.001 (nm) Last frame 2000 time 10000.000 Item #frames Timestep (ps) Step 2001 5 Time 2001 5 Lambda 0 Coords 2001 5 Velocities 0 Forces 0 Box 2001 5
The Simulation toked 11h32:38 ans the simulation speed was 20.790 ns/day. So, to reach 1 second of simulation, we had to wait around 131,780.954 years. But this is a bit to long for this Project, so we just used the results we got. The potential energy was fluctuating about -9.19218e+05 kJ/mol with a range of about 0.15e+04 kJ/mol. These information are given in the different log-files provided by the simulation. To create the pdb file for the visualization, we used the command trjconv -s 282a_md.tpr -f 282a_md.tpr.xtc -o protein.pdb -pbc nojump -dt 10.
To create the images, we saved each frame in PyMol as an image, which we converted into gif format by for file in *.png; do convert "$file" "$(basename $file .png).gif"; done. Than we were able to use gifsicle to create an animated gif with the command gifsicle *.gif -loop.
As one can see in Figure 2.1 and Figure 2.2, we have mostly a motion in space. The motion within the protein, like it is shown by the normal mode analysis is not identifiable. There is no movement of the beta-sheet or helical region against each other.
Quality assurance
For the quality assurance, we checked the convergence of different relevant parameters like the temperature, pressure or density. If they reach a thermal equilibrium, we can be sure our simulation has not to be extended any further.
Convergence of energy terms
To create the necessary files to check the convergence of energy terms we used the command:
g_energy -f 282a_md.tpr.edr -o xvg/xxxx.xvg.
Temperature
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Temperature 297.941 0.0052 1.11152 0.00591841 (K)
Figure 2.3 shows the temperature over the simulation time. The energy fluctuates between 294K (20.85°C) and 302K (28.85°C) around a mean of about 298K (24.85°C). So, the protein is simulated in a temperature range at which a hypothermia would be the consequence. As the temperature stays the whole simulation time in the same range, we can rate this as a stable region, but the fluctuation is about 8K. So the region is quite large and as there is no trend observable, so we can not see any convergence.
Pressure
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Pressure 1.00547 0.014 82.9753 0.0439014 (bar)
Figure 2.4 shows the pressure during the simulation. The pressure is not cnovereged as it ranges from -300 bar to 300 bar, which corresponds to about 600 atmospheres. The normal pressure in a human cell is about 0.37 bar, so the pressure during the simulation is not in the physiological range.
Energy
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Potential -919218 38 888.988 -238.954 (kJ/mol) Kinetic En. 164142 2.9 612.357 3.26073 (kJ/mol) Total Energy -755076 37 1093.34 -235.696 (kJ/mol)
Figure 2.5 shows the energy terms of the simulation. Included are the potential, the kinetic and the total energy. The total energy is calculated by substracting the potential energy from the kinetic energy. All energy terms reach the end of the simulation with about the same value as at the begining of the simulation. So, the fluctuation is just in a small range.
Volume
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Volume 694.169 0.02 0.527806 -0.0281068 (nm^3)
The simulation volume shown in Figure 2.6 moves around 694 nm^3. A convergence is not obersvable as the fluctuation is about 4 nm^3 during the whole simulation. But as the range of fluctuation is very small, we would classify this as a stable volume.
Density
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Density 999.333 0.029 0.759851 0.0403628 (kg/m^3)
The density of the simulation stays all the time in the same range (Figure 2.7), with just a small divergence. Therefor we can assume a stable density for the simulation.
Box
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Box-X 9.93864 9.7e-05 0.00251894 -0.000134026 (nm) Box-Y 9.93864 9.7e-05 0.00251894 -0.000134026 (nm) Box-Z 7.02768 6.9e-05 0.00178116 -9.47747e-05 (nm)
The size of the box containing the simulation volume is very constant (Figure 2.8).
Interaction Energy: Coulomb
Here we also take a look at the non covalent interaction energies between the protein and the solvent
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Coul-SR:Protein-non-Protein -21282.1 120 458.064 -645.984 (kJ/mol) Coul-14:Protein-non-Protein 0 0 0 0 (kJ/mol)
The coulomb energy is not converged at all like it is shown in Figure 2.9. So the high fluctuation means that we have a high fluctuation interaction with the surrounding solvet. The fluctuation is much higher than for the wild-type. That could mean, that the physiological environment is not matched at all. This could also be the reason for the strong movement in space. So the electrostatic interferences pushes the protein in one direction. As Coul-14:Protein-non-Protein is just 0, we did not include this term to the figure.
Interaction Energy: Van der Waals
We also take a look at the non covalent interaction energies between the protein and the solvent
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- LJ-SR:Protein-non-Protein -2230.49 11 135.734 -66.2062 (kJ/mol) LJ-14:Protein-non-Protein 0 0 0 0 (kJ/mol)
The Van der Waals energy describes the close range interaction of non covalent bound atoms. Figure 2.10 shows the interaction energy of the protein with the solvent. Here we see a smaller fluctuation compared to the coulomb energy, but also here, a convergence is not obervable. The energy level is a bit lower as the one of the wild-type. As LJ-14:Protein-non-Protein is just 0, we did not include this term to the figure.
Minimum distances between periodic images
To create the corresponding .xvg files, we used the command g_mindist -f 282a_md.tpr.xtc -s 282a_md.tpr -od xvg/mindist_xxxx.xvg -pi.
Figure 2.11 shows the minimal periodic distance on the protein calculated by the simulation. Figure 2.12 shows the minimal periodic distance on the C-alpha groups. In both cases, the boxes 1,2 and 3 have the same values over time, so just box 3 is visible in the plots. Figure 2.12 shows an oscillation for the minimal distance and a slightly increasing tendence of the maximum internal distance.
For the C-alpha groups, the shortest periodic distance is 2.32134 (nm) at time 9040 (ps), between atoms 2472 and 4307. For the whole protein, the shortest periodic distance is 1.55939 (nm) at time 9125 (ps), between atoms 2488 and 4355. As the shortest periodic distance is whithin the whole protein smaller than just in the C-alpha groups, the smallest distance must be placed in a side chain. Also the fluctuation interval becomes smaller if we look at the whole protein. If the distance becomes shorter than the cut-off distance for the electrostatic interactions, the energy level would be increase dramaticly.
Root mean square fluctuations
To get the root mean square (RMSF) fluctuations of the protein, we used the command:
g_rmsf -f 282a_md.tpr.xtc -s 282a_md.tpr -o xvg/rmsf_per_res.xvg -ox average.pdb -oq bfactors.pdb -res
Figure 2.13 shows a different RMS fluctuation as Figure 1.13. The region about residue 196 which is the most flexible one in the wild-type, is still flexible but not as much as in the wild-type. Therefor the region about residue 41 is mor flexible than in the wild-type. So, we see a change in flexibility mostly at the C-terminus of the protein.
Figure 2.14 shows the average structure superimposed with the b-factor structur. Both structures are collored according to there b-factors. Blue reagions are less flexible than red region. We can observe a higher difference in structure in regions with higher b-factors. In this case, the possible functional region is mostly affected. The underlying beta-sheet is now more flexible than in the wild-type. Also the alpha-helices are now more flexible. This is quite impressive as the mutation is not interacting with this region.
Convergence of RMSD
To create all necessary files for the root mean square deviation (RMSD) calculation , we used the commands:
trjconv -f 282a_md.tpr.xtc -o traj_nojump.xtc -pbc nojumpg_rms -f traj_nojump.xtc -s 282a_md.tpr -o xvg/rmsd_XXX_vs_start.xvgused option 1,1 and option 1,4.echo 1 | trjconv -f traj_nojump.xtc -s 282a_md.tpr -o protein.xtcg_rms -f protein.xtc -s average.pdb -o xvg/rmsd_XXX_vs_average.xvgused option 1,1 and option 1,4.
In all cases, no plateau is reached. The difference between the comparison of the starting structure and the average structure with the all atom and backbone model is clearly visible. Whereas the RMSD against the starting structure is continuously increasing with some peaks(Figure 2.15, 2.17), the RMSD against the average structure shows an oscillating behavior(Figure 2.16, 2.18), even more than the wild-type. But in both cases (starting and average structure) the comparisson with the all atom and backbone show the same behavior (Figure 2.15, 2.17 and Figure 2.16, 2.18).
The starting structure should be the better measure for convergence because, if the RMSD against the starting structure remains nearly constant, an equilibrium is reached, while an constant RMSD against the average structure could also mean an equal change in both directions over time, like becoming a structure closer to the average and a becoming a structure with a higher RMSD.
Convergence of radius of gyration
Gyration is a measure of the size of an object, a surface, or an ensemble of points.
g_gyrate -f 282a_md.tpr.xtc -s 282a_md.tpr -o xvg/radius_gyration.xvg
There is no convergens of the gyration observable (Figure 2.19). Just RgX and RgY reach a common point, but this also happend before about 3000ps. so, we can not be sure if these two realy converge. But like in the case of the wild-type, the Z direction shows a diffrent behavior. The gyration moves in a range from 2.25 to 2.5.
Structural analysis: properties derived from configurations
Solvent accessible surface area
To get an idea which residues are accessible by the solvent, we calculated the Solvent accessible surface area (SASA) by using the command:
g_sas -f traj_nojump.xtc -s 282a_md.tpr -o xvg/solvent_acc_surf.xvg -oa xvg/atomic_sas.xvg -or xvg/residue_sas.xvg
Figure 2.20 shows the total solvent accessable surface area over time. The area is not changing over time, so we have no strong movement within the protein. There is no element adjustment observable. Like in the wild-type, an equal amount of hydrophobic and hydrophilig residues are accessable. But here, the amount of hydrophobic accessable residues is slightly higher. Figure 2.21 show the average accessable area of all atoms. Figure 2.22 show the same for each residue. In both cases, a core region wich is lies protected within the protein is not observable.
Hydrogen bonds
We also looked at the hydrogen bonds within the protein and the one which connect the protein with the surrounding solvent. Therefor we used the commands:
echo 1 1 | g_hbond -f traj_nojump.xtc -s 282a_md.tpr -num xvg/hydro_bonds_intra_prot.xvgecho 1 12 | g_hbond -f traj_nojump.xtc -s 282a_md.tpr -num xvg/hydro_bonds_prot_water.xvg
Figure 2.23 show that the number of hydrogen binds is very stable, also the number of pairs which can form a hydrogen bond is fluctuating less that is most likely because there is no strong movement of the different secondary structure elements against each other. So, there are just less hydrogen bonds which brake up or form new. Figure 2.24 shows the number of hydrogen bonds between the protein and the surrounding solvent. The number of hydrogen bonds here is about 3 times larger than within the protein. This is not surprising as the most residues are accessable by the solvent. In this case, the fluctuation is much stronger than within the protein. Most likely because of the movement of the protein in space within the solvent.
In general, there is just a slightly obervable difference compared to the wild-type at 4500 ps and about 8000 ps in the amount of hydrogen bonds.
Ramachandran (phi/psi) plot
The Ramachandran plot visulizes the phi and psi angles of the amino acid back bone. To get a plot of the hfe protein, we used the command:
g_rama -f traj_nojump.xtc -s 282a_md.tpr -o xvg/ramachandran.xvg
We can observe a lot of angle combinations (Figure 2.25) which are not part of the standard ramachandran plot (Figure 2.26). Specialy in the positive phi range. Most residues are part of beta-sheets and alpha helices. Also, the HFE protein contains, according to Figure 2.25, many left handed alpha helices. Interessting is, that the large clusters do not spread out that much as in the case of the wild-type. We would call this plot more regular.
Analysis of dynamics and time-averaged properties
Root mean square deviations again
g_rms -s 282a_md.tpr -f traj_nojump.xtc -f2 traj_nojump.xtc -m rmsd-matrix.xpm -dt 10xpm2ps -f rmsd-matrix.xpm -o rmsd-matrix.eps -rainbow blue
Figure 2.27 shows the rmsd over time, we see a stronger decrease of the rmsd over time in both dimensions as in the case of the wild-type. This probably mean, the mutated structure reaches the equillibrium faster than the wild type.
Cluster analysis
Here we analysed cluster distribution and the size of the clusters. We used the command:
echo 6 6 | g_cluster -s 282a_md.tpr -f traj_nojump.xtc -dm rmsd-matrix.xpm -dist xvg/rmsd_disstribution.xvg -o clusters.xpm -sz xvg/cluster_sizes.xvg -tr cluster_transitions.xpm -ntr xvg/cluster_transitions.xvg -clid xvg/cluster_id_over_time.xvg -cl clusters.pdb -cutoff 0.1 -method gromos -dt 10
575 clusters were found by g_cluster. The most clusters consists of just one structure. Just 1 clusters constits of 6 structures (Figure 2.28). The RMSD distribution shown in Figure 2.29 shows that the most structures have an RMSD larger than 0.2. This indicates that the flexibility is not affected by the mutation.
Distance RMSD
To get the dRMSD, we used the command:
g_rmsdist -s 282a_md.tpr -f traj_nojump.xtc -o xvg/distance_rmsd.xvg
Figure 2.30 shows the same behavior as Figure 2.15. Also no convergence is observable. But, the dRMSD is not so spiky. The peaks about 1500 ps, 5000 ps and 7000 ps are not that high. So, Figure 2.30 is a smoother version of the data shown in Figure 2.15.
Mutation 8b [C282S]
First of all, we checked the resulting file with gmxcheck.
- gmxcheck -f 282b_md.tpr.xtc
Result
Reading frame 0 time 0.000 # Atoms 68600 Precision 0.001 (nm) Last frame 2000 time 10000.000 Item #frames Timestep (ps) Step 2001 5 Time 2001 5 Lambda 0 Coords 2001 5 Velocities 0 Forces 0 Box 2001 5
The Simulation toked 11h37:37 ans the simulation speed was 20.641 ns/day. So, to reach 1 second of simulation, we had to wait around 132996.409 years. But this is a bit to long for this Project, so we just used the results we got. The potential energy was fluctuating about -9.19215e+05 kJ/mol with a range of about 0.15e+04 kJ/mol. These information are given in the different log-files provided by the simulation. To create the pdb file for the visualization, we used the command trjconv -s 282b_md.tpr -f 282b_md.tpr.xtc -o protein.pdb -pbc nojump -dt 10.
To create the images, we saved each frame in PyMol as an image, which we converted into gif format by for file in *.png; do convert "$file" "$(basename $file .png).gif"; done. Than we were able to use gifsicle to create an animated gif with the command gifsicle *.gif -loop.
As one can see in Figure 3.1 and Figure 3.2, we have mostly a motion in space. The motion within the protein, like it is shown by the normal mode analysis is not identifiable. There is no movement of the beta-sheet or helical region against each other.
Quality assurance
For the quality assurance, we checked the convergence of different relevant parameters like the temperature, pressure or density. If they reach a thermal equilibrium, we can be sure our simulation has not to be extended any further.
Convergence of energy terms
To create the necessary files to check the convergence of energy terms we used the command:
g_energy -f 282b_md.tpr.edr -o xvg/xxxx.xvg.
Temperature
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Temperature 297.942 0.0062 1.11548 0.00209918 (K)
Figure 3.3 shows the temperature over the simulation time. The energy fluctuates between 294K (20.85°C) and 302K (28.85°C) around a mean of about 298K (24.85°C). So, the protein is simulated in a temperature range at which a hypothermia would be the consequence. As the temperature stays the whole simulation time in the same range, we can rate this as a stable region, but the fluctuation is about 8K. So the region is quite large and as there is no trend observable, so we can not see any convergence.
Pressure
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Pressure 1.00504 0.0098 82.9603 -0.0293434 (bar)
Figure 3.4 shows the pressure during the simulation. The pressure is not cnovereged as it ranges from -300 bar to 300 bar, which corresponds to about 600 atmospheres. The normal pressure in a human cell is about 0.37 bar, so the pressure during the simulation is not in the physiological range.
Energy
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Potential -919215 33 894.443 -178.792 (kJ/mol) Kinetic En. 164143 3.4 614.547 1.15673 (kJ/mol) Total Energy -755072 34 1102.52 -177.636 (kJ/mol)
Figure 3.5 shows the energy terms of the simulation. Included are the potential, the kinetic and the total energy. The total energy is calculated by substracting the potential energy from the kinetic energy. All energy terms reach the end of the simulation with about the same value as at the begining of the simulation. So, the fluctuation is just in a small range.
Volume
Energy Average Err.Est. RMSD Tot-Drift
Volume 694.06 0.026 0.544105 -0.162848 (nm^3)
The simulation volume shown in Figure 3.6 moves around 694 nm^3. A convergence is not obersvable as the fluctuation is about 4 nm^3 during the whole simulation. But as the range of fluctuation is very small, we would classify this as a stable volume.
Density
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Density 999.395 0.037 0.783482 0.234471 (kg/m^3)
The density of the simulation stays all the time in the same range (Figure 3.7), with just a small divergence. Therefor we can assume a stable density for the simulation.
Box
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Box-X 9.93812 0.00012 0.00259699 -0.000777245 (nm) Box-Y 9.93812 0.00012 0.00259699 -0.000777245 (nm) Box-Z 7.02731 8.8e-05 0.00183635 -0.000549606 (nm)
The size of the box containing the simulation volume is very constant (Figure 3.8).
Interaction Energy: Coulomb
Here we also take a look at the non covalent interaction energies between the protein and the solvent
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Coul-SR:Protein-non-Protein -20938.5 30 368.752 -132.83 (kJ/mol) Coul-14:Protein-non-Protein 0 0 0 0 (kJ/mol)
The coulomb energy is not converged at all like it is shown in Figure 3.9. So the high fluctuation means that we have a high fluctuation interaction with the surrounding solvet. The fluctuation is much higher than for the wild-type. That could mean, that the physiological environment is not matched at all. This could also be the reason for the strong movement in space. So the electrostatic interferences pushes the protein in one direction. As Coul-14:Protein-non-Protein is just 0, we did not include this term to the figure.
Interaction Energy: Van der Waals
We also take a look at the non covalent interaction energies between the protein and the solvent
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- LJ-SR:Protein-non-Protein -2086.53 9.1 137.417 9.01909 (kJ/mol) LJ-14:Protein-non-Protein 0 0 0 0 (kJ/mol)
The Van der Waals energy describes the close range interaction of non covalent bound atoms. Figure 3.10 shows the interaction energy of the protein with the solvent. Here we see a smaller fluctuation compared to the coulomb energy, but also here, a convergence is not obervable. The energy level is a bit lower as the one of the wild-type. As LJ-14:Protein-non-Protein is just 0, we did not include this term to the figure.
Minimum distances between periodic images
To create the corresponding .xvg files, we used the command g_mindist -f 282a_md.tpr.xtc -s 282a_md.tpr -od xvg/mindist_xxxx.xvg -pi.
Figure 2.11 shows the minimal periodic distance on the protein calculated by the simulation. Figure 2.12 shows the minimal periodic distance on the C-alpha groups. In both cases, the boxes 1,2 and 3 have the same values over time, so just box 3 is visible in the plots. Figure 2.12 shows an oscillation for the minimal distance and a slightly increasing tendence of the maximum internal distance.
For the C-alpha groups, the shortest periodic distance is 2.32134 (nm) at time 9040 (ps), between atoms 2472 and 4307. For the whole protein, the shortest periodic distance is 1.55939 (nm) at time 9125 (ps), between atoms 2488 and 4355. As the shortest periodic distance is whithin the whole protein smaller than just in the C-alpha groups, the smallest distance must be placed in a side chain. Also the fluctuation interval becomes smaller if we look at the whole protein. If the distance becomes shorter than the cut-off distance for the electrostatic interactions, the energy level would be increase dramaticly.
Root mean square fluctuations
To get the root mean square (RMSF) fluctuations of the protein, we used the command:
g_rmsf -f 282a_md.tpr.xtc -s 282a_md.tpr -o xvg/rmsf_per_res.xvg -ox average.pdb -oq bfactors.pdb -res
Figure 2.13 shows a different RMS fluctuation as Figure 1.13. The region about residue 196 which is the most flexible one in the wild-type, is still flexible but not as much as in the wild-type. Therefor the region about residue 41 is mor flexible than in the wild-type. So, we see a change in flexibility mostly at the C-terminus of the protein.
Figure 2.14 shows the average structure superimposed with the b-factor structur. Both structures are collored according to there b-factors. Blue reagions are less flexible than red region. We can observe a higher difference in structure in regions with higher b-factors. In this case, the possible functional region is mostly affected. The underlying beta-sheet is now more flexible than in the wild-type. Also the alpha-helices are now more flexible. This is quite impressive as the mutation is not interacting with this region.
Convergence of RMSD
To create all necessary files for the root mean square deviation (RMSD) calculation , we used the commands:
trjconv -f 282a_md.tpr.xtc -o traj_nojump.xtc -pbc nojumpg_rms -f traj_nojump.xtc -s 282a_md.tpr -o xvg/rmsd_XXX_vs_start.xvgused option 1,1 and option 1,4.echo 1 | trjconv -f traj_nojump.xtc -s 282a_md.tpr -o protein.xtcg_rms -f protein.xtc -s average.pdb -o xvg/rmsd_XXX_vs_average.xvgused option 1,1 and option 1,4.
In all cases, no plateau is reached. The difference between the comparison of the starting structure and the average structure with the all atom and backbone model is clearly visible. Whereas the RMSD against the starting structure is continuously increasing with some peaks(Figure 2.15, 2.17), the RMSD against the average structure shows an oscillating behavior(Figure 2.16, 2.18), even more than the wild-type. But in both cases (starting and average structure) the comparisson with the all atom and backbone show the same behavior (Figure 2.15, 2.17 and Figure 2.16, 2.18).
The starting structure should be the better measure for convergence because, if the RMSD against the starting structure remains nearly constant, an equilibrium is reached, while an constant RMSD against the average structure could also mean an equal change in both directions over time, like becoming a structure closer to the average and a becoming a structure with a higher RMSD.
Convergence of radius of gyration
Gyration is a measure of the size of an object, a surface, or an ensemble of points.
g_gyrate -f 282a_md.tpr.xtc -s 282a_md.tpr -o xvg/radius_gyration.xvg
There is no convergens of the gyration observable (Figure 2.19). Just RgX and RgY reach a common point, but this also happend before about 3000ps. so, we can not be sure if these two realy converge. But like in the case of the wild-type, the Z direction shows a diffrent behavior. The gyration moves in a range from 2.25 to 2.5.
Structural analysis: properties derived from configurations
Solvent accessible surface area
To get an idea which residues are accessible by the solvent, we calculated the Solvent accessible surface area (SASA) by using the command:
g_sas -f traj_nojump.xtc -s 282a_md.tpr -o xvg/solvent_acc_surf.xvg -oa xvg/atomic_sas.xvg -or xvg/residue_sas.xvg
Figure 2.20 shows the total solvent accessable surface area over time. The area is not changing over time, so we have no strong movement within the protein. There is no element adjustment observable. Like in the wild-type, an equal amount of hydrophobic and hydrophilig residues are accessable. But here, the amount of hydrophobic accessable residues is slightly higher. Figure 2.21 show the average accessable area of all atoms. Figure 2.22 show the same for each residue. In both cases, a core region wich is lies protected within the protein is not observable.
Hydrogen bonds
We also looked at the hydrogen bonds within the protein and the one which connect the protein with the surrounding solvent. Therefor we used the commands:
echo 1 1 | g_hbond -f traj_nojump.xtc -s 282a_md.tpr -num xvg/hydro_bonds_intra_prot.xvgecho 1 12 | g_hbond -f traj_nojump.xtc -s 282a_md.tpr -num xvg/hydro_bonds_prot_water.xvg
Figure 2.23 show that the number of hydrogen binds is very stable, also the number of pairs which can form a hydrogen bond is fluctuating less that is most likely because there is no strong movement of the different secondary structure elements against each other. So, there are just less hydrogen bonds which brake up or form new. Figure 2.24 shows the number of hydrogen bonds between the protein and the surrounding solvent. The number of hydrogen bonds here is about 3 times larger than within the protein. This is not surprising as the most residues are accessable by the solvent. In this case, the fluctuation is much stronger than within the protein. Most likely because of the movement of the protein in space within the solvent.
In general, there is just a slightly obervable difference compared to the wild-type at 4500 ps and about 8000 ps in the amount of hydrogen bonds.
Ramachandran (phi/psi) plot
The Ramachandran plot visulizes the phi and psi angles of the amino acid back bone. To get a plot of the hfe protein, we used the command:
g_rama -f traj_nojump.xtc -s 282a_md.tpr -o xvg/ramachandran.xvg
We can observe a lot of angle combinations (Figure 2.25) which are not part of the standard ramachandran plot (Figure 2.26). Specialy in the positive phi range. Most residues are part of beta-sheets and alpha helices. Also, the HFE protein contains, according to Figure 2.25, many left handed alpha helices. Interessting is, that the large clusters do not spread out that much as in the case of the wild-type. We would call this plot more regular.
Analysis of dynamics and time-averaged properties
Root mean square deviations again
g_rms -s 282a_md.tpr -f traj_nojump.xtc -f2 traj_nojump.xtc -m rmsd-matrix.xpm -dt 10xpm2ps -f rmsd-matrix.xpm -o rmsd-matrix.eps -rainbow blue
Figure 2.27 shows the rmsd over time, we see a stronger decrease of the rmsd over time in both dimensions as in the case of the wild-type. This probably mean, the mutated structure reaches the equillibrium faster than the wild type.
Cluster analysis
Here we analysed cluster distribution and the size of the clusters. We used the command:
echo 6 6 | g_cluster -s 282a_md.tpr -f traj_nojump.xtc -dm rmsd-matrix.xpm -dist xvg/rmsd_disstribution.xvg -o clusters.xpm -sz xvg/cluster_sizes.xvg -tr cluster_transitions.xpm -ntr xvg/cluster_transitions.xvg -clid xvg/cluster_id_over_time.xvg -cl clusters.pdb -cutoff 0.1 -method gromos -dt 10
575 clusters were found by g_cluster. The most clusters consists of just one structure. Just 1 clusters constits of 6 structures (Figure 2.28). The RMSD distribution shown in Figure 2.29 shows that the most structures have an RMSD larger than 0.2. This indicates that the flexibility is not affected by the mutation.
Distance RMSD
To get the dRMSD, we used the command:
g_rmsdist -s 282a_md.tpr -f traj_nojump.xtc -o xvg/distance_rmsd.xvg
Figure 2.30 shows the same behavior as Figure 2.15. Also no convergence is observable. But, the dRMSD is not so spiky. The peaks about 1500 ps, 5000 ps and 7000 ps are not that high. So, Figure 2.30 is a smoother version of the data shown in Figure 2.15.
References
<references />
