Molecular Dynamics Simulations Analysis Hemochromatosis

Hemochromatosis>>Task 10: Molecular dynamics simulations analysis

Short task description

Detailed description: Molecular dynamics simulations analysis

Protocol

A protocol with a description of the data acquisition and other scripts used for this task is available here.

Molecular dynamics simulations for HFE

Note: All pictures/graphs shown here are from the first run (in case of 1a6zC[wildtype]-pictures) or second run (in case of R224W- or C282S-mutation). The reason for this is described under LINKTOMINDISTTODO.

Calculation statistics

<figtable id="tab:simulation_stats"> Statistics of the MD simulations

Input	Calc. time	Calc. speed	time to reach 1 s
Wildtype	13h31:15	17.750 ns/day	154350,8 years
C282S	13h35:05	17.667 ns/day	155075,9 years
R224W	13h35:02	17.668 ns/day	155067,1 years

</figtable>

GMXcheck revealed for all calculations that all 2001 frames were calculated, resulting in a 10ns model.

Trajectory visualization

<figtable id="tab:traj_vis">

Table 2: Visualization of the molecular dynamics simulation trajectories. The calculated states (red) are superimposed on the PDB structure of 1a6zC (green). The timeframe of 10000ps was divided into 51 frames. The trajectories from left to right are: wildtype, R224W, and C282S.

	Wildtype	R224W	C282S

</figtable>

<xr id="tab:traj_vis"/> shows the trajectories for the three MD simulations (WT, R224W, and C282S). The wildtype and R224W mostly maintain the starting structure, though some bigger parts are shifted as a whole. In both cases the upper and lower helices of the MHC I domain seem to move closer together and the C1 domain (right part in the figure) shifts around, but retains its general structure. The trajectory for C282S exhibits much less movement as the previous two which are quite flexible. Towards the end of the simulation there is a deformation of the C1 domain which is most likely caused by the missing disulfide bond due to the mutation. This deformation presumably prevents the later (but crucial) binding of beta-2-microglobulin to HFE.

Energies

Pressure

<figtable id="tab:pressure">

Table 3: different pressures of the three calculated models over time. The red line denotes the average over 100 steps (500ps). From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutatuin at position 282 (C282S)

</figtable>

The plots in <xr id="tab:pressure"/> show the pressures of the calculated systems over time. These show that, although the pressures differ greatly in some cases, the average is still at about 0 (with minor fluctuations).

Temperature

<figtable id="tab:temperature">

Table 4: different temperature energies of the three calculated models over time. The red line denotes the average over 100 steps (500ps). From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

</figtable>

The next thing we calculated were the temperatures. For all three models they can be seen in <xr id="tab:temperature"/>. The maximal deviation from the average is about 4 degrees for all models.

Potential

<figtable id="tab:potential">

Table 5: different potential energies of the three calculated models over time. The red line denotes the average over 100 steps (500ps). From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

</figtable>

With gromacs we could also extract the potentials, as can be seen in <xr id="tab:potential"/>. As in the plots before the average fluctuates around the same value for all three models themselves. However, these points differ slightly:

The average potential of the wildtype and the C282S mutation tend to be around the same (~-9.195e+05) whereas the R224W mutation potential is slightly higher (~-9.19e+05).

Total energy

<figtable id="tab:total_energy">

Table 6: different total energies of the three calculated models over time. The red line denotes the average over 100 steps (500ps). From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

</figtable>

In <xr id="tab:total_energy"/> the values of the total energies are denoted over the different states in time. Again we get an average with minor fluctuation at around the same value for each model.

All these plots show the same behavior with one exception: average around the same value and look different between the three models. This can be expected as minor changes can introduce or eradicate bindings, therefore changing the overall energies which then influence all further steps. The exception is the potential of the R224W mutation which is slightly higher than the other two models' potentials.

Minimum distance between periodic boundary cells

<figtable id="minDistRun1">

Table 7: different total energies of the three calculated models over time. From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

</figtable>

The first calculations for the mutations resulted in the minimum distances of <xr id="minDistRun1"/>. As there should be at least 2 nm distance in between at all time one can see that the mutations show the opposite. Therefore it might be possible that the protein affects itself which is not desired. To see if this states were calculated just by chance (random fluctuations that built up over time into an undesired direction) we repeated the calculations for all three models.

<figtable id="minDistRun2">

Table 8: different total energies of the three calculated models over time. The calculations are from the second run. From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

</figtable>

