Unsupervised and supervised AI on molecular dynamics simulations reveals complex characteristics of HLA-A2-peptide immunogenicity.

Immunologic recognition of peptide antigens bound to class I major histocompatibility complex (MHC) molecules is essential to both novel immunotherapeutic development and human health at large. Current methods for predicting antigen peptide immunogenicity rely primarily on simple sequence representations, which allow for some understanding of immunogenic features but provide inadequate consideration of the full scale of molecular mechanisms tied to peptide recognition. We here characterize contributions that unsupervised and supervised artificial intelligence (AI) methods can make toward understanding and predicting MHC(HLA-A2)-peptide complex immunogenicity when applied to large ensembles of molecular dynamics simulations. We first show that an unsupervised AI method allows us to identify subtle features that drive immunogenicity differences between a cancer neoantigen and its wild-type peptide counterpart. Next, we demonstrate that a supervised AI method for class I MHC(HLA-A2)-peptide complex classification significantly outperforms a sequence model on small datasets corrected for trivial sequence correlations. Furthermore, we show that both unsupervised and supervised approaches reveal determinants of immunogenicity based on time-dependent molecular fluctuations and anchor position dynamics outside the MHC binding groove. We discuss implications of these structural and dynamic immunogenicity correlates for the induction of T cell responses and therapeutic T cell receptor design.

Analysis of Training and Seed Bias in Small Molecules Generated with a Conditional Graph-Based Variational Autoencoder─Insights for Practical AI-Driven Molecule Generation.

The application of deep learning to generative molecule design has shown early promise for accelerating lead series development. However, questions remain concerning how factors like training, data set, and seed bias impact the technology's utility to medicinal and computational chemists. In this work, we analyze the impact of seed and training bias on the output of an activity-conditioned graph-based variational autoencoder (VAE). Leveraging a massive, labeled data set corresponding to the dopamine D2 receptor, our graph-based generative model is shown to excel in producing desired conditioned activities and favorable unconditioned physical properties in generated molecules. We implement an activity-swapping method that allows for the activation, deactivation, or retention of activity of molecular seeds, and we apply independent deep learning classifiers to verify the generative results. Overall, we uncover relationships between noise, molecular seeds, and training set selection across a range of latent-space sampling procedures, providing important insights for practical AI-driven molecule generation.

Roughness of Molecular Property Landscapes and Its Impact on Modellability.

In molecular discovery and drug design, structure-property relationships and activity landscapes are often qualitatively or quantitatively analyzed to guide the navigation of chemical space. The roughness (or smoothness) of these molecular property landscapes is one of their most studied geometric attributes, as it can characterize the presence of activity cliffs, with rougher landscapes generally expected to pose tougher optimization challenges. Here, we introduce a general, quantitative measure for describing the roughness of molecular property landscapes. The proposed roughness index (ROGI) is loosely inspired by the concept of fractal dimension and strongly correlates with the out-of-sample error achieved by machine learning models on numerous regression tasks.

Self-Focusing Virtual Screening with Active Design Space Pruning.

High-throughput virtual screening is an indispensable technique utilized in the discovery of small molecules. In cases where the library of molecules is exceedingly large, the cost of an exhaustive virtual screen may be prohibitive. Model-guided optimization has been employed to lower these costs through dramatic increases in sample efficiency compared to random selection. However, these techniques introduce new costs to the workflow through the surrogate model training and inference steps. In this study, we propose an extension to the framework of model-guided optimization that mitigates inference costs using a technique we refer to as design space pruning (DSP), which irreversibly removes poor-performing candidates from consideration. We study the application of DSP to a variety of optimization tasks and observe significant reductions in overhead costs while exhibiting similar performance to the baseline optimization. DSP represents an attractive extension of model-guided optimization that can limit overhead costs in optimization settings where these costs are non-negligible relative to objective costs, such as docking.

Simplified, interpretable graph convolutional neural networks for small molecule activity prediction.

We here present a streamlined, explainable graph convolutional neural network (gCNN) architecture for small molecule activity prediction. We first conduct a hyperparameter optimization across nearly 800 protein targets that produces a simplified gCNN QSAR architecture, and we observe that such a model can yield performance improvements over both standard gCNN and RF methods on difficult-to-classify test sets. Additionally, we discuss how reductions in convolutional layer dimensions potentially speak to the "anatomical" needs of gCNNs with respect to radial coarse graining of molecular substructure. We augment this simplified architecture with saliency map technology that highlights molecular substructures relevant to activity, and we perform saliency analysis on nearly 100 data-rich protein targets. We show that resultant substructural clusters are useful visualization tools for understanding substructure-activity relationships. We go on to highlight connections between our models' saliency predictions and observations made in the medicinal chemistry literature, focusing on four case studies of past lead finding and lead optimization campaigns.

