Explicit Solvent Molecular Dynamics Research Articles

Significant Refinement of Protein Structure Models Using a Residue-Specific Force Field.

An important application of all-atom explicit-solvent molecular dynamics (MD) simulations is the refinement of protein structures from low-resolution experiments or template-based modeling. A critical requirement is that the native structure is stable with the force field. We have applied a recently developed residue-specific force field, RSFF1, to a set of 30 refinement targets from recent CASP experiments. Starting from their experimental structures, 1.0 μs unrestrained simulations at 298 K retain most of the native structures quite well except for a few flexible terminals and long internal loops. Starting from each homology model, a 150 ns MD simulation at 380 K generates the best RMSD improvement of 0.85 Å on average. The structural improvements roughly correlate with the RMSD of the initial homology models, indicating possible consistent structure refinement. Finally, targets TR614 and TR624 have been subjected to long-time replica-exchange MD simulations. Significant structural improvements are generated, with RMSD of 1.91 and 1.36 Å with respect to their crystal structures. Thus, it is possible to achieve realistic refinement of protein structure models to near-experimental accuracy, using accurate force field with sufficient conformational sampling.

Probing the Huntingtin 1-17 Membrane Anchor on a Phospholipid Bilayer by Using All-Atom Simulations

Mislocalization and aggregation of the huntingtin protein are related to Huntington’s disease. Its first exon—more specifically the first 17 amino acids (Htt17)—is crucial for the physiological and pathological functions of huntingtin. It regulates huntingtin’s activity through posttranslational modifications and serves as an anchor to membrane-containing organelles of the cell. Recently, structure and orientation of the Htt17 membrane anchor were determined using a combined solution and solid-state NMR approach. This prompted us to refine this model by investigating the dynamics and thermodynamics of this membrane anchor on a POPC bilayer using all-atom, explicit solvent molecular dynamics and Hamiltonian replica exchange. Our simulations are combined with various experimental measurements to generate a high-resolution atomistic model for the huntingtin Htt17 membrane anchor on a POPC bilayer. More precisely, we observe that the single α-helix structure is more stable in the phospholipid membrane than the NMR model obtained in the presence of dodecylphosphocholine detergent micelles. The resulting Htt17 monomer has its hydrophobic plane oriented parallel to the bilayer surface. Our results further unveil the key residues interacting with the membrane in terms of hydrogen bonds, salt-bridges, and nonpolar contributions. We also observe that Htt17 equilibrates at a well-defined insertion depth and that it perturbs the physical properties—order parameter, thickness, and area per lipid—of the bilayer in a manner that could favor its dimerization. Overall, our observations reinforce and refine the NMR measurements on the Htt17 membrane anchor segment of huntingtin that is of fundamental importance to its biological functions.

Stability of ion triplets in ionic liquid/lithium salt solutions: Insights from implicit and explicit solvent models and molecular dynamics simulations

Binding energies of ion triplets formed in ionic liquids by Li(+) with two anions have been studied using quantum-chemical calculations with implicit and explicit solvent supplemented by molecular dynamics (MD) simulations. Explicit solvent approach confirms variation of solute-ionic liquid interactions at distances up to 2 nm, resulting from structure of solvation shells induced by electric field of the solute. Binding energies computed in explicit solvent and from the polarizable continuum model approach differ largely, even in sign, but relative values generally agree between these two models. Stabilities of ion triplets obtained in quantum-chemical calculations for some systems disagree with MD results; the discrepancy is attributed to the difference between static optimized geometries used in quantum chemical modeling and dynamic structures of triplets in MD simulations.

PH-dependent transient conformational states control optical properties in cyan fluorescent protein.

