Challenges and Advances in Computational Chemistry and Physics, 2010, 247-266

Published by Springer

The original publication is available at www.springerlink.com

Volume 12, Kinetics and Dynamics: From Nano- to Bio-Scale

Eds. P. Paneth, and A. Dybala-Defratyka

Chapter 9. Quantum-Mechanical Molecular Dynamics of Charge Transfer

Victor M. Anisimov, Claudio N. Cavasotto

School of Health Information Sciences, University of Texas at Houston, Houston, 7000 Fannin St, TX 77030

Abstract

Computational studies of biological macromolecules are challenging due to large size of biomolecules, their conformational flexibility, and the need in explicit water solvation in order to simulate conditions close to experiment. Under these circumstances studying molecular systems via quantum-mechanical calculations becomes exceedingly difficult. Natural is the attempt to reduce the complex quantum-mechanical picture to a more tractable one by accommodating classical-mechanical principles. However, the simplified models may overlook important physics details of atomic interactions. To avoid such potential pitfalls higher level of theory methods should be available to conduct validation studies. Using semiempirical linear scaling quantum-mechanical LocalSCF method we performed molecular dynamics simulation of ubiquitin in explicit water. The simulation revealed various deviations from the classical mechanics picture. The average charge on amino acids varied depending on their environment. We observed charge transfer channels transmitting electric charge between amino acids in sync with protein motion. We also noticed that the excess charge transferred from protein to water creates a charge cloud around the protein. The observed global dynamic effects of charge transfer represent a new previously unaccounted degree of freedom of biomolecules which requires QM treatment in order to obtain more accurate dynamics of biomolecules at atomic resolution.

Keywords

NDDO method, PM5 Hamiltonian, ubiquitin, water droplet, spherical boundary potential, QM MD, charge transfer, VFL approximation, LocalSCF, linear scaling.

Abbreviations

AM1 – Austin model 1; DFT – density functional theory; HF/6-31G* - Hartree-Fock method using Pople 6-31G* basis set; LocalSCF – local self consistent field; MD – molecular dynamics; NDDO – neglect of diatomic differential overlap; PM3 – parametric method 3; PM5 – parametric method 5; QM – quantum mechanics; RM1 – Recife model 1; SCF – self-consistent field; VFL – variational finite localized molecular orbital approximation.

9.1 Introduction

Despite of their computational cost quantum-mechanical (QM) methods offer many important potential uses in computational biology. The data obtained from QM calculations can be used to validate and assist in improving empirical force fields. QM methods provide a possibility to study problems which are beyond the reach of classical models, particularly the processes accompanied by significant redistribution of electronic density e.g. enzymatic reactions and photochemical processes, to name a few. However the progress of such applications is limited by severe performance bottleneck of QM methods. Going to progressively larger systems the resource requirement of QM calculations quickly becomes prohibitive. The resource requirement scales as power of 4 of the number of basis functions for Hartree-Fock and Density Functional Theory (DFT) methods and as power of 5 for MP2 method. Coupled cluster calculations using singles and doubles scale as power of 6. Therefore depending on the level of theory the conventional QM methods are practically limited to systems containing 100 - 1000 atoms, whereas typical biological systems including solvent encounter tens of thousands of atoms minimum.

Performance of QM methods has recently been significantly improved due to the development of linear scaling methods.[1,2] In the linear scaling regimen the computation time increases linearly of the number of employed basis functions. Now when the scalability problem seems to be solved it might look like the road for QM methods to biological application is eventually going to be wide open. However, making such significant step forward was only to discover even more impenetrable wall erected by structural flexibility of biological systems.

Since the pioneering work of McCammon et al.[3] it has long been recognized that biomolecules are dynamic entities and that their biological structure and function are linked to their thermal motions. This makes insufficient the standard method of attack of quantum chemistry in the form of geometry optimization which was developed for and very efficiently applied to small molecule problems. Applying geometry optimization to large and flexible biomacromolecules is going to give information only about one of the million of energetically possible conformational local minima. Collecting such limited data about few conformations is obviously insufficient to make computational predictions about macroscopic properties of large biomolecular systems.

This problem is known as configurational sampling and it requires application of molecular dynamics (MD). According to statistical mechanics the complexity of free energy profiling of a flexible molecular system grows exponentially with the number of atoms in the system. After facing this sobering fact it becomes apparent that the further expansion of QM methods toward biological problems should go via development of computationally inexpensive approximate QM methods. Semiempirical QM methods based on neglect of diatomic differential overlap (NDDO) are particularly attractive in this regard. They are about 1000-fold faster than chemically accurate DFT methods. Correspondingly, same order of the resource saving is expected during MD simulation based on the given level of QM theory. Therefore semiempirical QM methods provide the first outpost on the long journey to bring QM platform to solution of biological problems.

QM/MM is another example of the approximate platform where a most critical part of the system is treated at QM level whereas the least important part of the system is treated by classical force fields. Typical applications of QM/MM include enzymatic reactions.[4] The QM/MM methods are not necessarily limited to semiempirical theory; quite contrarily selecting small QM area makes feasible using a higher level of theory, usually DFT methods. Although useful this approach cannot be a substitute for all-atom QM treatment. For instance, QM/MM cannot be used to validate and further improve force fields or to further advance our knowledge about structure and dynamics of biomacromolecules, the task which requires consistent application of a single level of theory to the entire molecular system.

