Molecular Dynamcis analysis

^{by Robert Greil and Cedric Landerer}

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.

Figure 1: Motion of the wild-type protein, Cartoon representation

Figure 2: Motion of the wild-type protein

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 and Figure 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

Figure 3: Temperature over time of the molecular dynamics simulation of the wild-type protein

Energy                      Average   Err.Est.       RMSD  Tot-Drift 
-------------------------------------------------------------------------------
Temperature                 297.942     0.0055    1.11471  0.0143589  (K)

Figure 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

Figure 4: Pressure over time of the molecular dynamics simulation of the wild-type protein

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Pressure                   0.995939      0.014     83.256  0.0113633  (bar)

Figure 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

Figure 5: Energy over time of the molecular dynamics simulation of the wild-type protein

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 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

Figure 6: Volume over time of the molecular dynamics simulation of the wild-type protein

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Volume                      694.068      0.012   0.512013 0.00522608  (nm^3)

The simulation volume shown in Figure 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

Figure 7: Densety over time of the molecular dynamics simulation of the wild-type protein

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 7), with just a small divergence. Therefor we can assume a stable density for the simulation.

Box

Figure 8: Box size in X,Y and Z direction over time of the molecular dynamics simulation of the wild-type protein. The X and Y values are equal over time.

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 8).

Interaction Energy: Coulomb

Figure 9: Coulomb energy of the simulation over time.

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 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

Figure 10: Van der Waals energy of the simulation over time.

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 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 11 shows the minimal periodic distance on the protein calculated by the simulation. Figure 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 11: Minimal periodic distance on the protein

Figure 12: Minimal periodic distance on the C-alpha groups

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 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 13: RMSF per residue of HFE in nm

Figure 14: HFE protein colored according to the b-factors. The average and the b-factor structure are superimposed.

Figure 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 nojump
• g_rms -f traj_nojump.xtc -s ref_md.tpr xvg/rmsd_XXX_vs_start.xvg used option 1,1 and option 1,4.
• echo 1 | trjconv -f traj_nojump.xtc -s ref_md.tpr -o protein.xtc
• g_rms -f protein.xtc -s average.pdb -o xvg/rmsd_all_vs_average.xvg used option 1,1 and option 1,4.

Figure 15: RMSD over time of all atoms versus the starting structure

Figure 16: RMSD over time of all atoms versus the average structure

Figure 17: RMSD over time of backbone versus the starting structure

Figure 18: RMSD over time of backbone versus the average structure

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 15, 17), the RMSD against the average structure shows an oscillating behavior(Figure 16, 18). But in both cases (starting and average structure) the comparisson with the all atom and backbone show the same behavior (Figure 15,17 and Figure 16,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

Figure 19: radius of gyration in nm

There is no convergens of the gyration observable (Figure 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: Area of solvent accessable surface over time according to hydrophilic and hydrophobic residues.

Figure 1.21: Average area of accessable surface per atom and standard deviation in nm².

Figure 1.22: Average area of accessable surface per residue and standard deviation in nm².

Figure 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 21 show the average accessable area of all atoms. Figure 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.xvg
• echo 1 12 | g_hbond -f traj_nojump.xtc -s ref_md.tpr -num xvg/hydro_bonds_prot_water.xvg

Figure 23: Hydrogen bonds within the protein in black and the number of pairs within 0.35 nm which are able to form a hydrogen bond in red.

Figure 24: Hydrogen bonds between the protein and the solvent in black and the number of pairs within 0.35 nm which are able to form a hydrogen bond in red.

Figure 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 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.

Salt bridges

Not working yet.. maybe later... maybe not

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 25: Ramachandran plot of the HFE protein.

Figure 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 XX) which are not part of the standard ramachandran plot (Figure XX). Specialy in the positive phi range. Most residues are part of beta-sheets and alpha helices. Also, the HFE protein contains, according to Figure XXX, 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 10
• xpm2ps -f rmsd-matrix.xpm -o rmsd-matrix.eps -rainbow blue

Figure 1.27: Dotplot of the RMSD of HFE against HFE over time.

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

• 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

Figure 1.28: Average area of accessable surface per atom and standard deviation in nm².

Figure 1.29: Average area of accessable surface per residue and standard deviation in nm².

Distance RMSD

• g_rmsdist -s ref_md.tpr -f traj_nojump.xtc -o xvg/distance_rmsd.xvg

Figure 1.27:

Mutation 8a [C282Y]

still waiting to finish // LRZ won't do anything because of LRZ cluster update

Mutation 8b [C282S]

still waiting to finish // LRZ won't do anything because of LRZ cluster update

References

<references />