Combining Docking Pose Rank and Structure with Deep Learning Improves Protein-Ligand Binding Mode Prediction over a Baseline Docking Approach.

We present a simple, modular graph-based convolutional neural network that takes structural information from protein-ligand complexes as input to generate models for activity and binding mode prediction. Complex structures are generated by a standard docking procedure and fed into a dual-graph architecture that includes separate subnetworks for the ligand bonded topology and the ligand-protein contact map. Recent work has indicated that data set bias drives many past promising results derived from combining deep learning and docking. Our dual-graph network allows contributions from ligand identity that give rise to such biases to be distinguished from effects of protein-ligand interactions on classification. We show that our neural network is capable of learning from protein structural information when, as in the case of binding mode prediction, an unbiased data set is constructed. We next develop a deep learning model for binding mode prediction that uses docking ranking as input in combination with docking structures. This strategy mirrors past consensus models and outperforms a baseline docking program (AutoDock Vina) in a variety of tests, including on cross-docking data sets that mimic real-world docking use cases. Furthermore, the magnitudes of network predictions serve as reliable measures of model confidence.

Role of water and steric constraints in the kinetics of cavity-ligand unbinding.

A key factor influencing a drug's efficacy is its residence time in the binding pocket of the host protein. Using atomistic computer simulation to predict this residence time and the associated dissociation process is a desirable but extremely difficult task due to the long timescales involved. This gets further complicated by the presence of biophysical factors such as steric and solvation effects. In this work, we perform molecular dynamics (MD) simulations of the unbinding of a popular prototypical hydrophobic cavity-ligand system using a metadynamics-based approach that allows direct assessment of kinetic pathways and parameters. When constrained to move in an axial manner, the unbinding time is found to be on the order of 4,000 s. In accordance with previous studies, we find that the cavity must pass through a region of sharp wetting transition manifested by sudden and high fluctuations in solvent density. When we remove the steric constraints on ligand, the unbinding happens predominantly by an alternate pathway, where the unbinding becomes 20 times faster, and the sharp wetting transition instead becomes continuous. We validate the unbinding timescales from metadynamics through a Poisson analysis, and by comparison through detailed balance to binding timescale estimates from unbiased MD. This work demonstrates that enhanced sampling can be used to perform explicit solvent MD studies at timescales previously unattainable, to our knowledge, obtaining direct and reliable pictures of the underlying physiochemical factors including free energies and rate constants.

How hydrophobic drying forces impact the kinetics of molecular recognition.

A model of protein-ligand binding kinetics, in which slow solvent dynamics results from hydrophobic drying transitions, is investigated. Molecular dynamics simulations show that solvent in the receptor pocket can fluctuate between wet and dry states with lifetimes in each state that are long enough for the extraction of a separable potential of mean force and wet-to-dry transitions. We present a diffusive surface hopping model that is represented by a 2D Markovian master equation. One dimension is the standard reaction coordinate, the ligand-pocket separation, and the other is the solvent state in the region between ligand and binding pocket which specifies whether it is wet or dry. In our model, the ligand diffuses on a dynamic free-energy surface which undergoes kinetic transitions between the wet and dry states. The model yields good agreement with results from explicit solvent molecular dynamics simulation and an improved description of the kinetics of hydrophobic assembly. Furthermore, it is consistent with a "non-Markovian Brownian theory" for the ligand-pocket separation coordinate alone.

Theory and simulations of quantum glass forming liquids.

A comprehensive microscopic dynamical theory is presented for the description of quantum fluids as they transform into glasses. The theory is based on a quantum extension of mode-coupling theory. Novel effects are predicted, such as reentrant behavior of dynamical relaxation times. These predictions are supported by path integral ring polymer molecular dynamics simulations. The simulations provide detailed insight into the factors that govern slow dynamics in glassy quantum fluids. Connection to other recent work on both quantum glasses as well as quantum optimization problems is presented.

Efficient multiple time scale molecular dynamics: Using colored noise thermostats to stabilize resonances.

Multiple time scale molecular dynamics enhances computational efficiency by updating slow motions less frequently than fast motions. However, in practice, the largest outer time step possible is limited not by the physical forces but by resonances between the fast and slow modes. In this paper we show that this problem can be alleviated by using a simple colored noise thermostatting scheme which selectively targets the high frequency modes in the system. For two sample problems, flexible water and solvated alanine dipeptide, we demonstrate that this allows the use of large outer time steps while still obtaining accurate sampling and minimizing the perturbation of the dynamics. Furthermore, this approach is shown to be comparable to constraining fast motions, thus providing an alternative to molecular dynamics with constraints.

Displaced path integral formulation for the momentum distribution of quantum particles.

