Simulating Water Exchange to Buried Binding Sites

Reverse Mapping of Coarse Grained Polyglutamine Conformations from PRIME20 Sampling

An inverse coarse-graining protocol is presented for generating and validating atomistic structures of large (bio-) molecules from conformations obtained via a coarse-grained sampling method. Specifically, the protocol is implemented and tested based on the (coarse-grained) PRIME20 protein model (P20/SAMC), and the resulting all-atom conformations are simulated using conventional biomolecular force fields. The phase space sampling at the coarse-grained level is performed with a stochastical approximation Monte Carlo approach. The method is applied to a series of polypeptides, specifically dimers of polyglutamine with varying chain length in aqueous solution. The majority (>70 %) of the conformations obtained from the coarse-grained peptide model can successfully be mapped back to atomistic structures that remain conformationally stable during 10 ns of molecular dynamics simulations. This work can be seen as the first step towards the overarching goal of improving our understanding of protein aggregation phenomena through simulation methods.

Free Energy Density of a Fluid and Its Role in Solvation and Binding

The concept that a fluid has a position-dependent free energy density appears in the literature but has not been fully developed or accepted. We set this concept on an unambiguous theoretical footing via the following strategy. First, we set forth four desiderata that should be satisfied by any definition of the position-dependent free energy density, f(R), in a system comprising only a fluid and a rigid solute: its volume integral, plus the fixed internal energy of the solute, should be the system free energy; it deviates from its bulk value, fbulk, near a solute but should asymptotically approach fbulk with increasing distance from the solute; it should go to zero where the solvent density goes to zero; and it should be well-defined in the most general case of a fluid made up of flexible molecules with an arbitrary interaction potential. Second, we use statistical thermodynamics to formulate a definition of the free energy density that satisfies these desiderata. Third, we show how any free energy density satisfying the desiderata may be used to analyze molecular processes in solution. In particular, because the spatial integral of f(R) equals the free energy of the system, it can be used to compute free energy changes that result from the rearrangement of solutes as well as the forces exerted on the solutes by the solvent. This enables the use of a thermodynamic analysis of water in protein binding sites to inform ligand design. Finally, we discuss related literature and address published concerns regarding the thermodynamic plausibility of a position-dependent free energy density. The theory presented here has applications in theoretical and computational chemistry and may be further generalizable beyond fluids, such as to solids and macromolecules.

Resolving binding pathways and solvation thermodynamics of plant hormone receptors

Plant hormones are small molecules that regulate plant growth, development, and responses to biotic and abiotic stresses. They are specifically recognized by the binding site of their receptors. In this work, we resolved the binding pathways for eight classes of phytohormones (auxin, jasmonate, gibberellin, strigolactone, brassinosteroid, cytokinin, salicylic acid, and abscisic acid) to their canonical receptors using extensive molecular dynamics simulations. Furthermore, we investigated the role of water displacement and reorganization at the binding site of the plant receptors through inhomogeneous solvation theory. Our findings predict that displacement of water molecules by phytohormones contributes to free energy of binding via entropy gain and is associated with significant free energy barriers for most systems analyzed. Also, our results indicate that displacement of unfavorable water molecules in the binding site can be exploited in rational agrochemical design. Overall, this study uncovers the mechanism of ligand binding and the role of water molecules in plant hormone perception, which creates new avenues for agrochemical design to target plant growth and development.

Modern Alchemical Free Energy Methods for Drug Discovery Explained

This Perspective provides a contextual explanation of the current state-of-the-art alchemical free energy methods and their role in drug discovery as well as highlights select emerging technologies. The narrative attempts to answer basic questions about what goes on "under the hood" in free energy simulations and provide general guidelines for how to run simulations and analyze the results. It is the hope that this work will provide a valuable introduction to students and scientists in the field.

Comparison of Receptor-Ligand Restraint Schemes for Alchemical Absolute Binding Free Energy Calculations

