Molecular Dynamcis analysis
by Robert Greil and Cedric Landerer
Contents
- 1 Wildtype
- 2 Mutation 8a [C282Y]
- 3 Mutation 8b [C282S]
- 4 Discussion
- 5 References
Wildtype
First of all, we checked the resultfile 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 took 6h33:50 and 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 converged 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 subtracting the potential energy from the kinetic energy. All energy terms reach the end of the simulation with about the same value as at the beginning 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 observable 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 solvent. 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 observable. 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 within 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 dramatically.
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 structure. Both structures are colored according to there b-factors. Blue regions 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 comparison 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 convergence 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 totally 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 accessible 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 hydrophilic residues are accessible. Figure 1.21 show the average accessible area of all atoms. Figure 1.22 show the same for each residue. In both cases, a core region which 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 accessible 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> visualizes 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). Specially 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 movement after 1200 ps is less strong then before. The protein converges after 1200 ps.
Cluster analysis
Here we analyzed 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 consits 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 took 11h32:38 and 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 converged 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 subtracting the potential energy from the kinetic energy. All energy terms reach the end of the simulation with about the same value as at the beginning 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 observable 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 solvent. 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 observable. 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 tendency 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 within 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 dramatically.
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 more 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 structure. Both structures are colored according to there b-factors. Blue regions 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 comparison 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 convergence of the gyration observable (Figure 2.19). Just RgX and RgY reach a common point, but this also happened before about 3000ps. so, we can not be sure if these two really converge. But like in the case of the wild-type, the Z direction shows a different 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 accessible 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 hydrophilic residues are accessible. But here, the amount of hydrophobic accessible residues is slightly higher. Figure 2.21 show the average accessible 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 accessible 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 observable 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 visualizes 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). Specially 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. Interesting 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 equilibrium faster than the wild type.
Cluster analysis
Here we analyzed 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 consists 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 took 11h37:37 and 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 converged 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 subtracting the potential energy from the kinetic energy. All energy terms reach the end of the simulation with about the same value as at the beginning 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 observable 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 solvent. 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 observable. 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 282b_md.tpr.xtc -s 282b_md.tpr -od xvg/mindist_xxxx.xvg -pi.
Figure 3.11 shows the minimal periodic distance on the protein calculated by the simulation. Figure 3.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 3.12 shows an oscillation for the minimal distance and a slightly increasing tendency of the maximum internal distance.
For the C-alpha groups, the shortest periodic distance is 2.94439 (nm) at time 3770 (ps), between atoms 2513 and 3324. For the whole protein, the shortest periodic distance is 2.17028 (nm) at time 2345 (ps), between atoms 2532 and 3330. As the shortest periodic distance is within 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 dramatically.
Root mean square fluctuations
To get the root mean square (RMSF) fluctuations of the protein, we used the command:
g_rmsf -f 282b_md.tpr.xtc -s 282b_md.tpr -o xvg/rmsf_per_res.xvg -ox average.pdb -oq bfactors.pdb -res
Figure 3.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, shows here the same flexibility like in the wild-type. The region about residue 86 is no less flexible than in the wild-type.
Figure 3.14 shows the average structure superimposed with the b-factor structure. Both structures are colored according to there b-factors. Blue regions are less flexible than red region. We can observe a higher difference in structure in regions with higher b-factors. This mutation shows less impact than the other one on the same position. The flexibility is all in all pretty much the same like in the wild-type.
Convergence of RMSD
To create all necessary files for the root mean square deviation (RMSD) calculation , we used the commands:
trjconv -f 282b_md.tpr.xtc -o traj_nojump.xtc -pbc nojumpg_rms -f traj_nojump.xtc -s 282b_md.tpr -o xvg/rmsd_XXX_vs_start.xvgused option 1,1 and option 1,4.echo 1 | trjconv -f traj_nojump.xtc -s 282b_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 3.15, 3.17), the RMSD against the average structure increases at the beginning, but drops about 800 ps to its initial conditions (Figure 3.16, 3.18). But in both cases (starting and average structure) the comparison with the all atom and backbone show the same behavior (Figure 3.15, 3.17 and Figure 3.16, 3.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 282b_md.tpr.xtc -s 282b_md.tpr -o xvg/radius_gyration.xvg
There is no convergence of the gyration observable (Figure 2.19). Just RgZ and RgY reach a common point. The gyration moves in a range from 2.25 to 2.5. The behavior of RgZ is not comparable with the other mutation and the wild-type.
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 282b_md.tpr -o xvg/solvent_acc_surf.xvg -oa xvg/atomic_sas.xvg -or xvg/residue_sas.xvg
Figure 3.20 shows the total solvent accessible 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 hydrophilic residues are accessible. So, it is not possible to differentiate between the wild-type and this mutation. Figure 3.21 show the average accessible area of all atoms. Figure 3.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 282b_md.tpr -num xvg/hydro_bonds_intra_prot.xvgecho 1 12 | g_hbond -f traj_nojump.xtc -s 282b_md.tpr -num xvg/hydro_bonds_prot_water.xvg
Figure 3.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 3.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 accessible 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. The results is nearly the same as for the wild-type.
Ramachandran (phi/psi) plot
The Ramachandran plot visualizes 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 282b_md.tpr -o xvg/ramachandran.xvg
We can observe a lot of angle combinations (Figure 3.25) which are not part of the standard ramachandran plot (Figure 3.26). Specially in the positive phi range. Most residues are part of beta-sheets and alpha helices. Also, the HFE protein contains, according to Figure 3.25, many left handed alpha helices. Interesting is, that, also, the large clusters do not spread out that much as in the case of the wild-type. We would call this plot more regular. But compared to the other mutation, the negativ phi area is more dense.
Analysis of dynamics and time-averaged properties
Root mean square deviations again
g_rms -s 282b_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 3.27 shows the rmsd over time, the RMSD decreases about 1500 ps. This is the behavior we can observe in Figure 3.16 and 3.18. Also this matrix is comparable with the wild-type.
Cluster analysis
Here we analyzed cluster distribution and the size of the clusters. We used the command:
echo 6 6 | g_cluster -s 282b_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
400 clusters were found by g_cluster. The most clusters consists of just one structure. Just 1 clusters consists of 9, 8 and 7 structures (Figure 3.28). The RMSD distribution shown in Figure 3.29 shows that the most structures have an RMSD larger than 0.2. This indicates that the flexibility is not affected by the mutation. This obersavtion holds for all three simulations.
Distance RMSD
To get the dRMSD, we used the command:
g_rmsdist -s 282b_md.tpr -f traj_nojump.xtc -o xvg/distance_rmsd.xvg
Figure 3.30 shows the same behavior as Figure 3.15. Also no convergence is observable. But a bit smoother at the end of the simulation.
Discussion
In all cases, the simulation parameter, like temperature, pressure and so on behave the same. So, we assume that all simulation reached the same state and are comparable.
The Solvent accessible surface remains the same in all simulations, just in the case of 282a, we see a slightly change towards the hydrophobic residues on the surface. For the wild-type and the 282b structure, an equal amount of hydrophilic and hydrophobic residues are accessible. The high hydrophilic fraction of the surface is probably due to the fact that the HFE protein is located in the plasma membrane, whereat the membrane region is not part of the simulated structure.
The Ramachandran plots of the simulations show, that the wild-type is more diverse than the two mutated structures. A larger difference is specially in the positive phi range observable (Figures 1.25 2.25 3.25). Here, many angle combinations present in the wild-type are missing in the mutated structures.
The two dimensional Root mean square deviation over time shows in all three cases the same picture as the one dimensional counterpart. For the wild-type and the 282b structure, we see a slow divergence than in the case of 282a. WT and 282b show a larger divergence in the early states of the simulation while 282a is changing just slightly. This indicates a less flexibility in this structure and a larger one in the other two. This is exactly what one can observe in the flexibility analysis.
The flexibility is the most interesting part. While the structure 282a is less flexible, the mutated structure 282b is more flexible than the wild-type. The structure 282a has a flexibility index (sum over the rmsf) of 49.57 and the structure has a flexibility index of 53.457, while the wild-type has one of 51.4624 which is very close to the mean of 51.4966. So, the wild-type can be called a mean structure in this case. We can observe four interesting points in the flexibility distribution shown in Figure 4.1. These regions are also marked in the structures shown in Figure 4.2. The first one is within a loop connecting two beta-strands. Here the structure 282a is more flexible than the other structures, so this could causes issues by forming a beta-sheet by the two beta-strands. Therefore the stability of the whole protein could be affected. In the second case, both mutated structures are less flexible than the wild type. The position is also whiting a loop, this time connecting a beta-sheet with an alpha helix. As this helix is assumed to be part of an active site, a loss in flexibility could cause issues at binding at the transferrin protein. So, a key-hole principle can no longer act. The third case is also a transfer from a beta-sheet into an alpha-helix. Here we have an interesting change in flexibility as first, the flexibility is increased at the structure 282b and 282a follows the wild-type, while just a few residues later the flexibility changes and the structure follows the wild-type. The other mutated structure (282a) behaves the opposite. The fourth case shows a difference at a very flexible region. The mutation 282a causes a high loss in flexibility at this point, while the flexibility is not affected by the other mutation. All cases have in common that they are not connected to the mutated side in sequence and in structure. Just in case 4, the mutated residue is interacting with the beta-sheet which is connected by the flexible region marked as 4. The mutated residue interacts with a cysteine (Figure 4.3). In the wild-type structure a tyrosine is placed. In case of 282a, the tyrosine is replaced by an serin, which are both polar and neutral, and both have an OH-group to form a hydrogen bond. In the case of 282b, the thyronine is replaced by an cystein, which also could form a hydrogen bond<ref>Gregoret LM, Rader SD, Fletterick RJ, Cohen FE.: Hydrogen bonds involving sulfur atoms in proteins.</ref> but as on the other site also a cystein is placed, so a disulfide bond can be formed. This disulfide bond is more stable as the hydrogen bonds since no water can attack to break up the secondary structure. Therefor it is surprisingly, that the 282b structure is flexible as the wild-type while the flexibility is decreased by the serin in structure 282a. To get a better idea what happened here exactly, the flexibility for all other possible mutations should be calculated.
 Figure 4.1: Comparison of the flexibility per residue in nm. Also the difference of the mutated structures compared to the wild-type is shown. a value close to zero indicates no change in flexibility while a value larger zero indicate a higher flexibility in the wild-type structure and a value smaller zero a higher flexibility in the mutated structure. |
All in all, we see two different mutation which changes the flexibility all over the protein, not just in the mutated region. As the region is not interaction directly with a binding partner, and also the biochemical properties if the amino acid is not changed, the loss of function is most likely due to the change in flexibility. A major question which remains unanswered is, why are so many regions affected which are placed not even in the same domain. As the HFE forms a complex with Beta-2-Microglobulin (B2M) to stabilize itself and avoid scissors movement like it is observable in the normal mode analysis, the beta-sheet close to the mutated residue could be affected (Figure 4.3) and the stabilizing effect can not take place. It also would be very interesting to see which of the other possible amino acids lead to a gain or loss in flexibility to compare this an analyze the effect.
References
<references />
