Article Outline

Article Information

PubMed

Related Articles

  • …more

Copyright © 2008 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 95, Issue 3, 1080-1088, 1 August 2008

doi:10.1529/biophysj.108.131110

Biophysical Theory and Modeling

Developing Polarized Protein-Specific Charges for Protein Dynamics: MD Free Energy Calculation of pKa Shifts for Asp26/Asp20 in Thioredoxin

Changge Ji*, Ye Mei* and John Z.H. Zhang*, †, Go To Corresponding Author,  

* Institute of Theoretical and Computational Chemistry, Key Laboratory of Mesoscopic Chemistry of Ministry of Education, School of Chemistry and Chemical Engineering, Nanjing University, Nanjing, China
† Department of Chemistry, New York University, New York, New York

Address reprint requests to John Z. H. Zhang.

Abstract

Ab initio quantum mechanical calculation of protein in solution is carried out to generate polarized protein-specific charge(s) (PPC) for molecular dynamics (MD) stimulation of protein. The quantum calculation of protein is made possible by developing a fragment-based quantum chemistry approach in combination with the implicit continuum solvent model. The computed electron density of protein is utilized to derive PPCs that represent the polarized electrostatic state of protein near the native structure. These PPCs are atom-centered like those in the standard force fields and are thus computationally attractive for molecular dynamics simulation of protein. Extensive MD simulations have been carried out to investigate the effect of electronic polarization on the structure and dynamics of thioredoxin. Our study shows that the dynamics of thioredoxin is stabilized by electronic polarization through explicit comparison between MD results using PPC and AMBER charges. In particular, MD free-energy calculation using PPCs accurately reproduced the experimental value of pKa shift for ionizable residue Asp26 buried inside thioredoxin, whereas previous calculations using standard force fields overestimated pKa shift by twice as much. Accurate prediction of pKa shifts by rigorous MD free energy simulation for ionizable residues buried inside protein has been a significant challenge in computational biology for decades. This study presented strong evidence that electronic polarization of protein plays an important role in protein dynamics.

Introduction

Molecular dynamics (MD) simulations of biological systems are now routine exercises in computational biology. Such computational advances have had a huge impact on our ability to understand protein structure and dynamics, protein-ligand binding, protein-protein interaction, drug design, etc. These large-scale MD simulations are made possible by the availability of force fields (FFs), which provide simple and fast evaluation of molecular interactions in biological systems. A great deal of progress has been made over the past decades in the development of molecular mechanics FFs for use in biological simulations [1,2,3,4,5,6](1–6). When parameters are extensively fitted to experimental data, accurate models can be produced. However, for the overwhelming majority of important biological applications, only a limited amount of relevant experimental data is available. Applications of the current molecular mechanics force fields have given decent pictures of the details of atomic motions and energetic interactions, but the level of quantitative accuracy is uncertain and in many cases, not adequate for reliable predictive studies.

An important component of the interaction in proteins is the electrostatic interaction, which plays a significant role in protein structure and function [7,8,9,10,11](7–11). Processes such as protein folding, protein-ligand binding, protein-protein interaction, electron transfer, proton binding and release, enzyme reaction, etc. are largely driven by electrostatic interactions. In particular, proton binding is a pH-dependent process and is strongly influenced by electrostatic interactions with the protein and the local environment [12,13,14](12–14). It is well known that the electrostatic interaction between protein side chains depends not only on their distances but also on their locations in the protein and their local solvent environment. In proton-binding, the protein and solvent undergo dielectric relaxation involving electronic polarization and displacement of atomic groups. As is understood, the local electrostatic environment inside a protein is inhomogeneous and hydrophobic, which is generally less favorable for ionization [15](15). Thus, ionizable (or charged) residues in the interior of proteins can have substantial pKa shifts compared to those on the surfaces of proteins or isolated amino acids in solution [16](16).

Current standard force fields, e.g., CHARMM and AMBER, are amino-acid-specific and therefore very portable. Despite great success in application of the standard FFs, there are fundamental limitations in their applications. Specifically, these FFs are amino-acid-specific or mean-field-like and therefore fail to give accurate representation of the electrostatics of the specific protein environment which is highly inhomogeneous and protein-specific. For example, two amino acids of the same type in the same or different proteins should have quite different charge status due to their different electrostatic environments. Thus, using the same set of charges for them is obviously inadequate. In summary, the current amino-acid-based FFs are incapable of describing the polarization state of a particular protein structure, e.g., native protein structure. To overcome this fundamental deficiency of the standard FFs, efforts have been made to develop polarizable force-field models for peptides and proteins [17,18](17,18). The polarizable FFs are theoretically attractive but their practical application is much more complicated than for the standard FFs. Consequently, simulations of biomolecular systems with polarizable FFs are still uncommon [19](19) and uncertainties regarding the accuracy and validity of the underlying theoretical models used to derive polarizable FFs remain [20](20).