Alchemical absolute binding free energy calculations are of increasing interest in drug discovery. These calculations require restraints between the receptor and ligand to restrict their relative positions and, optionally, orientations. Boresch restraints are commonly used, but they must be carefully selected in order to sufficiently restrain the ligand and to avoid inherent instabilities. Applying multiple distance restraints between anchor points in the receptor and ligand provides an alternative framework without inherent instabilities which may provide convergence benefits by more strongly restricting the relative movements of the receptor and ligand. However, there is no simple method to calculate the free energy of releasing these restraints due to the coupling of the internal and external degrees of freedom of the receptor and ligand. Here, a method to rigorously calculate free energies of binding with multiple distance restraints by imposing intramolecular restraints on the anchor points is proposed. Absolute binding free energies for the human macrophage migration inhibitory factor/MIF180, system obtained using a variety of Boresch restraints and rigorous and nonrigorous implementations of multiple distance restraints are compared. It is shown that several multiple distance restraint schemes produce estimates in good agreement with Boresch restraints. In contrast, calculations without orientational restraints produce erroneously favorable free energies of binding by up to approximately 4 kcal mol-1. These approaches offer new options for the deployment of alchemical absolute binding free energy calculations.

WaterKit: Thermodynamic Profiling of Protein Hydration Sites

Water desolvation is one of the key components of the free energy of binding of small molecules to their receptors. Thus, understanding the energetic balance of solvation and desolvation resulting from individual water molecules can be crucial when estimating ligand binding, especially when evaluating different molecules and poses as done in High-Throughput Virtual Screening (HTVS). Over the most recent decades, several methods were developed to tackle this problem, ranging from fast approximate methods (usually empirical functions using either discrete atom-atom pairwise interactions or continuum solvent models) to more computationally expensive and accurate ones, mostly based on Molecular Dynamics (MD) simulations, such as Grid Inhomogeneous Solvation Theory (GIST) or Double Decoupling. On one hand, MD-based methods are prohibitive to use in HTVS to estimate the role of waters on the fly for each ligand. On the other hand, fast and approximate methods show an unsatisfactory level of accuracy, with low agreement with results obtained with the more expensive methods. Here we introduce WaterKit, a new grid-based sampling method with explicit water molecules to calculate thermodynamic properties using the GIST method. Our results show that the discrete placement of water molecules is successful in reproducing the position of crystallographic waters with very high accuracy, as well as providing thermodynamic estimates with accuracy comparable to more expensive MD simulations. Unlike these methods, WaterKit can be used to analyze specific regions on the protein surface, (such as the binding site of a receptor), without having to hydrate and simulate the whole receptor structure. The results show the feasibility of a general and fast method to compute thermodynamic properties of water molecules, making it well-suited to be integrated in high-throughput pipelines such as molecular docking.

Accounting for Solvation Correlation Effects on the Thermodynamics of Water Networks in Protein Cavities

Macromolecular recognition and ligand binding are at the core of biological function and drug discovery efforts. Water molecules play a significant role in mediating the protein-ligand interaction, acting as more than just the surrounding medium by affecting the thermodynamics and thus the outcome of the binding process. As individual water contributions are impossible to measure experimentally, a range of computational methods have emerged to identify hydration sites in protein pockets and characterize their energetic contributions for drug discovery applications. Even though several methods model solvation effects explicitly, they focus on determining the stability of specific water sites independently and neglect solvation correlation effects upon replacement of clusters of water molecules, which typically happens in hit-to-lead optimization. In this work, we rigorously determine the conjoint effects of replacing all combinations of water molecules in protein binding pockets through the use of the RE-EDS multistate free-energy method, which combines Hamiltonian replica exchange (RE) and enveloping distribution sampling (EDS). Applications on the small bovine pancreatic trypsin inhibitor and four proteins of the bromodomain family illustrate the extent of solvation correlation effects on water thermodynamics, with the favorability of replacement of the water sites by pharmacophore probes highly dependent on the composition of the water network and the pocket environment. Given the ubiquity of water networks in biologically relevant protein targets, we believe our approach can be helpful for computer-aided drug discovery by providing a pocket-specific and a priori systematic consideration of solvation effects on ligand binding and selectivity.

Enhanced Grand Canonical Sampling of Occluded Water Sites Using Nonequilibrium Candidate Monte Carlo

Water molecules play a key role in many biomolecular systems, particularly when bound at protein-ligand interfaces. However, molecular simulation studies on such systems are hampered by the relatively long time scales over which water exchange between a protein and solvent takes place. Grand canonical Monte Carlo (GCMC) is a simulation technique that avoids this issue by attempting the insertion and deletion of water molecules within a given structure. The approach is constrained by low acceptance probabilities for insertions in congested systems, however. To address this issue, here, we combine GCMC with nonequilibium candidate Monte Carlo (NCMC) to yield a method that we refer to as grand canonical nonequilibrium candidate Monte Carlo (GCNCMC), in which the water insertions and deletions are carried out in a gradual, nonequilibrium fashion. We validate this new approach by comparing GCNCMC and GCMC simulations of bulk water and three protein binding sites. We find that not only is the efficiency of the water sampling improved by GCNCMC but that it also results in increased sampling of ligand conformations in a protein binding site, revealing new water-mediated ligand-binding geometries that are not observed using alternative enhanced sampling techniques.

