Characterizing the Free-Energy Landscape of MDM2 Protein-Ligand Interactions by Steered Molecular Dynamics Simulations

Insights into the Interaction Mechanisms of Peptide and Non-Peptide Inhibitors with MDM2 Using Gaussian-Accelerated Molecular Dynamics Simulations and Deep Learning

Inhibiting MDM2-p53 interaction is considered an efficient mode of cancer treatment. In our current study, Gaussian-accelerated molecular dynamics (GaMD), deep learning (DL), and binding free energy calculations were combined together to probe the binding mechanism of non-peptide inhibitors K23 and 0Y7 and peptide ones PDI6W and PDI to MDM2. The GaMD trajectory-based DL approach successfully identified significant functional domains, predominantly located at the helixes α2 and α2', as well as the β-strands and loops between α2 and α2'. The post-processing analysis of the GaMD simulations indicated that inhibitor binding highly influences the structural flexibility and collective motions of MDM2. Calculations of molecular mechanics-generalized Born surface area (MM-GBSA) and solvated interaction energy (SIE) not only suggest that the ranking of the calculated binding free energies is in agreement with that of the experimental results, but also verify that van der Walls interactions are the primary forces responsible for inhibitor-MDM2 binding. Our findings also indicate that peptide inhibitors yield more interaction contacts with MDM2 compared to non-peptide inhibitors. Principal component analysis (PCA) and free energy landscape (FEL) analysis indicated that the piperidinone inhibitor 0Y7 shows the most pronounced impact on the free energy profiles of MDM2, with the piperidinone inhibitor demonstrating higher fluctuation amplitudes along primary eigenvectors. The hot spots of MDM2 revealed by residue-based free energy estimation provide target sites for drug design toward MDM2. This study is expected to provide useful theoretical aid for the development of selective inhibitors of MDM2 family members.

Drug Discovery by Automated Adaptation of Chemical Structure and Identity

Computer-aided drug design offers the potential to dramatically reduce the cost and effort required for drug discovery. While screening-based methods are valuable in the early stages of hit identification, they are frequently succeeded by iterative, hypothesis-driven computations that require recurrent investment of human time and intuition. To increase automation, we introduce a computational method for lead refinement that combines concerted dynamics of the ligand/protein complex via molecular dynamics simulations with integrated Monte Carlo-based changes in the chemical formula of the ligand. This approach, which we refer to as ligand-exchange Monte Carlo molecular dynamics, accounts for solvent- and entropy-based contributions to competitive binding free energies by coupling the energetics of bound and unbound states during the ligand-exchange attempt. Quantitative comparison of relative binding free energies to reference values from free energy perturbation, conducted in vacuum, indicates that ligand-exchange Monte Carlo molecular dynamics simulations sample relevant conformational ensembles and are capable of identifying strongly binding compounds. Additional simulations demonstrate the use of an implicit solvent model. We speculate that the use of chemical graphs in which exchanges are only permitted between ligands with sufficient similarity may enable an automated search to capture some of the benefits provided by human intuition during hypothesis-guided lead refinement.

Thermodynamics and Kinetics of Drug-Target Binding by Molecular Simulation

Computational studies play an increasingly important role in chemistry and biophysics, mainly thanks to improvements in hardware and algorithms. In drug discovery and development, computational studies can reduce the costs and risks of bringing a new medicine to market. Computational simulations are mainly used to optimize promising new compounds by estimating their binding affinity to proteins. This is challenging due to the complexity of the simulated system. To assess the present and future value of simulation for drug discovery, we review key applications of advanced methods for sampling complex free-energy landscapes at near nonergodicity conditions and for estimating the rate coefficients of very slow processes of pharmacological interest. We outline the statistical mechanics and computational background behind this research, including methods such as steered molecular dynamics and metadynamics. We review recent applications to pharmacology and drug discovery and discuss possible guidelines for the practitioner. Recent trends in machine learning are also briefly discussed. Thanks to the rapid development of methods for characterizing and quantifying rare events, simulation's role in drug discovery is likely to expand, making it a valuable complement to experimental and clinical approaches.

