Methods for Efficiently and Accurately Computing Quantum Mechanical Free Energies for Enzyme Catalysis

Enzyme activity is inherently linked to free energies of transition states, ligand binding, protonation/deprotonation, etc.; these free energies, and thus enzyme function, can be affected by residue mutations, allosterically induced conformational changes, and much more. Therefore, being able to predict free energies associated with enzymatic processes is critical to understanding and predicting their function. Free energy simulation (FES) has historically been a computational challenge as it requires both the accurate description of inter- and intramolecular interactions and adequate sampling of all relevant conformational degrees of freedom. The hybrid quantum mechanical molecular mechanical (QM/MM) framework is the current tool of choice when accurate computations of macromolecular systems are essential. Unfortunately, robust and efficient approaches that employ the high levels of computational theory needed to accurately describe many reactive processes (ie, ab initio, DFT), while also including explicit solvation effects and accounting for extensive conformational sampling are essentially nonexistent. In this chapter, we will give a brief overview of two recently developed methods that mitigate several major challenges associated with QM/MM FES: the QM non-Boltzmann Bennett's acceptance ratio method and the QM nonequilibrium work method. We will also describe usage of these methods to calculate free energies associated with (1) relative properties and (2) along reaction paths, using simple test cases with relevance to enzymes examples.

CHARMM: the biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

Simulation studies of the protein-water interface. II. Properties at the mesoscopic resolution

We report molecular dynamics (MD) simulations of three protein-water systems (ubiquitin, apo-calbindin D(9K), and the C-terminal SH2 domain of phospholipase C-gamma1), from which we compute the dielectric properties of the solutions. Since two of the proteins studied have a net charge, we develop the necessary theory to account for the presence of charged species in a form suitable for computer simulations. In order to ensure convergence of the time correlation functions needed for the analysis, the minimum length of the MD simulations was 20 ns. The system sizes (box length, number of waters) were chosen so that the resulting protein concentrations are comparable to experimental conditions. A dielectric component analysis was carried out to analyze the contributions from protein and water to the frequency-dependent dielectric susceptibility chi(omega) of the solutions. Additionally, an even finer decomposition into protein, two solvation shells, and the remaining water (bulk water) was carried out. The results of these dielectric decompositions were used to study protein solvation at mesoscopic resolution, i.e., in terms of protein, first and second solvation layers, and bulk water. This study, therefore, complements the structural and dynamical analyses at molecular resolution that are presented in the companion paper. The dielectric component contributions from the second shell and bulk water are very similar in all three systems. We find that the proteins influence the dielectric properties of water even beyond the second solvation shell, in agreement with what was observed for the mean residence times of water molecules in protein solutions. By contrast, the protein contributions, as well as the contributions of the first solvation shell, are system specific. Most importantly, the protein and the first water shell around ubiquitin and apo-calbindin are anticorrelated, whereas the first water shell around the SH2 domain is positively correlated.

Simulation studies of the protein-water interface. I. Properties at the molecular resolution

We report molecular dynamics simulations of three globular proteins: ubiquitin, apo-calbindin D(9K), and the C-terminal SH2 domain of phospholipase C-gamma1 in explicit water. The proteins differ in their overall charge and fold type and were chosen to represent to some degree the structural variability found in medium-sized proteins. The length of each simulation was at least 15 ns, and larger than usual solvent boxes were used. We computed radial distribution functions, as well as orientational correlation functions about the surface residues. Two solvent shells could be clearly discerned about charged and polar amino acids. Near apolar amino acids the water density near such residues was almost devoid of structure. The mean residence time of water molecules was determined for water shells about the full protein, as well as for water layers about individual amino acids. In the dynamic properties, two solvent shells could be characterized as well. However, by comparison to simulations of pure water it could be shown that the influence of the protein reaches beyond 6 A, i.e., beyond the first two shells. In the first shell (r < or =3.5 A), the structural and dynamical properties of solvent waters varied considerably and depended primarily on the physicochemical properties of the closest amino acid side chain, with which the waters interact. By contrast, the solvent properties seem not to depend on the specifics of the protein studied (such as the net charge) or on the secondary structure element in which an amino acid is located. While differing considerably from the neat liquid, the properties of waters in the second solvation shell (3.5< r < or =6 A) are rather uniform; a direct influence from surface amino acids are already mostly shielded.