Resorting to semiempirical formalism does not exhaust the potential of approximate QM methods. Additional performance can be obtained by reformulating the large-scale problem. Such method was proposed by Anikin et al.[5] suggesting to seek the molecular density matrix in the basis of constrained molecular orbitals. The complete account of the theory of the variational finite localized molecular orbital (VFL) approximation and the linear scaling LocalSCF method was given by Anisimov et al.[6]

9.2 Theoretical part

In LocalSCF method and in the underlying VFL approximation the solution to large-scale molecular problem is sought via decomposition of the potential energy of the system, E, as function of density matrix, P, in Taylor series and neglecting the second and higher order terms δ⁽²⁾

     (9.1)

Such expansion given by eq. 9.1 basically implies that if we know an approximate density matrix, P(0), and the perturbation is small we can improve the density matrix toward the final solution. This is effectively utilized in the VFL approximation to maintain orthogonality of compact molecular orbitals (CMO). Considering the non-orthogonality as a small perturbation to the fully orthogonal case we can additionally expand the inverse overlap matrix S^-1 in Taylor series preserving only the linear terms

     (9.2)

This leads to the final expression for the total energy, where quadratic term represents a penalty function which helps minimizing the nonorthogonality between CMOs

     (9.3)

Here, W is penalty weight, I is diagonal unit matrix, C is matrix of linear coefficients of CMOs, and C^t is its transpose.

Starting from the initial atomic orbital (AO) expansion of CMOs the VFL approximation can variationally determine density matrix corresponding to this expansion. During the self consistent field (SCF) refinement of linear coefficients of CMOs the VFL method holds fixed the AO expansion of each CMO. Therefore VFL approximation is the method for variational determination of density matrix under the constraint of fixed CMO expansion. Correspondingly, besides pointing to the compact size of molecular orbitals the CMO abbreviation also stands for constrained molecular orbitals.

Since the VFL method assumes that the CMO expansion is known before starting the calculations and provides no mechanism to determine the optimal expansion, the VFL method alone is insufficient to perform QM calculations on large molecular systems. The complementary task of finding the optimal CMO expansion and the corresponding molecular density matrix is resolved in the linear scaling LocalSCF method, which is based on the VFL approximation. Similar to VFL the LocalSCF method also relies on the perturbation expansion (eq. 9.1) but uses it to determine the most energetically favorable direction for expansion of CMO i on trial AO μ based on the value of energy derivative

     (9.4)

The trial atomic orbital μ will be added to CMO i if the corresponding energy derivative defined by eq. 9.4 is the largest among the derivative values corresponding to other alternative AO expansion directions. The success of accurately determining the most energetically favorable expansion depends on the quality of the available density matrix. Having the density matrix P(0) fully converged assures determining best possible and minimally necessary AO expansion leading to rapidly converging total energy and density matrix via the series of successive expansions.

The localized character of CMOs offers a beneficially restricted choice of possible expansions onto only the nearest atomic centers. Further limiting the CMO expansion to the most energetically favorable AOs helps keeping CMOs maximally localized while providing the best possible improvement to the density matrix. In this way LocalSCF elegantly demonstrates the flexibility advantages of the MO-based linear scaling formalism vs density functional (DFT) approach in dealing with large-scale molecular problems.

Useful in interpretation of the results of quantum-chemical calculations is the concept of partial atomic charges, Q_A, where capital letter A denotes atom index, index i runs over CMOs, and index a corresponds to AOs centered on atom A.

     (9.5)

Throughout this work we will be using Mulliken algorithm as our method for calculation of partial atomic charges (eq. 9.5). In semiempirical methods Mulliken charges reduce to Coulson charges due to the NDDO approximation but referring to them as Mulliken charges is still formally correct and preferable for the sake of theoretical consistency. The Mulliken charges are most natural for the semiempirical framework and unlike other charge derivation schemes do not invoke additional parameters. Besides, since semiempirical methods are parametric ones they are free from the basis set superposition error. Consequentially, the troubles ab initio methods have with the Mulliken population analysis do not extend to semiempirical methods.

9.3 The notion of charge transfer

When considering molecular dynamics of biomolecules it is common to question the need in the QM level of theory when classical force fields are doing seemingly well. This is a rhetorical dilemma of any established field of science. However, it has long been recognized that the fixed-charge approximation staying behind the empirical force fields may not be optimal when considering environments of different polarity.[7] Correspondingly, force field parameters optimized to reproduce properties of a model compound in hydrophobic environment may poorly perform in hydrophilic environment and vice versa. Addressing such problems has long been the focus of polarizable force field development.[8-10] However, despite of the achieved significant progress,[11] the description of induced electrostatic effects at classical mechanics level turned out to be an exceedingly difficult problem, and after two decades of intense efforts additional development is still necessary.[12] Besides, the classical polarizable model captures only a portion of the large spectrum of induced electrostatic effects while the other part – intermolecular charge transfer, an inherently quantum mechanical phenomenon, has no classical mechanics formulation. This is a particularly important aspect to consider in the design of computational methods, since according to Nadig et al.[13] energetic consequence of charge transfer effects is greater than that of polarization effects.

Historically there was much uncertainty about the significance of charge transfer effects in molecular systems. Small systems were too small to show any sensible amount of charge transfer whereas studies on large systems were technically not feasible. However the progress in computer hardware and development of linear scaling methods provided the opportunity to reconsider the previous assumptions. The subject of charge transfer in ion solvation was extensively studied by theoretical[14-16] and experimental methods.[17] Charge transfer in DNA has received particularly strong attention.[18,19] Charge transfer in glycine dimer in water solution was confirmed via DFT MD by Peraro et al.[20] Merz and co-workers[13] predicted -2.0 electron units of charge transfer from a small Cold Shock protein to water by performing single-point PM3 calculations on 100 protein-water snapshots generated by force field MD. Despite of the later finding[21] that single-point calculations and particularly PM3 are not fully reliable for estimation of charge transfer, this work along with other studies[22] fueled the interest to the subject of charge transfer in biological systems.