A recently engineered mutant of cyan fluorescent protein (WasCFP) that exhibits pH-dependent absorption suggests that its tryptophan-based chromophore switches between neutral (protonated) and charged (deprotonated) states depending on external pH. At pH 8.1, the latter gives rise to green fluorescence as opposed to the cyan color of emission that is characteristic for the neutral form at low pH. Given the high energy cost of deprotonating the tryptophan at the indole nitrogen, this behavior is puzzling, even if the stabilizing effect of the V61K mutation in proximity to the protonation/deprotonation site is considered. Because of its potential to open new avenues for the development of optical sensors and photoconvertible fluorescent proteins, a mechanistic understanding of how the charged state in WasCFP can possibly be stabilized is thus important. Attributed to the dynamic nature of proteins, such understanding often requires knowledge of the various conformations adopted, including transiently populated conformational states. Transient conformational states triggered by pH are of emerging interest and have been shown to be important whenever ionizable groups interact with hydrophobic environments. Using a combination of the weighted-ensemble sampling method and explicit-solvent constant pH molecular dynamics (CPHMD(MSλD)) simulations, we have identified a solvated transient state, characterized by a partially open β-barrel where the chromophore pKa of 6.8 is shifted by over 20 units from that of the closed form (6.8 and 31.7, respectively). This state contributes a small population at low pH (12% at pH 6.1) but becomes dominant at mildly basic conditions, contributing as much as 53% at pH 8.1. This pH-dependent population shift between neutral (at pH 6.1) and charged (at pH 8.1) forms is thus responsible for the observed absorption behavior of WasCFP. Our findings demonstrate the conditions necessary to stabilize the charged state of the WasCFP chromophore (namely, local solvation at the deprotonation site and a partial flexibility of the protein β-barrel structure) and provide the first evidence that transient conformational states can control optical properties of fluorescent proteins.

Insights into the structure of intrastrand cross-link DNA lesion-containing oligonucleotides: G[8-5m]T and G[8-5]C from molecular dynamics simulations.

Oxidatively generated complex DNA lesions occur more rarely than single-nucleotide defects, yet they play an important role in carcinogenesis and aging diseases because they have proved to be more mutagenic than simple lesions. Whereas their formation pathways are rather well understood, the field suffers from the absence of structural data that are crucial for interpreting the lack of repair. No experimental structures are available for oligonucleotides featuring such a lesion. Hence, the detailed structural basis of such damaged duplexes has remained elusive. We propose the use of explicit solvent molecular dynamics simulations to build up damaged oligonucleotides containing two intrastrand cross-link defects, namely, the guanine-thymine and guanine-cytosine defects. Each of these lesions, G[8-5m]T and G[8-5]C, is placed in the middle of a dodecameric sequence, which undergoes an important structural rearrangement that we monitor and analyze. In both duplexes, the structural evolution is dictated by the more favorable stacking of guanine G6, which aims to restore π-stacking with the 3' purine nucleobase. Subsequently, transient formation of hydrogen bonds with a strand shifting is observed. Our simulations are combined with density functional theory to rationalize the structural evolution. We report converging computational evidence that the G[8-5m]T- and G[8-5]C-containing structures evolve toward "abasic-like" duplexes, with a stabilization of the interstrand pairing noncovalent interactions. Meanwhile, both lesions restore B-helicity within tens of nanoseconds. The identification of plausible structures characterizes the last hydrogen abstraction step toward the formation of such defects as a non-innocent chemical reaction.

Toward the correction of effective electrostatic forces in explicit-solvent molecular dynamics simulations: restraints on solvent-generated electrostatic potential and solvent polarization.