Water Networks in Complexes between Proteins and FDA-Approved Drugs

Water molecules at protein-ligand interfaces are often of significant pharmaceutical interest, owing in part to the entropy which can be released upon the displacement of an ordered water by a therapeutic compound. Protein structures may not, however, completely resolve all critical bound water molecules, or there may be no experimental data available. As such, predicting the location of water molecules in the absence of a crystal structure is important in the context of rational drug design. Grand canonical Monte Carlo (GCMC) is a computational technique that is gaining popularity for the simulation of buried water sites. In this work, we assess the ability of GCMC to accurately predict water binding locations, using a dataset that we have curated, containing 108 unique structures of complexes between proteins and Food and Drug Administration (FDA)-approved small-molecule drugs. We show that GCMC correctly predicts 81.4% of nonbulk crystallographic water sites to within 1.4 Å. However, our analysis demonstrates that the reported performance of water prediction methods is highly sensitive to the way in which the performance is measured. We also find that crystallographic water sites with more protein/ligand hydrogen bonds and stronger electron density are more reliably predicted by GCMC. An analysis of water networks revealed that more than half of the structures contain at least one ligand-contacting water network. In these cases, displacement of a water site by a ligand modification might yield unexpected results if the larger network is destabilized. Cooperative effects between waters should therefore be explicitly considered in structure-based drug design.

Functional Water Networks in Fully Hydrated Photosystem II

Water channels and networks within photosystem II (PSII) of oxygenic photosynthesis are critical for enzyme structure and function. They control substrate delivery to the oxygen-evolving center and mediate proton transfer at both the oxidative and reductive endpoints. Current views on PSII hydration are derived from protein crystallography, but structural information may be compromised by sample dehydration and technical limitations. Here, we simulate the physiological hydration structure of a cyanobacterial PSII model following a thorough hydration procedure and large-scale unconstrained all-atom molecular dynamics enabled by massively parallel simulations. We show that crystallographic models of PSII are moderately to severely dehydrated and that this problem is particularly acute for models derived from X-ray free electron laser (XFEL) serial femtosecond crystallography. We present a fully hydrated representation of cyanobacterial PSII and map all water channels, both static and dynamic, associated with the electron donor and acceptor sides. Among them, we describe a series of transient channels and the attendant conformational gating role of protein components. On the acceptor side, we characterize a channel system that is absent from existing crystallographic models but is likely functionally important for the reduction of the terminal electron acceptor plastoquinone Q_B. The results of the present work build a foundation for properly (re)evaluating crystallographic models and for eliciting new insights into PSII structure and function.

Meta-Analysis Reveals That Absolute Binding Free-Energy Calculations Approach Chemical Accuracy

Systematic and quantitative analysis of the reliability of formally exact methods that in silico calculate absolute protein-ligand binding free energies remains lacking. Here, we provide, for the first time, evidence-based information on the reliability of these methods by statistically studying 853 cases from 34 different research groups through meta-analysis. The results show that formally exact methods approach chemical accuracy (error = 1.58 kcal/mol), even if people are challenging difficult tasks such as blind drug screening in recent years. The geometrical-pathway-based methods prove to possess a better convergence ability than the alchemical ones, while the latter have a larger application range. We also reveal the importance of always using the latest force fields to guarantee reliability and discuss the pros and cons of turning to an implicit solvent model in absolute binding free-energy calculations. Moreover, based on the meta-analysis, an evidence-based guideline for in silico binding free-energy calculations is provided.

Evaluating the use of absolute binding free energy in the fragment optimisation process

