View source for Molecular Dynamics Simulations Analysis Gaucher Disease

Technical details can be found in our [[Gaucher_Task10_Protocol|protocol]].

== A brief check of results ==
We briefly checked the length of the MD simulations via <tt>gmxcheck</tt>. Since we used the AGroS paramter <tt>-lengthMD 10</tt>, we expected the MD simulation to be 10 ns long. <tt>gmxcheck</tt> reported for each simulation 2001 frames and a time step of 5 ps per frame which amounts to the expected length of 2001*5ps=10ns.

We further investigated the runtime of the MD simulation by looking at the log files. The simulation of the WT took approximately twice as long as the other two mutations since we used only 16 instead of 32 cores. Performing a MD simulation over 1s would take more than 100000 years if the same hardware were used than we did.

<figtable id="tab:runtime">
{|style="border-collapse: separate; border-style: solid; border-spacing: 0; border-width: 2px 0 2px 0; text-align:center" width="400px"
|- style="background-color: lightgrey" 
! Name !! Runtime !! ns/day !! years/1s
|- 
| colspan=5 style="border-collapse: separate; border-style: solid; border-spacing: 0; border-width: 1px 0 0 0"|
|-
| WT || 16h25:59 || 14.605 || 187536
|-
| L470P || 9h20:49 || 25.677 || 106699
|-
| W209R || 9h18:47 || 25.770 || 106314
|}
<caption>Runtime of the MD simulations and the estimated runtime for a simulation length of 1s.</caption>
</figtable>
The computation of the energy term <tt>Coul(T) + LJ [W3-W3]</tt> accounted for most (47%) of the floating point operations.

== Visualization ==
For getting a first impression of the MD simulations, we visualized the structures stored in the trajectory files. A course grained animation of the three simulations with 29 frames is shown in <xr id="fig:viz"/>. Altogether, the animations indicate that all three simulations terminated successfully without sever structural changes. Based on the visual inspection, we could not detect any rearrangements which would suggest a malign effect of one of the mutations. However, the active site residue E235 and E340 might be more flexible in L470P compared to the wildtype.

<figure id="fig:viz">
{|style="border-collapse: separate; border-style: solid; border-spacing: 0; border-width: 1px"
|<subfigure>[[File:wt.gif|thumb|300px|<caption>Wildtype</caption>]]</subfigure>
|<subfigure>[[File:L470P.gif|thumb|300px|<caption>L470P</caption>]]</subfigure>
|<subfigure>[[File:W209R.gif|thumb|300px|<caption>W209R</caption>]]</subfigure>
|}
<small><caption>Visualization of the MD simulations. Red: the active site residues E235 and E340; Blue: the mutated residue.</caption></small>
</figure>

== Quality assurance ==
We checked the quality of the MD simulations regarding the convergence of the pressure, temperature, potential energy, and total energy. Strongly fluctuating thermodynamic parameters and insufficiently converging parameters might hint at problems during the MD simulations or a lacking equilibration.

<xr id="fig:qual_pres"/> shows the development of the pressure along the MD simulations. In all three cases, the pressure fluctuates between -200 and 200 bar but the average pressure of 1 bar does not change. In contrast the the wildtype which exhibits a negative drift, the pressure slightly increases in case of L470P and W209R.

<figure id="fig:qual_pres">
{|style="border-collapse: separate; border-style: solid; border-spacing: 0; border-width: 1px"
|<subfigure>[[File:qual_pres_wt.png|thumb|300px|<caption>Wildtype (1.0069/-0.0423231)</caption>]]</subfigure>
|<subfigure>[[File:qual_pres_L470P.png|thumb|300px|<caption>L470P (0.996977/0.0144563)</caption>]]</subfigure>
|<subfigure>[[File:qual_pres_W209R.png|thumb|300px|<caption>W209R (1.00887/0.00535799)</caption>]]</subfigure>
|}
<small><caption>Development of the pressure during the MD simulation. In brackets: average/total-drift.</caption></small>
</figure>

The temperature of 297.9 K remains constant and varies by only +/- 2 K (<xr id="fig:qual_temp"/>).

