The future of molecular dynamics simulations in drug discovery

Molecular dynamics simulations can now track rapid processes—those occurring in less than about a millisecond—at atomic resolution for many biologically relevant systems. These simulations appear poised to exert a significant impact on how new drugs are found, perhaps even transforming the very process of drug discovery. We predict here future results we can expect from, and enhancements we need to make in, molecular dynamics simulations over the coming 25 years, and in so doing set out several Grand Challenges for the field. In the context of the problems now facing the pharmaceutical industry, we ask how we can best address drug discovery needs of the next quarter century using molecular dynamics simulations, and we suggest some possible approaches.     

Retinal counterion switch mechanism in vision evaluated by molecular simulations

Photoisomerization of the retinylidene chromophore of rhodopsin is the starting point in the vision cascade. A counterion switch mechanism that stabilizes the retinal protonated Schiff base (PSB) has been proposed to be an essential step in rhodopsin activation. On the basis of vibrational and UV-visible spectroscopy, two counterion switch models have emerged. In the first model, the PSB is stabilized by Glu181 in the meta I state, while in the most recent proposal, it is stabilized by Glu113 as well as Glu181. We assess these models by conducting a pair of microsecond scale, all-atom molecular dynamics simulations of rhodopsin embedded in a 99-lipid bilayer of SDPC, SDPE, and cholesterol (2:2:1 ratio) varying the starting protonation state of Glu181. Theoretical simulations gave different orientations of retinal for the two counterion switch mechanisms, which were used to simulate experimental 2H NMR spectra for the C5, C9, and C13 methyl groups. Comparison of the simulated 2H NMR spectra with experimental data supports the complex-counterion mechanism. Hence, our results indicate that Glu113 and Glu181 stabilize the retinal PSB in the meta I state prior to activation of rhodopsin.     

Understanding the nonfouling mechanism of surfaces through molecular simulations of sugar-based self-assembled monolayers

This paper presents a molecular simulation study of the interactions of a protein (lysozyme) with self-assembled monolayers (SAMs) of mannitol and sorbitol terminated alkanethiols in the presence of explicit water molecules and ions. The all-atom simulations were performed to calculate the force generated on the protein as a function of its distance above the SAM surfaces. The structural and dynamic properties of water molecules both above the SAM surfaces and around the SAM head groups were analyzed to provide a better understanding of the nonfouling behavior of the sugar-based SAM surfaces. Results from this work suggest that both mannitol and sorbitol SAMs generate a tightly bound, structured water layer around the SAM chains. This hydration layer creates a repulsive force on the protein when it approaches the surface, resulting in a nonfouling surface despite the presence of hydrogen-bond donor groups. This work demonstrates the importance of strong surface-water interactions for surface resistance to nonspecific protein adsorption.     

Specificity and molecular mechanism of abortificient action of prostaglandins

The paper reports comparison of electrostatic charge and energy distribution on the basis of the CNDO /2 method for six forms of prostaglandins–PGF_2α, PGF_1β, PGE₁, PGE₂, PGB₁, and PGA₁–having diverse physiological action. The isopotential mapping done in three dimension showed that the lower value of electrostatic potential and proximity of the two low energy regions around O₉ and O₁₁ in PGE₂ and PGF_2α is probably responsible for their higher abortificient activity. We also compare here the variation of the long- and short-range interaction between ring–chain and chain–chain portion of different forms and compared them with the variation in their action.     

The mechanism of vesicle fusion as revealed by molecular dynamics simulations

We describe molecular dynamics simulations elucidating the molecular details of the process of fusion for small lipid vesicles. The simulations are based on a coarse grained (CG) lipid model that accurately represents the lamellar state of a variety of phospholipids and enables us to observe intermediate stages during fusion at near atomic detail. Simulations were conducted on a variety of systems containing common phospholipids such as phosphatidylcholine (PC), phosphatidylethanolamine (PE), lysoPC, and mixtures of the above. The fusion intermediates found are in general agreement with the stalk-pore mechanism. Transient pores sometimes form adjacent to the stalk, however, resulting in the mixing of lipids from the outer and inner monolayers. The speed of stalk formation and the opening of the fusion pore can be modulated by altering the lipid composition in qualitative agreement with experimental observations.     