Key to the fragment optimisation process within drug design is the need to accurately capture the changes in affinity that are associated with a given set of chemical modifications. Due to the weakly binding nature of fragments, this has proven to be a challenging task, despite recent advancements in leveraging experimental and computational methods. In this work, we evaluate the use of Absolute Binding Free Energy (ABFE) calculations in guiding fragment optimisation decisions, retrospectively calculating binding free energies for 59 ligands across 4 fragment elaboration campaigns. We first demonstrate that ABFEs can be used to accurately rank fragment-sized binders with an overall Spearman's r of 0.89 and a Kendall τ of 0.67, although often deviating from experiment in absolute free energy values with an RMSE of 2.75 kcal/mol. We then also show that in several cases, retrospective fragment optimisation decisions can be supported by the ABFE calculations. Comparing against cheaper endpoint methods, namely N_wat-MM/GBSA, we find that ABFEs offer better ranking power and correlation metrics. Our results indicate that ABFE calculations can usefully guide fragment elaborations to maximise affinity.

py-MCMD: Python Software for Performing Hybrid Monte Carlo/Molecular Dynamics Simulations with GOMC and NAMD

py-MCMD, an open-source Python software, provides a robust workflow layer that manages communication of relevant system information between the simulation engines NAMD and GOMC and generates coherent thermodynamic properties and trajectories for analysis. To validate the workflow and highlight its capabilities, hybrid Monte Carlo/molecular dynamics (MC/MD) simulations are performed for SPC/E water in the isobaric-isothermal (NPT) and grand canonical (GC) ensembles as well as with Gibbs ensemble Monte Carlo (GEMC). The hybrid MC/MD approach shows close agreement with reference MC simulations and has a computational efficiency that is 2 to 136 times greater than traditional Monte Carlo simulations. MC/MD simulations performed for water in a graphene slit pore illustrate significant gains in sampling efficiency when the coupled-decoupled configurational-bias MC (CD-CBMC) algorithm is used compared with simulations using a single unbiased random trial position. Simulations using CD-CBMC reach equilibrium with 25 times fewer cycles than simulations using a single unbiased random trial position, with a small increase in computational cost. In a more challenging application, hybrid grand canonical Monte Carlo/molecular dynamics (GCMC/MD) simulations are used to hydrate a buried binding pocket in bovine pancreatic trypsin inhibitor. Water occupancies produced by GCMC/MD simulations are in close agreement with crystallographically identified positions, and GCMC/MD simulations have a computational efficiency that is 5 times better than MD simulations. py-MCMD is available on GitHub at https://github.com/GOMC-WSU/py-MCMD.

Comparison of Grand Canonical and Conventional Molecular Dynamics Simulation Methods for Protein-Bound Water Networks

Water molecules play important roles in all biochemical processes. Therefore, it is of key importance to obtain information of the structure, dynamics, and thermodynamics of water molecules around proteins. Numerous computational methods have been suggested with this aim. In this study, we compare the performance of conventional and grand-canonical Monte Carlo (GCMC) molecular dynamics (MD) simulations to sample the water structure, as well GCMC and grid-based inhomogeneous solvation theory (GIST) to describe the energetics of the water network. They are evaluated on two proteins: the buried ligand-binding site of a ferritin dimer and the solvent-exposed binding site of galectin-3. We show that GCMC/MD simulations significantly speed up the sampling and equilibration of water molecules in the buried binding site, thereby making the results more similar for simulations started from different states. Both GCMC/MD and conventional MD reproduce crystal-water molecules reasonably for the buried binding site. GIST analyses are normally based on restrained MD simulations. This improves the precision of the calculated energies, but the restraints also significantly affect both absolute and relative energies. Solvation free energies for individual water molecules calculated with and without restraints show a good correlation, but with large quantitative differences. Finally, we note that the solvation free energies calculated with GIST are ∼5 times larger than those estimated by GCMC owing to differences in the reference state.

Spatially Resolved Hydration Thermodynamics in Biomolecular Systems

Water is essential for the structure, dynamics, energetics, and thus the function of biomolecules. It is a formidable challenge to elicit, in microscopic detail, the role of the solvation-related driving forces of biomolecular processes, such as the enthalpy and entropy contributions to the underlying free-energy landscape. In this Perspective, we discuss recent developments and applications of computational methods that provide a spatially resolved map of hydration thermodynamics in biomolecular systems and thus yield atomic-level insights to guide the interpretation of experimental observations. An emphasis is on the challenge of quantifying the hydration entropy, which requires characterization of both the motions of the biomolecules and of the water molecules in their surrounding.

Accurate determination of protein:ligand standard binding free energies from molecular dynamics simulations