<figure id="fig:qual_temp">
{|style="border-collapse: separate; border-style: solid; border-spacing: 0; border-width: 1px"
|<subfigure>[[File:qual_temp_wt.png|thumb|300px|<caption>Wildtype (297.912/0.00443415)</caption>]]</subfigure>
|<subfigure>[[File:qual_temp_L470P.png|thumb|300px|<caption>L470P (297.907/-0.00147577)</caption>]]</subfigure>
|<subfigure>[[File:qual_temp_W209R.png|thumb|300px|<caption>W209R (297.915/0.0095005)</caption>]]</subfigure>
|}
<small><caption>Development of the temperature during the MD simulation. In brackets: average/total-drift.</caption></small>
</figure>

Both the potential energy and the total energy reduces at the beginning of the MD simulations (<xr id="fig:qual_pot"/>, <xr id="fig:qual_tot"/>). L470P shows the strongest shift towards the negative range, i.e. the stability of the structure increases most.

<figure id="fig:qual_pot">
{|style="border-collapse: separate; border-style: solid; border-spacing: 0; border-width: 1px"
|<subfigure>[[File:qual_pot_wt.png|thumb|300px|<caption>Wildtype (-939801/-450.018)</caption>]]</subfigure>
|<subfigure>[[File:qual_pot_L470P.png|thumb|300px|<caption>L470P (-940007/-705.14)</caption>]]</subfigure>
|<subfigure>[[File:qual_pot_W209R.png|thumb|300px|<caption>W209R (-940568/-355.406)</caption>]]</subfigure>
|}
<small><caption>Development of the potential energy during the MD simulation. In brackets: average/total-drift.</caption></small>
</figure>

<figure id="fig:qual_tot">
{|style="border-collapse: separate; border-style: solid; border-spacing: 0; border-width: 1px"
|<subfigure>[[File:qual_tot_wt.png|thumb|300px|<caption>Wildtype (-768901/-447.47)</caption>]]</subfigure>
|<subfigure>[[File:qual_tot_L470P.png|thumb|300px|<caption>L470P (-769094/-705.985)</caption>]]</subfigure>
|<subfigure>[[File:qual_tot_W209R.png|thumb|300px|<caption>W209R (-769675/-349.952)</caption>]]</subfigure>
|}
<small><caption>Development of the total energy during the MD simulation. In brackets: average/total-drift.</caption></small>
</figure>

Self-interactions between periodic images are unphysical and might disturb the MD simulations. We therefore checked the minimal distance to adjacent cells by taking (a) only C-alpha atoms and (b) also side-chain atoms into account (<xr id="fig:mindist"/>).

In all three simulations, the minimal distance slightly decreases but does not falls below the cut-off distance of 2 nm. If only C-alpha atoms are taking into account, the minimal distance is bout 0.5 nm higher compared to using all protein atoms. This corresponds to an average side-chain length of 5 A. For W209R, the minimal distance drops below 2 nm at 9180 ps but then increases again such that this event does not impact the MD simulation.

<figure id="fig:mindist">
{|style="border-collapse: separate; border-style: solid; border-spacing: 0; border-width: 1px"
|<subfigure>[[File:mindist_wt.png|thumb|300px|<caption>Wildtype (6860: 2.72/3820: 2.08)</caption>]]</subfigure>
|<subfigure>[[File:mindist_L470P.png|thumb|300px|<caption>L470P (365: 2.28 / 7395: 2.08)</caption>]]</subfigure>
|<subfigure>[[File:mindist_W209R.png|thumb|300px|<caption>W209R (9720: 2.61 / 9180: 1.95)</caption>]]</subfigure>
|}
<small><caption>Minimal distance to the adjacent cell with respect C-alpha atoms (black) and all protein atoms (red). In brackets: time frame and minimal distance for c-alpha/protein atoms.</caption></small>
</figure>

== Flexibility ==

After we had verified that our MD simulations terminated probably, we investigated the flexibility of the wildtype and the mutant structures. We further compared the flexibility of L470P and W209R to the wildtype for identifying changes in the flexibility which might indicate perturbations of glycosylceramidase activity. 

=== RMSF ===

For quantifying the flexibility we used the Root Mean Square Fluctation (RMSF). The RMSF of residue i is the mean square deviation from the position of i in the average structure. The average structures of all three simulations are shown in <xr id="fig:rmsf_avg"/> which were compared to the structure used for initializing the production run. Although the structures obviously changes during the simulation (<xr id="fig:viz"/>), the average structure hardly deviated from the initial structure in all three cases. Furthermore, we could not identify any significant changes between the average structure of the wildtype, L470P, and W209R.