Despite considerable advances in computing power, atomistic simulations under nonperiodic boundary conditions, with Coulombic electrostatic interactions and in systems large enough to reduce finite-size associated errors in thermodynamic quantities to within the thermal energy, are still not affordable. As a result, periodic boundary conditions, systems of microscopic size and effective electrostatic interaction functions are frequently resorted to. Ensuing artifacts in thermodynamic quantities are nowadays routinely corrected a posteriori, but the underlying configurational sampling still descends from spurious forces. The present study addresses this problem through the introduction of on-the-fly corrections to the physical forces during an atomistic molecular dynamics simulation. Two different approaches are suggested, where the force corrections are derived from special potential energy terms. In the first approach, the solvent-generated electrostatic potential sampled at a given atom site is restrained to a target value involving corrections for electrostatic artifacts. In the second approach, the long-range regime of the solvent polarization around a given atom site is restrained to the Born polarization, i.e., the solvent polarization corresponding to the ideal situation of a macroscopic system under nonperiodic boundary conditions and governed by Coulombic electrostatic interactions. The restraints are applied to the explicit-water simulation of a hydrated sodium ion, and the effect of the restraints on the structural and energetic properties of the solvent is illustrated. Furthermore, by means of the calculation of the charging free energy of a hydrated sodium ion, it is shown how the electrostatic potential restraint translates into the on-the-fly consideration of the corresponding free-energy correction terms. It is discussed how the restraints can be generalized to situations involving several solute particles. Although the present study considers a very simple system only, it is an important step toward the on-the-fly elimination of finite-size and approximate-electrostatic artifacts during atomistic molecular dynamics simulations.Electronic supplementary materialThe online version of this article (doi:10.1007/s00214-014-1600-8) contains supplementary material, which is available to authorized users.

DMD4B-HYDRA: Toward a Novel Discrete Molecular Dynamics Protein Model

The task of protein structure prediction is computationally demanding, particularly for intrinsically disordered proteins (IDPs), which by definition lack a well-defined folded structure. IDPs include amyloidogenic proteins associated with human disease, such as Alzheimer's and Parkinson's diseases, which are known to be critically involved in the corresponding pathologies. For most protein and peptide sequences with significantly more than 100 amino acids, it is not yet possible to examine the details of the folding landscape by fully atomistic molecular dynamics (MD) in explicit solvent. This problem can be solved through a two-stage MD, first exploring the protein dynamics using coarse-grain (CG) protein models to sample the phase space at a diminished computational effort and second by using a wide range of CG structures converted into fully atomistic representation as initial conformations for all-atom MD in explicit solvent. We here use the four-bead protein model in combination with discrete molecular dynamics (DMD), known as DMD4B-HYDRA force field, which has been previously successfully applied to studies of amyloid b-protein folding and assembly, to examine protein folding dynamics of five proteins with known native structures. We elucidate both the ability of DMD4B-HYDRA force field to capture the folded structure of these proteins and aspects that are less well reproduced by the current force field. We then demonstrate that the follow-up all-atom MD dynamics compensates for structure prediction difficulties of the DMD4B-HYDRA force field itself and provides feedback to be used in the future DMD4B-HYDRA force field development. Our preliminary results thus confirm that DMD4B-HYDRA combined with fully atomistic MD in explicit solvent presents a powerful computational approach to studying IDPs and elucidate structural dynamics of a broad range of amyloidogenic proteins.

Constant pH Molecular Dynamics in Explicit Solvent with Enveloping Distribution Sampling and Hamiltonian Exchange

We present a new computational approach for constant pH simulations in explicit solvent based on the combination of the enveloping distribution sampling (EDS) and Hamiltonian replica exchange (HREX) methods. Unlike constant pH methods based on variable and continuous charge models, our method is based on discrete protonation states. EDS generates a hybrid Hamiltonian of different protonation states. A smoothness parameter s is used to control the heights of energy barriers of the hybrid-state energy landscape. A small s value facilitates state transitions by lowering energy barriers. Replica exchange between EDS potentials with different s values allows us to readily obtain a thermodynamically accurate ensemble of multiple protonation states with frequent state transitions. The analysis is performed with an ensemble obtained from an EDS Hamiltonian without smoothing, s = ∞, which strictly follows the minimum energy surface of the end states. The accuracy and efficiency of this method is tested on aspartic...

Alternative Base-Pairing and Conformational Sampling in Loop a of the Hairpin Ribozyme

