What if molecular dynamics

For a multi-chain protein, the chains may not be held together when using the standard fitting procedure. A robust manner is a two-step fitting:. When prompted to select groups, select the backbone for the fitting and the protein for the output.

After all frames have loaded, you can show the frames by clicking the forward button in the main window. You will see the dynamics of the peptide. The default view is a stick representation. You can change the representation or actually view multiple representations through opening the menu "Graphics - Representation". To get a view of the conformational space visited, all frames can be shown at the same time. We want to show all of the atoms again. Now, click on this representation in the menu.

To show all frames at the same time, click on "Trajectory" tab in the Graphical Representations menu. Instead of the word "now near the bottom , type , which will show the first 11 frames, so if you want to see all of them, figure out how many frames there are For a good view, you may need to decrease the thickness of the bonds a bit - you can do that from the "Draw style" tab.

After a first visual inspection of the trajectory, it's time to perform some more thorough checks regarding quality of the simulations. This quality assurance QA involves tests for the convergence of thermodynamic parameters, such as the temperature, the pressure, the potential and the kinetic energy. More commonly phrased, QA tries to assess whether an equilibrium was reached.

Convergence is also checked in terms of the structure, through the root mean square deviation RMSD against the starting structure and the average structure. Next to that, it has to be checked that there have not been interactions between adjacent periodic images, as such interactions could lead to unphysical effects. A final QA test involves the calculation of the root mean square fluctuations of atoms, which can be compared to crystallographic b-factors.

We start off with the extraction of some thermodynamic parameters from the energy file. The following properties will be investigated: temperature, pressure, potential energy, kinetic energy, unit cell volume, density, and the box dimensions. Most of these have been checked already for some of the preparation steps. Energy analysis is performed with the tool gmx energy. This program reads in the energy file, which is produced during the simulation and has the extension.

Give the following command to run gmx energy: gmx energy -f ener. This will give a list of the energy and related terms which are stored in the. In our case, there are probably more than 60 terms in the energy file, each of which may be extracted and graphed.

The first set of entries correspond to different energy terms in the force field. Also note the group of entries which list the energy terms split over the groups Protein and Non-Protein, including the interactions between the two. To extract the temperature, you can type "Temperature" and press Enter twice, or enter the number for the temperature from the list, then enter 0 and press Enter.

Have a look at the graph with the program xmgrace and see how the temperature fluctuates around the value specified K. Have a look also at the kinetic and potential energies and see how the values converge. Also, check the volume. If the values have not converged, this indicates that the simulation has not yet reached thermal equilibrium and that it should be extended before doing further analysis.

Moreover, the period over which equilibration takes place should not be included in the analysis. Here, for the sake of simplicity, we will neglect these considerations and use the results from the simulations as they are. The equilibration of some terms takes longer than that of other terms. In particular, the temperature readily converges to its equilibrium value, whereas the interaction between different parts of the system may take much longer.

In general, you would need to investigate the properties of interest for their convergence. One of the most important things to check in terms of quality assurance is that there have not been direct interactions between periodic images. Since periodic images are identical, such interactions are unphysical self-interactions. Imagine that a peptide with a dipole would have direct interactions.

Then the attraction between the two ends of the same molecule over the periodic boundary would affect and invalidate the results relating to the native behaviour of the peptide.

To verify that such interactions have not taken place, we calculate the minimal distance between periodic images at each time. This is done using mindist. Not only direct interactions are of concern, but indirect effects, mediated through the water can also pose a problem.

For instance, the peptide can cause water ordening in the first four shells of water. That corresponds to about one nanometer. Ideally, the minimal distance should therefore not be less than two nanometer. Next to inspection of energies and such, convergence of the simulation towards equilibrium should also be inspected from the relaxation of the structure.

Usually, such relaxation is only considered in terms of the Euclidean distance from the structure to a reference structure, e. This distance is termed the root mean square deviation RMSD.

However, it is recommended also to investigate the relaxation towards the average structure, i. The reasons for this will be set out in the next paragraph. But to calculate the RMSD against the average structure, requires first to obtain the average.

This structure can be obtained as a side product from calculation of the root mean square fluctuations RMSF. The RMSF captures, for each atom, the fluctuation about its average position. This gives insight into the flexibility of regions of the peptide. The RMSF and the average structure are calculated with gmx rmsf. We are most interested in the fluctuations on a per residue basis, which is controlled by the flag -res.

Have a look at the graph of the RMSF with xmgrace and identify the flexible and rigid regions. Can you identify these regions also by looking at the superimposed structures? The RMSD is commonly used as an indicator of convergence of the structure towards an equilibrium state. As mentioned above, the RMSD is merely a distance measure and is most meaningful for low values. The RMSD is calculated using the program gmx rms.