The resulting minimal distances can be seen in <xr id="minDistRun2"/>. We therefore decided to use the model of the first calculation for the wildtype, and the models of the second calculation for the mutation types.

RMSF for protein and C-alpha

Protein based

<figtable id="tab:rmsf_prot">

Table 9: different RMS fluctuations (based on the whole protein) of the three calculated models over time. From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

</figtable>

In general the RMS fluctuations of all three models look similar (see <xr id="tab:rmsf_prot"/>). The most differing graph is the one of R224W which shows a major peak at residues 220-235 as well as only a small peak at around residue 20.

C-Alpha based

<figtable id="tab:rmsf_ca">

Table 10: different RMS fluctuations (based on the the C-alpha atoms of the protein) of the three calculated models over time. From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

</figtable>

The c-alpha based RMSF (see <xr id="tab:rmsf_ca"/>) shows the same behavior as the whole-protein-based ones: the major differences are at around residue 20 and 220-235 of the R224W mutation.

Statistical values

With a script (provided by Dr. Marc Offman) we calculated the p-values for the RMS fluctuations to identify significant differences between the wildtype and the mutants. Surprisingly all three pairwise comparisons had very low p-values:

• Wildtype and R224W: 8.079405e-57
• Wildtype and C282S: 2.87453e-56
• R224W and C282S: 5.45069e-25

These results are to be considered with care though as the significance does not have to be caused by the mutations, but could also be caused by random events during the simulation. In fact, the wildtype alone (simulation 1 compared to simulation 2) showed up as highly significant (different).

Pymol analysis of average and bfactor

<figtable id="tab:avg_bfactor">

Table 11: Pictures of the model averages (average over MD calculated states) colored by the b-factor. The range is from blue (bfactor value beneath threshold [500]) to red (high b-factor values). From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

</figtable>

In this part we evaluate the model averages and b-factors of each position. Because it is an average over all timesteps this averaged structure can be impossible in nature.

As one can see (compare: <xr id="tab:avg_bfactor"/>) there is a big change of b-factors when comparing the wildtype and both mutations we calculated. From both mutations the R224W one shows a bigger difference to the wildtype. As expected the positions "at the edges" and those not in a secondary structure tend to have higher fluctuations/higher b-factors. It is worth noting, that even the beta sheets on the right side of the R224W picture (position 180 and higher) have (compared to wildtype and C282S mutation) pretty high b-factors. Also one can see that a little helix is inserted into the average structure of the R224W average model (right part of the picture, high b-factor [red], position ~220-225).

Radius of gyration

<figtable id="tab:gyration_ca">

Table 12: different gyrations (based on the the C-alphas of the backbone of the protein) of the three calculated models over time. From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

</figtable>

<figtable id="tab:gyration_prot">

Table 13: different gyrations (based on the the whole protein) of the three calculated models over time. From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

</figtable>

These graphs (<xr id="tab:gyration_ca"/> and <xr id="tab:gyration_prot"/>) show the radii of gyration in general as well as in each dimension. From start to end there is a slight increase in general gyration-radius for all models. However the fluctuation of the radius is the weakest for the wildtype and the C282S mutation. There seem to be most fluctuation in the R224W model. Also the R224W mutation has in general the highes radius of gyration.

Another striking difference between the models is the low radius of gyration in the Z dimension for the mutations, whereas this type of gyration is in the wildtype only low for the first ~3000ps. In exchange both mutations gain more radius of gyration in the Y dimension compared to the wildtype.

Solvent accessible surface area

<figtable id="tab:sas">

Table 14: display of the different solvent accessible surface sizes of the three calculated models over time. From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

</figtable>
The preceding graphs (see <xr id="tab:sas"/>) show the Solvent accessible surface. As the hydrophilic, total and D Gsolv values are very similar we will just talk about the hydrophobic values.

The most similar ones between those three are the wildtype and the R224W mutation which differ only in the strength of the fluctuation, of which the mutation has the higher one. Both fluctuate around about the same value.

The less similar one is the C282S hydrophobic curve. It shows about the same fluctuation strength as the wildtype but in general a decrease towards the end can be seen. Therefore this might be an indicator, that the hydrophobic residues turn to the inside, but we cant be sure as there is no increase in the other values.

<figtable id="tab:sas_res">

Table 15: display of the different solvent accessible surface sizes (normalized to per residue values) of the three calculated models over time. From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

</figtable>

In the area per residue averaged over the trajectory no clear difference can be seen (<xr id="tab:sas_res"/>).

Hydrogen-bonds between protein and protein / protein and water