In another study Komeiji et al.[23] predicted protein-to-water charge transfer  of -0.8 electron units for protein ubiquitin by performing HF/6-31G* single-point calculations over five empirical snapshots. The limited number of configurations considered in their study was compensated by the use of ab initio Hartree-Fock level of theory, which is the highest level of theory so far applied to the problem in question. This result most definitely confirmed the presence of charge transfer in biological systems. Yet two major questions remained – first, what are the implications of charge transfer for biological simulations and, second, how to turn the knowledge about charge transfer into a predictive biomolecular simulation?

Considering the first question it is not entirely clear that the omission of charge transfer of -0.8 electron units by classical force fields is something to worry about. For a typical protein-water system consisting of 10,000 atoms the charge transfer would be at the order of 10^-4 electron units per atom. Such electrostatic perturbation will be entirely lost in the numerical noise of MD simulation. So the significance of charge transfer for biological applications still remains to be demonstrated.

Second, putting temporarily aside the proof and simply assuming the significance of charge transfer effects the derived question would be how to properly revise the computational protocol of biomolecules in order to make the practical use of charge transfer effects? Both above simulations of Merz and Komeiji were performed via single-point QM calculations on snapshots extracted from force field MD trajectories. Unfortunately, such two-step procedure has little or no use for typical biomolecular applications which are usually concerned with sampling a canonical statistical mechanics ensemble in order to relate the computational model to the actual experimental conditions. The QM total energies computed based on the empirical snapshots do not represent any statistical ensemble because of the expected differences between the QM and force field potential energy hypersurfaces. Therefore averaging the results of such QM calculations is not going to produce a value corresponding to any experimental condition. To generate a proper statistical ensemble one has to run a QM MD simulation.

Such QM MD simulation could also bring us closer to the solution of the problem of significance of charge transfer. There is some sort of inconsistency in studying charge transfer effects based on empirical trajectories which are generated under the convention of fixed charges. In computational biology a significance of a phenomenon is necessarily understood via its effect on structure and dynamics of a biomolecule. Putting it forward, if the effect in question does not influence the dynamics trajectory it has likely no significance at all for biomolecular applications. Hence, if charge transfer is important it must be present during the production of the MD trajectory. Since this condition is not satisfied when QM calculations are performed over empirical trajectories, it becomes impossible to determine physical significance of the computed charge transfer. This dilemma can only be resolved by performing QM MD simulation and analyzing the generated data.

9.3.1 QM MD of ubiquitin in explicit water

To clarify the issue of significance of charge transfer in protein-solvent systems we performed a series of QM MD simulations on 12199 atoms system representing ubiquitin in 30 Å water droplet for the duration of 20 ps[21] using AM1,[24] RM1,[25] PM3[26] and PM5[27] Hamiltonians. The size of the water droplet was selected to provide 10 Å water cushion between protein and vacuum. The system boundaries were restrained by application of a spherical boundary potential preventing water molecules from leaving the system. Constant temperature of the system was maintained at 310K by application of Nose-Hoover thermostat.[28] No restraints were applied to bonds involving hydrogen atoms. Integration time step was set to 1 fs. Numerical integration of dynamics equations was performed by using a constant volume – constant temperature algorithm.[29]

Despite of finding of Merz and co-workers that AM1 provides best agreement with correlated ab initio methods on the value of charge transfer[30] this Hamiltonian turned out to be unsuitable for MD simulations in presence of explicit water. AM1 resulted in unrealistically high water density of about 2 g/cm³ and showed an unphysical tendency to water dissociation at very first picosecond of the simulation. From this we conclude that AM1 is not applicable for QM MD simulations of biomolecules in presence of explicit water. This is kind of a sad finding considering the widespread popularity of AM1 in biomolecular applications.[31]

Another Hamiltonian, RM1, representing recent reparameterization of AM1, showed improved results. It produced correct water density and stronger OH bond than that in AM1. Unfortunately the original deficiency of, AM1 related to the underestimated strength of OH bond in water, is still somehow lurking in RM1 leading to water dissociation happening after 5 ps of simulation. Correspondingly, both AM1 and RM1 Hamiltonians were excluded from our further QM MD simulations of ubiquitin in water.

For information, considering the experimental strength of OH bond in water, a spontaneous dissociation of a given water molecule is expected on 11-hours time scale[32] which makes it unobservable in regular MD simulations. Contrarily to AM1 and RM1, the PM3 and PM5 MD simulations were entirely stable and predicted protein-to-water charge transfer of -1.60 and -0.45 electron units, respectively.[21]

Having the data of Komeiji as the reference point (-0.8) we can see now that PM3 overestimates charge transfer as previously noticed by Merz.[30] However, the reparameterized variant, PM5, shows quite promising result. The agreement between PM5 MD and HF/6-31G* is nearly quantitative. There are also reasons to believe that both semiempirical and ab initio levels of theory might come to even a closer agreement. Since the HF/6-31G* simulation was performed on a rigid force field geometry one may expect that the HF/6-31G* result would approach the PM5 MD value given the opportunity to relax the structure at the HF/6-31G* level. In addition to that, a smaller value for charge transfer seems more credible considering the multitude of factors increasing the appearance of charge transfer in ab initio calculations. These could be strained geometry, basis set superposition error, insufficient sampling, missing correlation effects and so on. 