In this article, we present an alternative force field containing polarized protein-specific charges for MD study of proteins. The polarized protein-specific charge (PPC) builds protein polarization into the atomic charges and is therefore computationally attractive for standard MD simulation without any dynamical complications. The PPC is not “polarizable” per se, but it is derived from first-principle quantum solvation calculation of protein in the native (or a given) structure. Thus, PPCs correctly represent the electronically polarized state of the protein and therefore provide accurate electrostatic interaction near the native structure. In our approach, the quantum chemistry calculation of protein solvation is made possible by combining a recently developed fragment-based scheme, molecular fractionation with conjugate caps (MFCC), with a continuum-solvent model. In this article, we employ linearized Poisson-Boltzmann method [21,22](21,22) to solve the self-consistent reaction-field equation coupled with quantum chemistry calculation of the solute using the MFCC scheme [23,24,25,26,27](23–27). The converged electron density of protein fragment (amino acid) is fitted to generate partial charges for every amino acid in the protein using the RESP method [28,29](28,29). The resulting fitted atomic partial charges are protein-specific and they correctly represent the polarized electronic state of the protein in the native (or other given) structure. Because the PPC is atom-centered and maintains the same simplicity as the standard charges in AMBER, it can be easily applied in MD simulation without any additional complication. Thus, we expect PPCs to provide much improved electrostatic interactions in MD simulation of protein near its native structure, both in structure and dynamics. For application of the PPC in MD simulation, one simply replaces the standard charges from the AMBER force field by PPC while keeping the rest of the force parameters intact.

To study the possible effect of PPCs in MD simulation of proteins, we carried out a series of MD studies to investigate various structural and dynamical properties of thioredoxin. The new MD results using PPCs are explicitly compared with those obtained from the corresponding MD studies but using the AMBER charges. In particular, MD free-energy simulation is carried out to predict the pKa shift of thioredoxin for buried Asp26. Previous MD free energy calculations using both AMBER and CHARMM give pKa shifts that are twice as large as the experimental value [30](30). To allow direct comparison with those results, we follow faithfully the computational methods and procedures in Simonson et al. [30](30) to calculate pKa shifts except that the standard AMBER charges are replaced by our calculated PPC for thioredoxin. This shall eliminate any possible dynamical uncertainties in the comparison of results.


Theoretical approach

In the continuum-solvent model, the solute (protein) is represented by a charge distribution ρ(r) embedded in a cavity surrounded by a polarizable medium with dielectric constant ϵ. The solute charge distribution ρ(r) polarizes the dielectric medium and creates a reaction field which acts back to polarize the solute until equilibrium is reached. The reaction field acting on the solute can be effectively represented by that of induced charges on the cavity surface according to the classical electrostatic theory. By discretizing the induced charges on the surface of the cavity and iteratively solving the quantum chemistry equation for the solute in an external reaction field created by the surface charges, one obtains the popular PCM method [31](31) as recently generalized to protein solvation [27,32](27,32). However, for large proteins, the PCM model requires many discrete surface charges and thus makes the solution of linear equation difficult computationally. In the current approach, we numerically solve the Poisson-Boltzmann (PB) equation to obtain the reaction field as was done before by other researchers [22,34](22,34). This avoids the solution of large linear equations and is computationally more attractive for large proteins.

The basic procedures in fitting atomic charges of protein in our approach can be described as follows: First, gas phase calculation of protein is performed with the MFCC approach to obtain initial electron density of the protein for the given structure as described in earlier publications [24](24). The calculated electron density is used to fit atomic charges using the RESP program [35](35). The charge-fitting philosophy used here was the same as that used in the AMBER force field and this makes the guarantee that the PPC was consistent with other parameters of the AMBER force field. Solution of the PB equation is then carried out to obtain the reaction field from which to generate discrete surface charges on the cavity surface.