Protein-Protein

<figtable id="tab:hbonds_pp">

Table 16: the number of hydrogen bonds inside the protein of the three calculated models over time. From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

</figtable>

The calculations show that in each of the three cases (see <xr id="tab:hbonds_pp"/>) the number of bonds within the protein as well as their distances tend to stay the same.

Protein-Water

<figtable id="tab:hbonds_pw">

Table 17: the number of hydrogen bonds of the protein with water of the three calculated models over time. From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

</figtable>

Although the inner-protein bonds seem to stay the same (compare section Protein-Protein), the bonds formed with hydrogen show a different behavior (compare <xr id="tab:hbonds_pw"/>).

In case of the wildtype, both the number of hydrogen bonds as well as the number of pairs within 0.35nm are (compared to the rest) fairly low at first, rising in the first ~600ps. This may be an indication that at first a dense protein state is existent. Also from ~8000-8600ps there is a drop in the number of hydrogen bonds, whereas the number of pairs within 0.35nm does not show an equal behavior. Overall the numbers tend to be at the same level each after the first 600ps rise.

The R224W model shows a different behavior over time:

The first steps show a similar "both numbers low and rising till 600ps" behavior. But instead of having a fluctuation around one constant value there seems to be a slight decrease over time as well as bigger fluctuations. This affects both, the number of bonds as well as pairs within 0.35nm.

The calculated model of the C282S mutation shows a different behavior at start than both preceding described, but a similar behavior to the R224W mutation:

For the number of hydrogen bonds there are very high fluctuations at the first 1000ps (rising till 100ps, drop to 300ps, rise till 1000ps) with a slight overall decrease afterwards like for the R224W mutation. The number of pairs within 0.35nm seem to show a similar behavior, however because of the fluctuations this is not very clear.

Ramachandran plots

<figtable id="tab:ramachandran">

Table 18: Ramachandran Plots of the three calculated models. From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

	Wildtype	R224W	C282S

</figtable>

The ramachandran plots (cf. <xr id="tab:ramachandran"/>) seem to follow the rules for the allowed and forbidden regions with a few exceptions. What is noticable though is that all of the allowed regions are almost filled to the maximum (i.e. no trend towards certain regions). Another difference are the regions above and below the left-handed alpha helix area (Psi: 110 to 180 and -180 to -150; Phi: 50 to 80) which are missing in the R224W mutation and almost non-existant in the C282S mutation. Though these areas have no significant structural element associated with them. Apart from that the plots appear almost the same, although they are hard to analyse as the dots are spread over a wide area and not clustered within distinct regions.

RMSD matrix

<figtable id="tab:rmsd_matrix_prot">

Table 19: RMSD matrices of the three calculated models over time (based on the whole protein) showing the RMSD between two models. From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

	Wildtype	R224W	C282S

</figtable>

<figtable id="tab:rmsd_matrix_mcb">

Table 20: RMSD matrices of the three calculated models over time (based on the mainchain and C-betas) showing the RMSD between two models. From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

	Wildtype	R224W	C282S

</figtable>

The RMSD matrices for the whole protein (cf. <xr id="tab:rmsd_matrix_prot"/>) and the mainchain and C-beta atoms (cf. <xr id="tab:rmsd_matrix_mcb"/>) are almost identical. The C-beta matrices have slighty lower RMSD values which is not surprising as only a subset of the atoms is taken into consideration. The small difference between the two matrix groups suggests that most of the structural changes are based on backbone rearrangements and not on orientational changes (rotations) of the residues.

The wildtype seems to periodically change between conformations as the RMSD goes up (green-yellow) and down (light blue) along the x-axis (time) for most of the different structural states (y-axis). The only noticable changes are in the at the beginning as the structure that is present in the first 1000ps seems to be quickly discarded (highest RMSD with the other states). Overall the changes are minor as the maximum RMSD is about 0.65 which is still considered quite low.

The R224W mutations appears to be very stable in the beginning, but has two structural changes towards the end. A moderate one between 6500-8000ps and a strong one between 9500-10000ps. For the rest of the time it has even less structural fluctuations than the wildtype.

The RMSD fluctuations for C282S exhibit an opposite behavior to R224W. There is a major structural change at 2000ps and some minor changes up to 3000ps. After 3000ps the structure remains almost the same until the end of the simulation. As the mutation is malign it can be assumed that it is a non-functional structure in which the proteins seems to be trapped.

Internal RMSD

<figtable id="tab:rmsd_vs_start">