Molecular mechanism of inhibitor bindings to bromodomain-containing protein 9 explored based on molecular dynamics simulations and calculations of binding free energies

Recently, bromodomain-containing protein 9 (BRD9) has been a prospective therapeutic target for anticancer drug design. Molecular dynamics (MD) simulations combined with molecular mechanics generalized Born surface area (MM-GBSA) method were adopted to explore binding modes of three inhibitors (5SW, 5U2, and 5U6) to BRD9 and identify the hot spot of the inhibitor-BRD9 binding. The results indicate that the inhibitor 5SW has the strongest binding ability to BRD9 among the current three inhibitors. Furthermore, the rank of the binding free energies predicted by MM-GBSA approach agrees with that determined by the experimental values. In addition, inhibitor-residue interactions were computed by using residue-based free-energy decomposition method and the results suggest that residue His42 produces the CH-H interactions, residues Asn100, Ile53 and Val49 produce the CH-[Formula: see text] interactions with three inhibitors and Tyr106, Phe45 and Phe44 generate the π-π interactions with inhibitors. Notably, the residue Asn140 forms hydrogen bonding interactions with three inhibitors. This research is expected to provide useful molecular basis and dynamics information at atomic levels for the design of potent inhibitors inhibiting the activity of BRD9.

Lessons learned about steered molecular dynamics simulations and free energy calculations

The calculation of free energy profiles is central in understanding differential enzymatic activity, for instance, involving chemical reactions that require QM-MM tools, ligand migration, and conformational rearrangements that can be modeled using classical potentials. The use of steered molecular dynamics (sMD) together with the Jarzynski equality is a popular approach in calculating free energy profiles. Here, we first briefly review the application of the Jarzynski equality to sMD simulations, then revisit the so-called stiff-spring approximation and the consequent expectation of Gaussian work distributions and, finally, reiterate the practical utility of the second-order cumulant expansion, as it coincides with the parametric maximum-likelihood estimator in this scenario. We illustrate this procedure using simulations of CO, both in aqueous solution and in a carbon nanotube as a model system for biologically relevant nanoheterogeneous environments. We conclude the use of the second-order cumulant expansion permits the use of faster pulling velocities in sMD simulations, without introducing bias due to large dispersion in the non-equilibrium work distribution.

Atomistic Analysis of ToxN and ToxI Complex Unbinding Mechanism

ToxIN is a triangular structure formed by three protein toxins (ToxNs) and three specific noncoding RNA antitoxins (ToxIs). To respond to stimuli, ToxI is preferentially degraded, releasing the ToxN. Thus, the dynamic character is essential in the normal function interactions between ToxN and ToxI. Here, equilibrated molecular dynamics (MD) simulations were performed to study the stability of ToxN and ToxI. The results indicate that ToxI adjusts the conformation of 3' and 5' termini to bind to ToxN. Steered molecular dynamics (SMD) simulations combined with the recently developed thermodynamic integration in 3nD (TI3nD) method were carried out to investigate ToxN unbinding from the ToxIN complex. The potentials of mean force (PMFs) and atomistic pictures suggest the unbinding mechanism as follows: (1) dissociation of the 5' terminus from ToxN, (2) missing the interactions involved in the 3' terminus of ToxI without three nucleotides (G31, A32, and A33), (3) starting to unfold for ToxI, (4) leaving the binding package of ToxN for three nucleotides of ToxI, (5) unfolding of ToxI. This work provides information on the structure-function relationship at the atomistic level, which is helpful for designing new potent antibacterial drugs in the future.

Exploring molecular mechanism of allosteric inhibitor to relieve drug resistance of multiple mutations in HIV-1 protease by enhanced conformational sampling