The catalytic cycle of the hairpin ribozyme begins with the ‘docking’ between loops A and B and the consideration of this step as the rate limiting step for the nucleolytic catalysis presents an intriguing question. The kinetic rates of docking previously reported by ourselves and others indicate that conformational sampling may be a factor that limits the docking transition rates. To assess conformational sampling in the pre-docking form of the hairpin ribozyme, we have run several explicit solvent Molecular Dynamics (MD) simulations of hairpin ribozyme loop A domain totaling 2.4 µs. We observed three dominant conformers and other minor states identified using hydrogen bonding and base stacking in the loop region. Targeted Molecular Dynamics (TMD) was used to model the pathway between a pair of conformations using CHARMM36/NAMD. The Molecular Mechanics Poisson-Boltzmann Surface Area (MMPB/SA) approximation was applied to predict conformational free energy and energy differences between various conformations. MMPB/SA correctly predicted the major conformations to be lower in free energy than the minor conformations. The observation of multiple conformations with different energies seem to underscore the formation of a rugged energy landscape with multiple accessible shallow energy minima for loop A of the hairpin ribozyme, consistent with a possible role for conformational sampling in the observed slow rates of docking.

Fluctuations within the Hydrogen Bond Network Modulate Cooperativity Across the Conformational Ensemble of Protein Structures

Proteins are dynamic molecules that exhibit conformational changes limited by cross-linking H-bonds. Here we monitor how fluctuations within the hydrogen bond network (HBN) alter protein flexibility. Specifically, we characterize fluctuations in the HBN within two TEM β-lactamase enzymes by employing the Distance Constraint Model (DCM) on 100 ns explicit solvent molecular dynamics (MD) trajectories. A key feature of the DCM is that it employs network rigidity explicitly as a long-range mechanical interaction that accounts for nonadditivity in conformational entropy. Taking advantage of the fundamental link between mechanical and thermodynamic properties, we quantify protein flexibility using a variety of metrics that are calculated within the native state ensemble, and elucidate how fluctuations in the HBN modulate protein flexibility. We observe that the overall flexibility of the proteins is weakly anti-correlated to the overall strength of the HBN. Interestingly, equilibrium fluctuations lead to dramatic loss/gain in correlated rigidity throughout the structure. This translates to conformations that are more flexible, characterized by redundant cross-linking H-bonds break more readily upon slight temperature increases. These conformations also have tendency towards reduced cooperativity within unfolding transitions. The essential property leading to cooperative and non-cooperative structural networks is H-bond independence, which strongly depends on how the H-bonds are distributed. Dense regions of H-bonds that form rigid clusters do so with reduced conformational entropy reductions, thus increasing cooperativity within the folding transition and increasing the co-rigid nature of the native structure. A less dense HBN reduces molecular cooperativity because H-bonds are more likely to be independent. While cooperativity is a hallmark of protein stability, the extent of this cooperativity is modulated by fluctuations within the HBN. Consequently, native state HBN fluctuations have important ramifications to function, suggesting a more complete paradigm of structure/flexibility/function.

Mechanism of Proton Transport of the M2 Proton Channel Studied by Constant pH Molecular Dynamics

The M2 protein of the influenza virus is a proton-selective homotetrameric channel. During virus entry, M2 is activated by the low pH of the endosome and transports protons into the virion, initiating viral uncoating. A histidine tetrad in the pore of the M2 transmembrane domain is responsible for pH activation and proton selectivity. A number of different mechanisms have been proposed for proton transport of M2 based on experimental and computational studies. To test these hypotheses, we applied the explicit-solvent continuous constant pH molecular dynamics method to study the pH-dependent conformational dynamics of M2 in explicit membrane. The calculated pKa values of the histidine tetrad are comparable to experimental values. The C-terminal openning of M2 became wider when pH was lowered. This work provides novel insight into the coupled protonation and conformational dynamics of the M2 proton channel.

Energetic contributions of residues to the formation of early amyloid-β oligomers.

