Difference between revisions of "Canavan Task 10 - Molecular Dynamics Simulations"

From Bioinformatikpedia
(Visualization)
(Visualization)
Line 448: Line 448:
 
==Visualization==
 
==Visualization==
   
[[File:movie.gif]]
 
   
[... wt_movie.avi]
+
[http://www.cip.ifi.lmu.de/~vorberg/wt_movie.avi wt_movie.avi]
   
 
==Radius of gyration ==
 
==Radius of gyration ==

Revision as of 15:42, 30 July 2012

Protocol

Further information and commands can be found in the protocol.

Initial Checks

For all three runs, there are 2000 time frames, with a resolution of 5 psec. Therefore the whole simulation ran for 10000 psec = 10 nsec. No errors occured during the run and the simulation finished properly.

<figtable id="gmxcheck">

<xr nolink id="gmxcheck"/> Results from gmxcheck and the logfile
WT K213E A305E
run time 5h 22:37 5h 08:54 5h 08:35
atoms outside of box 406, 408, 480, 482, 483, 484, 485, 486, 500, 501,.. 406, 480, 482, 483, 484, 485, 486, 500, 502, 503,.. 406, 408, 480, 482, 483, 484, 485, 486, 500, 501,..
Last frame 2000 time 10000.000 2000 time 10000.000 2000 time 10000.000

</figtable>

Energies

For all analysed thermodynamical parameters convergence could be observed. Though, for pressure the values vary enormously, but the average pressure is close to the specified value of 1 bar.

Temperature

The temperature fluctuates about the reference temperature of 298K. There are only small deviations from this value, which is reflected in the small error estimates.

<figtable id="energies_temp">

<xr nolink id="energies_temp"/> Quality check for convergence of energy values: temperature.
WT K213E A305E
Reference Value 298 K 298 K 298 K
Average 297.914 K 297.908 297.917
Err.Est. 0.0072 0.0059 0.0047
RMSD 1.38693 1.38991 1.3981
Total Drift -0.00042403 (K) 0.00867898 (K) 0.0154339 (K)
Plot <figure id="wt_temp">
<xr nolink id="wt_temp"/>
Temperature of the system for the wildtype smoothed over a window of 20.
</figure>
<figure id="k213_scwrl_temp">
<xr nolink id="k213_scwrl_temp"/>
Temperature of the system for K213E smoothed over a window of 20.
</figure>
<figure id="a305_scwrl_temp.png">
<xr nolink id="a305_scwrl_temp.png"/>
Temperature of the system for A305E smoothed over a window of 20.
</figure>

</figtable>

Pressure

The pressure varies enormously about the reference pressure of 1 bar. This results in large RMSD values. Yet, the average pressure over the whole simulation is very close to the specified pressure of 1 bar and therefore the error estimates are relatively small.


<figtable id="energies_temp">

<xr nolink id="energies_temp"/> Quality check for convergence of energy values: pressure.
WT K213E A305E
Reference Value 1.0 (Berendsen barostat) 1.0 (Berendsen barostat) 1.0 (Berendsen barostat)
Average 1.00763 1.00683 1.00797
Err.Est. 0.018 0.015 0.022
RMSD 111.943 112.761 112.987
Total Drift -0.0713928 (bar) -0.0585283 (bar) -0.100316 (bar)
Plot <figure id="wt_pressure">
<xr nolink id="wt_pressure"/>
Pressure of the system for the wildtype smoothed over a window of 20.
</figure>
<figure id="k213_scwrl_pressure">
<xr nolink id="k213_scwrl_pressure"/>
Pressure of the system for K213E smoothed over a window of 20.
</figure>
<figure id="a305_scwrl_pressure">
<xr nolink id="a305_scwrl_pressure"/>
Pressure of the system for A305E smoothed over a window of 20.
</figure>

</figtable>

Potential Energy

The potential energy decreases over simulation time, which leads to the conclusion, that the conformation of the structure is energetically optimized during the simulation. The decrease in energy is reflected in the large negative drift values.

<figtable id="energies_potential">

<xr nolink id="energies_potential"/> Quality check for convergence of energy values: potential energy.
WT K213E A305E
Average -592161 -585993 -583187
Err.Est. 55 65 46
RMSD 721.947 726.81 722.076
Total Drift -252.112 (kJ/mol) -425.823 (kJ/mol) -290.165 (kJ/mol)
Plot <figure id="CD_wt_potenergy">
<xr nolink id="CD_wt_potenergy"/>
Potential energy of the system for the wildtype smoothed over a window of 20
</figure>
<figure id="k213_scwrl_poten">
<xr nolink id="k213_scwrl_poten"/>
Potential energy of the system for K213E smoothed over a window of 20
</figure>
<figure id="a305_scwrl_potenergy">
<xr nolink id="a305_scwrl_potenergy"/>
Potential energy of the system for A305E smoothed over a window of 20
</figure>

</figtable>

Total energy

The total energy values resemble the potential energy values. In contrast to the potential energy, kinetic energies are included in the calculation of the total energy. Again the total energy decreases during the simulation.

<figtable id="energies_total">

<xr nolink id="energies_total"/> Quality check for convergence of energy values: total energy.
WT K213E A305E
Average -485680 -480905 -478506
Err.Est. 54 64 45
RMSD 886.726 888.522 887.724
Total Drift -252.262 (kJ/mol) -422.763 (kJ/mol) -284.743 (kJ/mol)
Plot <figure id="wt_tot_energy">
<xr nolink id="wt_tot_energy"/>
Total energy of the system for the wildtype, smoothed over a window of 20 timesteps.
</figure>
<figure id="k213_scwrl_toten">
<xr nolink id="k213_scwrl_toten"/>
Total energy of the system for K213E, smoothed over a window of 20 timesteps.
</figure>
<figure id="a305e_scwrl_totenergy">
<xr nolink id="a305e_scwrl_totenergy"/>
Total energy of the system for A305E smoothed over a window of 20
</figure>

</figtable>

distances between periodic boundaries

We calculated the minimum distance between periodic images for the whole protein (not only C-alpha atoms). The suggested distance limit of 2nm is undercut at some timesteps during the simulation. Especially for the mutant K213E the distance often is below 2nm. This might have caused undesired unphysical interactions.

<figtable id="periodic_boundary">

<xr nolink id="periodic_boundary"/> Minimum distances between periodic images for the simulation run calculated with g_mindist.
WT K213E A305E
shortest dist 1.6456 (nm) 1.47431 (nm) 1.59859 (nm)
at time step 7675 (ps) 7515 (ps) 2460 (ps)
between atoms 15 and 4507 598 and 4497 597 and 4507
Plot <figure id="wt_pi">
<xr nolink id="wt_pi"/>
Minimum distance of periodic images for the wildtype protein. Between timesteps 6000-9000 the minimal distance between periodic images is below 2nm.
</figure>
<figure id="k213e_pi">
<xr nolink id="k213e_pi"/>
Minimum distance of periodic images for K213E. Almost during the whole simulation the minimal distance is less than 2 nm.
</figure>
<figure id="a305e_pi">
<xr nolink id="a305e_pi"/>
Minimum distance of periodic images for A305E. The minimum distance between periodic images is larger than 2nm for most of the simulation.
</figure>

</figtable>

RMSF

Only small fluctuations can be observed for the residues of the three proteins.

For all three proteins there is a peak around residues 60-75, which defines this region as rather flexible. Especially for K213E, there is a strong peak of more than 0.35 nm. In the b-factor figures, this region is highlighted with a red circle in the lower right corner.

There is another flexible region formed by residues 220-230. This region is also highlighted by a red circle in the b-factor figures in the lower left corner.

For the wildtype, there is a region between residues 120 and 180 that is especially rigid. When looking at the structure, one finds that this regions defines the core of the protein.

The figures of the b-factors of the proteins visualize the fluctuation plots. Most parts of the protein are rather rigid and only some exposed loops have higher bfactors. These loopWhen looking at the bfactors one finds the same results. s correspong to the regions identified in the fluctuation plots.

Interestingly, K213E has rather low b-factor values, than one might expect from the fluctuations plot, whereas A305E has higher b-factors than the wildtype and seems to be more motile.


<figtable id="rmsf">

<xr nolink id="rmsf"/> Root Mean Square Fluctuations for the MD simulation run calculated per residue with g_rmsf.
WT K213E A305E
RMSF Plot <figure id="wt_rmsf_plot">
<xr nolink id="wt_rmsf_plot"/>
Per residue root-mean-square-fluctution plot for the wildtype protein. There is a rigid region around residues 120-180.
</figure>
<figure id="k213e_rmsf_plot">
<xr nolink id="k213e_rmsf_plot"/>
Per residue root-mean-square-fluctution plot for K213E.
</figure>
<figure id="a3o5e_rmsf_plot">
<xr nolink id="a3o5e_rmsf_plot"/>
Per residue root-mean-square-fluctution plot for A305E.
</figure>
B-Factors <figure id="wt_bfactors">
<xr nolink id="wt_bfactors"/>
The wildtype protein colored according to the calculated B-Factors. Highly flexible loops (according to residues from the fluctuations plot) are emphasized by red circles.
</figure>
<figure id="k213e_bfactors">
<xr nolink id="k213e_bfactors"/>
K213E colored according to the calculated B-Factors. Highly flexible loops (according to residues from the fluctuations plot) are emphasized by red circles.
</figure>
<figure id="a305e_bfactors">
<xr nolink id="a305e_bfactors"/>
A305E colored according to the calculated B-Factors. Highly flexible loops (according to residues from the fluctuations plot) are emphasized by red circles.
</figure>

</figtable>



g_rmsf also generates an unphysical average structure. One can see, that for regions with high b-facotrs, the averaged structure shows several possible residue conformers.

<figtable id="average">

<figure id="wt_average">
<xr nolink id="wt_average"/>
Close-up of the loop region formed by residues 60-75 and the average residue conformers for this loop. The higher the b-factor, the more possible residue conformers are found.
</figure>

</figtable>

Convergence of RMSD

As expected, the RMSD increases when using the starting structure as a reference: Over the simulation the structure changes and deviates more and more from the starting structure. Yet these changes are not tremendous (RMSD < 0.2 nm), as the starting structure is derived from the crystal structure and therefore should already have adopted a optimal conformation.

When using only C-alpha atoms for the RMSD calculation instead of the whole protein, the values are smaller. This is because sidechains are the most flexible elements in a structure and cause the higher RMSD compared to only C-alpha atoms.

When taking the average structure as reference, the RMSD is higher at the beginning of the simulation and finally converges as the structure reaches an equilibrium. Only for A305E the RMSD increases in the last 2000 timesteps and therefore does not show convergence! Another interesting point is, that for both mutants the RMSD values are much higher than for the wildtype: For the wildtype the RMSD towards the average structure is about 0.48 nm at the end of the simulation, whereas for K213E it is 0.69 nm and for A305E it is 0.61 nm.


<figtable id="rmsd">

<xr nolink id="rmsd"/> Root Mean Square Deviations of the structure over time towards different reference structures calculated with g_rmsd.
WT K213E A305E
RMSD whole protein towards starting structure <figure id="wt_rmsd_all_vs_first">
<xr nolink id="wt_rmsd_all_vs_first"/>
RMSD of wildtype protein towards the starting structure over the simulation time. RMSD increases fast and converges at the end of the simulation.
</figure>
<figure id="k213e_rmsd_all_vs_first">
<xr nolink id="k213e_rmsd_all_vs_first"/>
RMSD of K213E mutant towards the starting structure over the simulation time. RMSD increases fast and converges at the end of the simulation.
</figure>
<figure id="a3o5e_rmsd_all_vs_first">
<xr nolink id="a3o5e_rmsd_all_vs_first"/>
RMSD of A305E mutant towards the starting structure over the simulation time. RMSD increases fast and converges at the end of the simulation.
</figure>
RMSD C-alpha atoms towards starting structure <figure id="wt_rmsd_calpha_vs_first">
<xr nolink id="wt_rmsd_calpha_vs_first"/>
C-alpha RMSD of the wildtype protein towards the starting structure over the simulation time. At the beginning of the simulation, the RMSD increases a bit and converges towards the end.
</figure>
<figure id="k213_scwrl_rmsd_calpha_first">
<xr nolink id="k213_scwrl_rmsd_calpha_first"/>
C-alpha RMSD of K213E mutant towards the starting structure over the simulation time. At the beginning of the simulation, the RMSD increases a bit and converges towards the end.
</figure>
<figure id="a305e_scwrl_rmsd_calpha_first">
<xr nolink id="a305e_scwrl_rmsd_calpha_first"/>
C-alpha RMSD of A305E mutant towards the starting structure over the simulation time. At the beginning of the simulation, the RMSD increases a bit and converges towards the end.
</figure>
RMSD whole protein towards average structure <figure id="wt_rmsd_all_vs_average">
<xr nolink id="wt_rmsd_all_vs_average"/>
RMSD of wildtype protein towards the average structure over the simulation time. RMSD decreases as the protein structure converges towards an average structure.
</figure>
<figure id="k213e_rmsd_all_vs_average">
<xr nolink id="k213e_rmsd_all_vs_average"/>
RMSD of K213E mutant towards the average structure over the simulation time. RMSD decreases as the protein structure converges towards an average structure.
</figure>
<figure id="a305e_rmsd_all_vs_average">
<xr nolink id="a305e_rmsd_all_vs_average"/>
RMSD of A305E mutant towards the average structure over the simulation time. RMSD decreases at the beginning but increases for the last 2000 time steps. A305E does not show structural convergence.
</figure>
RMSD protein backbone towards average structure <figure id="wt_rmsd_calpha_vs_average">
<xr nolink id="wt_rmsd_calpha_vs_average"/>
RMSD of the wildtype backbone towards the average structure over the simulation time. RMSD decreases as the protein structure converges towards an average structure.
</figure>
<figure id="k213_scwrl_rmsd_calpha_average">
<xr nolink id="k213_scwrl_rmsd_calpha_average"/>
RMSD of the K213E mutant backbone towards the average structure over the simulation time. RMSD decreases as the protein structure converges towards an average structure.
</figure>
<figure id="a305e_scwrl_rmsd_calpha_average">
<xr nolink id="a305e_scwrl_rmsd_calpha_average"/>
RMSD of the A305E mutant backbone towards the average structure over the simulation time. RMSD varies a lot over simulation time, but in the end has the same value as at the beginning of the simulation.
</figure>

</figtable>

Visualization

wt_movie.avi

Radius of gyration

Against our expectations, the radius of gyration increases for the proteins. As the energy of the system decreases during the run, we would expect that the protein becomes more compact. One idea is, that we used the monomeric form of the protein for the simulation, whereas in the crystal structure it is a dimer. Therefore the monomer might have a different energetical optimal conformation than compared to its dimer bound form.

Yet, the changes are only minor and within a range of less than 0.05 nanometer. Thus, the increase in the radius of gyration is not of great impact.


<figtable id="rog">

<xr nolink id="rog"/> Radius of gyration of the proteins for the simulation run calculated with g_gyrate.
WT K213E A305E
Plot <figure id="wt_rg">
<xr nolink id="wt_rg"/>
Radius of gyration of the wildtype protein during the MD simulation.
</figure>
<figure id="k213_scwrl_rog">
<xr nolink id="k213_scwrl_rog"/>
Radius of gyration of the K213E mutant during the MD simulation.
</figure>
<figure id="a305e_scwrl_rog">
<xr nolink id="a305e_scwrl_rog"/>
Radius of gyration of the A305E mutant during the MD simulation.
</figure>

</figtable>