Pär Söderhjelm

On the Convergence of QM/MM Energies

We have studied the convergence of QM/MM calculations with respect to the size of the QM system. We study a proton transfer between a first-sphere cysteine ligand and a second-sphere histidine group in [Ni,Fe] hydrogenase and use a 446-atom model of the protein, treated purely with QM methods as a reference. We have tested 12 different ways to redistribute charges close to the junctions (to avoid overpolarization of the QM system), but once the junctions are moved away from the active site, there is little need to redistribute the charges. We have tested 13 different variants of QM/MM approaches, including two schemes to correct errors caused by the truncation of the QM system. However, we see little gain from such correction schemes; on the contrary, they are sensitive to the charge-redistribution scheme and may cause large errors if charges are close to the junctions. In fact, the best results were obtained with a mechanical embedding approach that does not employ any correction scheme and ignores polarization. It gives a mean unsigned error for 40 QM systems of different sizes of 7 kJ/mol with a maximum error of 28 kJ/mol. The errors can be significantly decreased if bonds between the QM and MM system (junctions) are moved one residue away from all active-site residues. Then, most QM/MM variants give mean unsigned errors of 5-9 kJ/mol, maximum errors of 16-35 kJ/mol, and only five to seven residues give an error of over 5 kJ/mol. In general, QM/MM calculations converge faster with system size than pure QM calculations.

How Accurate Can a Force Field Become? A Polarizable Multipole Model Combined with Fragment-wise Quantum-Mechanical Calculations

A new method to accurately estimate the interaction energy between a large molecule and a smaller ligand is presented. The method approximates the electrostatic and induction contributions classically by multipole and polarizability expansions, but uses explicit quantum-mechanical fragment calculations for the remaining (nonclassical) contributions, mainly dispersion and exchange repulsion. Thus, it represents a limit of how accurate a force field can ever become for interaction energies if pairwise additivity of the nonclassical term is assumed (e.g., all general-purpose force fields). The accuracy is tested by considering protein-ligand model systems for which the true MP2/6-31G* interaction energies can be computed. The method is shown to be more accurate than related fragmentation approaches. The remaining error (2-5 and approximately10 kJ/mol for neutral and charged ligands, respectively) can be decreased by including the polarizing effect from surrounding fragments in the quantum-mechanical calculations.

Ligand Affinities Estimated by Quantum Chemical Calculations

We present quantum chemical estimates of ligand-binding affinities performed, for the first time, at a level of theory for which there is a hope that dispersion and polarization effects are properly accounted for (MP2/cc-pVTZ) and at the same time effects of solvation, entropy, and sampling are included. We have studied the binding of seven biotin analogues to the avidin tetramer. The calculations have been performed by the recently developed PMISP approach (polarizable multipole interactions with supermolecular pairs), which treats electrostatic interactions by multipoles up to quadrupoles, induction by anisotropic polarizabilities, and nonclassical interactions (dispersion, exchange repulsion, etc.) by explicit quantum chemical calculations, using a fragmentation approach, except for long-range interactions that are treated by standard molecular-mechanics Lennard-Jones terms. In order to include effects of sampling, 10 snapshots from a molecular dynamics simulation are studied for each biotin analogue. Solvation energies are estimated by the polarized continuum model (PCM), coupled to the multipole-polarizability model. Entropy effects are estimated from vibrational frequencies, calculated at the molecular mechanics level. We encounter several problems, not previously discussed, illustrating that we are first to apply such a method. For example, the PCM model is, in the present implementation, questionable for large molecules, owing to the use of a surface definition that gives numerous small cavities in a protein.

Ligand-Binding Affinity Estimates Supported by Quantum-Mechanical Methods.

One of the largest challenges of computational chemistry is calculation of accurate free energies for the binding of a small molecule to a biological macromolecule, which has immense implications in drug development. It is well-known that standard molecular-mechanics force fields used in most such calculations have a limited accuracy. Therefore, there has been a great interest in improving the estimates using quantum-mechanical (QM) methods. We review here approaches involving explicit QM energies to calculate binding affinities, with an emphasis on the methods, rather than on specific applications. Many different QM methods have been employed, ranging from semiempirical QM calculations, via density-functional theory, to strict coupled-cluster calculations. Dispersion and other empirical corrections are mandatory for the approximate methods, as well as large basis sets for the stricter methods. QM has been used for the ligand, for a few crucial groups around the ligand, for all the closest atoms (200-1000 atoms), or for the full receptor-ligand complex, but it is likely that with a proper embedding it might be enough to include all groups within ∼6 Å of the ligand. Approaches involving minimized structures, simulations of the end states of the binding reaction, or full free-energy simulations have been tested.