9.3.2 Charge transfer inside protein

In our PM5 MD simulation of formally neutrally charged ubiquitin in water droplet[21] the instantaneous net protein charge fluctuates around the average value of 0.45 electron units (Fig. 9.1). Standard deviation of the protein charge is 0.07 electron units; minimal and maximal values of the protein charge are 0.25 and 0.76 electron units, respectively.

Fig. 9.1. Protein net charge fluctuation during MD simulation. Reproduced from Ref. 21 with permission. Copyright (c) Wiley-VCH Verlag GmbH & Co. KGaA

At this point the obtained results alone do not move us much further in answering the primary question of this study whether charge transfer has or has no any biological significance. Yet a critical step is achieved in making sure that our semiempirical QM MD prediction of protein-to-water charge transfer closely agrees with the reference ab initio value of Komeiji. This assures the predictive capability of our method and supports going into deeper analysis of the accumulated data. The obtained QM MD trajectory provides the necessary means to reconstruct the time evolution of charge transfer taking place inside and outside of the protein and obtaining expected values of physical quantities as the canonical ensemble trajectory average.

9.3.3 Charge transfer channel

When considering charge transfer in macromolecules first step coming to mind is to identify the location of charge donors and recipients. According to the convention adopted by polarizable force field the charge fluctuations are confined to single residues[33] and cannot go beyond their boundaries. Any attempt to permit the charge to freely travel between amino acids would be technically so difficult to implement and to optimize the necessary parameters so such opportunity is even not considered in the modern polarizable force fields. Of course, in reality it is hard to know beforehand the location of charge donors and charge recipients. This makes the charge transfer phenomenon beyond the rich of empirical force fields. Contrarily to that, quantum mechanics permits any imaginable charge distribution in the system and offers a ready mathematical apparatus to find the distribution corresponding to the minimum of potential energy of the system.

In the domain of electrostatic phenomena, protein-to-water charge transfer is only a tip of the iceberg. As our QM MD calculations reveal, all parts of the simulation system are actively involved in various charge transfer processes. Tracking these processes is not all easy. Since all electrons in the system are non-distinguishable we can only hope to find patterns in the way some regions of the system depleting the electron density while other areas simultaneously acquiring the excess density. Following this technique we observed a distinctly localized charge flow channel (Fig. 9.2) formed by amino acids E34 and K11 participating in a salt-bridge interaction.

Fig. 9.2. Synchronous charging and discharging of E34-K11 residues following the pattern of their motion. Horizontal lines across the charge fluctuation profiles denote average charge value. Reproduced from Ref. 21 with permission. Copyright (C) Wiley-VCH Verlag GmbH & Co. KGaA

The top panel in Fig. 9.2 shows net charge fluctuation from -1.11 to -0.61 electron units on glutamic acid residue E34 along the dynamics trajectory. The net charge on its salt-bridge partner lysine K11 varies exactly in anti-phase from +0.60 to +1.12 electron units. In any moment of time this pair carries nearly perfect zero net charge, which helps us to conclude that these two residues are involved in some sort of a space-confined dynamic charge exchange, operating almost independently from the rest of charge transfer processes happening around. During their thermal fluctuation these amino acids pump through their bodies about 0.5 units of electron charge, which is alone equal to the value of protein-to-water charge transfer. Since 0.5 units is a half of the effective charge associated with these titratable amino acids such magnitude of charge fluctuation strongly disagrees with the popular idea of force fields assigning a fixed unit charge to amino acid residues. This also shows limitations of confining charge transfer to the boundaries of individual amino acids.

The non-classical character of this salt bridge is revealed in the way the amino acid charges change. We already mentioned that during dynamics trajectory this pair maintains electric neutrality in the sea of charge fluctuation. Nearly all the electric charge lost by one amino acid goes to its partner. Most remarkably, the charge pump operates synchronously as the distance between amino acids changes (Fig. 9.2). When E34 and K11 come closer to each other the electron density travels from E34 to K11 discharging the so called capacitor formed by the amino acid pair. During the discharging process the energy is released to produce a mechanical work. When distance between the donor and acceptor sites increases the electron density flows in opposite direction from K11 to E34 thus charging up the capacitor and accumulating the potential energy. This complex picture of intermolecular charge exchange directly connects the charge transfer effects with protein structural dynamics making them functionally linked.

The effect discussed above is a distinct quantum-mechanical phenomenon, which is overlooked in typical biomolecular simulations. Yet some may argue that it also can be described by using classical mechanics methods. For instance enclosing the charge transfer channel in a single structural unit within a charge fluctuation polarizable force field scheme may do the work. However, since we do not know beforehand where such charge transfer channels will be located and how to handle the dynamic processes of creating and disrupting the charge transfer channels during biomolecular dynamics the empirical force field treatment of intermolecular charge transfer will not be practically feasible. Although similar problems are successfully handled by reactive force fields[34] they are typically limited to isotropic systems while biomolecules are highly anisotropic. Therefore, QM methods are necessary in order to treat intermolecular charge transfer in biological systems.