First calculate the RMSD for all peptide atoms, using the starting structure as a reference:. This time the RMSD settles at a lower value, which is the result of excluding the, often flexible, side chain atoms. In both cases the RMSD increases to a plateau value.

This means that the structure of the peptide reaches a certain distance from the reference structure and then keeps that distance more or less. However, with increasing distance, the amount of conformations available also increases. This means that, although the RMSD has reached a plateau value, the structure may still be progressing towards its equilibrium state.

For this reason, it is advisable also to check the convergence towards the average structure. Compare the resulting graphs to the previous ones. Note at which point the RMSD values level off. Which is a better measure for convergence? As a final part of the QA, we calculate the radius of gyration. Do a ls -lrt to see which files were written by mdrun. For this we will use VMD as a trajectory viewer. Other possibilities would be Pymol , or gOpenmol.

As standard the view on the system is shown as perspective. To change the speed of the movie you find a controller next to the arrow symbol. We have simulated Argon at K, which is above the boiling point, so we have simulated argon gas. After the NPT steps, a production simulation at conditions of constant volume NVT was performed in three replicates for each of the models and the crystal structure ns.

The equilibration and production steps used a time step of 2 and 4 fs, respectively. Selection of these time steps was possible due to the hydrogen mass repartitioning scheme being employed in ACEMD [ ]. A smooth switching function for the cut-off was applied, starting at 7. The size of the cell was set to prevent non-bonded interactions between the protein and its periodic boundary image.

The D 3 R crystal structure was used to perform redocking calculations of eticlopride. A conformational search of eticlopride was performed using the obconformer tool the Open Babel Package, version 2. The best conformer was further refined through conjugate gradient and steepest descent algorithms by using the obminimize tool the Open Babel Package, version 2.

AutoDock [ ] release 4. Docking input files were prepared through the AutoDockTools ADT package[ ] using a grid of 40x80x points with a spacing of 0. AutoGrid4 was used to generate grid maps. The Lamarckian genetic algorithm LGA was employed with a population size of individuals, a maximum number of 2,, energy evaluations, a maximum number of 27, generations, and runs. Ten diverse poses of eticlopride were selected based on visual inspection.

Three replicates of ns were generated for each complex. The five largest clusters from each model were selected for further analysis. All structures were first aligned to the initial crystal structure prepared using the OPLS protocol. Prior to the calculation, the membrane, solvent, salt, hydrogens, ligand, and terminus caps were removed.

The MolProbity scores were calculated using the MolProbity[ 28 ] module in Phenix[ ] with default settings. Each of the initial models and their corresponding five cluster centroids were aligned to the OPLS prepared crystal structure and prepared for docking with DOCK v3. Hydrogen positions from the MD preparations and simulations were maintained. Docking grids were generated using the program blastermaster with default parameters.

The binding energy was calculated as the sum of electrostatic and van der Waals interaction terms, corrected for ligand desolvation. These energy terms were calculated from precalculated grids with a default box size. The electrostatic grid was calculated using Qnifft [ ]. This library was docked with DOCK 3. The ligands were matched with the adaptive sampling routine with an initial distance tolerance of 0.

The best scoring conformation was minimized for a maximum of steps or until converged to 0. Table A. Protocols used in simulations of the D 3 R models and crystal structure. Table B.

Percentage of improved models after MD refinement based on analysis of side chains in the binding site. Table C. Table D. Percentage of improved models after MD refinement depending on the number of included cluster centroids 1—5. Table E. Table F. Fig A. The model M 5 was treated as an outlier and is not included in the comparison.

Fig B. Fluctuations in the TM region of D 3 R crystal structure and models. Fig C. Fig D. Effect of MD refinement on protein quality scores. Fig E. Fig F. Ligand poses from the simulations initiated from the crystal structure of D 3 R.

The receptor is shown as cartoons and the ligand in sticks. The best MD refined models and the crystal structure are colored green and grey, respectively. Fig G. Fig H. Fig I. Virtual screening performance of the initial and MD refined structures. EF1 results for the crystal structure and initial models are shown as bars.

The authors thank Ennys Gheyouche for his contributions to this work. Abstract The determination of G protein-coupled receptor GPCR structures at atomic resolution has improved understanding of cellular signaling and will accelerate the development of new drug candidates. Author summary The determination of detailed molecular structures of proteins is important to understand how they function and is an essential component of many drug discovery campaigns.