Dynamic mechanism for the autophosphorylation of CheA histidine kinase: molecular dynamics simulations

The two-component system (TCS) is an important signal transduction component for most bacteria. This signaling pathway is mediated by histidine kinases via autophosphorylation between P1 and P4 domains. Taking chemotaxis protein CheA as a model of TCS, the autophosphorylation mechanism of the TCS histidine kinases has been investigated in this study by using a computational approach integrated homology modeling, ligand-protein docking, protein-protein docking, and molecular dynamics (MD) simulations. Four nanosecond-scale MD simulations were performed on the free P4 domain, P4-ATP, P4-TNPATP, and P1-P4-ATP complexes, respectively. Upon its binding to the binding pocket of P4 with a folded conformation, ATP gradually extends to an open state with help from a water molecule. Meanwhile, ATP forms two hydrogen bonds with His413 and Lys494 at this state. Because of the lower energy of the folded conformations, ATP shrinks back to its folded conformations, leading to the rupture of the hydrogen bond between ATP and Lys494. Consequently, Lys494 moves away from the pocket entrance, resulting in an open of the ATP lid of P4. It is the open state of P4 that can bind tightly to P1, where the His45 of P1 occupies a favorable position for its autophosphorylation from ATP. This indicates that ATP is not only a phosphoryl group donor but also an activator for CheA phosphorylation. Accordingly, a mechanism of the autophosphorylation of CheA is proposed as that the ATP conformational switch triggers the opening of the ATP lid of P4, leading to P1 binding tightly, and subsequently autophosphorylation from ATP to P1.     

Steered molecular dynamics simulations for uncovering the molecular mechanisms of drug dissociation and for drug screening: A test on the focal adhesion kinase

Drug‐binding kinetics could play important roles in determining the efficacy of drugs and has caught the attention of more drug designers. Using the dissociation of 1H‐pyrrolo[2,3‐b]‐pyridines from the focal adhesion kinase as an example, this work finds that steered molecular dynamics simulations could help screen compounds with long‐residence times. It also reveals a two‐step mechanism of ligand dissociation resembling the release of ADP from protein kinase A reported earlier. A phenyl group attaching to the pyrrole prolongs residence time by creating a large activation barrier for transition from the bound to the intermediate state when it becomes exposed to the solvent. Principal component analysis shows that ligand dissociation does not couple with large‐scale collective motions of the protein involving many of its amino acids. Rather, a small subset of amino acids dominates. Some of these amino acids do not contact the ligands directly along the dissociation pathways and could exert long‐range allosteric effects. © 2018 Wiley Periodicals, Inc.     

Reaction mechanism of soluble epoxide hydrolase: insights from molecular dynamics simulations

Molecular dynamics simulations have been performed to gain insights into the catalytic mechanism of the hydrolysis of epoxides to vicinal diols by soluble epoxide hydrolase (sEH). The binding of a substrate, 1S,2S-trans-methylstyrene oxide, was studied in two conformations in the active site of the enzyme. It was found that only one is likely to be found in the active enzyme. In the preferred conformation the phenyl group of the substrate is pi-sandwiched between two aromatic residues, Tyr381 and His523, whereas the other conformation is pi-stacked with only one aromatic residue, Trp334. Two simulations were carried out to 1 ns for each conformation to evaluate the protonation state of active site residue His523. It was found that a protonated histidine is essential for keeping the active site from being disrupted. Long time scale, 4 ns, molecular dynamics simulation was done for the structure with the most likely combination of binding conformation and protonation state of His523. Near Attack Conformers (NACs) are present 5.3% of the time and nucleophilic attack on either epoxide carbon atom, approximately 75% on C(1) and approximately 25% on C(2), is found. A maximum of one hydrogen bond between the epoxide oxygen and either of the active site tyrosines, Tyr465 and Tyr381, is present, in agreement with experimental mutagenesis results that reveal a slight loss in activity if one tyrosine is mutated and essential loss of all activity upon double mutation of the two tyrosines in question. It was found that a hydrogen bond from Tyr465 to the substrate oxygen is essential for controlling the regioselectivity of the reaction. Furthermore, a relationship between the presence of this hydrogen bond and the separation of reactants was found. Two groups of amino acid segments were identified each as moving collectively. Furthermore, an overall anti-correlation was found between the movements of these two individually collectively moving groups, made up by parts of the cap-region, including the two tyrosines, and the site of the catalytic triad, respectively. This overall anti-correlated collective domain motion is, perhaps, involved in the conversion of E.NAC to E.TS.     