The proton momentum distribution, accessible by deep inelastic neutron scattering, is a very sensitive probe of the potential of mean force experienced by the protons in hydrogen-bonded systems. In this work we introduce a novel estimator for the end-to-end distribution of the Feynman paths, i.e., the Fourier transform of the momentum distribution. In this formulation, free particle and environmental contributions factorize. Moreover, the environmental contribution has a natural analogy to a free energy surface in statistical mechanics, facilitating the interpretation of experiments. The new formulation is not only conceptually but also computationally advantageous. We illustrate the method with applications to an empirical water model, ab initio ice, and one dimensional model systems.

Tunneling and delocalization effects in hydrogen bonded systems: a study in position and momentum space.

Novel experimental and computational studies have uncovered the proton momentum distribution in hydrogen bonded systems. In this work, we utilize recently developed open path integral Car-Parrinello molecular dynamics methodology in order to study the momentum distribution in phases of high pressure ice. Some of these phases exhibit symmetric hydrogen bonds and quantum tunneling. We find that the symmetric hydrogen bonded phase possesses a narrowed momentum distribution as compared with a covalently bonded phase, in agreement with recent experimental findings. The signatures of tunneling that we observe are a narrowed distribution in the low-to-intermediate momentum region, with a tail that extends to match the result of the covalently bonded state. The transition to tunneling behavior shows similarity to features observed in recent experiments performed on confined water. We corroborate our ice simulations with a study of a particle in a model one-dimensional double well potential that mimics some of the effects observed in bulk simulations. The temperature dependence of the momentum distribution in the one-dimensional model allows for the differentiation between ground state and mixed state tunneling effects.

Binding Specificity Determines the Cytochrome P450 3A4 Mediated Enantioselective Metabolism of Metconazole.

Cytochrome P450 3A4 (CYP3A4) is a promiscuous enzyme, mediating the biotransformations of ∼50% of clinically used drugs, many of which are chiral molecules. Probing the interactions between CYP3A4 and chiral chemicals is thus essential for the elucidation of molecular mechanisms of enantioselective metabolism. We developed a stepwise-restrained-molecular-dynamics (MD) method to model human CYP3A4 in a complex with cis-metconazole (MEZ) isomers and performed conventional MD simulations with a total simulation time of 2.2 µs to probe the molecular interactions. Our current study, which employs a combined experimental and theoretical approach, reports for the first time on the distinct conformational changes of CYP3A4 that are induced by the enantioselective binding of cis-MEZ enantiomers. CYP3A4 preferably metabolizes cis-RS MEZ over the cis-SR isomer, with the resultant enantiomer fraction for cis-MEZ increasing rapidly from 0.5 to 0.82. cis-RS MEZ adopts a more extended structure in the active pocket with its Cl atom exposed to the solvent, whereas cis-SR MEZ sits within the hydrophobic core of the active pocket. Free-energy-perturbation calculations indicate that unfavorable van der Waals interactions between the cis-MEZ isomers and the CYP3A4 binding pocket predominantly contribute to their binding-affinity differences. These results demonstrate that binding specificity determines the cytochrome P450 3A4 mediated enantioselective metabolism of cis-MEZ.

Accelerating physical simulations of proteins by leveraging external knowledge.

It is challenging to compute structure-function relationships of proteins using molecular physics. The problem arises from the exponential scaling of the computational searching and sampling of large conformational spaces. This scaling challenge is not met by today's methods, such as Monte Carlo, simulated annealing, genetic algorithms, or molecular dynamics (MD) or its variants such as replica exchange. Such methods of searching for optimal states on complex probabalistic landscapes are referred to more broadly as Explore-and-Exploit (EE), including in contexts such as computational learning, games, industrial planning and modeling military strategies. Here we describe a Bayesian method, called MELD, that 'melds' together explore-and-exploit approaches with externally added information that can be vague, combinatoric, noisy, intuitive, heuristic, or from experimental data. MELD is shown to accelerate physical MD simulations when using experimental data to determine protein structures; for predicting protein structures by using heuristic directives; and when predicting binding affinities of proteins from limited information about the binding site. Such Guided Explore-and-Exploit approaches might also be useful beyond proteins and beyond molecular science.

Computed Binding of Peptides to Proteins with MELD-Accelerated Molecular Dynamics.

It has been a challenge to compute the poses and affinities for binding of peptides to proteins by molecular dynamics (MD) simulations. Such computations would be valuable for capturing the physics and the conformational freedom of the molecules, but they are currently too computationally expensive. Here we describe using MELD (Modeling Employing Limited Data)-accelerated MD for finding the binding poses and approximate relative binding free energies for flexible-peptide/protein interactions. MELD uses only weak information about the binding motif and not the detailed binding mode that is typically required by other free-energy-based methods. We apply this technique to study binding of P53-derived peptides to MDM2 and MDMX. We find that MELD finds correct poses, that the binding induces the peptide into the correct helical conformation, and that it is capable of roughly estimating relative binding affinities. This method may be useful in peptide drug discovery.

