Force Field X – Publications

Tollefson, M. R., Gogal, R. A., Weaver, A. M., Schaefer, A. M., Marini, R. J., Azaiez, H., Kolbe, D. L., Wang, D., Weaver, A. E., Casavant, T. L., Braun, T. A., Smith, R. J. H., and Schnieders, M. J., Assessing Variants of Uncertain Significance Implicated in Hearing Loss Using a Comprehensive Deafness Proteome Submitted: Hearing loss is the leading sensory deficit, affecting ~5% of the population. It exhibits remarkable heterogeneity across 223 genes with 6,328 pathogenic missense variants, making deafness-specific expertise a prerequisite for ascribing phenotypic consequences to genetic variants. Deafness-implicated variants are curated in the Deafness Variation Database (DVD) after classification by a genetic hearing loss expert panel and thorough informatics pipeline. However, seventy percent of the 128,167 missense variants in the DVD are “variants of uncertain significance” (VUS) due to insufficient evidence for classification. Here, we use the deep learning protein prediction algorithm, AlphaFold2, to curate structures for all DVD genes. We refine these structures with global optimization and the AMOEBA force field and use DDGun3D to predict folding free energy differences (∆∆GFold) for all DVD missense variants. We find that 5,772 VUSs have a large, destabilizing ∆∆GFold that is consistent with pathogenic variants. When also filtered for CADD scores (>25.7), we determine 3,456 VUSs are likely pathogenic at a probability of 99.0%. These VUSs affect 119 patients (~3% of cases) sequenced by the OtoSCOPE targeted panel. Approximately half of these patients previously received an inconclusive report, and reclassification of these VUSs as pathogenic provides a new genetic diagnosis for six patients.

2022

Nessler, A. J., Okada, O., Hermon, M. J., Nagata, H., and Schnieders, M. J., Progressive alignment of crystals: reproducible and efficient assessment of crystal structure similarity. Journal of Applied Crystallography 2022, 55, 1528-1537.: During in silico crystal structure prediction of organic molecules, millions of candidate structures are often generated. These candidates must be compared to remove duplicates prior to further analysis (e.g. optimization with electronic structure methods) and ultimately compared with structures determined experimentally. The agreement of predicted and experimental structures forms the basis of evaluating the results from the Cambridge Crystallographic Data Centre (CCDC) blind assessment of crystal structure prediction, which further motivates the pursuit of rigorous alignments. Evaluating crystal structure packings using coordinate root-mean-square deviation (RMSD) for N molecules (or N asymmetric units) in a reproducible manner requires metrics to describe the shape of the compared molecular clusters to account for alternative approaches used to prioritize selection of molecules. Described here is a flexible algorithm called Progressive Alignment of Crystals (PAC) to evaluate crystal packing similarity using coordinate RMSD and introducing the radius of gyration (Rg) as a metric to quantify the shape of the superimposed clusters. It is shown that the absence of metrics to describe cluster shape adds ambiguity to the results of the CCDC blind assessments because it is not possible to determine whether the superposition algorithm has prioritized tightly packed molecular clusters (i.e. to minimize Rg) or prioritized reduced RMSD (i.e. via possibly elongated clusters with relatively larger Rg). For example, it is shown that when the PAC algorithm described here uses single linkage to prioritize molecules for inclusion in the superimposed clusters, the results are nearly identical to those calculated by the widely used program COMPACK. However, the lower Rg values obtained by the use of average linkage are favored for molecule prioritization because the resulting RMSDs more equally reflect the importance of packing along each dimension. It is shown that the PAC algorithm is faster than COMPACK when using a single process and its utility for biomolecular crystals is demonstrated. Finally, parallel scaling up to 64 processes in the open-source code Force Field X is presented.
Dybeck, E. C., Thiel, A., Schnieders, M. J., Pickard, F. C. I. V., Wood, G. P. F., Krzyzaniak, J. F., and Hancock, B. C., A Comparison of Methods for Computing Relative Anhydrous–Hydrate Stability with Molecular Simulation, Crystal Growth and Design 2022, 23, 142-167.: The transformation of a pharmaceutical solid from an anhydrous crystal into a hydrated form during drug development represents a risk to a product’s safety and efficacy due to the potential impact on stability, bioavailability, and manufacturability. In this work, we examine 10 classical free energy simulation protocols to evaluate the thermodynamic stability of hydrated crystals relative to their anhydrous forms. Molecular dynamics simulations are used to compute the Gibbs free energies of the crystals of three pharmaceutically relevant systems using two fixed-charge potentials, GAFF and OPLS, as well as the polarizable AMOEBA model. In addition, we explore a variety of water models, including TIP3P, TIP4P, and AMOEBA, for both the interstitial water and the effects of ambient humidity. The AMOEBA model predicts free energy values most consistent with experimental measurements among the models examined. The benefits of a fully polarizable water model relative to fixed-charged models appear to derive predominantly from a better treatment of water’s dipole moment in the crystalline phase. Despite this improved physical treatment, we find that no single model produces reliable predictions of the phase boundary between hydrated and anhydrous crystals from theory alone. However, we show that accurate phase diagrams can be constructed from the simulations by introducing a single experimentally determined coexistence point. With this single experimental data point as input, the phase boundary is correctly predicted within 10% relative humidity on the temperature range of 15 to 75 °C for all three systems examined. Furthermore, we demonstrate that with this known coexistence point as an input, the differences between the various potentials and the water models become insignificant, as all models yield accurate phase boundaries regardless of whether polarization is included due to significant temperature-dependent error cancellation between models.