On the performance of quantum chemical methods to predict solvatochromic effects: The case of acrolein in aqueous solution

The performance of the Hartree-Fock method and the three density functionals B3LYP, PBE0, and CAM-B3LYP is compared to results based on the coupled cluster singles and doubles model in predictions of the solvatochromic effects on the vertical n-->pi* and pi-->pi* electronic excitation energies of acrolein. All electronic structure methods employed the same solvent model, which is based on the combined quantum mechanics/molecular mechanics approach together with a dynamical averaging scheme. In addition to the predicted solvatochromic effects, we have also performed spectroscopic UV measurements of acrolein in vapor phase and aqueous solution. The gas-to-aqueous solution shift of the n-->pi* excitation energy is well reproduced by using all density functional methods considered. However, the B3LYP and PBE0 functionals completely fail to describe the pi-->pi* electronic transition in solution, whereas the recent CAM-B3LYP functional performs well also in this case. The pi-->pi* excitation energy of acrolein in water solution is found to be very dependent on intermolecular induction and nonelectrostatic interactions. The computed excitation energies of acrolein in vacuum and solution compare well to experimental data.

Accurate Predictions of Nonpolar Solvation Free Energies Require Explicit Consideration of Binding-Site Hydration

Continuum solvation methods are frequently used to increase the efficiency of computational methods to estimate free energies. In this paper, we have evaluated how well such methods estimate the nonpolar solvation free-energy change when a ligand binds to a protein. Three different continuum methods at various levels of approximation were considered, viz., the polarized continuum model (PCM), a method based on cavity and dispersion terms (CD), and a method based on a linear relation to the solvent-accessible surface area (SASA). Formally rigorous double-decoupling thermodynamic integration was used as a benchmark for the continuum methods. We have studied four protein-ligand complexes with binding sites of varying solvent exposure, namely the binding of phenol to ferritin, a biotin analogue to avidin, 2-aminobenzimidazole to trypsin, and a substituted galactoside to galectin-3. For ferritin and avidin, which have relatively hidden binding sites, rather accurate nonpolar solvation free energies could be obtained with the continuum methods if the binding site is prohibited to be filled by continuum water in the unbound state, even though the simulations and experiments show that the ligand replaces several water molecules upon binding. For the more solvent exposed binding sites of trypsin and galectin-3, no accurate continuum estimates could be obtained, even if the binding site was allowed or prohibited to be filled by continuum water. This shows that continuum methods fail to give accurate free energies on a wide range of systems with varying solvent exposure because they lack a microscopic picture of binding-site hydration as well as information about the entropy of water molecules that are in the binding site before the ligand binds. Consequently, binding affinity estimates based upon continuum solvation methods will give absolute binding energies that may differ by up to 200 kJ/mol depending on the method used. Moreover, even relative energies between ligands with the same scaffold may differ by up to 75 kJ/mol. We have tried to improve the continuum solvation methods by adding information about the solvent exposure of the binding site or the hydration of the binding site, and the results are promising at least for this small set of complexes.

Protein Influence on Electronic Spectra Modeled by Multipoles and Polarizabilities.

We have developed automatic methods to calculate multipoles and anisotropic polarizabilities for all atoms and bond centers in a protein and to include such a model in the calculation of electronic properties at any level of quantum mechanical theory. This approach is applied for the calculation of the electronic spectra of retinal in rhodopsin at the CASPT2//CASSCF level (second-order multiconfigurational perturbation theory) for the wild-type protein, as well as two mutants and isorhodopsin in QM/MM structures based on two crystal structures. We also perform a detailed investigation of the importance and distance dependence of the multipoles and the polarizabilities for both the absolute and the relative absorption energies. It is shown that the model of the surrounding protein strongly influences the spectrum and that different models give widely different results. For example, the Amber 1994 and 2003 force fields give excitation energies that differ by up to 16 kJ/mol. For accurate excitation energies, multipoles up to quadrupoles and anisotropic polarizabilities are needed. However, interactions with residues more than 10 Å from the chromophore can be treated with a standard polarizable force field without any dipoles or quadrupoles.

Accurate Reaction Energies in Proteins Obtained by Combining QM/MM and Large QM Calculations

