Spectral gap optimization of order parameters for sampling complex molecular systems.

In modern-day simulations of many-body systems, much of the computational complexity is shifted to the identification of slowly changing molecular order parameters called collective variables (CVs) or reaction coordinates. A vast array of enhanced-sampling methods are based on the identification and biasing of these low-dimensional order parameters, whose fluctuations are important in driving rare events of interest. Here, we describe a new algorithm for finding optimal low-dimensional CVs for use in enhanced-sampling biasing methods like umbrella sampling, metadynamics, and related methods, when limited prior static and dynamic information is known about the system, and a much larger set of candidate CVs is specified. The algorithm involves estimating the best combination of these candidate CVs, as quantified by a maximum path entropy estimate of the spectral gap for dynamics viewed as a function of that CV. The algorithm is called spectral gap optimization of order parameters (SGOOP). Through multiple practical examples, we show how this postprocessing procedure can lead to optimization of CV and several orders of magnitude improvement in the convergence of the free energy calculated through metadynamics, essentially giving the ability to extract useful information even from unsuccessful metadynamics runs.

Advancing Drug Discovery through Enhanced Free Energy Calculations.

A principal goal of drug discovery project is to design molecules that can tightly and selectively bind to the target protein receptor. Accurate prediction of protein-ligand binding free energies is therefore of central importance in computational chemistry and computer aided drug design. Multiple recent improvements in computing power, classical force field accuracy, enhanced sampling methods, and simulation setup have enabled accurate and reliable calculations of protein-ligands binding free energies, and position free energy calculations to play a guiding role in small molecule drug discovery. In this Account, we outline the relevant methodological advances, including the REST2 (Replica Exchange with Solute Temperting) enhanced sampling, the incorporation of REST2 sampling with convential FEP (Free Energy Perturbation) through FEP/REST, the OPLS3 force field, and the advanced simulation setup that constitute our FEP+ approach, followed by the presentation of extensive comparisons with experiment, demonstrating sufficient accuracy in potency prediction (better than 1 kcal/mol) to substantially impact lead optimization campaigns. The limitations of the current FEP+ implementation and best practices in drug discovery applications are also discussed followed by the future methodology development plans to address those limitations. We then report results from a recent drug discovery project, in which several thousand FEP+ calculations were successfully deployed to simultaneously optimize potency, selectivity, and solubility, illustrating the power of the approach to solve challenging drug design problems. The capabilities of free energy calculations to accurately predict potency and selectivity have led to the advance of ongoing drug discovery projects, in challenging situations where alternative approaches would have great difficulties. The ability to effectively carry out projects evaluating tens of thousands, or hundreds of thousands, of proposed drug candidates, is potentially transformative in enabling hard to drug targets to be attacked, and in facilitating the development of superior compounds, in various dimensions, for a wide range of targets. More effective integration of FEP+ calculations into the drug discovery process will ensure that the results are deployed in an optimal fashion for yielding the best possible compounds entering the clinic; this is where the greatest payoff is in the exploitation of computer driven design capabilities. A key conclusion from the work described is the surprisingly robust and accurate results that are attainable within the conventional classical simulation, fixed charge paradigm. No doubt there are individual cases that would benefit from a more sophisticated energy model or dynamical treatment, and properties other than protein-ligand binding energies may be more sensitive to these approximations. We conclude that an inflection point in the ability of MD simulations to impact drug discovery has now been attained, due to the confluence of hardware and software development along with the formulation of "good enough" theoretical methods and models.

Efficient sampling of puckering states of monosaccharides through replica exchange with solute tempering and bond softening.