Recently, allosteric regulations of HIV-1 protease (PR) are suggested as a promising approach to relieve drug resistance of mutations toward inhibitors targeting the active site of PR. Replica-exchange molecular dynamics (REMD) simulations and normal mode analysis (NMA) are integrated to enhance conformational sampling of PR. Molecular mechanics generalized Born surface area (MM-GBSA) method was applied to calculate binding free energies of three inhibitors APV, DRV, and NIT to the wild-type (WT) and multidrug resistance (MDR) PRs. The results suggest that binding free energies of APV and DRV are decreased in the MDR PR relative to the WT PR, suggesting drug resistance of mutations on these two inhibitors. However, the binding ability of the allosteric inhibitor NIT is not impaired in the MDR PR. In addition, internal dynamics analysis based on REMD simulations proves that mutations hardly produce obvious effect on the conformation of the MDR PR in comparison to the WT PR. Scanning of hydrophobic contacts and hydrogen bond contacts of inhibitors with residues of PRs on the concatenated trajectories of REMD demonstrates that mutations change the symmetric interaction networks of APV and DRV with PR, but do not generate obvious influence on the asymmetric interaction network of NIT with PR. In summary, allosteric inhibitor NIT can adapt the MDR PR better than those inhibitors toward the active site of PR, thus allosteric inhibitors of PR may be a possible channel to overcome drug resistance of PR.

Dynamics revelation of conformational changes and binding modes of heat shock protein 90 induced by inhibitor associations

Heat shock protein 90 (Hsp90) has been an attractive target of potential drug design for antitumor treatment. The current work integrates molecular dynamics (MD) simulations, calculations of binding free energy, and principal component (PC) analysis with scanning of inhibitor-residue interaction to probe the binding modes of inhibitors YK9, YKJ and YKI to Hsp90 and identify the hot spot of the inhibitor-Hsp90 binding. The results suggest that the introductions of two groups G1 and G2 into YKJ and YKI strengthen the binding ability of YKJ and YKI to Hsp90 compared to YK9. PC analysis based MD trajectories prove that inhibitor bindings exert significant effects on the conformational changes, internal dynamics and motion modes of Hsp90, especially for the helix α2 and the loops L1 and L2. The calculations of residue-based free energy decomposition and scanning of the inhibitor-Hsp90 interaction suggest that six residues L107, G108, F138, Y139, W162 and F170 construct the common hot spot of the inhibitor-residue interactions. Moreover the substitutions of the groups G1 and G2 in YKJ and YKI lead to two additional hydrogen bonding interactions and multiple hydrophobic interactions for bindings of YKJ and YKI to Hsp90. This work is also expected to contribute theoretical hints for the design of potent inhibitors toward Hsp90.

Zinc ion-induced conformational changes in new Delphi metallo-β-lactamase 1 probed by molecular dynamics simulations and umbrella sampling

The hydrolysis of a β-lactam core ring caused by new Delphi metallo-β-lactamase 1 (NDM-1) with the help of two zinc cofactors induces significant resistance toward β-lactam antibiotics. Molecular dynamics (MD) simulations and the umbrella sampling method are integrated to study the conformational change mechanism of NDM-1 mediated by zinc ion binding. The statistical analyses of interaction contacts of the antibiotic ampicillin (AMP) with residues based on MD trajectories suggest that two Zn ions are essential for maintaining the binding of AMP with NDM-1. Umbrella sampling simulations further reveal that double-Zn coordination exerts strong restriction on the motions of loop L10 relative to loops L3 and L4. Principal component (PC) analysis also demonstrates that zinc ion binding totally inhibits the motion extent of NDM-1 and changes internal motion modes in NDM-1. We expect that the current study can provide significant dynamical information involving conformational changes of NDM-1 for the development of efficient inhibitors to decrease drug resistance of NDM-1 toward antibiotics.

Mutation L1196M-induced conformational changes and the drug resistant mechanism of anaplastic lymphoma kinase studied by free energy perturbation and umbrella sampling

Anaplastic lymphoma kinase (ALK) has been regarded as a promising drug target in the treatment of tumors and the mutation L1196M induces different levels of drug resistance toward the existing inhibitors.

Computing the binding affinity of a ligand buried deep inside a protein with the hybrid steered molecular dynamics

