How Robust Is the Ligand Binding Transition State?

For many drug targets, it has been shown that the kinetics of drug binding (e.g., on rate and off rate) is more predictive of drug efficacy than thermodynamic quantities alone. This motivates the development of predictive computational models that can be used to optimize compounds on the basis of their kinetics. The structural details underpinning these computational models are found not only in the bound state but also in the short-lived ligand binding transition states. Although transition states cannot be directly observed experimentally due to their extremely short lifetimes, recent successes have demonstrated that modeling the ligand binding transition state is possible with the help of enhanced sampling molecular dynamics methods. Previously, we generated unbinding paths for an inhibitor of soluble epoxide hydrolase (sEH) with a residence time of 11 min. Here, we computationally modeled unbinding events with the weighted ensemble method REVO (resampling of ensembles by variation optimization) for five additional inhibitors of sEH with residence times ranging from 14.25 to 31.75 min, with average prediction accuracy within an order of magnitude. The unbinding ensembles are analyzed in detail, focusing on features of the ligand binding transition state ensembles (TSEs). We find that ligands with similar bound poses can show significant differences in their ligand binding TSEs, in terms of their spatial distribution and protein-ligand interactions. However, we also find similarities across the TSEs when examining more general features such as ligand degrees of freedom. Together these findings show significant challenges for rational, kinetics-based drug design.

Flexible Topology: A Dynamic Model of a Continuous Chemical Space.

Ligand design problems involve searching chemical space for a molecule with a set of desired properties. As chemical space is discrete, this search must be conducted in a pointwise manner, separately investigating one molecule at a time, which can be inefficient. We propose a method called "Flexible Topology", where a ligand is composed of a set of shapeshifting "ghost" atoms, whose atomic identities and connectivity can dynamically change over the course of a simulation. Ghost atoms are guided toward their target positions using a translation-, rotation-, and index-invariant restraint potential. This is the first step toward a continuous model of chemical space, where a dynamic simulation can move from one molecule to another by following gradients of a potential energy function. This builds on a substantial history of alchemy in the field of molecular dynamics simulation, including the Lambda dynamics method developed by Brooks and co-workers [X. Kong and C.L. Brooks III, J. Chem. Phys. 105, 2414 (1996)], but takes it to an extreme by associating a set of four dynamical attributes with each shapeshifting ghost atom that control not only its presence but also its atomic identity. Here, we outline the theoretical details of this method, its implementation using the OpenMM simulation package, and some preliminary studies of ghost particle assembly simulations in vacuum. We examine a set of 10 small molecules, ranging in size from 6 to 50 atoms, and show that Flexible Topology is able to consistently assemble all of these molecules to high accuracy, beginning from randomly initialized positions and attributes.

Support Vector Regression-Based Monte Carlo Simulation of Flexible Water Clusters.

Molecular simulations based on classical force fields are computationally efficient but lack accuracy due to the empirical formulation of non-bonded interactions. Quantum mechanical (QM) methods, albeit accurate, have inhibitory computational costs for large molecules and clusters. Hence, to overcome the bottleneck, machine learning (ML)-based methods have been employed in the recent years. We had earlier reported a combined scheme of many-body expansion (MBE) and ML to predict the interaction energies of rigid water clusters. In this work, we proceed toward building a flexible water model using the ML-MBE scheme. This ML-MBE scheme has an error of <1% for interaction energy prediction in comparison to the parent QM method for flexible water decamers. Machine learning-based Monte Carlo simulations (MLMC) are performed with this water model, and the structural properties of these configurations are compared with those obtained from ab initio molecular dynamics (AIMD) and the TIP3P classical force field. The radial distribution functions, tetrahedral order parameters, and number of hydrogen bonds in AIMD and MLMC have a similar qualitative and quantitative trend, whereas the classical force fields show a significant deviation.

Machine learning prediction of interaction energies in rigid water clusters.

Classical force fields form a computationally efficient avenue for calculating the energetics of large systems. However, due to the constraints of the underlying analytical form, it is sometimes not accurate enough. Quantum mechanical (QM) methods, although accurate, are computationally prohibitive for large systems. In order to circumvent the bottle-neck of interaction energy estimation of large systems, data driven approaches based on machine learning (ML) have been employed in recent years. In most of these studies, the method of choice is artificial neural networks (ANN). In this work, we have shown an alternative ML method, support vector regression (SVR), that provides comparable accuracy with better computational efficiency. We have further used many body expansion (MBE) along with SVR to predict interaction energies in water clusters (decamers). In the case of dimer and trimer interaction energies, the root mean square errors (RMSEs) of the SVR based scheme are 0.12 kcal mol-1 and 0.34 kcal mol-1, respectively. We show that the SVR and MBE based scheme has a RMSE of 2.78% in the estimation of decamer interaction energy against the parent QM method in a computationally efficient way.