A molecular-level understanding of the structure, dynamics, and reactivity of carbohydrates is fundamental to the understanding of a range of key biological processes. The six-membered pyranose ring, a central component of biological monosaccharides and carbohydrates, has many different puckering conformations, and the conformational free energy landscape of these biologically important monosaccharides remains elusive. The puckering conformations of monosaccharides are separated by high energy barriers, which pose a great challenge for the complete sampling of these important conformations and accurate modeling of these systems. While metadynamics or umbrella sampling methods have been used to study the conformational space of monosaccharides, these methods might be difficult to generalize to other complex ring systems with more degrees of freedom. In this paper, we introduce a new enhanced sampling method for the rapid sampling over high energy barriers that combines our previously developed enhanced sampling method REST (replica exchange with solute tempering) with a bond softening (BOS) scheme that makes a chemical bond in the ring weaker as one ascends the replica ladder. We call this new method replica exchange with solute tempering and bond softening (REST/BOS). We demonstrate the superior sampling efficiency of REST/BOS over other commonly used enhanced sampling methods, including temperature replica exchange method and REST. The conformational free energy landscape of four biologically important monosaccharides, namely, α-glucose, ß-glucose, ß-mannose, and ß-xylose, is studied using REST/BOS, and results are compared with previous experimental and theoretical studies.

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.

Predicting reaction coordinates in energy landscapes with diffusion anisotropy.

We consider a range of model potentials with metastable states undergoing molecular dynamics coupled to a thermal bath in the high friction regime and consider how the optimal reaction coordinate depends on the diffusion anisotropy. For this we use our recently proposed method "spectral gap optimization of order parameters (SGOOP)" [P. Tiwary and B. J. Berne, Proc. Natl. Acad. Sci. U. S. A. 113, 2839 (2016)]. We show how available information about dynamical observables in addition to static information can be incorporated into SGOOP, which can then be used to accurately determine the "best" reaction coordinate for arbitrary anisotropies. We compare our results with transmission coefficient calculations and published benchmarks wherever applicable or available, respectively.

How a Kinase Inhibitor Withstands Gatekeeper Residue Mutations.

Mutations in the gatekeeper residue of kinases have emerged as a key way through which cancer cells develop resistance to treatment. As such, the design of gatekeeper mutation resistant kinase inhibitors is a crucial way forward in increasing the efficacy of a broad range of anticancer drugs. In this work we use atomistic simulations to provide detailed thermodynamic and structural insight into how two inhibitors of cSrc kinase, namely, a commercial drug and type I kinase inhibitor Dasatinib and the type II inhibitor RL45, respectively fail and succeed in being effective against the T338M gatekeeper residue mutation in the kinase binding site. Given the well-known limitations of atomistic simulations in sampling biomolecular systems, we use an enhanced sampling technique called free energy perturbation with replica exchange solute tempering (FEP/REST). Our calculations find that the type I inhibitor Dasatinib binds favorably to the wild type but unfavorably to T338M mutated kinase, while RL45 binds favorably to both. The predicted relative binding free energies are well within 1 kcal/mol accuracy compared to experiments. We find that Dasatinib's impotency against gatekeeper residue mutations arises from a loss of ligand-kinase hydrogen bonding due to T338M mutation and from steric hindrance due to the presence of an inflexible phenyl ring close to the ligand. On the other hand, in the type II binding RL45, the central phenyl ring has very pronounced flexibility. This leads to the inhibitor overcoming effects of steric clashes on mutation and maintaining an electrostatically favorable "edge-to-face" orientation with a neighboring phenylalanine residue. Our work provides useful insight into the mechanisms of mutation resistant kinase inhibitors and demonstrates the usefulness of enhanced sampling techniques in computational drug design.

Kramers turnover: From energy diffusion to spatial diffusion using metadynamics.

We consider the rate of transition for a particle between two metastable states coupled to a thermal environment for various magnitudes of the coupling strength using the recently proposed infrequent metadynamics approach [P. Tiwary and M. Parrinello, Phys. Rev. Lett. 111, 230602 (2013)]. We are interested in understanding how this approach for obtaining rate constants performs as the dynamics regime changes from energy diffusion to spatial diffusion. Reassuringly, we find that the approach works remarkably well for various coupling strengths in the strong coupling regime, and to some extent even in the weak coupling regime.