Designing a reliable computational methodology to calculate protein:ligand standard binding free energies is extremely challenging. The large change in configurational enthalpy and entropy that accompanies the association of ligand and protein is notoriously difficult to capture in naive brute-force simulations. Addressing this issue, the present protocol rests upon a rigorous statistical mechanical framework for the determination of protein:ligand binding affinities together with the comprehensive Binding Free-Energy Estimator 2 (BFEE2) application software. With the knowledge of the bound state, available from experiments or docking, application of the BFEE2 protocol with a reliable force field supplies in a matter of days standard binding free energies within chemical accuracy, for a broad range of protein:ligand complexes. Limiting undesirable human intervention, BFEE2 assists the end user in preparing all the necessary input files and performing the post-treatment of the simulations towards the final estimate of the binding affinity.

Fast Equilibration of Water between Buried Sites and the Bulk by Molecular Dynamics with Parallel Monte Carlo Water Moves on Graphical Processing Units

Molecular dynamics (MD) simulations of proteins are commonly used to sample from the Boltzmann distribution of conformational states, with wide-ranging applications spanning chemistry, biophysics, and drug discovery. However, MD can be inefficient at equilibrating water occupancy for buried cavities in proteins that are inaccessible to the surrounding solvent. Indeed, the time needed for water molecules to equilibrate between the bulk solvent and the binding site can be well beyond what is practical with standard MD, which typically ranges from hundreds of nanoseconds to a few microseconds. We recently introduced a hybrid Monte Carlo/MD (MC/MD) method, which speeds up the equilibration of water between buried cavities and the surrounding solvent, while sampling from the thermodynamically correct distribution of states. While the initial implementation of the MC functionality led to considerable slowing of the overall simulations, here we address this problem with a parallel MC algorithm implemented on graphical processing units. This results in speed-ups of 10-fold to 1000-fold over the original MC/MD algorithm, depending on the system and simulation parameters. The present method is available for use in the AMBER simulation software.

Water molecules at protein-drug interfaces: computational prediction and analysis methods

The fundamental importance of water molecules at drug-protein interfaces is now widely recognised and a significant feature in structure-based drug design. Experimental methods for analysing the role of water in drug binding have many challenges, including the accurate location of bound water molecules in crystal structures, and problems in resolving specific water contributions to binding thermodynamics. Computational analyses of binding site water molecules provide an alternative, and in principle complete, structural and thermodynamic picture, and their use is now commonplace in the pharmaceutical industry. In this review, we describe the computational methodologies that are available and discuss their strengths and weaknesses. Additionally, we provide a critical analysis of the experimental data used to validate the methods, regarding the type and quality of experimental structural data. We also discuss some of the fundamental difficulties of each method and suggest directions for future study.

Challenges Encountered Applying Equilibrium and Nonequilibrium Binding Free Energy Calculations

Binding free energy calculations have become increasingly valuable to drive decision making in drug discovery projects. However, among other issues, inadequate sampling can reduce accuracy, limiting the value of the technique. In this paper, we apply absolute binding free energy calculations to ligands binding to T4 lysozyme L99A and HSP90 using equilibrium and nonequilibrium approaches. We highlight sampling problems encountered in these systems, such as slow side chain rearrangements and slow changes of water placement upon ligand binding. These same types of challenges are also likely to show up in other protein-ligand systems, and we propose some strategies to diagnose and test for such problems in alchemical free energy calculations. We also explore similarities and differences in how the equilibrium and the nonequilibrium approaches handle these problems. Our results show the large amount of work still to be done to make free energy calculations robust and reliable and provide insight for future research in this area.

Allosteric communication in DNA polymerase clamp loaders relies on a critical hydrogen-bonded junction

Clamp loaders are AAA+ ATPases that load sliding clamps onto DNA. We mapped the mutational sensitivity of the T4 bacteriophage sliding clamp and clamp loader by deep mutagenesis, and found that residues not involved in catalysis or binding display remarkable tolerance to mutation. An exception is a glutamine residue in the AAA+ module (Gln 118) that is not located at a catalytic or interfacial site. Gln 118 forms a hydrogen-bonded junction in a helical unit that we term the central coupler, because it connects the catalytic centers to DNA and the sliding clamp. A suppressor mutation indicates that hydrogen bonding in the junction is important, and molecular dynamics simulations reveal that it maintains rigidity in the central coupler. The glutamine-mediated junction is preserved in diverse AAA+ ATPases, suggesting that a connected network of hydrogen bonds that links ATP molecules is an essential aspect of allosteric communication in these proteins.