We here suggest and test a new method to obtain stable energies in proteins for charge-neutral reactions by running large quantum mechanical (QM) calculations on structures obtained by combined QM and molecular mechanics (QM/MM) geometry optimization on several snapshots from molecular dynamics simulations. As a test case, we use a proton transfer between a metal-bound cysteine residue and a second-sphere histidine residue in the active site of [Ni,Fe] hydrogenase, which has been shown to be very sensitive to the surroundings. We include in the QM calculations all residues within 4.5 Å of the active site, two capped residues on each side of the active-site residues, and all charged groups that are buried inside the protein, which for this enzyme includes three iron-sulfur clusters, in total, 930 atoms. These calculations are performed at the BP86/def2-SV(P) level, but the energies are then extrapolated to the B3LYP/def2-TZVP level with a smaller QM system, and zero-point energy, entropy, and thermal effects are added. We test three approaches to model the remaining atoms of the protein solvent, viz., by standard QM/MM approaches using either mechanical or electrostatic embedding or by using a continuum solvation model for the large QM systems. Quite encouragingly, the three approaches give the same results within 14 kJ/mol, and variations in the size of the QM system do not change the energies by more than 8 kJ/mol, provided that the QM/MM junctions are not moved closer to the QM system. The statistical precision for the average over 10 snapshots is 1-3 kJ/mol.

Calculation of Protein-Ligand Interaction Energies by a Fragmentation Approach Combining High-Level Quantum Chemistry with Classical Many-Body Effects

We have developed a method to estimate accurate interaction energies between a full protein and a bound ligand. It is based on the recently proposed PMISP (polarizable multipole interaction with supermolecular pairs) method (Soderhjelm, P.; Ryde, U. J. Phys. Chem. A 2009, 113, 617), which treats electrostatic interaction by multipoles up to quadrupoles, induction by anisotropic polarizabilities, and nonclassical interactions by explicit quantum mechanical (QM) calculations, using a fragmentation approach. For a whole protein, electrostatics and induction are treated the same way, but for the nonclassical interactions, a Lennard-Jones term from a standard molecular mechanics (MM) force field (e.g., Amber) is used outside a certain distance from the ligand (4-7 A). This QM/MM variant of the PMISP method is carefully tested by varying this distance. Several approximations related to the classical interactions are also evaluated. It is found that one can speed up the calculation by using density functional theory to compute multipoles and polarizabilities but that a proper treatment of polarization is important. As a demonstration of the method, the interaction energies of two ligands bound to avidin are calculated at the MP2/aug-cc-pVTZ level, with an expected relative error of 1-2%.

Effect of geometry optimisations on QM-cluster and QM/MM studies of reaction energies in proteins

We have examined the effect of geometry optimization on energies calculated with the quantum mechanical (QM) cluster, combined QM and molecular mechanics (QM/MM), and big-QM approaches (very large single-point QM calculations taken from QM/MM-optimized structures, including all atoms within 4.5 Å of the minimal active site, all buried charged groups in the protein, and truncations moved at least three residues away from the active site). We studied a simple proton-transfer reaction between His-79 and Cys-546 in the active site of [Ni,Fe] hydrogenase and optimize QM systems of 50 different sizes (56-362 atoms). Geometries optimized with QM/MM are stable and reliable, whereas QM-cluster optimizations give larger changes in the structures and sometimes lead to large distortions in the active site if some hydrogen-bond partners to the metal ligands are omitted. Keeping 2-3 atoms for each truncated residue (rather than one) fixed in the optimization improves the results but does not solve all problems for the QM-cluster optimizations. QM-cluster energies in vacuum and a continuum solvent are insensitive to the geometry optimizations with a mean absolute change upon the optimizations of only 4-7 kJ/mol. This shows that geometry optimizations do not decrease the dependence of QM-cluster energies on how the QM system is selected; there is still a ∼60 kJ/mol difference between calculations in which groups have been added to the QM system according to their distance to the active site or based on QM/MM free-energy components. QM/MM energies do not show such a difference, but they converge rather slowly with respect to the size of the QM system, although the convergence is improved by moving truncations away from the active site. The big-QM energies are stable over the 50 different optimized structures, 57 ± 1 kJ/mol, although some smaller trends can be discerned. This shows that both QM-cluster geometries and energies should be interpreted with caution. Instead, we recommend QM/MM for geometry optimizations and energies calculated by the big-QM approach.