How wet should be the reaction coordinate for ligand unbinding?

We use a recently proposed method called Spectral Gap Optimization of Order Parameters (SGOOP) [P. Tiwary and B. J. Berne, Proc. Natl. Acad. Sci. U. S. A. 113, 2839 (2016)], to determine an optimal 1-dimensional reaction coordinate (RC) for the unbinding of a bucky-ball from a pocket in explicit water. This RC is estimated as a linear combination of the multiple available order parameters that collectively can be used to distinguish the various stable states relevant for unbinding. We pay special attention to determining and quantifying the degree to which water molecules should be included in the RC. Using SGOOP with under-sampled biased simulations, we predict that water plays a distinct role in the reaction coordinate for unbinding in the case when the ligand is sterically constrained to move along an axis of symmetry. This prediction is validated through extensive calculations of the unbinding times through metadynamics and by comparison through detailed balance with unbiased molecular dynamics estimate of the binding time. However when the steric constraint is removed, we find that the role of water in the reaction coordinate diminishes. Here instead SGOOP identifies a good one-dimensional RC involving various motional degrees of freedom.

Elasticity, structure, and relaxation of extended proteins under force.

Force spectroscopies have emerged as a powerful and unprecedented tool to study and manipulate biomolecules directly at a molecular level. Usually, protein and DNA behavior under force is described within the framework of the worm-like chain (WLC) model for polymer elasticity. Although it has been surprisingly successful for the interpretation of experimental data, especially at high forces, the WLC model lacks structural and dynamical molecular details associated with protein relaxation under force that are key to the understanding of how force affects protein flexibility and reactivity. We use molecular dynamics simulations of ubiquitin to provide a deeper understanding of protein relaxation under force. We find that the WLC model successfully describes the simulations of ubiquitin, especially at higher forces, and we show how protein flexibility and persistence length, probed in the force regime of the experiments, are related to how specific classes of backbone dihedral angles respond to applied force. Although the WLC model is an average, backbone model, we show how the protein side chains affect the persistence length. Finally, we find that the diffusion coefficient of the protein's end-to-end distance is on the order of 10(8) nm(2)/s, is position and side-chain dependent, but is independent of the length and independent of the applied force, in contrast with other descriptions.

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.

Unraveling quantum mechanical effects in water using isotopic fractionation.

When two phases of water are at equilibrium, the ratio of hydrogen isotopes in each is slightly altered because of their different phase affinities. This isotopic fractionation process can be utilized to analyze water's movement in the world's climate. Here we show that equilibrium fractionation ratios, an entirely quantum mechanical property, also provide a sensitive probe to assess the magnitude of nuclear quantum fluctuations in water. By comparing the predictions of a series of water models, we show that those describing the OH chemical bond as rigid or harmonic greatly overpredict the magnitude of isotope fractionation. Models that account for anharmonicity in this coordinate are shown to provide much more accurate results because of their ability to give partial cancellation between inter- and intramolecular quantum effects. These results give evidence of the existence of competing quantum effects in water and allow us to identify how this cancellation varies across a wide-range of temperatures. In addition, this work demonstrates that simulation can provide accurate predictions and insights into hydrogen fractionation.

On achieving high accuracy and reliability in the calculation of relative protein-ligand binding affinities.

We apply a free energy perturbation simulation method, free energy perturbation/replica exchange with solute tempering, to two modifications of protein-ligand complexes that lead to significant conformational changes, the first in the protein and the second in the ligand. The approach is shown to facilitate sampling in these challenging cases where high free energy barriers separate the initial and final conformations and leads to superior convergence of the free energy as demonstrated both by consistency of the results (independence from the starting conformation) and agreement with experimental binding affinity data. The second case, consisting of two neutral thrombin ligands that are taken from a recent medicinal chemistry program for this interesting pharmaceutical target, is of particular significance in that it demonstrates that good results can be obtained for large, complex ligands, as opposed to relatively simple model systems. To achieve quantitative agreement with experiment in the thrombin case, a next generation force field, Optimized Potentials for Liquid Simulations 2.0, is required, which provides superior charges and torsional parameters as compared to earlier alternatives.