The analysis of net charge associated with E34 and K11 residues shows additional deviations from the classical picture. The average charge on E34 and K11 is -0.92 and +0.92 units, respectively, which is by 8% smaller from the unit charges of -1 and +1 assigned to these amino acids by classical force fields. The charge on E34 and K11 amino acids fluctuates around the average with the standard deviation of 0.08 and 0.09, respectively. Similar in magnitude charge fluctuation was earlier observed by Yang and co-workers[35] in their QM/MM treatment of crambin in water where protein was treated quantum mechanically whereas solvent was considered at force field level.

It is apparent that fixing the residue charges to unit values as it is done in empirical force fields would result in a different interaction energy profile between amino acids and would ultimately lead to a different trajectory. It would also result in predicting different mechanical properties of macromolecules. Since mechanical flexibility of proteins is an inherent component of their biological function the correct handling of charge transfer is necessary in order to further improve the accuracy of computer simulations of biomolecules.

Looking at the smooth work of the charge transfer channel E34-K11 it is tempting to suggest a potential biological function of this unit. The mechanic motion of the protein generates an alternating electric current at the E34-K11 salt bridge which in turn produces an infra-red frequency electromagnetic wave. Such electromagnetic wave can be sensed by other charge transfer channels and through that obstruct or facilitate a specific protein motion. In fact the E34-K11 salt bridge in ubiquitin may itself serve as an electromagnetic sensor by converting the absorbed electromagnetic energy into mechanic motion. Such signaling mechanism may control protein-protein binding recognition. This in turn may be a part of a signal transduction mechanism consisting of the cascade of protein-protein binding events.[36] These are of course only abstract thoughts at this moment. Additional studies are necessary in order to prove the feasibility of such mechanism.

9.3.4 Inequality among same-type amino acids

It is common in empirical force fields that same-type amino acids carry same set of parameters - charges and polarizabilities. Therefore despite of saddle differences in their environment same-type amino acids usually exhibit very similar physical response. Such model is ideally suited to describe bulk properties of materials where average response is what matters. In such models large fluctuations of properties at the level of individual blocks are even unwelcomed due to the accompanied computational noise. Such problem already becomes visible in polarizable force fields vs. non-polarizable ones. It takes longer time to obtain same level of convergent properties in polarizable force fields, particularly in weak electrostatic field environments, than in non-polarizable ones.[37] Therefore, because of adopting a simplified electrostatic model empirical force fields are not capable to fully substitute quantum-mechanical treatment when atomic level resolution is necessary.[38] Next we will see what kind of unusual effects one can observe at atomic level.

While performing charge transfer analysis along dynamics trajectory it is reasonable to expect significant changes in charge distribution taking place at electron-rich or electron-deficient sites. To our great surprise, ubiquitin showed strong charge fluctuations happening at the sites of neutral amino acids too. Large charge fluctuation in the range of 0.50 electron units took place at formally neutral amino acid tyrosine Y59 (Fig. 9.3). During the course of the simulation its net charge varied from -0.39 to +0.11. This charge fluctuation could happen due to hydrogen-bond interaction of Y59 with negatively charged glutamic acid E51. The question is can we predict this outcome without performing QM calculation? As a neutrally charged amino acid the tyrosine is not expected to exhibit such large charge fluctuation. Definitely, typical tyrosine in a protein behaves exactly as expected. What we see in ubiquitin is a unique environment of Y59 which forces it to deviate from the average. The ability to pinpoint such details makes QM calculations irreplaceable when considering the ways of improving the predictive ability of computational methods in biology.

Fig. 9.3. Ubiquitin in water droplet. Denoted are amino acids showing the largest charge transfer effects.

Even more unusual is large charge fluctuation of 0.43 electron units happening on amino acid glycine G10. This simple amino acid assumes extreme net charge in the range from -0.20 to +0.23 electron units during the simulation. G10 does not participate in hydrogen-bond interactions with other amino acids but its neighbor K11 could be responsible for inducing the charge fluctuations on G10. Glycine G10 is fully exposed to solvent and its carbonyl oxygen is in dynamic charge exchange with water serving as the charge transfer partner to G10. Since glycine is not known for high polarization response abilities and because force fields treat all glycine amino acids equally, it is likely that the unique electrostatic properties of G10 in ubiquitin would be overlooked by polarizable force fields. According to our QM MD data, other glycine amino acids in ubiquitin behave nothing but ordinarily. This is another example of quantum-mechanical effects playing a significant role in the place where they are least expected.

Besides of showing deviations from the classical-mechanical picture the analysis of QM data also provides some evidences supporting classical model. For the majority of amino acids their average net charge agrees quite well with the classical assignment of unit charges. Same fact was earlier noticed by Yang and co-workers in their QM/MM simulation of crambin.[35] Excellent agreement between average QM and classical charges can be seen on amino acids G10, E16, and K48 in Table 9.1. Particularly interesting is the case of glycine G10. In the light of the above discussion we know that this amino acids exhibits quantum-mechanical effects in the form of unusually large instantaneous charge fluctuations whereas on average G10 looks like a classical neutrally charged amino acid. This is an illustrative example of a difference between the average response and actual atomic level resolution processes. Force fields are definitely good in describing the average properties whereas QM is necessary when simulations are performed to obtain atomic level resolution.

Table 9.1. Average^a and extreme net charges on selected amino acids.

^aData in brackets represent standard deviation.

^bCharge assigned by Force Field.

^cΔq = qQM_ave - qFF