Introduction Three-dimensional structures of proteins have contributed to detailed knowledge about biological processes at the molecular level.

Results Benchmarking set, simulation protocols, and assessment criteria A set of 30 models of the D 3 R in complex with the antagonist eticlopride was obtained from the GPCR Dock assessment [ 14 ]. Download: PPT. Fig 1. Structure of the D 3 R in complex with eticlopride.

Fig 2. Table 2. Fig 3. Fig 4. What is the influence of the quality of the receptor model on refinement of the ligand pose? Fig 7. Fig 8. MD refinement of D 3 R crystal structure with binding modes of eticlopride generated by molecular docking. Restraints on the TM region improve loop and ligand refinement The TM region tended to drift away from the reference during the unrestrained MD simulations.

Fig 9. Comparison of MD refinement with restrained and unrestrained TM region. Fig Discussion The goal of this study was to assess if models of GPCRs could be refined by MD simulations to improve predictions of receptor-ligand complexes. Docking of eticlopride to the D 3 R crystal structure The D 3 R crystal structure was used to perform redocking calculations of eticlopride. Analysis of MD simulations Clustering of the trajectories.

Molecular docking and ligand enrichment. Supporting information. S1 Appendix. Supporting information for MD simulations and molecular docking calculations. Acknowledgments The authors thank Ennys Gheyouche for his contributions to this work. References 1. The Protein Data Bank. Nucleic Acids Res. UniProt C. UniProt: a worldwide hub of protein knowledge. Kuhlman B, Bradley P.

Advances in protein structure prediction and design. Nat Rev Mol Cell Biol. Assessment of the protein-structure refinement category in CASP8. J Chem Theory Comput. Terashi G, Kihara D. Protein structure model refinement in CASP12 using short and long molecular dynamics simulations in implicit solvent.

Protein structure modeling and refinement by global optimization in CASP Assessment of the model refinement category in CASP A comprehensive map of molecular drug targets. Nat Rev Drug Discov. Trends in GPCR drug discovery: new agents, targets and indications. J Phys Chem B. View Article Google Scholar Nat Methods. J Comput Chem. Four-dimensional docking: a fast and accurate account of discrete receptor flexibility in ligand docking.

J Med Chem. Simulations of a protein crystal with a high resolution X-ray structure: evaluation of force fields and water models. Dynameomics: a comprehensive database of protein dynamics. Concerted interconversion between ionic lock substates of the beta 2 adrenergic receptor revealed by microsecond timescale molecular dynamics. Biophys J. Structure and function of an irreversible agonist-beta 2 adrenoceptor complex. MolProbity: all-atom structure validation for macromolecular crystallography.

Acta Crystallogr D Biol Crystallogr. Shen MY, Sali A. Statistical potential for assessment and prediction of protein structures. Protein Sci. How well can the accuracy of comparative protein structure models be predicted? Directory of useful decoys, enhanced DUD-E : better ligands and decoys for better benchmarking. Ligand pose and orientational sampling in molecular docking. PLoS One. Rapid context-dependent ligand desolvation in molecular docking.

J Chem Inf Model. Systematic validation of protein force fields against experimental data. How fast-folding proteins fold.

Refinement of protein structure homology models via long, all-atom molecular dynamics simulations. Evaluation of model refinement in CASP Structure refinement of membrane proteins via molecular dynamics simulations. Feig M, Mirjalili V. Protein structure refinement via molecular-dynamics simulations: What works and what does not?

Physics-based protein structure refinement through multiple molecular dynamics trajectories and structure averaging. Heo L, Feig M. What makes it difficult to refine protein models further via molecular dynamics simulations? Childers MC, Daggett V. The reliability of molecular dynamics simulations of the multidrug transporter P-glycoprotein in a membrane environment.

Int J Mol Sci. Assessment of the utility of contact-based restraints in accelerating the prediction of protein structure using molecular dynamics simulations. Zhou H, Zhou Y. Distance-scaled, finite ideal-gas reference state improves structure-derived potentials of mean force for structure selection and stability prediction. A composite score for predicting errors in protein structure models. Performance of virtual screening against GPCR homology models: Impact of template selection and treatment of binding site plasticity.

PLoS Comput Biol. Structure-based ligand discovery for the protein-protein interface of chemokine receptor CXCR4. Strategies for improved modeling of GPCR-drug complexes: blind predictions of serotonin receptors bound to ergotamine. Computational prediction of alanine scanning and ligand binding energetics in G-protein coupled receptors.

Pathway and mechanism of drug binding to G-protein-coupled receptors. Advances in free-energy-based simulations of protein folding and ligand binding.