Ligand binding to protein-binding pockets with wet and dry regions.

Biological processes often depend on protein-ligand binding events, yet accurate calculation of the associated energetics remains as a significant challenge of central importance to structure-based drug design. Recently, we have proposed that the displacement of unfavorable waters by the ligand, replacing them with groups complementary to the protein surface, is the principal driving force for protein-ligand binding, and we have introduced the WaterMap method to account this effect. However, in spite of the adage "nature abhors vacuum," one can occasionally observe situations in which a portion of the receptor active site is so unfavorable for water molecules that a void is formed there. In this paper, we demonstrate that the presence of dry regions in the receptor has a nontrivial effect on ligand binding affinity, and suggest that such regions may represent a general motif for molecular recognition between the dry region in the receptor and the hydrophobic groups in the ligands. With the introduction of a term attributable to the occupation of the dry regions by ligand atoms, combined with the WaterMap calculation, we obtain excellent agreement with experiment for the prediction of relative binding affinities for a number of congeneric ligand series binding to the major urinary protein receptor. In addition, WaterMap when combined with the cavity contribution is more predictive than at least one specific implementation [Abel R, Young T, Farid R, Berne BJ, Friesner RA (2008) J Am Chem Soc 130:2817-2831] of the popular MM-GBSA approach to binding affinity calculation.

Probing the chemistry of thioredoxin catalysis with force.

Thioredoxins are enzymes that catalyse disulphide bond reduction in all living organisms. Although catalysis is thought to proceed through a substitution nucleophilic bimolecular (S(N)2) reaction, the role of the enzyme in modulating this chemical reaction is unknown. Here, using single-molecule force-clamp spectroscopy, we investigate the catalytic mechanism of Escherichia coli thioredoxin (Trx). We applied mechanical force in the range of 25-600 pN to a disulphide bond substrate and monitored the reduction of these bonds by individual enzymes. We detected two alternative forms of the catalytic reaction, the first requiring a reorientation of the substrate disulphide bond, causing a shortening of the substrate polypeptide by 0.79 +/- 0.09 A (+/- s.e.m.), and the second elongating the substrate disulphide bond by 0.17 +/- 0.02 A (+/- s.e.m.). These results support the view that the Trx active site regulates the geometry of the participating sulphur atoms with sub-ångström precision to achieve efficient catalysis. Our results indicate that substrate conformational changes may be important in the regulation of Trx activity under conditions of oxidative stress and mechanical injury, such as those experienced in cardiovascular disease. Furthermore, single-molecule atomic force microscopy techniques, as shown here, can probe dynamic rearrangements within an enzyme's active site during catalysis that cannot be resolved with any other current structural biological technique.

Water's role in the force-induced unfolding of ubiquitin.

In atomic force spectroscopic studies of the elastomeric protein ubiquitin, the ß-strands 1-5 serve as the force clamp. Simulations show how the rupture force in the force-induced unfolding depends on the kinetics of water molecule insertion into positions where they can eventually form hydrogen bonding bridges with the backbone hydrogen bonds in the force-clamp region. The intrusion of water into this region is slowed down by the hydrophobic shielding effect of carbonaceous groups on the surface residues of ß-strands 1-5, which thereby regulates water insertion prior to hydrogen bond breakage. The experiments show that the unfolding of the mechanically stressed protein is nonexponential due to static disorder. Our simulations show that different numbers and/or locations of bridging water molecules give rise to a long-lived distribution of transition states and static disorder. We find that slowing down the translational (not rotational) motions of the water molecules by increasing the mass of their oxygen atoms, which leaves the force field and thereby the equilibrium structure of the solvent unchanged, increases the average rupture force; however, the early stages of the force versus time behavior are very similar for our "normal" and fictitious "heavy" water models. Finally, we construct six mutant systems to regulate the hydrophobic shielding effect of the surface residues in the force-clamp region. The mutations in the two termini of ß-sheets 1-5 are found to determine a preference for different unfolding pathways and change mutant's average rupture force.

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.