<figure id="fig:rmsf_avg">
{|style="border-collapse: separate; border-style: solid; border-spacing: 0; border-width: 1px"
|<subfigure>[[File:rsmf_avg_wt.png|thumb|300px|<caption>Wildtype</caption>]]</subfigure>
|<subfigure>[[File:rsmf_avg_L470P.png|thumb|300px|<caption>L470P</caption>]]</subfigure>
|<subfigure>[[File:rsmf_avg_W209R.png|thumb|300px|<caption>W209R</caption>]]</subfigure>
|}
<small><caption>Superposition of the average g_rmsf (green) structure onto the initial structure (blue).</caption></small>
</figure>

<xr id="fig:rmsf_c"/> shows the RMSF of the C-alpha atoms which were used to color the residues in <xr id="fig:rmsf"/>. The fluctuation of all atoms per residue is depicted in <xr id="fig:rmsf_p"/>. The fluctuation patterns of <xr id="fig:rmsf_p"/> coincide with <xr id="fig:rmsf_c"/> but the amplitudes are more pronounced since side chains can rotate without changing the position of the C-alpha atom.

In all three simulations, the region 15-75 is fluctuating most which covers flexible loop regions of the hydrolase domain shown in the top left corner of <xr id="fig:rmsf"/>. Furthermore, the loop 341-350 at the center is moving dynamically, in particular residue F347 which belongs to the residues with the highest RMSF in all three simulations. The flexibility of loop 341-350 might be in fact functionally relevant since it is at the entrance of the TIM-beta/alpha-barrel domain which includes the active site residues E235 and E340. In W209R, the adjacent loop 393-396 is also fluctuating conspicuously. 
 
<figure id="fig:rmsf">
{|style="border-collapse: separate; border-style: solid; border-spacing: 0; border-width: 1px"
|<subfigure>[[File:rsmf_wt.gif|thumb|300px|<caption>Wildtype. Most flexible residues: 169, 347, 441.</caption>]]</subfigure>
|<subfigure>[[File:rsmf_L470P.gif|thumb|300px|<caption>L470P. Most flexible residues: 28, 30, 49, 197, 347, 313.</caption>]]</subfigure>
|<subfigure>[[File:rsmf_W209R.gif|thumb|300px|<caption>W209R. Most flexible residues: 61, 63, 347, 348, 395, 397.</caption>]]</subfigure>
|}
<small><caption>MD simulations with residues colored by the RMSF per residue and C-alpha atom. Red: high flexibility; Green: intermediate flexibility; Blue: low flexibility; Purple: mutant residue; Sticks: residues with the highest RMSF.</caption></small>
</figure>

<figure id="fig:rmsf_c">
{|style="border-collapse: separate; border-style: solid; border-spacing: 0; border-width: 1px"
|<subfigure>[[File:rsmf_c_wt.png|thumb|300px|<caption>Wildtype</caption>]]</subfigure>
|<subfigure>[[File:rsmf_c_L470P.png|thumb|300px|<caption>L470P</caption>]]</subfigure>
|<subfigure>[[File:rsmf_c_W209R.png|thumb|300px|<caption>W209R</caption>]]</subfigure>
|}
<small><caption>RMSF per residue with respect to the C-alpha atom.</caption></small>
</figure>

<figure id="fig:rmsf_p">
{|style="border-collapse: separate; border-style: solid; border-spacing: 0; border-width: 1px"
|<subfigure>[[File:rsmf_p_wt.png|thumb|300px|<caption>Wildtype</caption>]]</subfigure>
|<subfigure>[[File:rsmf_p_L470P.png|thumb|300px|<caption>L470P</caption>]]</subfigure>
|<subfigure>[[File:rsmf_p_W209R.png|thumb|300px|<caption>W209R</caption>]]</subfigure>
|}
<small><caption>RMSF per residue with respect to all atoms.</caption></small>
</figure>

=== RMSF L470P vs. RMSF wildtype ===
For testing to which extent the motion of certain parts of L470P differ from the wildtype, we substracted the RMSF values of the wildtype from the RMSF values of L470P, i.e deltaRMSF = RMSFmut-RMSFwt. The difference of <xr id="fig:rmsf_c#2"/> and <xr id="fig:rmsf_c#1"/> is shown in <xr id="fig:rmsf_L470P_diff#1"/>: positive and negative values indicate regions in L470P which are more flexible and less flexible compared to the wildtype, respectively. 