2021

Corrigan, R. A., Qi, G., Thiel, A. C., Lynn, J. R., Walker, B. D., Casavant, T. L., Lagardere, L., Piquemal, J.-P., Ponder, J. W., Ren, P., and Schnieders, M. J., Implicit Solvents for the Polarizable Atomic Multipole AMOEBA Force Field. Journal of Chemical Theory and Computation 2021, 17 (4), 2323-2341.: Computational protein design, ab initio protein/RNA folding, and protein–ligand screening can be too computationally demanding for explicit treatment of solvent. For these applications, implicit solvent offers a compelling alternative, which we describe here for the polarizable atomic multipole AMOEBA force field based on three treatments of continuum electrostatics: numerical solutions to the nonlinear and linearized versions of the Poisson–Boltzmann equation (PBE), the domain-decomposition conductor-like screening model (ddCOSMO) approximation to the PBE, and the analytic generalized Kirkwood (GK) approximation. The continuum electrostatics models are combined with a nonpolar estimator based on novel cavitation and dispersion terms. Electrostatic model parameters are numerically optimized using a least-squares style target function based on a library of 103 small-molecule solvation free energy differences. Mean signed errors for the adaptive Poisson–Boltzmann solver (APBS), ddCOSMO, and GK models are 0.05, 0.00, and 0.00 kcal/mol, respectively, while the mean unsigned errors are 0.70, 0.63, and 0.58 kcal/mol, respectively. Validation of the electrostatic response of the resulting implicit solvents, which are available in the Tinker (or Tinker-HP), OpenMM, and Force Field X software packages, is based on comparisons to explicit solvent simulations for a series of proteins and nucleic acids. Overall, the emergence of performative implicit solvent models for polarizable force fields opens the door to their use for folding and design applications.

2019

Tollefson, M. R., Litman, J. M., Qi, G., O’Connell, C. E., Wipfler, M. J., Marini, R. J., Bernabe, H. V., Tollefson, W. T. A., Braun, T. A., Casavant, T. L., Smith, R. J. H. and M. J. Schnieders, Structural Insights into Hearing Loss Genetics from Polarizable Protein Repacking. Biophysical Journal 2019, 117 (3), 602-612.: Hearing loss is associated with ?8100 mutations in 152 genes, and within the coding regions of these genes are over 60,000 missense variants. The majority of these variants are classified as variants of uncertain significance to reflect our inability to ascribe a phenotypic effect to the observed amino acid change. A promising source of pathogenicity information is biophysical simulation, although input protein structures often contain defects because of limitations in experimental data and/or only distant homology to a template. Here, we combine the polarizable atomic multipole optimized energetics for biomolecular applications force field, many-body optimization theory, and graphical processing unit acceleration to repack all deafness-associated proteins and thereby improve average structure MolProbity score from 2.2 to 1.0. We then used these optimized wild-type models to create over 60,000 structures for missense variants in the Deafness Variation Database, which are being incorporated into the Deafness Variation Database to inform deafness pathogenicity prediction. Finally, this work demonstrates that advanced polarizable atomic multipole force fields are efficient enough to repack the entire human proteome.
Litman, J. M., Thiel, A. C., and M. J. Schnieders, Scalable Indirect Free Energy Method Applied to Divalent Cation-Metalloprotein Binding. Journal of Chemical Theory and Computation 2019, 15 (8), 4602-4614.: Many biological processes are based on molecular recognition between highly charged molecules such as nucleic acids, inorganic ions, charged amino acids, etc. For such cases, it has been demonstrated that molecular simulations with fixed partial charges often fail to achieve experimental accuracy. Although incorporation of more advanced electrostatic models (such as multipoles, mutual polarization, etc.) can significantly improve simulation accuracy, it increases computational expense by a factor of 5–20×. Indirect free energy (IFE) methods can mitigate this cost by modeling intermediate states at fixed-charge resolution. For example, an efficient “reference” model such as a pairwise Amber, CHARMM, or OPLS-AA force field can be used to derive an initial estimate, followed by thermodynamic corrections to a more advanced “target” potential such as the polarizable AMOEBA model. Unfortunately, all currently described IFE methods encounter difficulties reweighting more than ∼50 atoms between resolutions due to extensive scaling of both the magnitude of the thermodynamic corrections and their statistical uncertainty. We present an approach called “simultaneous bookending” (SB) that is fundamentally different from existing IFE methods based on a tunable sampling approximation, which permits scaling to thousands of atoms. SB is demonstrated on the relative binding affinity of Mg²⁺/Ca²⁺ to a set of metalloproteins with up to 2972 atoms, finding no statistically significant difference between direct AMOEBA results and those from correcting Amber to AMOEBA. The ability to change the resolution of thousands of atoms during reweighting suggests the approach may be applicable in the future to protein–protein binding affinities or nucleic acid thermodynamics.

2016