Concerning the differences between the quantum- and classical-mechanical models, according to our QM MD data for 76 ubiquitin amino acids their average charge deviates from the classical values on overall by 3% only. It explains the general success of empirical force fields in computational biology. However, the differences between classical and quantum mechanical pictures increase when looking at the details of individual amino acids. Most noticeable is charge difference on terminal residues C-ter and N-ter, which participate in charge transfer to water. Since these residues are fully exposed to solvent they act almost like free ions.

Even though for the majority of amino acids the classical units charge is a good match, it is not that good for K11, D32, E34, G51, and Y59. These are the examples where average charge of amino acids strongly depends on their environment. Most interesting is to compare average charge on same-type amino acids. The classical expected charge of +1 agrees with our QM MD prediction for K48 (Q_ave=1.00) but it is not applicable to K11 (Q_ave=0.92). The most unusual average charge of -0.84 electron units is predicted for glutamic acid E51, whereas charge on same-type amino acid E16 is quite classical (-0.99; see Table 9.1). From this comparison it becomes apparent that the assignment of a unified unit charge on amino acids in anisotropic protein environment is not fully correct. Hence, this is a factor reducing accuracy of computer simulation of biomolecules. Classical force fields can in principle manipulate with fractional amino acid charges. However in order to derive such charges one would need to run QM MD simulation first. Such computational protocol may indeed be of certain utility.

Fig. 9.4. Extreme range of amino acid net charge fluctuation. Rectangles, circles, and triangles denote positively charged, neutral, and negatively charged amino acids, respectively.

Besides of the non-equivalent average charge determined for same-type amino acids by our QM MD calculation, we can also see global-scale dynamic charge fluctuations in ubiquitin as displayed in Fig. 9.4. According to these data every amino acid including neutral ones experiences large charge fluctuation of 0.25 electron units or more. Supporting the obtained numerical values is our previous observation that PM5 MD is even less prone to exaggerate charge transfer than ab initio HF/6-31G* level. In contrast to empirical polarizable models artificially confining charge fluctuation to individual amino acid fragments, in QM picture the charge can travel from and to any distant part of the system. Since this effect cannot be described by polarizable force fields this makes the current findings even more compelling. By averaging out charge fluctuation and resorting to fixed charges the classical model makes itself computationally more attractive but also less sensitive to perturbations. The charge density freely floating in the system, as in QM picture, is a sensitive probe to the environment changes. Subtle changes in the environment alter protein charge flow and the latter instantaneously transmits this information to the level of macromolecular motion. Therefore integrating global charge transfer effects into computational scheme is the necessary way to improve accuracy of biomolecular simulations.

9.3.5 Protein-solvent charge transfer

So far we limited the analysis of charge transfer effects solely to the processes happening inside the protein. However, since the original charge transfer effects pertaining to biological systems were observed on protein-water interactions[13] it is reasonable to expect large charge transfer effects taking place at the protein-water interface. Therefore it is necessary to extend the charge transfer analysis to solvent molecules. Among water molecules the largest charge fluctuation of 0.24 electron units happens in vicinity of protein N-terminus. Residue K48, which is fully exposed to solvent, induces maximal charge fluctuation of 0.21 electron units to the nearest water molecules. Similar charge fluctuation on water molecules happens in vicinity of the other solvent exposed ionized amino acids.

Now we can go back to the original question of significance of protein-to-water charge transfer. It was presumed that the charge transfer value divided by the total number of atoms will be entirely lost. However, as we could see above, large charge fluctuations happen at the level of individual solvent molecules. It helps identifying a flaw in the initial reasoning involving division by the total number of atoms, which is unacceptable. The charge transfer effect would be really negligible when dealing with isotropic systems like neat liquids where the bulk properties dominate over microscopic details. But proteins are highly anisotropic systems and averaging over the number of atoms is not physical here. Such averaging compromises the main purpose of biomolecular simulations in deducing the atomic-level details which requires the precise accounting of actual physics at every single atomic position. We could also see that water molecules are not all equivalent as presumed in empirical force fields. The water molecules close to protein participate in dynamic charge transfer and carry large excess charge. This makes them different from the bulk water molecules.

The charge transfer effects go somewhat further than just invalidating the charge equity between individual solvent molecules at microscopic level. A deviation from the classical picture can even be seen on a mesoscopic scale of large group of solvent molecules. The charge transfer between protein and solvent molecules induces excess charge on solvent molecules close to protein and the charge value depends on the distance between water molecules and protein surface as depicted on Fig. 9.5. To obtain these data the charge on water molecule is computed as the sum of partial atomic Mulliken charges on hydrogen and oxygen atoms. The distance is measured between oxygen atom of water molecule and the nearest protein atom, which may be hydrogen; next the computed charge is averaged over other water molecules located at the same distance and over the QM dynamics trajectory.

According to our data the charge on water molecules changes its sign as the function of distance from protein (Fig. 9.5). Due to water geometry and our choice of measurement of the protein-water distance the nearest water molecules are those involved in solvation of positively charged protein sites (protein-H^+…OHH) with the charge intensity peak happening at 2.2 Å. Second characteristic peak takes place at 2.6 Å. It presumably describes water molecules solvating electron-donor sites of the protein (protein-O^¯…HHO). After the distance of 3.5 Å the charge on individual water molecules becomes close to zero.

Fig. 9.5. Excess charge on individual water molecules as a function of distance to protein, averaged over PM5 MD trajectory.