The coupling of MFCC/RESP charge-fitting method for protein with a Poison-Boltzmann continuum solvent should provide accurate estimates of the solvation free energies of proteins in water. Partial charges on each atom of the protein generated from the MFCC/RESP procedure were passed to the PB solver DELPHI [36](36) to determine the self-consistent reaction field. A set of induced surface charges qind on the dielectric boundary was derived that represent reaction field effects of solvent molecules. The dielectric solute/solvent boundary was defined by AMBER van der Waals radii [37](37) for atoms of solute molecule with a probe radius of 1.4Å. The internal dielectric constant, denoted as ϵsolute, was set to unity, for molecular polarizability is explicitly included in quantum mechanical calculation. The solvent dielectric constant, denoted as ϵsolvent, was set to 80. Grid density was set to 4.0grids/Å. Surface charges were then added as background charges in the next QM calculation for each capped fragment (CF). As mentioned above, partial charges of the protein and the screened surface charges polarize each other until converge was reached. The cycle stops only when the dipole of the protein and the surface charges both converge to within a certain numerical accuracy.

The electrostatic solvation free energy can be divided into two terms: solute polarization energy (Gpol) and polarized reaction field energy (Ges),

(1)
In the DELPHI program, the latter term can be expressed as
(2)
where is the ith induced charge on the dielectric boundary, and ϕi is the electrostatic potential on position of ith induced charge. In the MFCC approximation, the solute polarization energy is given by
(3)
where the polarization energies of individual fragment ΔEk and conjugate cap are given, respectively, by Mei et al. [27](27), as
(4)
(5)
where and ρk stand for electron density distribution of each fragment or conjugate cap in vacuum and in solvent, respectively.

The procedure to generate PPCs can be understood easily from the flow chart in Fig. 1. The quantum chemistry calculation of electron density of individual protein fragment is performed at the level of B3LYP/6-31G* for results reported in this article.

Display large version of this figure
Figure 1
The flow chart of the MFCC-PB computational protocol.

Results and discussion

Solvation and polarization of protein

The MFCC-PB computational protocol described in the previous section needs to be extensively tested numerically to validate its efficacy. In particular, we need to examine the effect of protein polarization as represented by PPC on protein solvation. We performed benchmark studies to examine the contribution of electronic polarization on solvation energy for a number of protein systems. The calculated contribution of electrostatic solvation energy for a number of protein systems is listed in.Table 1 For comparison, results calculated from the previous calculation using the MFCC-CPCM method were also listed in the table. Different from classical PB or generalized Born (GB) methods, solute polarization energy is quite significant in the calculated electrostatic solvation free energy. Our calculation shows that this term contributes 7∼19% to the total electrostatic solvation energy as shown in Table 1. This big part of polarization energy indicates that inclusion of the polarization effect in the current force field is needed.

Table 1 Comparison of electrostatic solvation free energies from MFCC-FDPB calculation and MFCC-CPCM (kcal/mol)
MFCC-FDPBMFCC-CPCM*
PeptidePDB IDGesGpolGele†Gpol/GeleGesGpolGele†
Amyloid1AMC−1013.13163.42−849.7019.23%−991.16105.15−886.01
BPTI1BPI−1405.71101.85−1303.857.81%−1420.2387.52−1332.71
Calbindin1CDN−2533.75352.35−2181.4116.15%−2402.94143.32−2259.62
Crambin1CBN−392.7247.59−345.1313.79%−406.1044.58−361.52
Lysozyme2BLX−2112.25224.06−1888.1911.87%−2034.35146.64−1887.71
RP719551RPB−289.2532.94−256.3012.85%−303.7436.14−267.60
Ubiquitin1UBQ−1129.00138.57−990.4413.99%−1091.4594.43−997.02
* See Mei et al. 27(27).
† Gele=Ges+Gpol.

To investigate the effect of protein polarization on various protein properties, we also compare calculated dipole moments of protein in solution calculated from the current PPC with those from AMBER03 force field. Since most proteins have a nonzero net charge, and the calculated dipole moment depends on the placement of a charged entity within the coordinate system, so the protein is first centered to its charge center. The charge center is defined as in which ri is the original coordinate of the ith atom in the protein, and qi is its partial charge. The dipole of the entire protein is then given by To examine more detailed local electrostatic environment, we calculated dipole moments for individual residues of proteins. Fig. 2 shows plots of dipole moment of individual residues calculated from PPC versus that from AMBER03 for a number of protein systems. The comparison shows that dipoles of most residues calculated from PPCs are generally larger than those from the AMBER charges. This may indicate the fact that more polarization of the residues is accounted for in the MFCC-PB charges.