Enhancing water sampling of buried binding sites using nonequilibrium candidate Monte Carlo

Water molecules can be found interacting with the surface and within cavities in proteins. However, water exchange between bulk and buried hydration sites can be slow compared to simulation timescales, thus leading to the inefficient sampling of the locations of water. This can pose problems for free energy calculations for computer-aided drug design. Here, we apply a hybrid method that combines nonequilibrium candidate Monte Carlo (NCMC) simulations and molecular dynamics (MD) to enhance sampling of water in specific areas of a system, such as the binding site of a protein. Our approach uses NCMC to gradually remove interactions between a selected water molecule and its environment, then translates the water to a new region, before turning the interactions back on. This approach of gradual removal of interactions, followed by a move and then reintroduction of interactions, allows the environment to relax in response to the proposed water translation, improving acceptance of moves and thereby accelerating water exchange and sampling. We validate this approach on several test systems including the ligand-bound MUP-1 and HSP90 proteins with buried crystallographic waters removed. We show that our BLUES (NCMC/MD) method enhances water sampling relative to normal MD when applied to these systems. Thus, this approach provides a strategy to improve water sampling in molecular simulations which may be useful in practical applications in drug discovery and biomolecular design.

Prediction of Drug Potencies of BACE1 Inhibitors: A Molecular Dynamics Simulation and MM_GB(PB)SA Scoring

Alzheimer’s disease (AD) is a progressive neurodegenerative brain disorder. One of the important therapeutic approaches of AD is the inhibition of β-site APP cleaving enzyme-1 (BACE1). This enzyme plays a central role in the synthesis of the pathogenic β-amyloid peptides (Aβ) in Alzheimer’s disease. A group of potent BACE1 inhibitors with known X-ray structures (PDB ID 5i3X, 5i3Y, 5iE1, 5i3V, 5i3W, 4LC7, 3TPP) were studied by molecular dynamics simulation and binding energy calculation employing MM_GB(PB)SA. The calculated binding energies gave Kd values of 0.139 µM, 1.39 nM, 4.39 mM, 24.3 nM, 1.39 mM, 29.13 mM, and 193.07 nM, respectively. These inhibitors showed potent inhibitory activities in enzymatic and cell assays. The Kd values are compared with experimental values and the structures are discussed in view of the energy contributions to binding. Drug likeness of these inhibitors is also discussed. Accommodation of ligands in the catalytic site of BACE1 is discussed depending on the type of fragment involved in each structure. Molecular dynamics (MD) simulations and energy studies were used to explore the recognition of the selected BACE1 inhibitors by Asp32, Asp228, and the hydrophobic flap. The results show that selective BACE1 inhibition may be due to the formation of strong electrostatic interactions with Asp32 and Asp228 and a large number of hydrogen bonds, in addition to π–π and van der Waals interactions with the amino acid residues located inside the catalytic cavity. Interactions with the ligands show a similar binding mode with BACE1. These results help to rationalize the design of selective BACE1 inhibitors.

Accounting for the Central Role of Interfacial Water in Protein-Ligand Binding Free Energy Calculations

Rigorous binding free energy methods in drug discovery are growing in popularity because of a combination of methodological advances, improvements in computer hardware, and workflow automation. These calculations typically use molecular dynamics (MD) to sample from the Boltzmann distribution of conformational states. However, when part or all of the binding sites is inaccessible to the bulk solvent, the time needed for water molecules to equilibrate between bulk solvent and the binding site can be well beyond what is practical with standard MD. This sampling limitation is problematic in relative binding free energy calculations, which compute the reversible work of converting ligand 1 to ligand 2 within the binding site. Thus, if ligand 1 is smaller and/or more polar than ligand 2, the perturbation may allow additional water molecules to occupy a region of the binding site. However, this change in hydration may not be captured by standard MD simulations and may therefore lead to errors in the computed free energy. We recently developed a hybrid Monte Carlo/MD (MC/MD) method, which speeds up the equilibration of water between bulk solvent and buried cavities, while sampling from the intended distribution of states. Here, we report on the use of this approach in the context of alchemical binding free energy calculations. We find that using MC/MD markedly improves the accuracy of the calculations and also reduces hysteresis between the forward and reverse perturbations, relative to matched calculations using only MD with or without the crystallographic water molecules. The present method is available for use in AMBER simulation software.

