Motivated by the long-term goal of understanding vectorial biological processes such as proton transport (PT) in biomolecular ion pumps, a number of developments were made to establish combined quantum mechanical/molecular mechanical (QM/MM) methods suitable for studying chemical reactions involving significant charge separation in the condensed phase. These developments were summarized and discussed with representative problems. Specifically, free energy perturbation and boundary potential methods for treating long-range electrostatics were implemented to test the robustness of QM/MM results for protein systems. It was shown that consistent models with sufficient sampling were able to produce quantitatively satisfactory results, such as pK a for titritable groups in the interior of T4-lysozyme, while an inconsistent treatment of electrostatics or lack of sufficient sampling may produce incorrect results. Modifications were made to an approximate density functional theory (SCC-DFTB) to improve the description of proton affinity and hydrogen-bonding, which are crucial for the treatment of PT in polar systems. Test calculations on water autoionization showed clearly that both improvements are necessary for quantitatively reliable results. Finally, the newly established SCC-DFTB/MM-GSBP protocol was used to explore mechanistic issues in carbonic anhydrase (CA). Preliminary results suggest that PT in CA occurs mainly through short water wires containing two water molecules in a thermally activated fashion. Although longer water wires occur with similar frequencies, PT along those pathways, on average, has substantially higher barriers, a result not expected based on previous studies. The fluctuations of water molecules peripheral to the water wire were found to make a larger impact on the PT energetics compared to polar protein residues in the active site, which are largely pre-organized and therefore have less tendency to reorganize during the reaction.

The accuracy of biological simulations depends, in large part, on the treatment of electrostatics. Due to the availability of accurate experimental values, calculation of pKa provides stringent evaluation of computational methods. The generalized solvent boundary potential (GSBP) and Ewald summation electrostatic treatments were recently implemented for combined quantum mechanical and molecular mechanics (QM/MM) simulations by our group. These approaches were tested by calculating pKa shifts due to differences in electronic structure and electrostatic environment; the shifts were determined for a series of small molecules in solution, using various electrostatic treatments, and two residues (His 31, Lys 102) in the M102K T4-lysozyme mutant with large pKa shifts, using the GSBP approach. The calculations utilized a free energy perturbation scheme with the QM/MM potential function involving the self-consistent charge density functional tight binding (SCC-DFTB) and CHARMM as the QM and MM methods, respectively. The study of small molecules demonstrated that inconsistent electrostatic models produced results that were difficult to correct in a robust manner; by contrast, extended electrostatics, GSBP, and Ewald simulations produced consistent results once a bulk solvation contribution was carefully chosen. In addition to the electrostatic treatment, the pKa shifts were also sensitive to the level of the QM method and the scheme of treating QM/MM Coulombic interactions; however, simple perturbative corrections based on SCC-DFTB/CHARMM trajectories and higher level single point energy calculations were found to give satisfactory results. Combining all factors gave a root-mean-square difference of 0.7 pKa units for the relative pKa values of the small molecules compared to experiment. For the residues in the lysozyme, an accurate pKa shift was obtained for His 31 with multiple nanosecond simulations. For Lys 102, however, the pKa shift was estimated to be too large, even after more than 10 nanosecond simulations for each lambda window; the difficulty was due to the significant, but slow, reorganization of the protein and water structure when Lys 102 was protonated. The simulations support that Lys 102 is deprotonated in the X-ray structure and the protein is highly destabilized when this residue is protonated.

A robust approach for dealing with electrostatic interactions for spherical boundary conditions has been implemented in the QM/MM framework. The development was based on the generalized solvent boundary potential (GSBP) method proposed by Im et al. [J. Chem. Phys. 114, 2924 (2001)], and the specific implementation was applied to the self-consistent-charge density-functional tight-binding approach as the quantum mechanics (QM) level, although extension to other QM methods is straightforward. Compared to the popular stochastic boundary-condition scheme, the new protocol offers a balanced treatment between quantum mechanics/molecular mechanics (QM/MM) and MM/MM interactions; it also includes the effect of the bulk solvent and macromolecule atoms outside of the microscopic region at the Poisson-Boltzmann level. The new method was illustrated with application to the enzyme human carbonic anhydrase II and compared to stochastic boundary-condition simulations using different electrostatic treatments. The GSBP-based QM/MM simulations were most consistent with available experimental data, while conventional stochastic boundary simulations yielded various artifacts depending on different electrostatic models. The results highlight the importance of carefully treating electrostatics in QM/MM simulations of biomolecules and suggest that the commonly used truncation schemes should be avoided in QM/MM simulations, especially in simulations that involve extensive conformational samplings. The development of the GSBP-based QM/MM protocol has opened up the exciting possibility of studying chemical events in very complex biomolecular systems in a multiscale framework.

No abstract