Table 21: RMSD of the calculated models over time against the beginning structure. From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

	Wildtype	R224W	C282S

</figtable>

<figtable id="tab:rmsd_vs_avg">

Table 22: RMSD of the calculated models over time against the average structure (average based on all models over time). From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

	Wildtype	R224W	C282S

</figtable>

<xr id="tab:rmsd_vs_start"/> and <xr id="tab:rmsd_vs_avg"/> show the RMSD values for the different structural states compared to the starting and average structure respectively.

The first (against starting structure) in particular reflects what could already be seen in the RMSD matrix section. The wildtype quickly assumes what seems to be its average structure and then periodically fluctuates around it. R224W slowly deviates from the starting structure until 6500ps where it has a moderate structural change and a strong one later at around 9500ps. C282S has a strong deviation at 2000ps and then remains quite stable till the end.

When compared to the average structure the wildtype reaches its favored conformation at around 2000ps. There are some peaks around 4300ps, 5100ps, 7900ps, and 9100ps, but the RMSD alsways goes down afterwards. At the very end, though, is the highest peak which suggests a disruption of the structural equilibrium. This might be caused by the drop in the MinDist between the boundaries below 2nm at the end of the simulation and therefore be a result of unwanted interactions between the simulated proteins.

In contrast to the wildtype R224W has a higher RMSD to its average structure (0.3 compared to 0.2). This suggests a higher number of different states which are still quite similar to the starting structure (RMSD ~0.4). At the end of the simulation the equilibrium is disrupted (like in the wildtype simulation), but this time the reason is unkown as the minimum distance between boundaries is one of the highes values in this time period.

Of the three simulations C282S shows the most stable average structure. It is accomplished by two major structural changes at around 2000ps and 3000ps. After that C282S exhibits the lowest RMSD fluctuations of the trio.

Cluster analysis

<figtable id="tab:cluster_size_prot">

Table 23: Graphs showing the cluster sizes of the three models. The clustering was based on the whole protein. From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

</figtable>

<figtable id="tab:cluster_size_ca">

Table 24: Graphs showing the cluster sizes of the three models. The clustering was based on the C-alpha atoms of the protein. From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

</figtable>

<xr id="tab:cluster_size_prot"/> and <xr id="tab:cluster_size_ca"/> show a clustering of the different structural states during the simulation. A cutoff of 0.18 was used for these clusters and they were separately calculated based on the whole protein (cf. <xr id="tab:cluster_size_prot"/>) and the C-alpha atoms (cf. <xr id="tab:cluster_size_ca"/>) only. As expected there are fewer clusters (conformations) for the C-alpha atoms only due to the exclusion of residue variations (rotations).

The wildtype and C282S show a similar behavior in that the first two to three clusters represent the the majority conformations (the ones that are present the most during th simulation). There are 33 (14 C-alpha only) clusters for the wildtype and 25 (13 C-alpha only) clusters for C282S. R224W in contrast has 42 (20 C-alpha only) different clusters and the majority of the conformations is spread throughout more clusters. This suggests that there are not only more conformations, but they are also maintained for a longer time or assumed more often by the HFE mutant (i.e. it is more flexible).

<figtable id="tab:top_clusters">

Table 25: Comparison between the two most abundant clusters based on a cutoff of 0.18. The bigger cluster is colored green, the smaller red. From left to right: 1a6zC (wildtype), mutation at position 224 (R224W) and mutation at position 282 (C282S)

	Wildtype	R224W	C282S

</figtable>

The two biggest clusters for each simulation are shown in <xr id="tab:top_clusters"/>. The clusters are aligned by the helices only which explains the good match for the MHC I domain, but it also helps to see the relative motions between both domains. The wildtype and C282S again are quite similar in that the both clusters have about the same structure. R224W in contrast shows quite strong inter-domain movement for the two most favored conformations. This further supports the suspected increased flexibility of the R224W mutant.

Conclusion

When comparing the mutants' simulations to the wildtype, R224W seems to show a bit more differences than C282S. Especially the whole protein's flexibility seems to be increased, but if this is a blessing or a curse is yet unkown. C282S would seem to be the less harmfull mutation of the two based on the MD results, yet it is known to be one of the more malicious ones.

This task showed the power of molecular dynamics simulations, but it also showed its weakness. While it provides a lot of data one could easily be deluded by the results of a single run (as the differences between our two runs showed). In order to get reliable result multiple simulations (5+) should be performed to get some kind of average or consensus and to eliminate outliers.

References

<references/>