Nessler, I. J., Litman, J. M., and M. J. Schnieders, Toward polarizable AMOEBA thermodynamics at fixed charge efficiency using a dual force field approach. Physical Chemistry Chemical Physics 2016, 18 (44), 30313-30322.: First principles prediction of the structure, thermodynamics and solubility of organic molecular crystals, which play a central role in chemical, material, pharmaceutical and engineering sciences, challenges both potential energy functions and sampling methodologies. Here we calculate absolute crystal deposition thermodynamics using a novel dual force field approach whose goal is to maintain the accuracy of advanced multipole force fields (e.g. the polarizable AMOEBA model) while performing more than 95% of the sampling in an inexpensive fixed charge (FC) force field (e.g. OPLS-AA). Absolute crystal sublimation/deposition phase transition free energies were determined using an alchemical path that grows the crystalline state from a vapor reference state based on sampling with the OPLS-AA force field, followed by dual force field thermodynamic corrections to change between FC and AMOEBA resolutions at both end states (we denote the three step path as AMOEBA/FC). Importantly, whereas the phase transition requires on the order of 200 ns of sampling per compound, only 5 ns of sampling was needed for the dual force field thermodynamic corrections to reach a mean statistical uncertainty of 0.05 kcal mol⁻¹. For five organic compounds, the mean unsigned error between direct use of AMOEBA and the AMOEBA/FC dual force field path was only 0.2 kcal mol⁻¹ and not statistically significant. Compared to experimental deposition thermodynamics, the mean unsigned error for AMOEBA/FC (1.4 kcal mol⁻¹) was more than a factor of two smaller than uncorrected OPLS-AA (3.2 kcal mol⁻¹). Overall, the dual force field thermodynamic corrections reduced condensed phase sampling in the expensive force field by a factor of 40, and may prove useful for protein stability or binding thermodynamics in the future.
Bell, D., Qi, R., Jing, Z., Xiang, J. Y., Mejias, C., Schnieders, M., Ponder, J., and P. Ren, Calculating binding free energies of host-guest systems using AMOEBA polarizable force field. Physical Chemistry Chemical Physics 2016, 18 (44), 30261-30269.: Molecular recognition is of paramount interest in many applications. Here we investigate a series of host–guest systems previously used in the SAMPL4 blind challenge by using molecular simulations and the AMOEBA polarizable force field. The free energy results computed by Bennett's acceptance ratio (BAR) method using the AMOEBA polarizable force field ranked favorably among the entries submitted to the SAMPL4 host–guest competition. In this work we conduct an in-depth analysis of the AMOEBA force field host–guest binding thermodynamics by using both BAR and the orthogonal space random walk (OSRW) methods. The binding entropy–enthalpy contributions are analyzed for each host–guest system. For systems of inordinate binding entropy–enthalpy values, we further examine the hydrogen bonding patterns and configurational entropy contribution. The binding mechanism of this series of host–guest systems varies from ligand to ligand, driven by enthalpy and/or entropy changes. Convergence of BAR and OSRW binding free energy methods is discussed. Ultimately, this work illustrates the value of molecular modelling and advanced force fields for the exploration and interpretation of binding thermodynamics.

2015

LuCore, S. D., Litman, J. M., Powers, K. T., Gao, S., Lynn, A. M., Tollefson, W. T. A., Fenn, T. D., Washington, M. T. and M. J. Schnieders, Dead-End Elimination with a Polarizable Force Field Repacks PCNA Models. Biophysical Journal 2015, 109 (4), 816-826.: A balance of van der Waals, electrostatic, and hydrophobic forces drive the folding and packing of protein side-chains. Although residue interactions are often approximated as being pairwise additive, in reality higher order ‘many-body’ contributions that depend on environment drive hydrophobic collapse and cooperative electrostatics. Beginning from dead-end elimination (DEE), we derive the first algorithm capable of deterministic repacking of side-chains compatible with ‘many-body’ energy functions. The approach is applied to seven PCNA X-ray crystallographic data sets with resolutions 2.5-3.8 Å (mean 3.0 Å). While PDB_Redo models average an R_free value of 29.5% and MolProbity score of 2.71 Å (77^th percentile), DEE with the polarizable AMOEBA force field lowered R_free by 2.8% to 26.7% and improved mean MolProbity score to atomic resolution at 1.25 Å (100^th percentile). For applications that depend on side-chain repacking, including structure refinement and protein design, the accuracy limitations of pairwise additivity can now be eliminated via polarizable or quantum mechanical potentials.
Shi, Y., Schnieders, M. J., Piquemal, J.-P., and P. Ren, Polarizable Force Fields for Biomolecular Modeling. Reviews in Computational Chemistry 2015 Edition: 28, Publisher: Springer, Editors: Kenny B. Lipkowitz: Molecular mechanics based modeling has been widely used in the study of chemical and biological systems. The classical potential energy functions and their parameters are referred to as force fields. Empirical force fields for biomolecules emerged in the early 1970's, followed by the first molecular dynamics simulations of the bovine pancreatic trypsin inhibitors (BPTI). Over the past 30 years, a great number of empirical molecular mechanics force fields, including AMBER, CHARMM, GROMOS, OPLS, and many others, have been developed. These force fields share similar functional forms, including valence interactions represented by harmonic oscillators, point dispersionrepulsion for van der Waals (vdW) interactions, and an electrostatic contribution based on fixed atomic partial charges. This generation of molecular mechanics force fields has been widely used in the study of molecular structures, dynamics, interactions, design and engineering.
Lipparini, F., Lagardère, L., Raynaud, C., Stamm, B., Cancès, E., Mennucci, B., Schnieders, M. J., Ren P., Maday, Y. and J.-P. Piquemal, Polarizable Molecular Dynamics in a Polarizable Continuum Solvent. Journal of Chemical Theory and Computation 2015 11 (2), 623-634.: We present for the first time polarizable molecular dynamics (MD) simulations within a polarizable continuum solvent with molecular shape cavities and exact solution of the mutual polarization. The key ingredients are a very efficient algorithm for solving the equations associated with the polarizable continuum, in particular, the domain decomposition Conductor-like Screening Model (ddCOSMO), a rigorous coupling of the continuum with the polarizable force field achieved through a robust variational formulation and an effective strategy to solve the coupled equations. The coupling of ddCOSMO with non variational force fields, including AMOEBA, is also addressed. The MD simulations are feasible, for real life systems, on standard cluster nodes; a scalable parallel implementation allows for further speed up in the context of a newly developed module in Tinker, named Tinker-HP. NVE simulations are stable and long term energy conservation can be achieved. This paper is focused on the methodological developments, on the analysis of the algorithm and on the stability of the simulations; a proof-of-concept application is also presented to attest the possibilities of this newly developed technique.