We further applied the two-sided Student's t-test for testing whether the mean of <xr id="fig:rmsf_L470P_diff#1"/> significantly deviates from zero, i.e. whether there is a significant difference in the fluctuation between the wildtype and L470P. For this, we first computed the test-statistic t=sqrt(n)*avg(deltaRMSF)/sd(deltaRMSF), where sqrt is the square root, avg the average, and sd the sampling standard deviation. Next, we computed the p-value P(|T| > t) using the t-distribution with n-1 degrees of freedom. And in deed, L470P is significantly (p-value 5e-13) more flexible compared to the wildtype.

In order to identify the regions which contribute most to the change in flexibility, we colored residues of L470P according to deltaRMSF. <xr id="fig:rmsf_L470P_diff#2"/> reveals that there are two loop regions (23-33, 42-53) underneath the hydrolase domain which are moving conspicuously. <xr id="fig:rmsf_L470P_diff#2"/> depicts that both loop regions are exposed and might be interacting with other constituents of the lysosomal membrane. These interactions might be disturbed by the increased flexibility due to L470P. 

<figure id="fig:rmsf_L470P_diff">
{|style="border-collapse: separate; border-style: solid; border-spacing: 0; border-width: 1px"
|-align="top"
|<subfigure>[[File:rmsf_L470P_diff.png|thumb|300px|<caption>Difference RMSF mutant - RMSF wildtype. The most flexible regions (23-33) and (42-53) are bordered by the red and orange vertical lines, respectively.</caption>]]</subfigure>
|<subfigure>[[File:rmsf_L470P_diff.gif|thumb|300px|<caption>MD simulation of L470P. The red loop regions are moving most.</caption>]]</subfigure>
|<subfigure>[[File:rmsf_L470P_diff_surface.png|thumb|300px|<caption>Surface of L470P. The flexible loops are on the surface of the protein.</caption>]]</subfigure>
|}
<small><caption>RMSF L470P compared to RMSF wildtype. The RMSF difference of Figure a were used for coloring residues. Red: higher flexibility; Green: same flexibility; Blue: lower flexibility; Purple: L470. Loop region 23-33 and 42-53 is significantly more flexible in L470P than in the wildtype.</caption></small>
</figure>

=== RMSF W209R vs. RMSF wildtype ===
We investigated the difference of W209R and the wildtype in the same way as described above. <xr id="fig:rmsf_W209R_diff#1"/> illustrates that there are single residues whose mobility changed in W209R relative to the wildtype. Specifically, E395 (<xr id="fig:rmsf#3"/>) is more mobile, and T180 is less mobile. However, the overall fluctuation does not significantly (p-value 0.34) differ from the wildtype. We found that the positively charged mutant residue R209 forms a hydrogen bond to T180 (<xr id="fig:rmsf_W209R_diff#2"/>, <xr id="fig:rmsf_W209R_diff#3"/>) and therefore reduces the flexibility of T180. This might be one reason for the malign effect of W209R.

<figure id="fig:rmsf_W209R_diff">
{|style="border-collapse: separate; border-style: solid; border-spacing: 0; border-width: 1px"
|-align="top"
|<subfigure>[[File:rmsf_W209R_diff.png|thumb|300px|<caption>Difference RMSF mutant - RMSF wildtype. T180 and S181 are less flexible due to interactions with R209.</caption>]]</subfigure>
|<subfigure>[[File:rmsf_W209R_diff.gif|thumb|300px|<caption>MD simulation of W209R. The blue loop adjacent to W209 (purple) is bordered red in Figure a.</caption>]]</subfigure>
|<subfigure>[[File:rmsf_W209R_diff_stab.png|thumb|300px|<caption>W209 forms a hydrogen bond to T180 which becomes therefore less flexible.</caption>]]</subfigure>
|}
<small><caption>RMSF W209R compared to RMSF wildtype. The RMSF difference of Figure a were used for coloring residues. Red: higher flexibility; Green: same flexibility; Blue: lower flexibility; Purple: W209. The flexibility does not significantly differ from the wildtype.</caption></small>
</figure>