Display large version of this figure
Figure 2
Comparison of dipole moments of individual residues calculated from PPC and those from AMBER03 charge for several protein systems.

Protein-specific atomic charges

Since the MFCC-PB charges correctly describe the electrostatic polarization of the protein, partial charges on the same type of amino acid but in different locations of the protein are generally different as expected. Partial charges on a particular residue are determined by its specific conformation and chemical environment due to other residues of the protein. To illustrate this feature, we listed some partial charges of lysozyme in Table 2. As in lysozyme, we find that charges on the same type of amino acid differ a lot from each other. For examples, the value of atomic charge ranges from −0.083 to −0.642 for backbone atom N in ASP, from −0.727 to −0.556 for backbone atom O in Ser, and from 0.025 to 0.304 for atom CG1 in Ile. This is the most important feature of the PPC, which is different from the standard amino-acid-based force fields with the same fixed atomic charges for a given residue.

Table 2 Some partial charges on lysozyme calculated by MFCC-FDPB
Atom*ResidueChargeAtom*ResidueChargeAtom*ResidueCharge
NAsp18−0.414OSer24−0.603CG1Ile550.024
NAsp48−0.082OSer50−0.727CG1Ile580.165
NAsp52−0.641OSer72−0.588CG1Ile780.304
NAsp66−0.433OSer85−0.582CG1Ile880.024
NAsp87−0.309OSer86−0.604CG1Ile980.137
NAsp101−0.510OSer91−0.657CG1Ile1240.063
NAsp119−0.484OSer100−0.604
−0.558† −0.580† 0.022†
* Atom name in PDB nomenclature.
† AMBER03 charge.

Calculations of pKa shifts for Asp26/Asp20 in thioredoxin

Although standard amino-acid-based partial charge model has worked well in modeling many of the macroscopic properties of proteins in MD simulation, it is expected to have difficulties when simulating properties that are more sensitive to local electrostatic environment. This is because the standard force-field charges are mean-field-like and they do not contain protein polarization and other protein-specific electrostatic information described above. For example, accurate prediction of pKa shifts for buried residues in proteins is a challenging problem in current MD simulation [30](30). Accurate calculation of pKa shifts is crucial for us to understand the mechanism of acid-base catalysis in enzymes. Activity of such enzymes requires that the catalytic residues exist in the appropriate protonation state. A popular approach used in pKa calculation is based on solving the PB equation in which the solvent is treated as a continuum dielectric medium [38,39,40](38–40). The PB approach is a mean-field theory and can give useful insight for some simple cases. However, the PB approach fails in more complicated situations because it is not a microscopic model and therefore could not account for detailed molecular processes. To accurately predict the pKa shift, one needs to correctly account for molecular factors affecting the proton-binding process. Thus, microscopic methods in which water molecules are explicitly included are needed to correctly predict pKa shift from first principles. However, some of the previous attempts using microscopic methods failed to give accurate prediction of pKa shift [30,41,42,43](30,41–43).

A recent molecular dynamics free energy (MDFE) study by Simonson et al. to calculate pKa shift for Asp26 in Thioredoxin demonstrated the problems with the current standard AMBER and CHARMM force fields. The results from the rigorous MDFE simulations with both force fields and from different runs all overestimate the pKa shift of Asp26 by∼4-5kcal/mol [30](30). The MDFE calculation of pKa by Simonson et al. is perhaps the most rigorous and extensive study so far and the results can serve as a benchmark for the system being studied. Simonson et al. concluded that the problem must arise in part from systematic errors in the force fields employed [30](30). In this article, the same rigorous MDFE procedures of Simonson et al. [30](30) are followed to calculate the pKa shifts for buried Asp26 and surface-exposed Asp20 but using the new PPC described previously in this article.

In MDFE calculation, the thermodynamic integration approach is employed in which the interaction potential is defined by