Computing the ligand-protein binding affinity (or the Gibbs free energy) with chemical accuracy has long been a challenge for which many methods/approaches have been developed and refined with various successful applications. False positives and, even more harmful, false negatives have been and still are a common occurrence in practical applications. Inevitable in all approaches are the errors in the force field parameters we obtain from quantum mechanical computation and/or empirical fittings for the intra- and inter-molecular interactions. These errors propagate to the final results of the computed binding affinities even if we were able to perfectly implement the statistical mechanics of all the processes relevant to a given problem. And they are actually amplified to various degrees even in the mature, sophisticated computational approaches. In particular, the free energy perturbation (alchemical) approaches amplify the errors in the force field parameters because they rely on extracting the small differences between similarly large numbers. In this paper, we develop a hybrid steered molecular dynamics (hSMD) approach to the difficult binding problems of a ligand buried deep inside a protein. Sampling the transition along a physical (not alchemical) dissociation path of opening up the binding cavity---pulling out the ligand---closing back the cavity, we can avoid the problem of error amplifications by not relying on small differences between similar numbers. We tested this new form of hSMD on retinol inside cellular retinol-binding protein 1 and three cases of a ligand (a benzylacetate, a 2-nitrothiophene, and a benzene) inside a T4 lysozyme L99A/M102Q(H) double mutant. In all cases, we obtained binding free energies in close agreement with the experimentally measured values. This indicates that the force field parameters we employed are accurate and that hSMD (a brute force, unsophisticated approach) is free from the problem of error amplification suffered by many sophisticated approaches in the literature.

Clarifying binding difference of ATP and ADP to extracellular signal-regulated kinase 2 by using molecular dynamics simulations

Extracellular signal-regulated kinase 2 is a promising target for designs and development of anticancer drugs. Molecular dynamics simulations and molecular mechanics Poisson-Boltzmann method were applied to study binding difference of ADP and ATP to extracellular signal-regulated kinase 2. The results prove that the binding ability of ATP to extracellular signal-regulated kinase 2 is stronger than that of ADP. Principal component analysis performed by using molecular dynamics trajectories suggests that binding of ADP and ATP to extracellular signal-regulated kinase 2 change motion directions of two helices α1 and α2. Residue-based free energy decomposition method was adopted to calculate contributions of separate residues to associations of ADP and ATP with extracellular signal-regulated kinase 2. The results show that ADP and ATP produce strong CH-π interactions with five residues Ile29, Val37, Ala50, Leu105, and Leu154. In addition, five hydrogen bonding interactions of ADP and ATP with residues Lys52, Gln103, Asp104, and Met106 also stabilize bindings of ADP and ATP to extracellular signal-regulated kinase 2. Overall, the CH-π interactions of ATP with five residues Ile29, Val37, Ala50, Leu105, and Leu154 are stronger than ADP. This study is expected to contribute a significant theoretical hint for designs of anticancer drugs targeting extracellular signal-regulated kinase 2.

Factors affecting the interactions between beta-lactoglobulin and fatty acids as revealed in molecular dynamics simulations

Beta-lactoglobulin (BLG), a bovine dairy protein, is a promiscuously interacting protein that can bind multiple hydrophobic ligands. Fatty acids (FAs), common hydrophobic molecules bound to BLG, are important sources of fuel for life because they yield large quantities of ATP when metabolized. The binding affinity increases with the length of the ligands, indicating the importance of the van der Waals (vdW) interactions between the hydrocarbon tail and the hydrophobic calyx of BLG. An exception to this rule is caprylic acid (OCA) which is two-carbon shorter but has a stronger binding affinity than capric acid. Theoretical calculations in the current literature are not accurate enough to shed light on the underlying physics of this exception. The computed affinity values are greater for longer fatty acids without respect for the caprylic exception and those values are generally several orders of magnitude away from the experimental data. In this work, we used hybrid steered molecular dynamics to accurately compute the binding free energies between BLG and the five saturated FAs of 8 to 16 carbon atoms. The computed binding free energies agree well with experimental data not only in rank but also in absolute values. We gained insights into the exceptional behavior of caprylic acid in the computed values of entropy and electrostatic interactions. We found that the electrostatic interaction between the carboxyl group of caprylic acid and the two amino groups of K60/69 in BLG is much stronger than the vdW force between the OCA's hydrophobic tail and the BLG calyx. This pulls OCA to the top of the beta barrel where it is easier to fluctuate, giving rise to greater entropy of OCA at the binding site.