Jennifer L. Knight

λ‐Dynamics free energy simulation methods

Free energy calculations are fundamental to obtaining accurate theoretical estimates of many important biological phenomena including hydration energies, protein‐ligand binding affinities and energetics of conformational changes. Unlike traditional free energy perturbation and thermodynamic integration methods, λ‐dynamics treats the conventional “λ” as a dynamic variable in free energy simulations and simultaneously evaluates thermodynamic properties for multiple states in a single simulation. In the present article, we provide an overview of the theory of λ‐dynamics, including the use of biasing and restraining potentials to facilitate conformational sampling. We review how λ‐dynamics has been used to rapidly and reliably compute relative hydration free energies and binding affinities for series of ligands, to accurately identify crystallographically observed binding modes starting from incorrect orientations, and to model the effects of mutations upon protein stability. Finally, we suggest how λ‐dynamics may be extended to facilitate modeling efforts in structure‐based drug design.

MATCH: An Atom- Typing Toolset for Molecular Mechanics Force Fields

We introduce a toolset of program libraries collectively titled multipurpose atom‐typer for CHARMM (MATCH) for the automated assignment of atom types and force field parameters for molecular mechanics simulation of organic molecules. The toolset includes utilities for the conversion of multiple chemical structure file formats into a molecular graph. A general chemical pattern‐matching engine using this graph has been implemented whereby assignment of molecular mechanics atom types, charges, and force field parameters are achieved by comparison against a customizable list of chemical fragments. While initially designed to complement the CHARMM simulation package and force fields by generating the necessary input topology and atom‐type data files, MATCH can be expanded to any force field and program, and has core functionality that makes it extendable to other applications such as fragment‐based property prediction. In this work, we demonstrate the accurate construction of atomic parameters of molecules within each force field included in CHARMM36 through exhaustive cross validation studies illustrating that bond charge increment rules derived from one force field can be transferred to another. In addition, using leave‐one‐out substitution it is shown that it is also possible to substitute missing intra and intermolecular parameters with ones included in a force field to complete the parameterization of novel molecules. Finally, to demonstrate the robustness of MATCH and the coverage of chemical space offered by the recent CHARMM general force field (Vanommeslaeghe, et al., J Comput Chem 2010, 31, 671), one million molecules from the PubChem database of small molecules are typed, parameterized, and minimized.

Surveying implicit solvent models for estimating small molecule absolute hydration free energies

Implicit solvent models are powerful tools in accounting for the aqueous environment at a fraction of the computational expense of explicit solvent representations. Here, we compare the ability of common implicit solvent models (TC, OBC, OBC2, GBMV, GBMV2, GBSW, GBSW/MS, GBSW/MS2 and FACTS) to reproduce experimental absolute hydration free energies for a series of 499 small neutral molecules that are modeled using AMBER/GAFF parameters and AM1‐BCC charges. Given optimized surface tension coefficients for scaling the surface area term in the nonpolar contribution, most implicit solvent models demonstrate reasonable agreement with extensive explicit solvent simulations (average difference 1.0–1.7 kcal/mol and R2 = 0.81–0.91) and with experimental hydration free energies (average unsigned errors = 1.1–1.4 kcal/mol and R2 = 0.66–0.81). Chemical classes of compounds are identified that need further optimization of their ligand force field parameters and others that require improvement in the physical parameters of the implicit solvent models themselves. More sophisticated nonpolar models are also likely necessary to more effectively represent the underlying physics of solvation and take the quality of hydration free energies estimated from implicit solvent models to the next level.

On the mechanism of crystalline polymorph selection by polymer heteronuclei.

The phase-selective crystallization of acetaminophen (ACM) using insoluble polymers as heteronuclei was investigated in a combined experimental and computational effort to elucidate the mechanism of polymer-induced heteronucleation (PIHn). ACM heteronucleates from supersaturated aqueous solution in its most thermodynamically stable monoclinic form on poly(n-butyl methacrylate), whereas the metastable orthorhombic form is observed on poly(methyl methacrylate). When ACM crystals were grown through vapor deposition, only the monoclinic polymorph was observed on each polymer. Each crystallization condition leads to a unique powder X-ray diffraction pattern with the major preferred orientation corresponding to the crystallographic faces in which these crystal phases nucleate from surfaces of the polymers. The molecular recognition events leading to these outcomes are elucidated with the aid of computed polymer-crystal binding energies using docking simulations. This investigation illuminates the mechanism by which phase selection occurs during the crystallization of ACM using polymers as heteronuclei, paving the way for the improvement of methods for polymorph selection and discovery based on heterogeneous nucleation promoters.