(6)
where U0 and U1 are energy functions corresponding, respectively, to the protonated and deprotonated states. The parameter λ varies from 0 to 1 representing the system going from the initial (protonated) to the final (deprotonated) state. During the perturbation process, only charges on the residue undergoing titration are changed. The free-energy derivative with respect to parameter λ is given by
(7)
Simulations for a series of values of λ are carried out and the free-energy change is constructed from integration over λ by Gaussian quadrature,
(8)
The free energy shift is given by
(9)
and the corresponding shift of pKa value is
(10)
The numerical details of our MDFE simulation closely follow that of Simonson et al. [30](30). The NMR structure of thioredoxin (Protein Data Bank, i.e., PDB, code 1XOA; and see Fig. 3) is used as the initial structure. The simulations were performed with the AMBER program and explicit water model (TIP3P) with periodic boundary conditions is used. Firstly, the system was equilibrated under NPT ensembles with 500-ps runs at 298 K. The simulation is then continued in the NVE ensemble. SHAKE [44](44) was used to constrain covalent bonds connected to hydrogens. Long-range electrostatic interactions were treated by the particle-mesh Ewald method. The time step is one femtosecond. The model compound we used is aspartic acid with n-acetyl and n-methylamide blocking groups. Charge distributions in the protonated and ionized states are shown in Fig. 4.

Display large version of this figure
Figure 3
The NMR structure of thioredoxin (PDB code 1XOA) showing buried Asp26 and surface-exposed Asp20.
Display large version of this figure
Figure 4
Partial charges of the model compound (2n-acetyl 1-1 n-methyl-aspartic acid 1-amide) representing change between protonated (upper values) and deprotonated (lower values) states.

The calculation results for buried Asp26 as well as the surface-exposed Asp20 and the corresponding model system are given in Table 3. Our calculated free-energy shift for deprotonation of Asp26 is 5.0kcal/mol, which is very close to the experimental value of 4.8kcal/mol. The simulation went all the way to 18 ns to provide adequate samplings for all the intermediate states of the protein. Fig. 5 shows the energy derivative of the potential with respect to the system parameter λ=0.5 as a function of simulation time. This corresponds to a fictitious intermediate state between the protonated and d-protonated states of the Asp26. For purpose of verification of our calculation, we also calculated the pKa shift of Asp20, a surface-exposed residue. Our simulation result gives a free-energy difference of 1.0kcal/mol, which is in quite reasonable agreement with the experimental value of zero considering the statistical uncertainties of the simulation. Our result shows that using PPC, which correctly describes the polarized electrostatic environment surrounding the charged residue in the protein, can give accurate prediction of pKa shift from rigorous MDFE simulations from first principles. In contrast, the same numerical simulation but using AMBER and CHARMM charges produced pKa shifts that are twice as large as the experimental value in Simonson et al. [30](30).

Table 3 Results from molecular dynamics free energy simulation for shifts of pKa values of Asp26 and Asp20 in thioredoxin
ModelAsp26Asp20
Run length (ns)61818
∂G/∂λ(λ=0.11270)10.13.512.3
∂G/∂λ(λ=0.5)−43.6−35.5−42.5
∂G/∂λ(λ=0.88279)−101.8−83.2−102.3
ΔG−44.9−39.8−43.9
ΔΔG 5.01.0
ΔΔG* 9.11.7
ΔΔG(expt) 4.80.0
All energies are in kcal/mol.
* Result from Simonson et al. 30(30) using the AMBER force field.

To estimate the statistical error, we performed block averaging analysis in which each trajectory of 6-ns duration is divided into four blocks. The statistical uncertainty of the derivatives is estimated as twice the standard deviation of block averages. As we can see from Table 3, the statistical error of the free energy is ∼1kcal/mol and the calculated theoretical value is within the statistical error bar from the experimental data for Asp26 and Asp20. Since our MDFE calculation is exactly the same as those in Simonson et al. [30](30) except that we replaced AMBER atomic charges by PPC while keeping all the other parameters in the force field, the different result must come from the effect of charges. In other words, the polarization effect of the protein contained in the PPC is responsible for the accurate prediction of pKa shift in Asp26. This result really demonstrated the important role that electronic polarization plays in delicate electrostatic interaction of protein.

The above result implies that the lack of electronic polarization may be a major culprit in some earlier MD free-energy calculations of pKa shifts for ionizable residues buried in the interiors of proteins. Of course, more such calculations for different protein systems are needed to confirm this conclusion. There are of course other important dynamical effects that have been proposed to explain the failure of MD simulations of pKa shifts. For example, Kato and Warshel have proposed that a major factor that determines the pKa values of buried groups are the structural changes that take place when an internal group is ionized [42,45](42,45). However, the possibility of large structural reorganization is not easy to study computationally because of problems associated with limited MD samplings. For example, the calculation of pKa shift for Glu66 in staphylococcal nuclease has presented a huge challenge so far. Kato and Warshel developed a special overcharging approach to allow for more structural reorganization and obtained encouraging pKa shift for Glu66 in staphylococcal nuclease [45](45). In reality, we expect that both electronic polarization and structural reorganization should play some roles in determining pKa shifts for buried residues and they are likely intertwined with each other. We believe that the actual level of contribution from a particular effect will be case-dependent. Clearly, further computational study will be needed to clarify the situation.