Low-weight amyloid-β (Aβ) oligomers formed at early stages of oligomerization rather than fibril assemblies seem to be the toxic components that drive neurodegeneration in Alzheimer's disease. Unfortunately, detailed knowledge of the structure of these early oligomers at the residue level is not yet available. In this study, we performed all-atom explicit solvent molecular dynamics simulations to examine the oligomerization process of Aβ10-35 monomers when forming dimers, trimers, tetramers and octamers, with four independent simulations of a total simulated time of 3 μs for each oligomer system. The decomposition of the stability free energy by MM-GBSA methodology allowed us to unravel the network of energetic interactions that stabilize such oligomers. The contribution of the intermonomeric van der Waals term is the most significant energy feature of the oligomerization process, consistent with the so-called hydrophobic effect. Furthermore, the decomposition of the stability free energy into residues and residue-pairwise terms revealed that it is mainly apolar interactions between the three specific hydrophobic fragments 31-35 (C-terminal region), 17-20 (central hydrophobic core) and 12-14 (N-terminal region) that are responsible for such a favourable effect. The conformation in which the hydrophobic cthr-chc interaction is oriented perpendicularly is particularly important. We propose three other model substructures that favour the oligomerization process and can thus be considered as molecular targets for future inhibitors. Understanding Aβ oligomerization at the residue level could lead to more efficient design of inhibitors of this process.

Interpretation of Solution X-Ray Scattering Data by Molecular Dynamics

Small and wide angle X-ray scattering (SWAXS) and molecular dynamics (MD) simulations are complementary approaches that probe conformational transitions of biomolecules in solution, even in a time-resolved manner. However, the structural interpretation of the scattering signals is challenging, while MD simulations frequently suffer from incomplete sampling or from a force field bias.To combine the advantages of both techniques, we present a method that incorporates solution scattering data into explicit-solvent MD simulations, termed SWAXS-driven MD, with the aim to direct the simulation into conformations satisfying the experimental data [1]. Because the calculations fully rely on explicit solvent, no fitting parameters associated with the solvation layer or excluded solvent are required, and the calculations remain valid at wide angles [2].The complementarity of SWAXS and MD is illustrated using three biological examples, namely a periplasmic binding protein, aspartate carbamoyltransferase, and a nuclear exportin. The examples suggest that SWAXS-driven MD is capable of refining structures against SWAXS data without foreknowledge of possible reaction paths. In turn, the SWAXS data accelerates conformational transitions in MD simulations and reduces the force field bias. Extensions towards the interpretation of anisotropic X-ray scattering data are discussed [3].

Site Identification by Ligand Competitive Saturation (SILCS) simulations for fragment-based drug design.

Fragment-based drug design (FBDD) involves screening low molecular weight molecules ("fragments") that correspond to functional groups found in larger drug-like molecules to determine their binding to target proteins or nucleic acids. Based on the principle of thermodynamic additivity, two fragments that bind nonoverlapping nearby sites on the target can be combined to yield a new molecule whose binding free energy is the sum of those of the fragments. Experimental FBDD approaches, like NMR and X-ray crystallography, have proven very useful but can be expensive in terms of time, materials, and labor. Accordingly, a variety of computational FBDD approaches have been developed that provide different levels of detail and accuracy.The Site Identification by Ligand Competitive Saturation (SILCS) method of computational FBDD uses all-atom explicit-solvent molecular dynamics (MD) simulations to identify fragment binding. The target is "soaked" in an aqueous solution with multiple fragments having different identities. The resulting computational competition assay reveals what small molecule types are most likely to bind which regions of the target. From SILCS simulations, 3D probability maps of fragment binding called "FragMaps" can be produced. Based on the probabilities relative to bulk, SILCS FragMaps can be used to determine "Grid Free Energies (GFEs)," which provide per-atom contributions to fragment binding affinities. For essentially no additional computational overhead relative to the production of the FragMaps, GFEs can be used to compute Ligand Grid Free Energies (LGFEs) for arbitrarily complex molecules, and these LGFEs can be used to rank-order the molecules in accordance with binding affinities.