Assessing the quality of absolute hydration free energies among CHARMM‐compatible ligand parameterization schemes

Multipurpose atom‐typer for CHARMM (MATCH), an atom‐typing toolset for molecular mechanics force fields, was recently developed in our laboratory. Here, we assess the ability of MATCH‐generated parameters and partial atomic charges to reproduce experimental absolute hydration free energies for a series of 457 small neutral molecules in GBMV2, Generalized Born with a smooth SWitching (GBSW), and fast analytical continuum treatment of solvation (FACTS) implicit solvent models. The quality of hydration free energies associated with small molecule parameters obtained from ParamChem, SwissParam, and Antechamber are compared. Given optimized surface tension coefficients for scaling the surface area term in the nonpolar contribution, these automated parameterization schemes with GBMV2 and GBSW demonstrate reasonable agreement with experimental hydration free energies (average unsigned errors of 0.9–1.5 kcal/mol and R2 of 0.63–0.87). GBMV2 and GBSW consistently provide slightly more accurate estimates than FACTS, whereas Antechamber parameters yield marginally more accurate estimates than the current generation of MATCH, ParamChem, and SwissParam parameterization strategies. Modeling with MATCH libraries that are derived from different CHARMM topology and parameter files highlights the importance of having sufficient coverage of chemical space within the underlying databases of these automated schemes and the benefit of targeting specific functional groups for parameterization efforts to maximize both the breadth and the depth of the parameterized space.

Applying efficient implicit nongeometric constraints in alchemical free energy simulations

Several strategies have been developed for satisfying bond lengths, angle, and other geometric constraints in molecular dynamics simulations. Advanced variations of alchemical free energy perturbation simulations, however, also require nongeometric constraints. In our recently developed multisite λ‐dynamics simulation method, the conventional λ parameters that are associated with the progress variables in alchemical transformations are treated as dynamic variables and are constrained such that: 0 ≤ λi ≤ 1 and ∑  i = 1N λi = 1. Here, we present four functional forms of λ that implicitly satisfy these nongeometric constraints, whose values and forces are facile to compute and that yield stable simulations using a 2 fs integration timestep. Using model systems, we present the sampling characteristics of these functional forms and demonstrate the enhanced sampling profiles and improved convergence rates that are achieved by the functional form:

Deconstructing activation events in rhodopsin.

\lambda _i = {{e^{c\sin \theta _i } } \over {\sum\nolimits_{j = 1}^N {e^{c\sin \theta _j } } }}

Constant pH Molecular Dynamics Simulations of Nucleic Acids in Explicit Solvent

that oscillates between λi = 0 and λi = 1 and has relatively steep transitions between these endpoints.

pH-Dependent Dynamics of Complex RNA Macromolecules

Activation of class-A G-protein-coupled receptors (GPCRs) involves large-scale reorganization of the H3/H6 interhelical network. In rhodopsin (Rh), this process is coupled to a change in the protonation state of a key residue, E134, whose exact role in activation is not well understood. Capturing this millisecond pH-dependent process is a well-appreciated challenge. We have developed a scheme combining the harmonic Fourier beads (HFB) method and constant-pH molecular dynamics with pH-based replica exchange (pH-REX) to gain insight into the structural changes that occur along the activation pathway as a function of the protonation state of E134. Our results indicate that E134 is protonated as a consequence of tilting of H6 by ca. 4.0° with respect to its initial position and simultaneous rotation by ca. 23° along its principal axis. The movement of H6 is associated with breakage of the E247-R135 and R135-E134 salt bridges and concomitant release of the E134 side chain, which results in an increase in its pKa value above physiological pH. An increase in the hydrophobicity of the environment surrounding E134 leads to further tilting and rotation of H6 and upshift of the E134 pKa. Such atomic-level information, which is not accessible through experiments, refines the earlier proposed sequential model of Rh activation (see: Zaitseva, E.; et al. Sequential Rearrangement of Interhelical Networks Upon Rhodopsin Activation in Membranes: The Meta IIa Conformational Substate . J. Am. Chem. Soc. 2010, 132, 4815) and argues that the E134 protonation switch is both a cause and a consequence of the H6 motion.