Cucurbituril-resisted acylation of the anti-tuberculosis drug isoniazid via a supramolecular strategy

A chemical investigation reveals that the resistance to acylation of an anti-tuberculosis drug, isoniazid is a consequent result of the inclusion or exclusion of cucurbit[n]urils (n = 6 or 7). The (1)H NMR spectra analysis shows that the different interaction models of the isoniazid with the two cucurbiturils are dependent on the cavity size of the hosts. Quantum chemistry calculations with density functional theory method indicate that the interaction of the isoniazid with both cucurbiturils is through thermodynamic stabilization in both the gas phase and aqueous solution through hydrogen bonding on the portal carbonyls of the cucurbiturils. Electronic absorption titration spectra suggest the hosts and guest interact in a ratio of 1 : 1 with moderate binding constants. Acylation kinetics of isoniazid with various acylating agents in the presence of the cucurbiturils revealed that resistance is only dependent on the host-isoniazid ratio, and independent on the size of the cucurbiturils and the species of acylating agents.     

Wax inhibition by comb-like polymers: support of the incorporation-perturbation mechanism from molecular dynamics simulations

Deposition of wax on a cold surface is a serious problem in oil production. Progress in developing more effective wax inhibitors has been impeded by the lack of an established mechanism connecting the molecular structure to inhibitor efficiency. Some comb-like polymers having long alkyl side chains are known to decrease the rate of wax formation. Among several possible mechanisms, we investigate here the incorporation-perturbation mechanism. According to this mechanism, the inhibitor molecules in oil are preferentially partitioned (incorporation) toward the wax-rich (amorphous) wax deposits (soft wax), which then serves as a perturbation to slow down the ordering transition of soft amorphous wax into more stable but problematic hard wax crystals. Indeed, molecular dynamics simulations on an effective inhibitor molecule in both the oil phase and in the amorphous wax phase support the idea that the oil-to-wax partition of the inhibitor is energetically favorable. With the inhibitor molecule embedded, the structure of wax crystal is disturbed, significantly decreasing the order and significantly lowering the cohesive energy density relative to that of the pure wax crystal, supporting the slower transition from soft wax to hard wax. Thus, in the presence of an effective wax inhibitor, crystallization (formation of hard wax) is slowed dramatically, so that there is time to flush out the soft wax with a high-pressure flow inside the pipeline. This suggests design principles for developing improved wax inhibitors.     

Investigation of chiral recognition by molecular micelles with molecular dynamics simulations

Molecular dynamics simulations were used to characterize the binding of the chiral drugs chlorthalidone and lorazepam to the molecular micelle poly-(sodium undecyl-(L)-leucine-valine). The project’s goal was to characterize the nature of chiral recognition in capillary electrophoresis separations that use molecular micelles as the chiral selector. The shapes and charge distributions of the chiral molecules investigated, their orientations within the molecular micelle chiral binding pockets, and the formation of stereoselective intermolecular hydrogen bonds with the molecular micelle were all found to play key roles in determining where and how lorazepam and chlorthalidone enantiomers interacted with the molecular micelle.     

Polycaprolactone nanocomposite systems used to deliver ifosfamide anticancer drug: molecular dynamics simulations

