Molecular Dynamics Simulations Analysis Gaucher Disease
Technical details can be found in our protocol.
A brief check of results
We briefly checked the length of the MD simulations via gmxcheck. Since we used the AGroS paramter -lengthMD 10, we expected the MD simulation to be 10 ns long. gmxcheck 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">
| Name | Runtime | ns/day | years/1s | |
|---|---|---|---|---|
| WT | 16h25:59 | 14.605 | 187536 | |
| L470P | 9h20:49 | 25.677 | 106699 | |
| W209R | 9h18:47 | 25.770 | 106314 | |
Runtime of the MD simulations and the estimated runtime for a simulation length of 1s. </figtable> The computation of the energy term Coul(T) + LJ [W3-W3] 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">
</subfigure> </subfigure> </subfigure>| <subfigure> |
| <subfigure> |
| <subfigure> |
Visualization of the MD simulations. Red: the active site residues E235 and E340; Blue: the mutated residue. </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">
</subfigure> </subfigure> </subfigure>| <subfigure> |
| <subfigure> |
| <subfigure> |
Development of the pressure during the MD simulation. In brackets: average/total-drift. </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">
</subfigure> </subfigure> </subfigure>| <subfigure> |
| <subfigure> |
| <subfigure> |
Development of the temperature during the MD simulation. In brackets: average/total-drift. </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">
</subfigure> </subfigure> </subfigure>| <subfigure> |
| <subfigure> |
| <subfigure> |
Development of the potential energy during the MD simulation. In brackets: average/total-drift. </figure>
<figure id="fig:qual_tot">
</subfigure> </subfigure> </subfigure>| <subfigure> |
| <subfigure> |
| <subfigure> |
Development of the total energy during the MD simulation. In brackets: average/total-drift. </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">
</subfigure> </subfigure> </subfigure>| <subfigure> |
| <subfigure> |
| <subfigure> |
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. </figure>
RMSF
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.
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:rsmf_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">
</subfigure> </subfigure> </subfigure>| <subfigure> |
| <subfigure> |
| <subfigure> |
Superposition of the average g_rmsf (green) structure onto the initial structure (blue). </figure>
<xr id="fig:rmsf_c"/> shows the RMSF of the C-alpha atoms. In all three simulations, the backbone region 15-75 is fluctuating most. This region covers the loops of the hydrolase domain shown in the top left corner of <xr id="fig:viz"/>. The fluctuation of all atoms per residue is depicted in Figure <xr id="fig:rmsf_p"/>. The fluctuation patterns 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.
<figure id="fig:rmsf_c">
</subfigure> </subfigure> </subfigure>| <subfigure> |
| <subfigure> |
| <subfigure> |
RMSF per residue with respect to the C-alpha atom. </figure>
<figure id="fig:rmsf_p">
</subfigure> </subfigure> </subfigure>| <subfigure> |
| <subfigure> |
| <subfigure> |
RMSF per residue with respect to all atoms. </figure>
The B-factor of an atom describes the
<figure id="fig:rmsf">
</subfigure> </subfigure> </subfigure>| <subfigure> |
| <subfigure> |
| <subfigure> |
MD simulations with residues colored by the g_rsmf B-factor (rainbow spectrum). Sticks: residues with a B-factor > 120. Purple: the mutated residue. </figure>
<figure id="fig:rmsf_L470P_diff">
</subfigure> </subfigure> </subfigure>| <subfigure> |
| <subfigure> |
| <subfigure> |
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. </figure>
<figure id="fig:rmsf_W209R_diff">
| <subfigure> |
| <subfigure> |
| <subfigure> |
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. </figure>