2014

Park, J., Nessler, I., McClain, B., Macikenas D., Baltrusaitis, J. and M. J. Schnieders, Absolute Organic Crystal Thermodynamics: Growth of the Asymmetric Unit into a Crystal via Alchemy. Journal of Chemical Theory and Computation 2014 10, (7), 2781-2791.: The solubility of organic molecules is of critical importance to the pharmaceutical industry, however, robust computational methods to predict this quantity from first principles are lacking. Solubility can be computed from a thermodynamic cycle that decomposes standard state solubility into the sum of solid-vapor sublimation and vapor-liquid solvation free energies ΔG° _solubility = ΔG° _sub + ΔG° _solv. Over the past few decades, alchemical simulation methods to compute solvation free energy using classical force fields have become widely used. However, analogous methods for determining the free energy of the sublimation/deposition phase transition are currently limited by the necessity of a priori knowledge of the atomic coordinates of the crystal. Here we describe progress toward an alternative scheme based on growth of the asymmetric unit into a crystal via alchemy (GAUCHE). GAUCHE computes deposition free energy ΔG°_dep=-ΔG°_sub= k_BTln(V_AU⁄1 M) + ΔG°_AU + ΔG°_(AU→UC) as the sum of an entropic term to account for compressing a 1 M vapor into the molar volume of the crystal asymmetric unit (VAU), where kB is Boltzmann’s constant and T is temperature in degrees Kelvin, plus two simulation steps. In the first simulation step, the deposition free energy ΔG°_AU for a system composed of only NAU asymmetric unit (AU) molecule(s) is computed beginning from an arbitrary conformation in vacuum. In the second simulation step, the change in free energy ΔG°_(AU→UC) to expand the asymmetric unit degrees of freedom into a unit cell (UC) composed of NUC independent molecules is computed. This latter step accounts for the favorable free energy of removing the constraint that every symmetry mate of the asymmetric unit has an identical conformation and intermolecular interactions. The current work is based on NVT simulations, which requires knowledge of the crystal space group and unit cell parameters from experiment, but not a priori knowledge of crystalline atomic coordinates. GAUCHE was applied to 5 organic molecules whose sublimation free energy has been measured experimentally, based on the polarizable atomic multipole AMOEBA force field and more than a microsecond of sampling per compound. The mean signed and unsigned errors were only -0.4 and 1.2 kcal/mol, respectively, which indicates that GAUCHE is capable of accurately predicts sublimation thermodynamics.
Lipparini, F., Lagardère, L., Stamm, B., Cancès, E., Schnieders, M. J., Ren P., Maday, Y. and J.-P. Piquemal, Scalable Evaluation of Polarization Energy and Associated Forces in Polarizable Molecular Dynamics: I. Toward Massively Parallel Direct Space Computations. Journal of Chemical Theory and Computation 2014 10 (4), 1638-1651.: In this paper, we investigate various numerical strategies to compute the direct space polarization energy and associated forces in the context of the point dipole approximation (including damping) used in polarizable molecular dynamics. We present a careful mathematical analysis of the algorithms that have been implemented in popular production packages and applied to large test systems. We show that the classical Jacobi Over-Relaxation method (JOR) should not be used as its convergence requires a proper value of the relaxation parameter, whereas other strategies should be preferred. On a single node, Preconditioned Conjugate Gradient methods (PCG) and Jacobi algorithm coupled with the Direct Inversion in the Iterative Subspace (JI/DIIS) provide reliable stability/convergence and are roughly twice as fast as JOR. Moreover, both algorithms are suitable for massively parallel implementations. The lower requirements in terms of processes communications make JI/DIIS the method of choice for MPI and hybrid OpenMP/MPI paradigms for real life tests. Furthermore, using a predictor step as a guess along a molecular dynamics simulation provides another inexpensive, yet very effective, form of convergence acceleration. Overall, two to three orders of magnitude in time can be gained compared to the initial JOR single node approach to the final PGC or JI/DIIS parallel one combined with the predictors MD refinements. Such a speedup traces a new route for the high performance implementation of polarizable molecular dynamics and therefore extends the applicability of the technique as it will facilitate future multiscale QM/MM/continuum computations.