Molecular dynamics (MD) simulations were accomplished on polycaprolactone (PCL) nanocomposite systems composed of hydroxyapatite (HA) nanoparticles (0–4%) to deliver ifosfamide (IF) anticancer drug in order to achieve the most suitable drug delivery system (DDS). It was shown that the free volume (FV) was the greatest for the PCL-HA2-IF whereas the lowest value was measured for the PCL-HA0-IF. The fractional free volume (FFV) values varied similar to the FV values so that the PCL-HA2-IF had the maximum FFV (22.48%) but PCL-HA0-IF illustrated the minimum FFV (17.89%). The smallest interchain distances measured for the PCL-HA2-IF established that the greatest intermolecular interactions occurred in the PCL-HA2-IF. The highest diffusion coefficient (0.1267?×?10^?4 cm²/s) was obtained for the PCL-HA0-IF whereas the lowest one was achieved for the PCL-HA2-IF (0.0688?×?10^?4 cm²/s) that confirmed the drug diffusion was the slowest/most controlled in the PCL-HA2-IF which would bring about the most effective drug delivery.

     

Insights into the acylation mechanism of class A beta-lactamases from molecular dynamics simulations of the TEM-1 enzyme complexed with benzylpenicillin

Herein, we present results from molecular dynamics MD simulations ( approximately 1 ns) of the TEM-1 beta-lactamase in aqueous solution. Both the free form of the enzyme and its complex with benzylpenicillin were studied. During the simulation of the free enzyme, the conformation of the Omega loop and the interresidue contacts defining the complex H-bond network in the active site were quite stable. Most interestingly, the water molecule connecting Glu166 and Ser70 does not exchange with bulk solvent, emphasizing its structural and catalytic relevance. In the presence of the substrate, Ser130, Ser235, and Arg244 directly interact with the beta-lactam carboxylate via H-bonds, whereas the Lys234 ammonium group has only an electrostatic influence. These interactions together with other specific contacts result in a very short distance ( approximately 3 A) between the attacking hydroxyl group of Ser70 and the beta-lactam ring carbonyl group, which is a favorable orientation for nucleophilic attack. Our simulations also gave insight into the possible pathways for proton abstraction from the Ser70 hydroxyl group. We propose that either the Glu166 carboxylate-Wat1 or the substrate carboxylate-Ser130 moieties could abstract a proton from the nucleophilic Ser70.     

Insight into the C-F bond mechanism of molecular analogs for antibacterial drug design

The activities of biological molecules usually rely on both of intra-molecular and intermolecular interactions between their function groups. These interactions include interonic attraction theory, Van der Waal’s forces and the function of geometry on the individual molecules, whether they are naturally or synthetic. The purpose of this study was to evaluate the antibacterial activity of C-F bond compound using combination of experiments verification and theoretical calculation. We target on the insect natural products from the maggots of Chrysomyis megacephala Fabricius. Based on density functional theory(DFT) and B3LYP method, a theoretical study of the C-F bond on fluoride was designed to explore compounds 2 and 4 antibacterial structure–activity relationship. With the progress in DFT, first-principle calculation based on DFT has gradually become a routine method for drug design, quantum chemistry and other science fields.     

Probing the kinetics of membrane-mediated helix folding

The kinetics of peptide-membrane association have been studied previously using stopped-flow tryptophan fluorescence; however, such experiments do not directly report the coil-to-helix transition process, which is a hallmark of peptide-membrane interaction. Herein, we report a new method for directly assessing the kinetics of the helix formation accompanied by the peptide-membrane association. This method is based on the technique of fluorescence resonance energy transfer (FRET) and an amino acid FRET pair, p-cyano-L-phenylalanine and tryptophan. To demonstrate the utility of this method, we have studied the membrane-mediated helix folding dynamics of a mutant of magainin 2, an antibiotic peptide found in the skin of the African clawed frog, Xenopus laevis. Our results indicate that the coil-to-helix transition occurs during the binding of the peptide to the lipid vesicle (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine/1-palmitoyl-2-oleoyl-sn-glycero-3-[phospho-rac-(1-glycerol)], 3:1, wt/wt) but prior to the full insertion of the peptide into the hydrophobic region of the lipid bilayers.     

Coarse-grained molecular dynamics simulations of nanopatterning with multivalent inks