Comparative modeling of GABA(A) receptors: limits, insights, future developments

GABA(A) receptors are chloride ion channels that mediate fast synaptic transmission and belong to a superfamily of pentameric ligand-gated ion channels. The recently published crystal structure of the acetylcholine binding protein can be used as a template for comparative modeling of the extracellular domain of GABA(A) receptors. In this commentary, difficulties with comparative modeling at low sequence identity are discussed, the degree of structural conservation to be expected within the superfamily is analyzed and numerical estimates of model uncertainties in functional regions are provided. Topography of the binding sites at subunit-interfaces is examined and possible targets for rational mutagenesis studies are suggested. Allosteric motions are considered and a mechanism for mediation of positive cooperativity at the benzodiazepine site is proposed.

Towards a better description and understanding of biomolecular solvation

We introduce a flexible framework for the correct description of the solvation of biological macromolecules, the dielectric field equation (DFE). The formalism permits the use of any combination of quantum mechanical (QM), molecular mechanical (MM) and continuum electrostatic (CE) based techniques. For the CE region a method that yields the electric field rather than the potential is outlined. The DFE formalism makes clear the need to consider and to calibrate a dielectric boundary region surrounding the simulation system. The details of how to do this are presented for the case of the Ewald summation method; the effects are demonstrated by calculations of the dielectric properties and the spatially resolved Kirkwood G-factor, G(K)(r), of TIP3P water. Computing the dielectric properties of a multi-component system provides a sensitive method to better understand the solvation of biological macromolecules. Towards this goal a rigorous analysis of the dielectric properties of solvated biomolecules based on a decomposition of the frequency-dependent dielectric constant (or susceptibility) of the full system is presented. The meaning of our approach is investigated, and the results of a first application are reported. Using the method of Voronoi polyhedra, the dielectric properties of the first two solvation shells and bulk water are compared by re-analyzing a 12-ns trajectory of a zinc finger peptide in water [Löffler et al. J. Mol. Biol. 270 (1997) 520]. It is found that the first shell behaves considerably different; in addition, there is a non-negligible contribution to the total susceptibility of the system from coupling between the protein and the bulk water phase, i.e. the water molecules not in the immediate vicinity of the solute.

The meaning of component analysis: decomposition of the free energy in terms of specific interactions

Free energy simulations are of particular interest for the interpretation of macroscopic data in terms of microscopic interactions. This can be done by expressing calculated free energies as a sum of components that correspond to the contributions of different energy terms or different parts of the system. Since the resulting components depend on the integration path, care is required for their use. We show that a linear coupling scheme for the alchemical creation of a chemical identity corresponds to a particularly useful path because it leads to a symmetric decoupling of the free energy components. The path dependence also provides an additional degree of freedom that can be used to study different processes. This latter point is illustrated by a reinterpretation of a recent simulation on wild-type and mutant azurin by Mark and van Gunsteren.

Free energy simulations: the meaning of the individual contributions from a component analysis

A theoretical analysis is made of the decomposition into contributions from individual interactions of the free energy calculated by thermodynamic integration. It is demonstrated that such a decomposition, often referred to as "component analysis," is meaningful, even though it is a function of the integration path. Moreover, it is shown that the path dependence can be used to determine the relation of the contribution of a given interaction to the state of the system. To illustrate these conclusions, a simple transformation (Cl- to Br- in aqueous solution) is analyzed by use of the Reference Interaction Site Model-Hypernetted Chain Closure integral equation approach; it avoids the calculational difficulties of macromolecular simulation while retaining their conceptual complexity. The difference in the solvation free energy between chloride and bromide is calculated, and the contributions of the Lennard-Jones and electrostatic terms in the potential function are analyzed by the use of suitably chosen integration paths. The model is also used to examine the path dependence of individual contributions to the double free energy differences (delta delta G or delta delta A) that are often employed in free energy simulations of biological systems. The alchemical path, as contrasted with the experimental path, is shown to be appropriate for interpreting the effects of mutations on ligand binding and protein stability. The formulation is used to obtain a better understanding of the success of the Poisson-Boltzmann continuum approach for determining the solvation properties of polar and ionic systems.