2012

Ren, P., Chun, J., Thomas, D. G., Schnieders, M. J., Marucho, M., Zhang, J., and N. A. Baker, Biomolecular electrostatics and solvation: a computational perspective. Quarterly Reviews of Biophysics 2012 45, (4), 427-91.: An understanding of molecular interactions is essential for insight into biological systems at the molecular scale. Among the various components of molecular interactions, electrostatics are of special importance because of their long-range nature and their influence on polar or charged molecules, including water, aqueous ions, proteins, nucleic acids, carbohydrates, and membrane lipids. In particular, robust models of electrostatic interactions are essential for understanding the solvation properties of biomolecules and the effects of solvation upon biomolecular folding, binding, enzyme catalysis, and dynamics. Electrostatics, therefore, are of central importance to understanding biomolecular structure and modeling interactions within and among biological molecules. This review discusses the solvation of biomolecules with a computational biophysics view toward describing the phenomenon. While our main focus lies on the computational aspect of the models, we provide an overview of the basic elements of biomolecular solvation (e.g. solvent structure, polarization, ion binding, and non-polar behavior) in order to provide a background to understand the different types of solvation models.
Schnieders, M. J., J. Baltrusaitis, Y. Shi, G. Chattree, L. Zheng, W. Yang and P. Ren The Structure, Thermodynamics, and Solubility of Organic Crystals from Simulation with a Polarizable Force Field. Journal of Chemical Theory and Computation 2012 8, (5), 1721–36.: An important unsolved problem in materials science is prediction of the thermodynamic stability of organic crystals and their solubility from first principles. Solubility can be defined as the saturating concentration of a molecule within a liquid solvent, where the physical picture is of solvated molecules in equilibrium with their solid phase. Despite the importance of solubility in determining the oral bioavailability of pharmaceuticals, prediction tools are currently limited to quantitative structure–property relationships that are fit to experimental solubility measurements. For the first time, we describe a consistent procedure for the prediction of the structure, thermodynamic stability, and solubility of organic crystals from molecular dynamics simulations using the polarizable multipole AMOEBA force field. Our approach is based on a thermodynamic cycle that decomposes standard state solubility into the sum of solid–vapor sublimation and vapor–liquid solvation free energies ΔG° _solubility = ΔG° _sub + ΔG° _solv, which are computed via the orthogonal space random walk (OSRW) sampling strategy. Application to the n-alkylamides series from acetamide through octanamide was selected due to the dependence of their solubility on both amide hydrogen bonding and the hydrophobic effect, which are each fundamental to protein structure and solubility. On average, the calculated absolute standard state solubility free energies are accurate to within 1.1 kcal/mol. The experimental trend of decreasing solubility as a function of n-alkylamide chain length is recapitulated by the increasing stability of the crystalline state and to a lesser degree by decreasing favorability of solvation (i.e., the hydrophobic effect). Our results suggest that coupling the polarizable AMOEBA force field with an orthogonal space based free energy algorithm, as implemented in the program Force Field X , is a consistent procedure for predicting the structure, thermodynamic stability, and solubility of organic crystals.
Schnieders, M. J., T. S. Kaoud, C. Yan; K. N. Dalby and P. Ren Computational Insights for the Discovery of Non-ATP Competitive Inhibitors of MAP Kinases. Current Pharmaceutical Design 2012 18, (9), 1173-85.: Due to their role in cellular signaling mitogen activated protein (MAP) kinases represent targets of pharmaceutical interest. However, the majority of known MAP kinase inhibitors compete with cellular ATP and target an ATP binding pocket that is highly conserved in the 500 plus representatives of the human protein kinase family. Here we review progress toward the development of non-ATP competitive MAP kinase inhibitors for the extracellular signal regulated kinases (ERK1/2), the c-jun N-terminal kinases (JNK1/2/3) and the p38 MAPKs (α, β, γ, and δ). Special emphasis is placed on the role of computational methods in the drug discovery process for MAP kinases. Topics include recent advances in X-ray crystallography theory that improve the MAP kinase structures essential to structurebased drug discovery, the use of molecular dynamics to understand the conformational heterogeneity of the activation loop and inhibitors discovered by virtual screening. The impact of an advanced polarizable force field such as AMOEBA used in conjunction with sophisticated kinetic and thermodynamic simulation methods is also discussed.

2011