A Computational Study of Amyloid β-Protein Assembly in Crowded Environments

Alzheimer's disease is strongly associated with aberrant amyloid β-protein (Aβ) assembly into heterogeneous, metastable oligomeric assemblies with structures that have not been experimentally characterized yet. The 40 and 42 amino acids long Aβ40 and Aβ42 are the two predominant Aβ alloforms in the brain. Whereas Aβ40 and Aβ42 oligomer formation from monomeric state is still inaccessible to fully atomistic explicit-solvent molecular dynamics, Aβ40 and Aβ42 oligomers were structurally characterized using discrete molecular dynamics (DMD) and an intermediate-resolution protein model within the DMD4B-HYDRA implicit solvent force field, and the corresponding oligomer size distributions well matched the available in vitro data. In vivo, however, Aβ coexists with other biomolecules in a rather crowded environment. To understand the effect of crowding on Aβ oligomer formation, we used the DMD4B-HYDRA force field and added to an ensemble of 32 monomeric Aβ40 or Aβ42 peptides inert spherical “crowders” with a diameter of 0.5 nm at various concentrations to examine their effect on Aβ40 and Aβ42 oligomerization pathways. Our results show that crowding shifts oligomer size distributions towards smaller oligomer sizes and increases solubility of both peptides in a concentration-dependent way. The effect is stronger for Aβ42, where crowding abolishes the multimodal character of the oligomer size distribution. Our structural analysis revealed that the stability of larger oligomers is compromised by effective osmotic pressure exerted by the crowders, resulting in an increased rate of assembly breakage. While in vivo crowding agents are not inert as the crowders in our study, we here reveal that crowding-induced osmotic pressure strongly affects protein assembly dynamics, which is of significance to the disease.

Conformational Preferences Underlying Reduced Activity of a Thermophilic Ribonuclease H

The conformational basis for reduced activity of the thermophilic ribonuclease HI enzyme from Thermus thermophilus, compared to its mesophilic homolog from Escherichia coli, is elucidated using a combination of NMR spectroscopy and molecular dynamics (MD) simulations. Explicit-solvent all-atom MD simulations of the two wild-type proteins and an E. coli mutant in which a glycine residue is inserted after position 80 to mimic the T. thermophilus protein reproduce the differences in conformational dynamics determined from 15N spin-relaxation NMR spectroscopy of three loop regions that surround the active site and contain functionally important residues: the glycine-rich region, the handle region, and the β5/αE loop. Examination of the MD trajectories indicates that the thermophilic protein samples conformations productive for substrate binding and activity less frequently than the mesophilic enzyme, although these differences may manifest as either increased or decreased relative flexibility of the different regions. Additional MD simulations indicate that mutations increasing activity of the T. thermophilus enzyme at mesophilic temperatures do so by reconfiguring the local environments of the mutated sites to more closely resemble active conformations. Taken together, the results show that both locally increased and decreased flexibility contribute to an overall reduction in activity of T. thermophilus ribonuclease H compared to its mesophilic E. coli homolog.

COFFDROP: A Coarse-Grained Nonbonded Force Field for Proteins Derived from All-Atom Explicit-Solvent Molecular Dynamics Simulations of Amino Acids.