Fig. 9.5 reveals an asymmetry in the ability of water molecules to hold excess charge of particular sign with the apparent preference to acquire excess negative charge. This indicates that greater solute-to-solvent charge transfer might be expected for highly negatively charged biomolecules, e.g. DNA and RNA. No such charge asymmetry is assumed in the design of empirical water models, i.e., the instantaneous charge fluctuations on water molecule averaged due to mean field approximation are expected to be zero. Moreover, since empirical water models do not take into account intermolecular charge transfer their radial charge distribution function will be a flat horizontal line at value zero rather than the shape shown in Fig. 9.5.

Having obtained a non-classical water radial charge distribution function, next obvious step is to integrate this function in order to obtain a picture of cumulative excess charge transferred from protein to water as a function of distance. Such integral plot is displayed in Fig. 9.6. It basically goes over the main conclusions derived from Fig. 9.5. The nearest water shell of thickness of 2.2 Å carries total net charge of +0.04 electron units. The net positive charge quickly becomes outweighed by the excess negative charge and the water shell of thickness of 3.3 Å already encloses 95% of the entire electron charge transferred from the protein to water. The second minimum at 5.7 Å encloses almost 100% of the excess charge with the value reaching plateau at 8 Å.

Fig. 9.6. Integral plot of cumulative charge on water molecules as a function of distance to protein, averaged over PM5 MD trajectory.

Somewhat more diffuse excess charge localization as a function of thickness of water shell was suggested by HF/6-31G* computations of ubiquitin-water system by Komeiji et al.[23] According to their studies the protein-to-water charge transfer value converges at about 6 Å water shell whereas in our case it converged at 3.3 Å. Although the superiority of the ab initio level over semiempirical one is unquestionable it is possible that the ab initio data were biased toward more diffused picture because of geometric strain in the studied structures. Unlike our PM5 MD simulation the HF/6-31G* computations were performed on the structures which do not correspond to the corresponding QM energy minimum. The energy strain present in the molecular structure can artificially increase the appearance of charge transfer. To test this assumption we performed single-point (rigid-structure) PM5 calculations over 1000 ubiquitin-water snapshots generated by force field MD and got a prediction of -1.6 electron units for protein-to-water charge transfer. It clearly showed that performing QM calculations over not self-consistent geometries leads to overestimated charge transfer predictions.

In cases when QM MD is not available the accuracy in handling charge transfer can be significantly improved by performing geometry relaxation at the corresponding QM level. Performing PM5 geometry optimization on 100 empirical snapshots changed the prediction from -1.6 to -0.6 electron units,[21] which is in excellent agreement with our QM MD data and the data of Komeiji. Therefore taking into account unavailability of ab initio MD data on charge transfer for ubiquitin-water system and based on the variety of indirect evidences discussed above there are reasons to believe that our PM5 MD predictions on charge transfer in ubiquitin-water system are the most accurate to date.

9.4 Implications of charge transfer

In the conclusion, the undertaken QM MD simulation revealed complex physics of dynamic charge transfer in the ubiquitin-water system. We observed a distinctly localized charge transfer channel pumping electron density between interacting amino acids in sync with their thermal motion. Delocalized effects of charge transfer are responsible for global instantaneous charge fluctuations exceeding 0.25 electron units and affecting every amino acid in the system including neutrally charged and weakly polarizable ones. Same-type amino acids which are treated as equal by force fields show different average net charges in QM picture due to their interaction with environment. It means that the assignment of the same unit charge to same-type amino acids is not actually correct. As it follows from our QM MD data water molecules of the solvent are also not equal. Those water molecules located in close proximity to the protein carry instantaneous excess charge of up to 0.2 electron units per water molecule resulting in accumulation of a charge cloud around the protein.

This atomic-level picture is largely overlooked by classical mechanics treatment of electrostatic interactions. Due to significant charge transfer effects taking place at every part of the protein-water system and because of the link between charge transfer and protein dynamics the biomolecular systems can no longer be considered as merely classical mechanical entities unless at the cost of omission of important atomistic details. Although empirical force fields dealing with fixed atomic (or amino acid) charges are computationally attractive, the oversimplified electrostatic model is the major factor limiting accuracy of computer simulation of biological systems. Therefore QM MD represents an important physics-based approach for further improvement of accuracy of biomolecular simulations.

Acknowledgments