Fenn, T. D. and M. J. Schnieders, Polarizable atomic multipole X-ray refinement: Weighting schemes for macromolecular diffraction. Acta Crystallographica Section D 2011 67, (11), 957-65.: In the past, weighting between the sum of chemical and data-based targets in macromolecular crystallographic refinement was based on comparing the gradients or Hessian diagonal terms of the two potential functions. Here we demonstrate limitations of this scheme, especially in the context of a maximum likelihood target that is inherently weighted by the model and data errors. In fact, the congruence between the maximum likelihood target and a chemical potential based on polarizable atomic multipole electrostatics evaluated with Ewald summation has opened the door to a transferable, static weight. We derive from first principles an optimal static weight and demonstrate that it is transferable across a broad range of data resolutions in the context of our recent implementation of X-ray crystallographic refinement using the polarizable AMOEBA force field and show the resulting models are balanced with respect to optimizing both R _free and Molprobity scores. Conversely, the classical automatic weighting scheme is shown to lead to under or overfitting of the data and poor model geometry. We also highlight the benefits of this approach for low resolution diffraction data, where the need for prior chemical information is of particular importance. We demonstrate that this method is transferable between low and high resolution maximum likelihood based crystallographic refinement, which proves for the first time that resolution dependent parameterization of either the weight or chemical potential is unnecessary.
Fenn, T. D., M. J. Schnieders, M. Mustyakimov, C. Wu, V. S. Pande, P. Langan and A. T. Brunger, Reintroducing electrostatics into macromolecular crystallographic refinement: Application to neutron crystallography and DNA hydration. Structure 2011 19, (4), 523-33.: Most current crystallographic structure refinements augment the diffraction data with a priori information consisting of bond, angle, dihedral, planarity restraints, and atomic repulsion based on the Pauli exclusion principle. Yet, electrostatics and van der Waals attraction are physical forces that provide additional a priori information. Here, we assess the inclusion of electrostatics for the force field used for all-atom (including hydrogen) joint neutron/X-ray refinement. Two DNA and a protein crystal structure were refined against joint neutron/X-ray diffraction data sets using force fields without electrostatics or with electrostatics. Hydrogen-bond orientation/geometry favors the inclusion of electrostatics. Refinement of Z-DNA with electrostatics leads to a hypothesis for the entropic stabilization of Z-DNA that may partly explain the thermodynamics of converting the B form of DNA to its Z form. Thus, inclusion of electrostatics assists joint neutron/X-ray refinements, especially for placing and orienting hydrogen atoms.
Schnieders, M. J., T. D. Fenn and V. S. Pande, Polarizable atomic multipole X-ray refinement: Particle-mesh Ewald electrostatics for macromolecular crystals. Journal of Chemical Theory and Computation 2011 7, (4), 1141-56.: Refinement of macromolecular models from X-ray crystallography experiments benefits from prior chemical knowledge at all resolutions. As the quality of the prior chemical knowledge from quantum or classical molecular physics improves, in principle so will resulting structural models. Due to limitations in computer performance and electrostatic algorithms, commonly used macromolecules X-ray crystallography refinement protocols have had limited support for rigorous molecular physics in the past. For example, electrostatics is often neglected in favor of nonbonded interactions based on a purely repulsive van der Waals potential. In this work we present advanced algorithms for desktop workstations that open the door to X-ray refinement of even the most challenging macromolecular data sets using state-of-the-art classical molecular physics. First we describe theory for particle mesh Ewald (PME) summation that consistently handles the symmetry of all 230 space groups, replicates of the unit cell such that the minimum image convention can be used with a real space cutoff of any size and the combination of space group symmetry with replicates. An implementation of symmetry accelerated PME for the polarizable atomic multipole optimized energetics for biomolecular applications (AMOEBA) force field is presented. Relative to a single CPU core performing calculations on a P1 unit cell, our AMOEBA engine called Force Field X (FFX) accelerates energy evaluations by more than a factor of 24 on an 8-core workstation with a Tesla GPU coprocessor for 30 structures that contain 240000 atoms on average in the unit cell. The benefit of AMOEBA electrostatics evaluated with PME for macromolecular X-ray crystallography refinement is demonstrated via rerefinement of 10 crystallographic data sets that range in resolution from 1.7 to 4.5 Å. Beginning from structures obtained by local optimization without electrostatics, further optimization using AMOEBA with PME electrostatics improved agreement of the model with the data (R _free was lowered by 0.5%), improved geometric features such as favorable (φ,ψ) backbone conformations, and lowered the average potential energy per residue by over 10 kcal/mol. Furthermore, the MolProbity structure validation tool indicates that the geometry of these rerefined structures is consistent with X-ray crystallographic data collected up to 2.2 Å, which is 0.9 Å better than the actual mean quality (3.1 Å). We conclude that polarizable AMOEBA-assisted X-ray refinement offers advantages to methods that neglect electrostatics and is now efficient enough for routine use.

2010