Effect of solvation on the ionization of guanine nucleotide: A hybrid QM/EFP study.

Ionization of nucleobases is affected by their biological environment, which includes both the effect of adjacent nucleotides as well as the presence of water around it. Guanine and its nucleotide have the lowest ionization potentials among the various DNA bases. Therefore, the threshold of ionization is dependent on that of guanine and its characterization is crucial to the prediction of interaction of light with DNA. We investigate the effect of solvation on the vertical ionization energies (VIEs) of guanine and its nucleotide. In this work, we have used hybrid quantum mechanics/molecular mechanics (QM/MM) approach with effective fragment potential as the MM method of choice and equation-of-motion coupled-cluster for ionization potential with singles and doubles (EOM-IP-CCSD) as the QM method. The performance of the hybrid scheme with respect to the full QM method shows an accuracy of ≤ 0.02-0.04 eV. The lowest few ionizations of the nucleotide are found to be from different parts of the moiety, that is, the nucleic acid base, phosphate, or sugar, and these ionization energies are very closely spaced giving rise to a very complicated spectrum. Furthermore, microsolvation has large effects on these ionizations and can lead to red or blue shift depending on the position of the water molecule. Even a single water molecule can change the order of ionized states in the nucleotide. The VIEs of the bulk solvated chromophores are predicted and compared to existing experimental spectra. The predominant role of polarization in the solvatochromic shift is noticed. © 2017 Wiley Periodicals, Inc.

An interaction energy driven biased sampling technique: A faster route to ionization spectra in condensed phase.

We introduce a computationally efficient approach for calculating spectroscopic properties, such as ionization energies (IEs) in the condensed phase. Discrete quantum mechanical/molecular mechanical (QM/MM) approaches for spectroscopic properties in a dynamic system, such as aqueous solution, need a large sample space to obtain converged estimates, especially for the cases where particle (electron) number is not conserved, such as IEs or electron affinities (EAs). We devise a biased sampling technique based on an approximate estimate of interaction energy between the solute and solvent, that accelerates the convergence and therefore, reduces the computational cost significantly. The approximate interaction energy also provides a good measure of the spectral width of the chromophores in the condensed phase. This technique has been tested and benchmarked for (i) phenol, (ii) HBDI anion (hydroxybenzylidene dimethyl imidazolinone), and (iii) thymine in water. © 2017 Wiley Periodicals, Inc.

Electrostatic Origin of the Red Solvatochromic Shift of DFHBDI in RNA Spinach.

Interactions with the environment tune the spectral properties of biological chromophores, e.g., fluorescent proteins. Understanding the relative contribution of the various types of noncovalent interactions in the spectral shifts can provide rational design principles toward developing new fluorescent probes. In this work, we investigate the origin of the red shift in the absorption spectra of the difluoro hydroxybenzylidene dimethyl imidazolinone (DFHBDI) chromophore in RNA spinach as compared to the aqueous solution. We systematically decompose the effects of various components of interactions, namely, stacking, hydrogen bonding, and long-range electrostatics, in order to elucidate the relative role of these interactions in the observed spectral behavior. We find that the absorption peak of DFHBDI is red-shifted by ∼0.35 eV in RNA relative to the aqueous solution. Earlier proposals from Huang and co-workers have implicated the stacking interactions between DFHBDI and nucleic acid bases to be the driving force behind the observed red shift. In contrast, our findings reveal that the long-range electrostatic interactions between DFHBDI and negatively charged RNA make the most significant contribution. Moreover, we notice that the opposing electrostatic fields due to the RNA backbone and the polarized water molecules around the RNA give rise to the resultant red shift. Our results emphasize the effect of strong heterogeneity in the various environmental factors that might be competing with each other.

Effect of Solvation on Electron Detachment and Excitation Energies of a Green Fluorescent Protein Chromophore Variant.

Hybrid quantum mechanics/molecular mechanics (QM/MM) is applied to the fluorinated green fluorescent protein (GFP) chromophore (DFHBDI) in its deprotonated form to understand the solvatochromic shifts in its vertical detachment energy (VDE) and vertical excitation energy (VEE). This variant of the GFP chromophore becomes fluorescent in an RNA environment and has a wide range of applications in biomedical and biochemical fields. From microsolvation studies, we benchmark (with respect to full QM) the accuracy of our QM/MM calculations with effective fragment potential (EFP) as the MM method of choice. We show that while the solvatochromic shift in the VEE is minimal (0.1 eV blue shift) and its polarization component is only 0.03 eV, the effect of the solvent on the VDE is quite large (3.85 eV). We also show by accurate calculations on the solvatochromic shift of the VDE that polarization accounts for ∼0.23 eV and therefore cannot be neglected. The effect of the counterions on the VDE of the deprotonated chromophore in solvation is studied in detail, and a charge-smearing scheme is suggested for charged chromophores.