Because the PPC represents the polarized electrostatic state of the protein in a given (native) structure, we expect there will be some difference in structural details of the protein from MD simulations using PPC in comparison to mean-field-like charges as in AMBER. We performed further MD calculations of thioredoxin using both PPC and AMBER charges to compare structural details. First, we examine the characteristics of RMSD from two sets of MD simulations using, respectively, PPC and AMBER03. They obviously exhibit different characteristics as shown in Fig. 6. Since all the other parameters in the force field used in both calculations are the same except the charges, the result in Fig. 6 clearly shows the important effect of protein polarization on detailed protein dynamics, which is electrostatic in nature. Since the formation of internal hydrogen bonds in protein polarize the donors and acceptors, this polarization effect is better represented by PPCs than by AMBER charges. As a result, the hydrogen bonds in protein should generally be more stable during MD calculations using PPC and this may explain the different characteristics of RMSD in MD simulation shown in Fig. 6 for thioredoxin.

Display large version of this figure
Figure 6
RMSD of backbone atoms of thioredoxin as a function of MD simulation time using AMBER03 (upper red curve) and the PPC (lower black curve).


Summary

We presented a new computational protocol for quantum mechanical calculation of protein in solution by combining a linear scaling fragmentation scheme for electronic structure of the protein and continuum dielectric model for the solvent in a self-consistent treatment. The computed electron density of protein in the native structure is employed to produce the polarized protein force field (charges). The derived protein-specific PPC correctly represents the polarized electrostatic state of protein near its native structure and can be easily employed in MD simulation of protein motion. Our computational study from MD simulation for thioredoxin revealed the following properties of PPC.

1. The PPC is protein-specific and thus correctly describes the polarized electrostatic state of the protein near the vicinity of native structure. The charges on the same type of residue are generally different depending on their specific local electrostatic environment. This is in contrast to the amino-acid-based charges that are mean-field-like and remain the same irrespective of the specific position in the protein.
2. By employing PPCs, MD free-energy simulation using thermodynamic integration approach accurately reproduces the experimental pKa shift for buried Asp26 in thioredoxin while the same calculations employing standard mean-field force fields of CHARMM and AMBER overestimate the pKa shift by approximately twice as much.

It is also important to point out the limitation of the PPC. Since PPC is based on a given structure (native structure in most cases) of a protein, it may not be appropriate or advantageous to use it to describe structures far off from the given structure. Fortunately, most MD studies of proteins are meant to simulate protein motions near the vicinity of native structures and PPC should be very attractive in such computational studies of proteins.


Acknowledgments

This work is partially supported by the National Science Foundation of China (grant No. 20773060) and National Basic Research Program of China (grant No. 2004CB719901). J.Z.H.Z. acknowledges partial financial support from the New York University Research Challenge Fund. C.J. thanks the Scientific Research Foundation of Graduate School of Nanjing University (grant No. 2006CL06).

References

1. Lifson, S., and Warshel, A. (1968). Consistent force field calculations of conformations vibrational spectra and enthalpies of cycloalkane and n-alkane molecules. J. Chem. Phys. 49, 5116–5129. CrossRef | PubMed

2. Brooks, B.R., Bruccoleri, R.E., Olafson, B.D., States, D.J., Swaminathan, S., and Karplus, M. (1983). CHARMM: a program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 4, 187–217. CrossRef | PubMed

3. Pearlman, D.A., Case, D.A., Caldwell, J.W., Ross, W.S., Cheatham, T.E., DelBolt, S., Ferguson, D., Seibel, G., and Kollman, P. (1995). AMBER, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules. Comput. Phys. Commun. 91, 1–41. PubMed

4. Hagler, A.T., Huler, E., and Lifson, S. (1977). Energy functions for peptides and proteins. I. Derivation of a consistent force field including the hydrogen bond from amide crystals. J. Am. Chem. Soc. 96, 5319–5327. CrossRef | PubMed