Fenn, T. D., M. J. Schnieders and A. T. Brunger, A smooth and differentiable bulk-solvent model for macromolecular diffraction. Acta Crystallographica Section D 2010, 66, (9), 1024-31.: Inclusion of low-resolution data in macromolecular crystallography requires a model for the bulk solvent. Previous methods have used a binary mask to accomplish this, which has proven to be very effective, but the mask is discontinuous at the solute–solvent boundary (i.e. the mask value jumps from zero to one) and is not differentiable with respect to atomic parameters. Here, two algorithms are introduced for computing bulk-solvent models using either a polynomial switch or a smoothly thresholded product of Gaussians, and both models are shown to be efficient and differentiable with respect to atomic coordinates. These alternative bulk-solvent models offer algorithmic improvements, while showing similar agreement of the model with the observed amplitudes relative to the binary model as monitored using R, R _free and differences between experimental and model phases. As with the standard solvent models, the alternative models improve the agreement primarily with lower resolution (>6 Å) data versus no bulk solvent. The models are easily implemented into crystallographic software packages and can be used as a general method for bulk-solvent correction in macromolecular crystallography. \[\textbf{F}_{{t}}=\textbf{F}_{{c}}-k_{{s}}\textbf{F}_{{m}}{{\rm e}^{-1/4\,B_{{s}}\left|\textbf{s}\right|^{2}}}\]
Fenn, T. D.*, M. J. Schnieders*, A. T. Brunger and V. S. Pande, Polarizable atomic multipole X-ray refinement: Hydration geometry and application to macromolecules. Biophysical Journal 2010, 98, (12), 2984-92 (*joint first authors).: We recently developed a polarizable atomic multipole refinement method assisted by the AMOEBA force field for macromolecular crystallography. Compared to standard refinement procedures, the method uses a more rigorous treatment of X-ray scattering and electrostatics that can significantly improve the resultant information contained in an atomic model. We applied this method to high-resolution lysozyme and trypsin data sets, and validated its utility for precisely describing biomolecular electron density, as indicated by a 0.4–0.6% decrease in the R- and R _free-values, and a corresponding decrease in the relative energy of 0.4–0.8 Kcal/mol/residue. The re-refinements illustrate the ability of force field electrostatics to orient water networks and catalytically relevant hydrogens, which can be used to make predictions regarding active site function, activity, and protein-ligand interaction energies. Re-refinement of a DNA crystal structure generates the zigzag spine pattern of hydrogen bonding in the minor groove without manual intervention. The polarizable atomic multipole electrostatics model implemented in the AMOEBA force field is applicable and informative for crystal structures solved at any resolution.
Ponder, J. W., C. Wu, P. Ren, V. S. Pande, J. D. Chodera, M. J. Schnieders, I. Haque, D. L. Mobley, D. S. Lambrecht, R. A. DiStasio, M. Head-Gordon, G. N. I. Clark, M. E. Johnson and T. Head-Gordon, Current Status of the AMOEBA Polarizable Force Field. The Journal of Physical Chemistry B 2010, 114, 2549-64.: Molecular force fields have been approaching a generational transition over the past several years, moving away from well-established and well-tuned, but intrinsically limited, fixed point charge models toward more intricate and expensive polarizable models that should allow more accurate description of molecular properties. The recently introduced AMOEBA force field is a leading publicly available example of this next generation of theoretical model, but to date, it has only received relatively limited validation, which we address here. We show that the AMOEBA force field is in fact a significant improvement over fixed charge models for small molecule structural and thermodynamic observables in particular, although further fine-tuning is necessary to describe solvation free energies of drug-like small molecules, dynamical properties away from ambient conditions, and possible improvements in aromatic interactions. State of the art electronic structure calculations reveal generally very good agreement with AMOEBA for demanding problems such as relative conformational energies of the alanine tetrapeptide and isomers of water sulfate complexes. AMOEBA is shown to be especially successful on protein−ligand binding and computational X-ray crystallography where polarization and accurate electrostatics are critical.

2009

Schnieders, M. J., T. D. Fenn, V. S. Pande and A. T. Brunger, Polarizable atomic multipole X-ray refinement: Application to peptide crystals. Acta Crystallographica Section D 2009, 65, (9), 952-65.: Recent advances in computational chemistry have produced force fields based on a polarizable atomic multipole description of biomolecular electrostatics. In this work, the Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA) force field is applied to restrained refinement of molecular models against X-ray diffraction data from peptide crystals. A new formalism is also developed to compute anisotropic and aspherical structure factors using fast Fourier transformation (FFT) of Cartesian Gaussian multipoles. Relative to direct summation, the FFT approach can give a speedup of more than an order of magnitude for aspherical refinement of ultrahigh-resolution data sets. Use of a sublattice formalism makes the method highly parallelizable. Application of the Cartesian Gaussian multipole scattering model to a series of four peptide crystals using multipole coefficients from the AMOEBA force field demonstrates that AMOEBA systematically underestimates electron density at bond centers. For the trigonal and tetrahedral bonding geometries common in organic chemistry, an atomic multipole expansion through hexadecapole order is required to explain bond electron density. Alternatively, the addition of interatomic scattering (IAS) sites to the AMOEBA-based density captured bonding effects with fewer parameters. For a series of four peptide crystals, the AMOEBA-IAS model lowered R _free by 20-40% relative to the original spherically symmetric scattering model.
Yue, S. J. Dian, M. J. Schnieders and P. Ren, Trypsin-ligand binding free energy calculation with AMOEBA. Engineering in Medicine and Biology Society. EMBC. Annual International Conference of the IEEE 2009, 2328-31.: The binding free energies of several benzamidine-like inhibitors to trypsin were examined using a polarizable molecular mechanics potential. All the computed binding free energies are in good agreement with the experimental data. From free energy decomposition, electrostatic interaction was indicated to be the driving force for the binding. MD simulations show that the ligands form hydrogen bonds with trypsin and water molecules nearby in a competitive fashion. While the binding free energy is independent of the ligand dipole moment, it shows a strong correlation with the ligand molecular polarizability.
Jiao, D., J. J. Zhang, R. E. Duke, G. H. Li, M. J. Schnieders and P. Ren, Trypsin-ligand binding free energies from explicit and implicit solvent simulations with polarizable potential. Journal of Computational Chemistry 2009, 30, (11), 1701-11.: We have calculated the binding free energies of a series of benzamidine-like inhibitors to trypsin with a polarizable force field using both explicit and implicit solvent approaches. Free energy perturbation has been performed for the ligands in bulk water and in protein complex with molecular dynamics simulations. The binding free energies calculated from explicit solvent simulations are well within the accuracy of experimental measurement and the direction of change is predicted correctly in all cases. We analyzed the molecular dipole moments of the ligands in gas, water and protein environments. Neither binding affinity nor ligand solvation free energy in bulk water shows much dependence on the molecular dipole moments of the ligands. Substitution of the aromatic or the charged group in the ligand results in considerable change in the solvation energy in bulk water and protein whereas the binding affinity varies insignificantly due to cancellation. The effect of chemical modification on ligand charge distribution is mostly local. Replacing benzene with diazine has minimal impact on the atomic multipoles at the amidinium group. We have also utilized an implicit solvent based end-state approach to evaluate the binding free energies of these inhibitors. In this approach, the polarizable multipole model combined with Poisson-Boltzmann/surface area (PMPB/SA) provides the electrostatic interaction energy and the polar solvation free energy. Overall the relative binding free energies obtained from the MM-PMPB/SA model are in good agreement with the experimental data.