Accounting for the Central Role of Interfacial Water in Protein-Ligand Binding Free Energy Calculations

Rigorous binding free energy methods in drug discovery are growing in popularity because of a combination of methodological advances, improvements in computer hardware, and workflow automation. These calculations typically use molecular dynamics (MD) to sample from the Boltzmann distribution of conformational states. However, when part or all of the binding sites is inaccessible to the bulk solvent, the time needed for water molecules to equilibrate between bulk solvent and the binding site can be well beyond what is practical with standard MD. This sampling limitation is problematic in relative binding free energy calculations, which compute the reversible work of converting ligand 1 to ligand 2 within the binding site. Thus, if ligand 1 is smaller and/or more polar than ligand 2, the perturbation may allow additional water molecules to occupy a region of the binding site. However, this change in hydration may not be captured by standard MD simulations and may therefore lead to errors in the computed free energy. We recently developed a hybrid Monte Carlo/MD (MC/MD) method, which speeds up the equilibration of water between bulk solvent and buried cavities, while sampling from the intended distribution of states. Here, we report on the use of this approach in the context of alchemical binding free energy calculations. We find that using MC/MD markedly improves the accuracy of the calculations and also reduces hysteresis between the forward and reverse perturbations, relative to matched calculations using only MD with or without the crystallographic water molecules. The present method is available for use in AMBER simulation software.

Alchemical Binding Free Energy Calculations in AMBER20: Advances and Best Practices for Drug Discovery

Predicting protein-ligand binding affinities and the associated thermodynamics of biomolecular recognition is a primary objective of structure-based drug design. Alchemical free energy simulations offer a highly accurate and computationally efficient route to achieving this goal. While the AMBER molecular dynamics package has successfully been used for alchemical free energy simulations in academic research groups for decades, widespread impact in industrial drug discovery settings has been minimal because of the previous limitations within the AMBER alchemical code, coupled with challenges in system setup and postprocessing workflows. Through a close academia-industry collaboration we have addressed many of the previous limitations with an aim to improve accuracy, efficiency, and robustness of alchemical binding free energy simulations in industrial drug discovery applications. Here, we highlight some of the recent advances in AMBER20 with a focus on alchemical binding free energy (BFE) calculations, which are less computationally intensive than alternative binding free energy methods where full binding/unbinding paths are explored. In addition to scientific and technical advances in AMBER20, we also describe the essential practical aspects associated with running relative alchemical BFE calculations, along with recommendations for best practices, highlighting the importance not only of the alchemical simulation code but also the auxiliary functionalities and expertise required to obtain accurate and reliable results. This work is intended to provide a contemporary overview of the scientific, technical, and practical issues associated with running relative BFE simulations in AMBER20, with a focus on real-world drug discovery applications.

Protein Matrix Control of Reaction Center Excitation in Photosystem II

Photosystem II (PSII) is a multisubunit pigment-protein complex that uses light-induced charge separation to power oxygenic photosynthesis. Its reaction center chromophores, where the charge transfer cascade is initiated, are arranged symmetrically along the D1 and D2 core polypeptides and comprise four chlorophyll (PD1, PD2, ChlD1, ChlD2) and two pheophytin molecules (PheoD1 and PheoD2). Evolution favored productive electron transfer only via the D1 branch, with the precise nature of primary excitation and the factors that control asymmetric charge transfer remaining under investigation. Here we present a detailed atomistic description for both. We combine large-scale simulations of membrane-embedded PSII with high-level quantum-mechanics/molecular-mechanics (QM/MM) calculations of individual and coupled reaction center chromophores to describe reaction center excited states. We employ both range-separated time-dependent density functional theory and the recently developed domain based local pair natural orbital (DLPNO) implementation of the similarity transformed equation of motion coupled cluster theory with single and double excitations (STEOM-CCSD), the first coupled cluster QM/MM calculations of the reaction center. We find that the protein matrix is exclusively responsible for both transverse (chlorophylls versus pheophytins) and lateral (D1 versus D2 branch) excitation asymmetry, making ChlD1 the chromophore with the lowest site energy. Multipigment calculations show that the protein matrix renders the ChlD1 → PheoD1 charge-transfer the lowest energy excitation globally within the reaction center, lower than any pigment-centered local excitation. Remarkably, no low-energy charge transfer states are located within the "special pair" PD1-PD2, which is therefore excluded as the site of initial charge separation in PSII. Finally, molecular dynamics simulations suggest that modulation of the electrostatic environment due to protein conformational flexibility enables direct excitation of low-lying charge transfer states by far-red light.