5. Nemethy, G., Pottle, M.S., and Scheraga, H.A. (1984). Energy parameters in polypeptides. 9. Updating of geometrical parameters, nonbonded interactions, and hydrogen bond interactions for the naturally occurring amino acids. J. Phys. Chem. 87, 1883–1887. CrossRef | PubMed

6. Berendsen, H.J.C., Postma, J.P.M, van Gunsteren, W.F., diNola, A., and Haak, J.R. (1984). Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81, 3684–3690. CrossRef | PubMed

7. Perutz, M. (1978). Electrostatic effects in proteins. Science 201, 1187–1191. PubMed

8. Warshel, A., and Russell, S. (1984). Calculations of electrostatic effects in biological systems and in solutions. Q. Rev. Biophys. 17, 283–342. PubMed

9. Matthew, J.B. (1985). Electrostatic effects in proteins. Annu. Rev. Biophys. Biophys. Chem. 14, 387–417. PubMed

10. Honig, B., and Nicholls, A. (1995). Classical electrostatics in biology and chemistry. Science 268, 1144–1149. PubMed

11. Davis, M.E., and McCammon, J.A. (1990). Electrostatics in biomolecular structure and dynamics. Chem. Rev. 90, 509–521. CrossRef | PubMed

12. Warshel, A. (1981). Calculations of enzymatic reactions—calculations of pKa, proton-transfer reactions, and general acid catalysis reactions in enzymes. Biochemistry 20, 3167–3177. PubMed

13. Schaefer, M., Vlijmen, H.V., and Karplus, M. (1998). Electrostatic contributions to molecular free energies in solution. Adv. Protein Chem. 51, 1–57. CrossRef | PubMed

14. Simonson, T. (2003). Electrostatics and dynamics of proteins. Rep. Prog. Phys. 66, 737–787. PubMed

15. Warshel, A. (1978). Energetics of enzyme catalysis. Proc. Natl. Acad. Sci. USA 75, 5250–5254. CrossRef | PubMed

16. Sternberg, M., Hayes, F., Russell, A., Thomas, P., and Fersht, A. (1987). Prediction of electrostatic effects of engineering of protein charges. Nature 330, 86–88. CrossRef | PubMed

17. Halgren, T.A. (2001). Polarizable force fields. Curr. Opin. Struct. Biol. 11, 236–242. CrossRef | PubMed

18. Kaminski, G.A., Stern, H.A., Berne, B.J., Friesner, R.A., Cao, Y.X.X., Murphy, R.B., Zhou, R.H., and Halgren, T.A. (2002). Development of a polarizable force field for proteins via ab initio quantum chemistry: first-generation model and gas phase tests. J. Comput. Chem. 23, 1515–1531. CrossRef | PubMed

19. Jorgensen, W.L. (2007). Special issue on polarization. J. Chem. Theory Comput. 3, 1877. PubMed

20. Warshel, A., Kato, M., and Pisliakov, A.V. (2007). Polarizable force fields: history, test cases, and prospects. J. Chem. Comput. 3, 2034–2045. PubMed

21. Honig, B., Sharp, K.A., and Yang, A.-S. (1993). Macroscopic models of aqueous solutions: biological and chemical applications. J. Phys. Chem. 97, 1101–1109. CrossRef | PubMed

22. Tannor, D.J., Marten, B., Murphy, R., Friesner, R.A., Sitkoff, D., Nicholls, D.A., Ringnalda, M., Goddard, W.A., and Honig, B. (1994). Accurate first principles calculation of molecular charge distributions and solvation energies from ab initio quantum mechanics and continuum dielectric theory. J. Am. Chem. Soc. 116, 11875–11882. CrossRef | PubMed

23. Zhang, D.W., and Zhang, J.Z.H. (2003). Molecular fractionation with conjugate caps for full quantum mechanical calculation of protein-molecule interaction energy. J. Chem. Phys. 119, 3599–3605. CrossRef | PubMed

24. Gao, A.M., Zhang, D.W., Zhang, J.Z.H., and Zhang, Y.K. (2004). An efficient linear scaling method for ab initio calculation of electron density of proteins. Chem. Phys. Lett. 394, 293–297. PubMed

25. Mei, Y., Zhang, D.W., and Zhang, J.Z.H. (2005). A new method for direct linear scaling calculation of electron density of protein. J. Phys. Chem. A 109, 2;, Letter.. PubMed

26. Mei, Y., Wu, E.L., Han, K.L., and Zhang, J.Z.H. (2006). Treating hydrogen bonding in ab initio calculation of biopolymers. Int. J. Quantum Chem. 106, 1267–1276. PubMed