2007

Schnieders, M. J. and J. W. Ponder, Polarizable atomic multipole solutes in a generalized Kirkwood continuum. Journal of Chemical Theory and Computation 2007, 3, (6), 2083-97.: The generalized Born (GB) model of continuum electrostatics is an analytic approximation to the Poisson equation useful for predicting the electrostatic component of the solvation free energy for solutes ranging in size from small organic molecules to large macromolecular complexes. This work presents a new continuum electrostatics model based on Kirkwood's analytic result for the electrostatic component of the solvation free energy for a solute with arbitrary charge distribution. Unlike GB, which is limited to monopoles, our generalized Kirkwood (GK) model can treat solute electrostatics represented by any combination of permanent and induced atomic multipole moments of arbitrary degree. Here we apply the GK model to the newly developed Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA) force field, which includes permanent atomic multipoles through the quadrupole and treats polarization via induced dipoles. A derivation of the GK gradient is presented, which enables energy minimization or molecular dynamics of an AMOEBA solute within a GK continuum. For a series of 55 proteins, GK electrostatic solvation free energies are compared to the Polarizable Multipole Poisson−Boltzmann (PMPB) model and yield a mean unsigned relative difference of 0.9%. Additionally, the reaction field of GK compares well to that of the PMPB model, as shown by a mean unsigned relative difference of 2.7% in predicting the total solvated dipole moment for each protein in this test set. The CPU time needed for GK relative to vacuum AMOEBA calculations is approximately a factor of 3, making it suitable for applications that require significant sampling of configuration space.
Schnieders, M. J., N. A. Baker, P. Ren and J. W. Ponder, Polarizable atomic multipole solutes in a Poisson-Boltzmann continuum. Journal of Chemical Physics 2007, 126, (12), 124114.: Modeling the change in the electrostatics of organic molecules upon moving from vacuum into solvent, due to polarization, has long been an interesting problem. In vacuum, experimental values for the dipole moments and polarizabilities of small, rigid molecules are known to high accuracy; however, it has generally been difficult to determine these quantities for a polar molecule in water. A theoretical approach introduced by Onsager [ J. Am. Chem. Soc. 58, 1486 (1936) ] used vacuum properties of small molecules, including polarizability, dipole moment, and size, to predict experimentally known permittivities of neat liquids via the Poisson equation. Since this important advance in understanding the condensed phase, a large number of computational methods have been developed to study solutes embedded in a continuum via numerical solutions to the Poisson-Boltzmann equation. Only recently have the classical force fields used for studying biomolecules begun to include explicit polarization in their functional forms. Here the authors describe the theory underlying a newly developed polarizable multipole Poisson-Boltzmann (PMPB) continuum electrostatics model, which builds on the atomic multipole optimized energetics for biomolecular applications (AMOEBA) force field. As an application of the PMPB methodology, results are presented for several small folded proteins studied by molecular dynamics in explicit water as well as embedded in the PMPB continuum. The dipole moment of each protein increased on average by a factor of 1.27 in explicit AMOEBA water and 1.26 in continuum solvent. The essentially identical electrostatic response in both models suggests that PMPB electrostatics offers an efficient alternative to sampling explicit solvent molecules for a variety of interesting applications, including binding energies, conformational analysis, and pKa prediction. Introduction of 150 mM salt lowered the electrostatic solvation energy between 2 and 13 kcal/mole, depending on the formal charge of the protein, but had only a small influence on dipole moments.