Molecular Simulations Identify Binding Poses and Approximate Affinities of Stapled α-Helical Peptides to MDM2 and MDMX.

Traditionally, computing the binding affinities of proteins to even relatively small and rigid ligands by free-energy methods has been challenging due to large computational costs and significant errors. Here, we apply a new molecular simulation acceleration method called MELD (Modeling by Employing Limited Data) to study the binding of stapled α-helical peptides to the MDM2 and MDMX proteins. We employ free-energy-based molecular dynamics simulations (MELD-MD) to identify binding poses and calculate binding affinities. Even though stapled peptides are larger and more complex than most protein ligands, the MELD-MD simulations can identify relevant binding poses and compute relative binding affinities. MELD-MD appears to be a promising method for computing the binding properties of peptide ligands with proteins.

Ab initio molecular dynamics simulation of the structure and proton transport dynamics of methanol-water solutions.

Ab initio molecular dynamics simulations are employed to study the structural and proton transport properties of methanol-water mixtures. Structural characteristics analyzed at two different methanol mole fractions (X(M) = 0.25 and X(M) = 0.5) reveal enhanced structuring of water as the methanol mole fraction increases in agreement with recent neutron diffraction experiments. The simulations reveal the existence of separate hydrogen-bonded water and methanol networks, also in agreement with the neutron diffraction data. The addition of a single proton to the X(M) = 0.5 mixture leads to an anomalous structural or Grotthuss-type diffusion mechanism of the charge defect in which water-to-water, methanol-to-water, and water-to-methanol proton transfer reactions play the dominant role with methanol-to-methanol transfers being much less significant. Unlike in bulk water, where coordination number fluctuations drive the proton transport process, suppression of the coordination number of waters in the first solvation shell of the defect diminish the importance of coordination number fluctuations as a driving force in the structural diffusion process. The charge defect is found to reside preferentially at the interface between water and methanol networks. The length of the ab initio molecular dynamics run (approximately 120 ps), allowed the diffusion constant of the charge defect to be computed, yielding a value of D = 4.2 x 10(-5) cm2/s when deuterium masses are assigned to all protons in the system. The relation of this value to excess proton diffusion in bulk water is discussed. Finally, a kinetic theory is introduced to identify the relevant time scales in the proton transfer/transport process.

Advances in free-energy-based simulations of protein folding and ligand binding.

Free-energy-based simulations are increasingly providing the narratives about the structures, dynamics and biological mechanisms that constitute the fabric of protein science. Here, we review two recent successes. It is becoming practical: first, to fold small proteins with free-energy methods without knowing substructures and second, to compute ligand-protein binding affinities, not just their binding poses. Over the past 40 years, the timescales that can be simulated by atomistic MD are doubling every 1.3 years--which is faster than Moore's law. Thus, these advances are not simply due to the availability of faster computers. Force fields, solvation models and simulation methodology have kept pace with computing advancements, and are now quite good. At the tip of the spear recently are GPU-based computing, improved fast-solvation methods, continued advances in force fields, and conformational sampling methods that harness external information.

Blind protein structure prediction using accelerated free-energy simulations.

We report a key proof of principle of a new acceleration method [Modeling Employing Limited Data (MELD)] for predicting protein structures by molecular dynamics simulation. It shows that such Boltzmann-satisfying techniques are now sufficiently fast and accurate to predict native protein structures in a limited test within the Critical Assessment of Structure Prediction (CASP) community-wide blind competition.

Are hydrodynamic interactions important in the kinetics of hydrophobic collapse?

We study the kinetics of assembly of two plates of varying hydrophobicity, including cases where drying occurs and water strongly solvates the plate surfaces. The potential of mean force and molecular-scale hydrodynamics are computed from molecular dynamics simulations in explicit solvent as a function of particle separation. In agreement with our recent work on nanospheres [J. Phys. Chem. B 2012,116, 378-389], regions of high friction are found to be engendered by large and slow solvent fluctuations. These slow fluctuations can be due to either drying or confinement. The mean first passage times for assembly are computed by means of molecular dynamics simulations in explicit solvent and by Brownian dynamics simulations along the reaction path. Brownian dynamics makes use of the potential of mean force and hydrodynamic profile that we determined. Surprisingly, we find reasonable agreement between full-scale molecular dynamics and Brownian dynamics, despite the role of slow solvent relaxation in the assembly process. We found that molecular-scale hydrodynamic interactions are essential in describing the kinetics of assembly.