27. Mei, Y., Ji, C.G., and Zhang, J.Z.H. (2006). A new quantum method for electrostatic solvation energy of protein. J. Chem. Phys. 125, 94906–94912. PubMed

28. Bayly, C.I., Cieplak, P., Cornell, W., and Kollman, P.A. (1993). A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. J. Phys. Chem. 97, 10269–10280. CrossRef | PubMed

29. Cornell, W.D., Cieplak, P., Bayly, C.I., and Kollman, P.A. (1993). Application of RESP charges to calculate conformational energies, hydrogen bond energies, and free energies of solvation. J. Am. Chem. Soc. 115, 9620–9631. CrossRef | PubMed

30. Simonson, T., Carlsson, J., and Case, D.A. (2004). Proton binding to proteins: pKa calculations with explicit and implicit solvent models. J. Am. Chem. Soc. 126, 4167–4180. CrossRef | PubMed

31. Miertuš, S., Scrocco, E., and Tomasi, J. (1981). Electrostatic interaction of a solute with a continuum. A direct utilization of ab initio molecular potentials for the prevision of solvent effects. Chem. Phys. 55, 117–129. PubMed

32. Fedorov, D.G., Kitaura, K., Li, H., Jensen, J.H., and Gordon, M.S. (2006). The polarizable continuum model (PCM) interfaced with the fragment molecular orbital method (FMO). J. Comput. Chem. 27, 976–985. CrossRef | PubMed

33. Reference deleted in proof.

34. Gogonea, V., and Merz, K.M. (1999). Fully quantum mechanical description of proteins in solution. Combining linear scaling quantum mechanical methodologies with the Poisson-Boltzmann equation. J. Phys. Chem. A 103, 5171–5188. PubMed

35. Case, D.A., Darden, T.A., Cheatham, T.E., Simmerling, C.L., Wang, J., Duke, R.E., Luo, R., Merz, K.M., Wang, B., Pearlman, D.A., et al. (2004). AMBER 8. (San Francisco: University of California). PubMed

36. Rocchia, W., Sridharan, S., Nicholls, A., Alexov, E., Chiabrera, A., and Honig, B. (2002). Rapid grid-based construction of the molecular surface and the use of induced surface charge to calculate reaction field energies: applications to the molecular systems and geometric objects. J. Comput. Chem. 23, 128–137. CrossRef | PubMed

37. Cornell, W.D., Cieplak, P., Bayly, C.I., Gould, I.R., Merz, K.M., Ferguson, D.M., Spellmeyer, D.C., Fox, T., Caldwell, J.W., and Kollman, P.A. (1995). A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J. Am. Chem. Soc. 117, 5179–5197. CrossRef | PubMed

38. Bashford, D., and Karplus, M. (1990). pKas of ionizable groups in proteins: atomic detail from a continuum electrostatic model. Biochemistry 29, 10219–10225. PubMed

39. Antosiewicz, J., McCammon, J., and Gilson, M. (1994). Prediction of pH-dependent properties of proteins. J. Mol. Biol. 238, 415–436. CrossRef | PubMed

40. Demchuk, E., and Wade, R. (1996). Improving the continuum dielectric approach to calculating pKas of ionizable groups in proteins. J. Phys. Chem. 100, 17373–17387. CrossRef | PubMed

41. Warshel, A., Sussman, F., and King, G. (1986). Free energy of charges in solvated proteins: microscopic calculations using a reversible charging process. Biochemistry 25, 8368–8372. PubMed

42. Warshel, A., Sharma, P.K., Kato, M., and Parson, W.W. (2006). Modeling electrostatic effects in proteins. Biochim. Biophys. Acta 1764, 1647–1676. PubMed

43. Riccardi, D., Schaefer, P., and Cui, Q. (2005). pKa calculations in solution and proteins with QM/MM free energy perturbation simulations: a quantitative test of QM/MM protocols. J. Phys. Chem. B 109, 17715–17733. PubMed

44. Ryckaert, J., Ciccotti, G., and Berendsen, H. (1977). Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 23, 327–341. PubMed

45. Kato, M., and Warshel, A. (2006). Using a charging coordinate in studies of ionization induced partial unfolding. J. Phys. Chem. B 110, 11566–11570. PubMed

Publication Information


Received: February 7, 2008
Accepted: March 18, 2008