Single homopolypeptide chains collapse into mechanically rigid conformations.

Huntington's disease is linked to the insertion of glutamine (Q) in the protein huntingtin, resulting in polyglutamine (polyQ) expansions that self-associate to form aggregates. While polyQ aggregation has been the subject of intense study, a correspondingly thorough understanding of individual polyQ chains is lacking. Here we demonstrate a single molecule force-clamp technique that directly probes the mechanical properties of single polyQ chains. We have made polyQ constructs of varying lengths that span the length range of normal and diseased polyQ expansions. Each polyQ construct is flanked by the I27 titin module, providing a clear mechanical fingerprint of the molecule being pulled. Remarkably, under the application of force, no extension is observed for any of the polyQ constructs. This is in direct contrast with the random coil protein PEVK of titin, which readily extends under force. Our measurements suggest that polyQ chains form mechanically stable collapsed structures. We test this hypothesis by disrupting polyQ chains with insertions of proline residues and find that their mechanical extensibility is sensitive to the position of the proline interruption. These experiments demonstrate that polyQ chains collapse to form a heterogeneous ensemble of conformations that are mechanically resilient. We further use a heat-annealing molecular dynamics protocol to extensively search the conformation space and find that polyQ can exist in highly mechanically stable compact globular conformations. The mechanical rigidity of these collapsed structures may exceed the functional ability of eukaryotic proteasomes, resulting in the accumulation of undigested polyQ sequences in vivo.

Observation of a dewetting transition in the collapse of the melittin tetramer.

Marked hydration changes occur during the self-assembly of the melittin protein tetramer in water. Hydrophobicity induces a drying transition in the gap between simple sufficiently large (more than 1 nm(2)) strongly hydrophobic surfaces as they approach each other, resulting in the subsequent collapse of the system, as well as a depletion of water next to single surfaces. Here we investigate whether the hydrophobic induced collapse of multidomain proteins or the formation of protein oligimers exhibits a similar drying transition. We performed computer simulations to study the collapse of the tetramer of melittin in water, and observed a marked water drying transition inside a nanoscale channel of the tetramer (with a channel size of up to two or three water-molecule diameters). This transition, although occurring on a microscopic length scale, is analogous to a first-order phase transition from liquid to vapour. We find that this drying is very sensitive to single mutations of the three isoleucines to less hydrophobic residues and that such mutations in the right locations can switch the channel from being dry to being wet. Thus, quite subtle changes in hydrophobic surface topology can profoundly influence the drying transition. We show that, even in the presence of the polar protein backbone, sufficiently hydrophobic protein surfaces can induce a liquid-vapour transition providing an enormous driving force towards further collapse. This behaviour was unexpected because of the absence of drying in the collapse of the multidomain protein 2,3-dihydroxybiphenyl dioxygenase (BphC).

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.

Urea denaturation by stronger dispersion interactions with proteins than water implies a 2-stage unfolding.

The mechanism of denaturation of proteins by urea is explored by using all-atom microseconds molecular dynamics simulations of hen lysozyme generated on BlueGene/L. Accumulation of urea around lysozyme shows that water molecules are expelled from the first hydration shell of the protein. We observe a 2-stage penetration of the protein, with urea penetrating the hydrophobic core before water, forming a "dry globule." The direct dispersion interaction between urea and the protein backbone and side chains is stronger than for water, which gives rise to the intrusion of urea into the protein interior and to urea's preferential binding to all regions of the protein. This is augmented by preferential hydrogen bond formation between the urea carbonyl and the backbone amides that contributes to the breaking of intrabackbone hydrogen bonds. Our study supports the "direct interaction mechanism" whereby urea has a stronger dispersion interaction with protein than water.