The Authors are grateful to the Texas Advanced Computing Center (TACC) (http://www.tacc.utexas.edu) for providing high-performance computing resources for this project.

9.5 References

[1] Goedecker S (1999) Linear scaling electronic structure methods. Rev Mod Phys 71:1085-1123.

[2] Goedecker S, Scuseria GE (2003) Linear scaling electronic structure methods in chemistry and physics. Comp Sci Eng 5:14-21.

[3] McCammon JA, Gelin BR, Karplus M (1977) Dynamics of folded proteins. Nature 267:585-590.

[4] Warshel A, Levitt M (1976) Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J Mol Biol 103:227-249.

[5] Anikin NA, Anisimov VM, Bugaenko VL et al. (2004) LocalSCF method for semiempirical quantum-chemical calculation of ultralarge biomolecules. J Chem Phys 121:1266-1270.

[6] Anisimov VM, Bugaenko VL (2009) QM/QM docking method based on the variational finite localized molecular orbital approximation. J Comp Chem 30:784-798.

[7] Halgren TA, Damm W (2001) Polarizable force fields. Curr Opin Struct Biol 11:236-242.

[8] Sprik M, Klein ML (1988) A polarizable model for water using distributed charge sites. J Chem Phys 89:7556-7560.

[9] Caldwell J, Dang LX, Kollman PA (1990) Implementation of nonadditive intermolecular potentials by use of molecular dynamics: Development of a water-water potential and water-ion cluster interactions. J Amer Chem Soc 112:9144-9147.

[10] Rick SW, Stuart SJ, Berne BJ (1994) Dynamical fluctuating charge force fields: Application to liquid water. J Chem Phys 101:6141.

[11] Lopes PEM, Roux B, MacKerell AD (2009) Molecular modeling and dynamics studies with explicit inclusion of electronic polarizability: Theory and applications. Theor Chem Acc 124:11-28.

[12] Anisimov VM, Vorobyov IV, Roux B et al. (2007) Polarizable empirical force field for the primary and secondary alcohol series based on the classical Drude model. J Chem Theory  Comput 3:1927-1946.

[13] Nadig G, Van Zant LC, Dixon SL et al. (1998) Charge-transfer interactions in macromolecular systems: A new view of the protein/water interface. J Am Chem Soc 120:5593-5594.

[14] Thompson WH, Hynes JT (2000) Frequency shifts in the hydrogen-bonded oh stretch in halide-water clusters. The importance of charge transfer. J Am Chem Soc 122:6278-6286.

[15] Tanaka M, Aida M (2004) An ab initio MO study on orbital interaction and charge distribution in alkali metal aqueous solution: Li+, Na+, and K+. J Sol Chem 33:887-901.

[16] Ikeda T, Hirata M, Kimura T (2005) Hydration of Y3+ ion: A Car-Parrinello molecular dynamics study. J Chem Phys 122:024510.

[17] Cappa CD, Smith JD, Wilson KR et al. (2005) Effects of alkali metal halide salts on the hydrogen bond network of liquid water. J Phys Chem B 109:7046-7052.

[18] Voityuk A, A., Siriwong K, Rösch N (2004) Environmental fluctuations facilitate electron-hole transfer from guanine to adenine in DNA pi stacks. Angew Chem Int Ed 43:624-627.

[19] Kubar T, Kleinekathofer U, Elstner M (2009) Solvent fluctuations drive the hole transfer in DNA: A mixed quantum-classical study. J Phys Chem B 113:13107-13117.

[20] Peraro MD, Raugei S, Carloni P et al. (2005) Solute-solvent charge transfer in aqueous solution. ChemPhysChem 6:1715-1718.

[21] Anisimov VM, Bugaenko VL, Cavasotto CN (2009) Quantum mechanical dynamics of charge transfer in ubiquitin in aqueous solution. ChemPhysChem 10:3194-3196.

[22] Balabin IA, Beratan DN, Skourtis SS (2008) Persistence of structure over fluctuations in biological electron-transfer reactions. Phys Rev Lett 101:158102.

[23] Komeiji Y, Ishida T, Fedorov DG et al. (2007) Change in a protein's electronic structure induced by an explicit solvent: An ab initio fragment molecular orbital study of ubiquitin. J Comp Chem 28:1750-1762.

[24] Dewar MJS, Zoebisch EG, Healy EF et al. (1985) Development and use of quantum mechanical molecular models. 76. AM1: A new general purpose quantum mechanical molecular model. J Am Chem Soc 107:3902-3909.

[25] Rocha GB, Freire RO, Simas AM et al. (2006) Rm1: A reparameterization of AM1 for H, C, N, O, P, S, F, Cl, Br, and I. J Comp Chem 27:1101-1111.

[26] Stewart JJP (1989) Optimization of parameters for semiempirical methods i. Method. J Comp Chem 10:209-220.

[27] Stewart JJP (2002) Mopac 2002. Fujitsu Ltd, Tokyo, Japan.

[28] Hoover WG (1985) Canonical dynamics: Equilibrium phase-space distributions. Phys Rev A 31:1695-1697.

[29] Ferrario M, Fionino A, Ciccotti G (1997) Long-time tails in two-dimensional fluids by molecular dynamics. Physica A 240:268-276.

[30] van der Vaart A, Merz KM (2002) Charge transfer in small hydrogen bonded clusters. J Chem Phys 116:7380-7388.

[31] Cavalli A, Carloni P, Recanatini M (2006) Target-related applications of first principles quantum chemical methods in drug design. Chem Rev 106:3497-3519.

[32] Geissler PL, Dellago C, Chandler D et al. (2001) Autoionization in liquid water. Science 291:2121-2124.

[33] Patel S, Mackerell AD, Jr, Brooks CL, III (2004) CHARMM fluctuating charge force field for proteins: Ii protein/solvent properties from molecular dynamics simulations using a nonadditive electrostatic model. J Comp Chem 25:1504-1514.

[34] van Duin ACT, Dasgupta S, Lorant F et al. (2001) ReaxFF: A reactive force field for hydrocarbons. J Phys Chem A 105:9396-9409.

[35] Liu H, Elstner M, Kaxiras E et al. (2001) Quantum mechanics simulation of protein dynamics on long timescale. Proteins 44:484-489.

[36] Tong L, Warren TC, King J et al. (1996) Crystal structures of the human p56lcksh2 domain in complex with two short phosphotyrosyl peptides at 1.0 å and 1.8 å resolution. J Mol Biol 256:601-610.

[37] Vorobyov IV, Anisimov VM, MacKerell ADJ (2005) Polarizable empirical force field for alkanes based on the classical drude oscillator model. J Phys Chem B 109:18988-18999.

[38] Giese TJ, York DM (2004) Many-body force field models based solely on pairwise coulomb screening do not simultaneously reproduce correct gas-phase and condensed-phase polarizability limits. J Chem Phys 120:9903-9906.