Enhancing Water Sampling in Free Energy Calculations with Grand Canonical Monte Carlo

The prediction of protein-ligand binding affinities using free energy perturbation (FEP) is becoming increasingly routine in structure-based drug discovery. Most FEP packages use molecular dynamics (MD) to sample the configurations of proteins and ligands, as MD is well-suited to capturing coupled motion. However, MD can be prohibitively inefficient at sampling water molecules that are buried within binding sites, which has severely limited the domain of applicability of FEP and its prospective usage in drug discovery. In this paper, we present an advancement of FEP that augments MD with grand canonical Monte Carlo (GCMC), an enhanced sampling method, to overcome the problem of sampling water. We accomplished this without degrading computational performance. On both old and newly assembled data sets of protein-ligand complexes, we show that the use of GCMC in FEP is essential for accurate and robust predictions for ligand perturbations that disrupt buried water.

Water Distribution and Clustering on the Lyophilized IgG1 Surface: Insight from Molecular Dynamics Simulations

Water has a critical role in the stability of the higher-order structure of proteins. In addition, it is considered to be a major destabilization factor for the physical and chemical stability of freeze-dried proteins and peptides. Physical and chemical aspects of protein/water relationships are commonly studied with the use of water vapor sorption isotherms for amorphous lyophilized proteins, which, in turn, are commonly analyzed using the Brunauer-Emmett-Teller (BET) equation to obtain the parameters, W_m and C_B. The parameter W_m is generally referred to as the "monolayer limit of adsorption" and has a narrow range of 6-8% for most proteins. In this study, the water distribution on an IgG1 surface is investigated by molecular dynamics (MD) simulations at different water contents. The monolayer of water molecules was found to have limited coverage of the protein surface, and the true monolayer coverage of the protein globule actually occurs at a hydration level above 30%. The distribution of water molecules on the IgG1 surface is also highly heterogeneous, and the heterogeneity is not considered in the BET theory. In this study, a mechanistic model has been developed to describe the water vapor sorption isotherm. This model is based on the analysis of the hydrogen bonding network extracted from the MD simulations. The model is consistent with the experimental Type-II isotherm, which is usually observed for proteins. The physical meaning of the BET monolayer was redefined as the onset of water cluster formation. A simple model to calculate the onset water level, W_m, is proposed based on the hydration of different amino acids, as determined from the MD simulations.

Utilizing Grand Canonical Monte Carlo Methods in Drug Discovery

The concepts behind targeting waters for potency and selectivity gains have been well documented and explored, although maximizing such potential gains can prove to be challenging. This problem is exacerbated in cases where there are multiple interacting waters, wherein perturbation of one water can affect the free energy landscape of the remaining waters. Knowing the right modification a priori is challenging, and computational approaches are ideally suited to help answer the key question of which substitution is best to try. Here, we use Grand Canonical Monte Carlo and the recent Grand Canonical Alchemical Perturbation methods to both understand and predict the effect of ligand-mediated water displacement when more than one water molecule is involved, as well as to understand how exploiting water networks can help govern selectivity.

Dissecting the Energetics of Intrinsically Disordered Proteins via a Hybrid Experimental and Computational Approach

Intrinsically disordered proteins (IDPs) play important roles in biology, but little is known about the energetics of their inter-residue interactions. Methods that have been successfully applied to analyze the energetics of globular proteins are not applicable to the fluctuating partially ordered ensembles populated by IDPs. A combined computational experimental strategy is introduced for analyzing the energetic role of individual residues in the free state of IDPs. The approach combines experimental measurements of the binding of wild-type and mutant IDPs to their partners with alchemical free energy calculations of the structured complexes. These data allow quantitative information to be deduced about the free state via a thermodynamic cycle. The approach is validated by the analysis of the effects of mutations upon the binding free energy of the ovomucoid inhibitor third binding domain to its partners and is applied to the C-terminal domain of the measles virus nucleoprotein, a 125-residue IDP involved in the RNA transcription and replication of measles virus. The analysis reveals significant inter-residue interactions in the unbound IDP and suggests a biological role for them. The work demonstrates that advances in force fields and computational hardware have now led to the point where it is possible to develop methods, which integrate experimental and computational techniques to reveal insights that cannot be studied using either technique alone.