A coarse-grained molecular dynamics (MD) model is developed to study the multivalent, or multisite, binding of small functionalized dendrimer molecules to beta-cyclodextrin-terminated self-assembled monolayers, the so-called "molecular printboards" used to print "ink" molecules on surfaces with a high degree of positional control and specificity. Some current and future bionanotechnology applications are in the creation of nanoparticle assemblies, directed protein assembly, platforms for biosensing, and cell:surface attachment. The coarse-grained model allows us to probe up to microsecond timescales and model ink diffusion, crucial for the application of the printboard in, for example, medical diagnostics. Recent all-atom MD simulations identified and quantified the molecular strain limiting the stability of nanopatterns created with small dendrimer inks, and explained the different patterns obtained experimentally with different dendrimer inks. In the present work, the all-atom simulations are "scaled up" to longer timescales via coarse graining, without incurring significant additional computational expense, and, crucially, without significant loss in atom-scale detail, the coarse-grained MD simulations yielding properties similar to those obtained from the all-atom simulations. The anchoring of the ink molecules to the monolayer is of multivalent nature and the degree of multivalency shows a sharp dependence on temperature, control of temperature thus providing a further operational "switch" for directed molecular assembly. The computational protocol developed can, in principle, be extended to model any multivalent assembly, for example, virus-cell complexation.     

Probing molecular mechanism of inhibitor bindings to bromodomain-containing protein 4 based on molecular dynamics simulations and principal component analysis

ABSTRACT

It is well known that bromodomain-containing protein 4 (BRD4) has been thought as a promising target utilized for treating various human diseases, such as inflammatory disorders, malignant tumours, acute myelogenous leukaemia (AML), bone diseases, etc. For this study, molecular dynamics (MD) simulations, binding free energy calculations, and principal component analysis (PCA) were integrated together to uncover binding modes of inhibitors 8P9, 8PU, and 8PX to BRD4(1). The results obtained from binding free energy calculations show that van der Waals interactions act as the main regulator in bindings of inhibitors to BRD4(1). The information stemming from PCA reveals that inhibitor associations extremely affect conformational changes, internal dynamics, and movement patterns of BRD4(1). Residue-based free energy decomposition method was wielded to unveil contributions of independent residues to inhibitor bindings and the data signify that hydrogen bonding interactions and hydrophobic interactions are decisive factors affecting bindings of inhibitors to BRD4(1). Meanwhile, eight residues Trp81, Pro82, Val87, Leu92, Leu94, Cys136, Asn140, and Ile146 are recognized as the common hot interaction spots of three inhibitors with BRD4(1). The results from this work are expected to provide a meaningfully theoretical guidance for design and development of effective inhibitors inhibiting of the activity of BRD4.     

A detailed reaction mechanism for hexamethyldisiloxane combustion via experiments and ReaxFF molecular dynamics simulations

Hexamethyldisiloxane (HMDSO) is one of the main impurities in the syngas produced from sewage and landfill plants. In order to utilize this syngas or control the characteristics of the generated silica particles, it is crucial to understand the chemical kinetics of HMDSO combustion. This study investigated the process of HMDSO combustion using synchrotron radiation mass spectrometry (SRMS), gas chromatography (GC), and ReaxFF molecular dynamics simulations. First, the force field used for ReaxFF simulation was validated by comparing the energies of different bond lengths, bond angles, and dihedral angles with the ones from DFT calculations. Good agreements were found. Then, ReaxFF simulations of HMDSO combustion with this force field were conducted under various conditions, which include different equivalence ratios (0.67, 1.0, and 1.5) and temperatures ranging from 2000 to 3500 K. The oxidation characteristics of HMDSO were analyzed, including the evolution of gas products and particle formation. Finally, based on the results from experiments and ReaxFF simulations, the reaction pathways, reaction lists, and reaction kinetics data during HMDSO combustion were obtained. A detailed reaction mechanism was proposed and validated by applying it in modeling the H₂/HMDSO/O₂ combustion systems. The temperature and part of the gas products such as CO and CO₂ as well as SiO could be well predicted.