We describe the derivation of a set of bonded and nonbonded coarse-grained (CG) potential functions for use in implicit-solvent Brownian dynamics (BD) simulations of proteins derived from all-atom explicit-solvent molecular dynamics (MD) simulations of amino acids. Bonded potential functions were derived from 1 μs MD simulations of each of the 20 canonical amino acids, with histidine modeled in both its protonated and neutral forms; nonbonded potential functions were derived from 1 μs MD simulations of every possible pairing of the amino acids (231 different systems). The angle and dihedral probability distributions and radial distribution functions sampled during MD were used to optimize a set of CG potential functions through use of the iterative Boltzmann inversion (IBI) method. The optimized set of potential functions—which we term COFFDROP (COarse-grained Force Field for Dynamic Representation Of Proteins)—quantitatively reproduced all of the “target” MD distributions. In a first test of the force field, it was used to predict the clustering behavior of concentrated amino acid solutions; the predictions were directly compared with the results of corresponding all-atom explicit-solvent MD simulations and found to be in excellent agreement. In a second test, BD simulations of the small protein villin headpiece were carried out at concentrations that have recently been studied in all-atom explicit-solvent MD simulations by Petrov and Zagrovic (PLoS Comput. Biol.2014, 5, e1003638). The anomalously strong intermolecular interactions seen in the MD study were reproduced in the COFFDROP simulations; a simple scaling of COFFDROP’s nonbonded parameters, however, produced results in better accordance with experiment. Overall, our results suggest that potential functions derived from simulations of pairwise amino acid interactions might be of quite broad applicability, with COFFDROP likely to be especially useful for modeling unfolded or intrinsically disordered proteins.

Mechanical Model of DNA Allostery.

The importance of allosteric effects in DNA is becoming increasingly appreciated, but the underlying mechanisms remain poorly understood. In this work, we propose a general modeling framework to study DNA allostery. We describe DNA in a coarse-grained manner by intra-base pair and base pair step coordinates, complemented by groove widths. Quadratic deformation energy is assumed, yielding linear relations between the constraints and their effect. Model parameters are inferred from standard unrestrained, explicit-solvent molecular dynamics simulations of naked DNA. We applied the approach to study minor groove binding of diamidines and pyrrole-imidazole polyamides. The predicted DNA bending is in quantitative agreement with experiment and suggests that diamidine binding to the alternating TA sequence brings the DNA closer to the A-tract conformation, with potentially important functional consequences. The approach can be readily applied to other allosteric effects in DNA and generalized to model allostery in various molecular systems.

MoD-QM/MM Structural Refinement Method: Characterization of Hydrogen Bonding in the Oxytricha nova G-Quadruplex.

A generalization of the Moving-Domain Quantum Mechanics/Molecular Mechanics (MoD-QM/MM) hybrid method [Gascon, J. A.; Leung, S. S. F.; Batista, E. R.; Batista, V. S. J. Chem. Theory Comput. 2006, 2, 175-186] is introduced to provide a self-consistent computational protocol for structural refinement of extended systems. The method partitions the system into molecular domains that are iteratively optimized as quantum mechanical (QM) layers embedded in their surrounding molecular environment to obtain an ab initio quality description of the geometry and the molecular electrostatic potential of the extended system composed of those constituent fragments. The resulting methodology is benchmarked as applied to model systems that allow for full QM optimization as well as through refinement of the hydrogen bonding geometry in Oxytricha nova guanine quadruplex for which several studies have been reported, including the X-ray structure and NMR data. Calculations of (1)H NMR chemical shifts based on the gauge independent atomic orbital (GIAO) method and direct comparisons with experiments show that solvated MoD-QM/MM structures, sampled from explicit solvent molecular dynamics simulations, allow for NMR simulations in much improved agreement with experimental data than models based on the X-ray structure or those optimized using classical molecular mechanics force fields.

Molecular dynamics in drug design.

Molecular dynamics (MD) simulations are useful tools for structure-based drug design. We review recent publications in which explicit solvent MD was used at the initial or final stages of high-throughput docking campaigns. In some cases, MD simulations of the protein target have been carried out before docking to generate a conformer of the protein which differs from the available crystal structure(s). Furthermore, MD runs have been performed after docking to assess the predicted binding modes of the top ranking compounds as final filter in silico or to guide chemical synthesis for hit optimization. We present examples of in silico discoveries of tyrosine kinase inhibitors and bromodomain antagonists whose binding mode was predicted by automated docking and further corroborated by MD simulations with final validation by X-ray crystallography.