SCIENCE CHINA Chemistry • ARTICLES •

December 2011 Vol.54 No.12: 1974–1981 doi: 10.1007/s11426-011-4399-3

Folding of EK peptide and its dependence on salt concentration and pH: A computational study MEI Ye1,2*, ZHANG DaWei3*, DUAN LiLi4, ZHANG QingGang4 & ZHANG John ZengHui1,2,5* 1

2

State Key Laboratory of Precision Spectroscopy; East China Normal University, Shanghai 200062, China Institute of Theoretical and Computational Science, East China Normal University, Shanghai 200062, China 3 Division of Chemistry and Biological Chemistry, School of Physical and Mathematical Sciences, Nanyang Technological University, Singapore 637371, Singapore 4 College of Physics and Electronics, Shandong Normal University, Jinan 250014, China 5 Department of Chemistry, New York University, New York 10003, USA Received July 5, 2011; accepted August 12, 2011; published online November 4, 2011

In this study, we apply both pairwise AMBER03 force field and the recently developed polarized force field to study the folding process of EK peptide under various ion strength and pH conditions. The polarized force field is based on our newly proposed adaptive hydrogen bond-specific charge (AHBC) scheme. These two force fields differ only by the atomic charges. Solvent effect is described with generalized Born models (IGB5 in AMBER 10 package). The result shows that although when applying AMBER03 charge, the helical structure is preferred, its dependence on salt concentration and pH is qualitatively wrong. While using AHBC the peptide finds its native structure within 10 ns, and then fluctuates around this folded state. Under high salt concentration or extreme pH conditions the calculated helical structure probability drops, which is in qualitative agreement with the experiment. Analysis of the atomic charges and the interaction between the donor-acceptor pair in main hydrogen bonds shows that the helical structure is stabilized when polarization effect is counted. It again shows that polarization effect is a very important improvement over traditional force field and is essential for protein folding. We also prove that the salt bridge interaction between 4-residue apart GLU and LYS residues is not critical to the stability of helical structure of EK peptide, but is merely an auxiliary factor, also in agreement with the experiment. EK peptide, polarization, salt concentration, pH

1 Introduction How protein folds to its precisely defined structure has been one of the most challenging problems for contemporary research. Understanding the mechanism of folding is thus important from both experimental and theoretical perspectives. Much effort has been devoted to study the folding process, but an atomic and femtosecond (fs) spatiotemporal resolution is still an unachieved goal. Besides the advances

*Corresponding author (email: [email protected]; [email protected]; [email protected]) © Science China Press and Springer-Verlag Berlin Heidelberg 2011

in the experimental aspect, molecular dynamics simulation plays a more and more important role in studying protein folding mechanism. However, the ability of computer simulation is limited by two factors, i.e., the accuracy of force field employed and the efficiency of phase space sampling. And often it is hard to tell the failure of folding simulation is due to the deficiency of the force field employed or the limited time scale of simulation [1]. Protein folds in microsecond generally, which is much longer than the capability of current computer simulation. Reducing the number of the degrees of freedom (DOF), such as implicit solvent model [2] and coarse-grained model [3], may extend the attainable time scale of simulation and strengthen the ability of phase chem.scichina.com

www.springerlink.com

Mei Y, et al.

Sci China Chem

space sampling. Other enhancing methods, such as metadynamics [4], may also improve the ability in finding the minimal free energy structure. As for the force field, much progress has been made in the last two decades [5–10]. However, most of the force fields, such as AMBER, CHARMM, OPLS, are based on the mean field approximation, i.e. the polarization effect is taken into account implicitly and homogeneously. This casts uncertainties on the accuracy and reliability of computational study of biomolecules, especially in long time molecular dynamics (MD) simulation studies. Although, polarizable force fields have been improved over recent years, their validities are still limited due to some difficulties in tuning of parameters [11]. Progress in molecular tailoring method has opened a new avenue for the study of biological molecules on quantum mechanical (QM) level [12–24]. Among various linear scaling quantum mechanical methods, the molecular fractionation with conjugate caps (MFCC) proposed in Zhang's group has been proved to be quite efficient for the electronic structure calculation for protein [15–17]. Recently, based on the MFCC method, a new charge scheme for protein termed the adaptive hydrogen bond-specific charge (AHBC) has been developed [25, 26]. AHBC is obtained by periodic restrained fitting to the electrostatic potential (RESP) calculated usually at density functional theory (DFT) level in continuum solvent. Thus AHBC is consistent with the AMBER force fields, and it can cooperate well with other AMBER force field parameters, like bond, angle and dihedral terms. Its applicability and advantage over traditional AMBER force field have been shown in several studies, including pKa prediction [26], structural stability [27, 28], NMR coupling [29, 30], protein/ligand binding [31] and peptide folding [25]. EK peptides are a series of 17-residue long polyalanines with glutamic and lysine residues embedded, designed by Marqusee and Baldwin [32]. In one of these peptides, the glutamic and lysine residues are spaced 4 residues or 1 residue apart (sequence: Ac-AEAAAKEAAAKEAAAKANH2), therefore named (i+4)EK peptide (EK peptide for short hereafter). They detected the intensity of circular dichroism to measure the proportion of helical residues in aqueous solution with different salt concentration and pH. And they found that the helical structure of EK peptide was very stable even at extreme pH conditions. It had 80% helical residues in 0.01 M NaCl at pH 7 and 0 °C, and increasing the ion strength or changing the pH lowered the helical probability. In their subsequent study of EK peptide and some other charge embedded polyalanine, they claimed that the folding of EK peptide was not due to the interaction of charged residues with the helix dipole nor to the hydrophobic interaction of sidechain, but to the high helical propensity of alanine residues [33, 34]. Due to its small size and simple secondary structure, it, as well as some of its related peptides, attracted several computational studies [35–43]. Wang et al. [37] successfully folded EK peptide within 10

December (2011) Vol.54 No.12

1975

ns under 0.01 M salt and pH 7 using AMBER03 force field and generalized Born (GB) solvation model proposed by Hawkins, Cramer and Truhlar and then parametrized by Tsui and Case (IGB1 in AMBER package) [44–47]. Jang et al. studied the folding of EK peptide with modified AMBER force field using replica exchange molecular dynamics (REMD) [38]. However, they did not study its dependence on salt concentration and pH. Dzubiella [40, 41] studied the salt-specific stability of this ALA rich peptide in explicit solvent, and found that Na+ had strong affinity to side chain carboxylates and backbone carbonyls, thereby weakening salt bridges and the main chain hydrogen bonds. In this work, we apply AHBC scheme to study the folding process of EK peptide under various ion strength and pH conditions. By comparing with the study using pairwise AMBER03 force field, we try to explain the failure of traditional force field in some studies of the dynamic process such as folding and large conformational change, and also to show the importance of explicit treatment of polarization effect into the MD simulations.

2 Method 2.1

Adaptive hydrogen bond-specific charge

Electron distribution varies with the conformation change. To accurately describe the electron distribution, on-the-fly charge update is necessary. Polarized atomic charges are obtained by periodic restrained fitting to electrostatic potential calculated at quantum mechanical level in continuum solvent [25]. As in our previous work [25], charge update is based on the main chain hydrogen bond variation to save the computational expense. In this regard, this charge scheme is termed AHBC. Quantum mechanical calculation is carried out using the MFCC scheme [48, 49] incorporating the Poisson-Boltzmann solvation model (MFCC-PB) [26]. Intra-protein polarization effect and solvation effect are both included in the quantum mechanical calculations, by taking them as background charges that polarize the electron distribution. Thus, AHBC should give a more accurate description of charge distribution than traditional force fields based on pairwise Coulomb potential. The whole protein is decomposed into residues with proper caps added to saturate the covalent bonds and to mimic the chemical environment. The initial electron density of each fragment is obtained through quantum mechanical calculation in gas phase, and then assembled to build the electron density of the whole protein [50]. In the QM calculation of each fragment, other residues are taken as background charges to include the intraprotein polarization effect. The RESP [51–53] program in Amber 10 package [54] is applied to atomic charge fitting to the calculated electrostatic potential for each residue. Charge flow between residues is not allowed in this fitting scheme. Thus, each residue has fixed and integer total charge. The atomic charges of the protein are

1976

Mei Y, et al.

Sci China Chem

then passed to the PB solver Delphi program [55] to generate discrete induced charges on cavity surface. The dielectric solute/solvent boundary is defined by AMBER van der Waals (vdW) radii [56] of solute atoms with a probe radius of 1.4 Å. The internal dielectric constant is set to unity, for molecular polarizability is explicitly included. The solvent dielectric constant is set to 80. Grid density is set to 4.0 grids/Å. Surface charges mimicking the solvation effect are then taken as additional background charges into the next cycle of QM calculations to fit new atomic charges. The solute and solvent polarize each other until convergence is reached when the dipole of the protein and the surface charges both converge within a pre-assigned numerical accuracy. This procedure quite follows the self-consistent reaction field theory (SCRF), except that the generation of ESP on solute-solvent surface is based on our MFCC scheme. Also, many body effect has been naturally introduced into the final atomic charges. All quantum mechanical calculations are performed at B3LYP/6-31G* level with Gaussian 03 software [57]. In this study, hydrogen bonds are checked periodically using HBplus [58]. If any main chain hydrogen bond between two residues is formed or broken, these two residues will be passed to charge update calculations, while atomic charges of other residues are kept intact. In this study, only the main chain hydrogen bond is considered, because it is accordant with the variation of secondary structure. Other criteria, like amplitude of conformation change or radius of gyration, may also be taken into account to further refine this scheme for charge update. To balance efficiency and accuracy, hydrogen bonds are checked every 1 ps. 2.2 Molecular dynamics simulation The initial linear structures of EK peptides are built by the sequence using LEaP module in AMBER 10. Titratable residues are protonated or deprotonated according to the pH (pH 2, 7 and 12). Under pH 7, glutamic acid and lysine possess 1 and +1 charge, respectively (denoted as E/K+). Under extreme pH, there is no special treatment of the system other than the protonation state of GLU or LYS. Under pH 2, both glutamic acid and lysine are protonated (E0/K+). While under pH 12, both glutamic acid and lysine are deprotonated (E/K0). Thus the salt bridges formed between 4-residue apart GLU and LYS residues under pH 7 switch to chargedipole interaction under extreme pH condition. More accurate treatment of protonation state of titrable residues such as constant pH MD [59] is ready to be applied. However the exact calculation of the protonation state along the simulation time is not a major concern. We are only interested in the structure of EK peptide at extreme pH. So we fix the protonation state of titrable residues during each MD simulation. The initial structure is first optimized with steepest descent minimization method for 10000 steps, and then further relaxed with conjugate gradient method until conver-

December (2011) Vol.54 No.12

gence. The minimized structure is then heated up to 300 K in 100 ps. Temperature is regulated using Langevin dynamics with the collision frequency set to 4 ps1. Solvation effect include using generalized Born (GB) models, denoted as IGB5 [60, 61] in AMBER package. The dielectric constants of the protein interior and of the solvent are respectively set to 1.0 and 78.5. Surface area terms, which may cause nonnative ion-pairing interactions [62], are neglected. Nonbond interactions are all counted without any cutoff. All the bonds involving hydrogen atoms are restrained by applying SHAKE algorithm [63]. Integration time step is 2 fs. Three salt concentrations are studied for pH 7, which are 0.01, 1.0 and 3.0 M. Only the screening effect of salt is taken into consideration, using the Debye-Huckel model [64]. Other impact such as steric effect or charge transfer for salt in the first solvation shell of protein cannot be included using GB model. And for pH 2 and 12, only 0.01 M salt concentration is studied. Trajectories are saved every 1 ps. Two force fields are employed in this study. The first one is AMBER03 force field and the second one is AHBC. In AHBC simulations, force field parameters other than atomic charges are kept the same as those in AMBER03 simulations. Totally, there are 10 simulations carried out (as listed in Table 1) to study the folding of EK peptide and its dependence on salt concentration and pH. And we also want to compare the results between AMBER03 charge and AHBC to see whether polarization effect plays an important role in folding of EK peptide. Due to large expense of CPU time for QM calculation, simulation time under AHBC is limited to 20 ns. While that under AMBER03 charge is much lengthier and extends to 260 ns each. Then we count the occurrence of helical residues to calculate the proportion of folded state. If any residue, as well as both of its covalent-bonded neighboring residues, is in helical conformation, this residue will be counted as helical residue. And the helical conformation is defined by two main chain dihedrals  and , if they fall in the range (90.0, 30.0) and (77.0, 17.0). The contribution of salt bridge interaction between the 4-residue spaced GLU and LYS residues to the total helical structure can be measured by the correlation function for the equilibrium trajectory as in ref. [65]: ( Di  D )( Pi  P ) CDP  , (1) ( Di  D )2 ( Pi  P )2 in which Di and Pi are the salt bridge distance and the number of helical residues at time i, and D and P are their average value respectively.

3

Results and discussion

In simulations using AMBER charge, the root mean square deviation (RMSD) of main chain heavy atoms from fully helical structure is shown in Figure 1. As can be seen, the folded state with RMSD below 2 Å can be reached very

Mei Y, et al.

Sci China Chem

1977

December (2011) Vol.54 No.12

Table 1 Averaged number of helical residues over the last 250 ns (for AMEBR charge) and 20 ns (for AHBC) of MD simulation at 300 K

Figure 1 Root mean square deviation from the fully helical structure for main chain heavy atoms in simulations under (a) 0.01 M, pH 7; (b) 1.0 M, pH 7; (c) 3.0 M, pH 7; (d) 0.01 M, pH 2 and (e) 0.01 M, pH 12 using AMBER charge.

soon (within 10 ns) under various solution conditions, and the folded state is the preferred structure, which is in agreement with Sun et al. [37]. They ran 40 trajectories using AMBER03 and IGB1 model. Although all the trajectories led to the folded state, they did not study the influence of salt and pH on the population of helical structure, which is the major concern of this study. It should be emphasized that the time for the first folding event taking place is not correlated with the experimental observation of folding speed, because atom diffusion is much faster in GB model than that in real solvent and there is only one trajectory for each condition, which is not adequate for the determination of folding speed. To obtain a more reasonable folding speed, more rigorous solvation models such as SPC/E [66], TIP4P-EW [67] or some others explicit water models should be implemented. The peptides travel between folded and unfolded states many times. So the simulation time is long enough to show the converged result. We also notice that folded state is more populated under pH 12. To give a more qualitative result, we calculate the helical content averaged over the last 250 ns, and list the result in Table 1. It can be read from the table that the averaged helical content does not show any dependence on ionic strength. The helical content does not change too much, while the ionic strength changes from 0.01 to 3.0 M. The peptide has even higher helical content in 3.0 M salt than that in 1.0 M one.

Salt conc. (M)

pH

Protonation

0.01 1.0 3.0 0.01 0.01

7 7 7 2 12

E/K+ E/K+ E/K+ E0/K+ E/K0

Helical content A03/IGB5 AHBC/IGB5 6.58 8.69 6.48 8.32 6.54 7.63 6.48 8.22 7.21 7.30

And under pH 12, the peptide has the most populated helical structure. The result is surely against the experimental observation, and it makes AMBER force field and generalized Born solvation model less reliable for the folding simulation of such a simple case. The failure may come from two sources of error, which are the mean-field fashion of AMBER force field and approximation of implicit solvent representation in generalized Born solvation model. It is of great interest to identify the main source of error, and it is valuable to include the polarization effect into MD simulation and justify its impact on folding of EK peptide. The RMSD of main chain heavy atoms from fully helical structure in simulations using AHBC scheme is depicted in Figure 2. It indicates that under all the conditions EK peptide finds the folded structure within 10 ns, which is the same as in the simulation with AMBER charge. Then, the structures fluctuate around the folded state, and show relatively strong stability. As also can be seen in Table 1, when AHBC is employed, the folded structure has larger population as using AMBER charge under all conditions. Although it is not quantitatively identical to the experiment, the result is at least qualitatively correct. High concentration of salt and extreme pH destabilizes the helical structure. Comparing the result from simulations using AMBER charge and AHBC, the difference lies only in the atomic charges, and it directly proves the importance of electronic polarization effect. In forming the helix, the accessible phase space by each residue is narrowed, i.e. folding into helix leads to entropy loss. The enthalpy benefit by the newly formed hydrogen bond may fully or partially compensate the entropy loss, depending on its strength. Without the explicit polarization effect, the enthalpy benefit by AMBER charge cannot fully compensate the entropy loss, which makes the folded structure less populated. When AHBC is employed, we noticed that the main chain amide hydrogen atom and carbonyl oxygen atom carry more charges. Thus, the main chain dipoles are enhanced and the enthalpy decrease is more beneficial, which make helical structure more favored. With increasing salt concentration, the screening effect between polar groups is more and more dominant, and it alleviates the electrostatic polarization. However, when the salt concentration increases from 0.01 to 1.0 M, the reduction of helical structure is not that pronounced, which does not quantitatively agree with the experiment. This may be due to the insufficient description of salt’s influence on the heli-

1978

Mei Y, et al.

Sci China Chem

Figure 2 Root mean square deviation from the fully helical structure for main chain heavy atoms in simulations under (a) 0.01 M, pH 7; (b) 1.0 M, pH 7; (c) 3.0 M, pH 7; (d) 0.01 M, pH 2 and (e) 0.01 M, pH 12 using AHBC scheme.

cal structure only by screening effect. Recent studies by Dzubiella show that salt has strong binding affinity with the peptide backbone, which causes large deformation of the folded structure [40, 41]. But this strong binding cannot be described by induced charges on solute-solvent interface in implicit solvent model. A more sophisticated calculation should be carried out in periodic water box with explicit salt ions. However, this calculation will be very time consuming, which hinders the application of AHBC scheme currently. At extreme pH, the helical structure probability is also lower than that in pH 7. Under pH 7, the salt bridge is very strong and may help stabilize the helical structure. But under extreme pH conditions, the salt bridge switches to charge-dipole interaction, which may loosen the side chain interaction and eventually decrease the stability of helical structure. During the equilibrium simulation (last 10 ns), the helical structure of central residues is more stable than that of the two terminals, which can be seen in the probabilities of helical structure for all the residues in various simulations as shown in Figure 3. The probabilities drop very fast from the center residues to the terminals. In all the simulations, the probabilities of helical structure differ mainly at the center residues. Along with the decreasing of the helical probability (from 0.01 M/pH 7 to 0.01 M/pH 12), the distribution changes from nearly unimodal distribution to bimodal distribution, showing that the center residues are more difficult to fold into helical structure under unfavorable conditions.

December (2011) Vol.54 No.12

This can be explained by the unfavorable enthalpy decrease caused by these destabilization effects, which cannot compensate the entropy loss in forming a long helix. Some representative snapshots in the simulation under 0.01 M salt and pH 7 are shown in Figure 4. It gives a clear illustration of the folding pathway. A small turn near the N-terminal forms at the very beginning of the simulation, and grows towards the C-terminal in the next 10 ns. In 5 ns, nearly 75 percent residues are in helical structure, and three turns are built. Small structural fluctuation undergoes in the last 10 ns, and deviation from the perfect long helix structure can be detected, like bending of helix, breaking of main chain hydrogen bonds and elongation of the distance between adjacent turns. This can be seen in Figure 5 from clustering analysis of the final 10 ns trajectory based on k-means algorithm [68]. In about 12.2% of the trajectory, EK peptide forms a bended helix. In other part of the trajectory, the conformations are very close, and it is in a good helical structure except for some residues near the C-terminal. The average salt bridge distances in equilibrium MD simulation employing AHBC scheme are shown in Table 2. When the solution becomes more and more unfavorable for helix formation, the salt bridge distance also increases, except for that in 3.0 M salt, which is smaller than that under 1 M salt. Because in 3 M salt, the peptide is much less ordered, and it undergoes more random dynamics. At extreme pH conditions, the salt bridge interaction changes to

Figure 3 Helical structure probabilities of residues in all the simulations using AHBC scheme.

Figure 4 Some representative snapshots in the simulation under 0.01 M salt and pH 7 using AHBC scheme.

Mei Y, et al.

Sci China Chem

ion-dipole interaction, which is much weaker, and the average distance increases by around 2.5 to 3 Å. At these distances, the ion-dipole interactions are weak, however helical structures remain. And the correlations between the number of helical residues and the average distance of native salt bridges are around –0.20, except for that in 3 M salt, which are not strong. These observations indicate that the main driving force of folding is on EK peptide itself, or say the intrinsic property of alanine residue, and the salt bridge interaction is merely an auxiliary stabilization potential, which agrees well with the conclusion in ref. [33]. Comparison of the result from simulations using AMBER charge and AHBC may give a clearer view of the difference in the simulations using AMBER and AHBC charge. The structure of EK peptide is a compromise among various interactions, especially the native main chain hydrogen bonds, the native salt bridges, and nonnative salt bridges (e.g. the interaction between GLU3 and LYS17 and that between GLU13 and LYS7). The first two interactions tend to hold the helical structure, while the last one tends to distort it. When nonpolarized force field is applied as AMBER charge, the native main chain hydrogen bonds are not strong enough, so the salt bridges may play more important roles. And the peptide may have large population for unfolded structure. But under extreme pH, the salt bridge interaction changes to much weaker charge-dipole interaction. Breaking the native interaction and building nonnative “salt bridge” interaction is not enthalpically favorable. Thus, under extreme pH conditions, the folded state has larger pop-

December (2011) Vol.54 No.12

1979

ulation than that in the simulation under pH 7 using AMBER charge. If the main chain hydrogen bonds are stabilized, as in the AHBC simulation, the folded structure is expected to be more populated. To prove this conjecture, we re-studied the last 10 ns trajectory of AHBC simulation under 0.01 M salt and pH 7, applying both AMBER03 and AHBC charges to calculate the “fake” interaction energy of the main chain hydrogen bonds. In AHBC simulation, the atomic charges change all the time along the simulation. However, when EK peptide is near its folded state, its atomic charges should not vary too much. For a qualitative comparison, we assign the atomic charge fitted at 5530 ps to calculate the interaction energy representing the AHBC simulation. In Figure 6, we find that using AHBC the total interaction energy of 15 main chain hydrogen bonds are about 20 kcal/mol lower than that using AMBER charge for the same trajectory. It shows that when polarization effect is included, the secondary structure is more stable than that under pairwise potential by about 1.3 kcal/mol for each of the native main chain hydrogen bond. This can explain the preference of native state in AHBC simulation.

Figure 6 Interaction energy of 15 native main chain hydrogen bonds in the last 10 ns of the AHBC simulation under 0.01 M and pH 7. (a) Applying AMBER charge; (b) applying AHBC at 5530 ps.

4 Conclusions Figure 5 Representative structures of five clusters from clustering analysis based on K-means algorithm for the final 10 ns simulation under 0.01 M salt and pH 7.

Table 2 Average salt bridge distances and the correlation between percentage of helical content and average salt bridge distance in equilibrium MD simulation using AHBC scheme Salt conc. (M) 0.01 1.0 3.0 0.01 0.01

pH 7 7 7 2 12

Distance (Å) 6.15 6.91 6.64 9.12 8.73

Correlation 0.21 0.20 0.11 0.20 0.24

In this work, we applied both AMBER03 charge and adaptive hydrogen-bond specific charge schemes to the folding simulation of EK peptide at various conditions to study the influence of salt concentration and pH. The result shows that although using AMBER03 charge, helical structure is preferred, its dependence on salt concentration and pH is qualitatively wrong. When polarization effect is included as in the application of AHBC while keeping other parameters intact, the calculated population of helical structure agrees qualitatively with the measured intensity of absorptions. And under all conditions, folding begins at the N-terminal and extends toward the C-terminal. When at equilibrium, the stability of center residues is more pronounced than that

1980

Mei Y, et al.

Sci China Chem

of the terminals. The result also shows that the salt bridge interaction between 4-residue apart GLU and LYS residues is not critical to the stability of helical structure of EK peptide, but is merely an auxiliary factor, which is also in agreement with the experiment. We also observe that the destabilization effect from salt ions is not that sensitive to the ion strength in this study, and this may be due to the incomplete description of salt by only its screening effect. More accurate study with explicit solvent and ions can be carried out later to study the salt effect. This work again shows that explicit treatment of polarization effect is essential in molecular dynamics simulation. This method can be further refined by considering the internal energy change due to charge variation, which is often confronted by molecular mechanics using discrete point charges to represent continuous distribution of electrons. MEI Ye is supported by the National Natural Science Foundation of China (20803034) and the Shanghai Rising-Star Program. ZHANG DaWei is supported by the Nanyang Technological University Start-up (M58110043). ZHANG ZH John acknowledges financial support from the National Natural Science Foundation of China (20933002) and Shanghai Pujiang Program (09PJ1404000). We also thank the High Performance Computer Center of East China Normal University and Shanghai Supercomputer Center for CPU time support.

December (2011) Vol.54 No.12

13 14

15

16 17

18 19

20

21

22 23

24 1 2 3 4

5

6

7

8

9

10

11 12

Freddolino PL, Park S, Roux B, Schulten K. Force field bias in protein folding simulations. Biophys J, 2009, 96: 3772–3780 Cramer CJ, Truhlar DG. Implicit solvation models: Equilibria, structure, spectra, and dynamics. Chem Rev, 1999, 99: 2161–2200 Voth GA. Coarse-Graining of Condensed Phase and Biomolecular Systems. Boca Raton: CRC Press, 2009 Laio A, Gervasio FL. Metadynamics: A method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science. Rep Prog Phys, 2008, 71: 126601 Lifson S, Warshel A. Consistent force field for calculations of conformations, vibrational spectra, and enthalpies of cycloalkane and n-alkane molecules, J Chem Phys, 1968, 49: 5116–5129 Brooks BR, Bruccoleri RE, Olafson BD, States DJ, Swaminathan S, Karplus M. CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. J Comput Chem, 1983, 4: 187–217 Pearlman DA, Case DA, Caldwell JW, Ross WS, Cheatham III TE, DeBolt S, Ferguson D, Seibel G, Kollman PA. 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, 1995, 91: 1–41 Hagler AT, Huler E, Lifson S. Energy functions for peptides and proteins. I. Derivation of a consistent force field including the hydrogen bond from amide crystals. J Am Chem Soc, 1974, 96: 5319–5327 Nemethy G, Pottle MS, Scheraga HA. Energy parameters in polypeptides. 9. Updating of geometrical parameters, non-bonding interactions and hydrogen bonding interactions for naturally occurring amino acids. J Phys Chem, 1983, 87: 1883–1887 Berendsen HJC, Postma JPM, van Gunsteren WF, DiNola A, Haak JR. Molecular dynamics with coupling to an external bath. J Chem Phys, 1984, 81: 3684–3690 Jorgensen WL. Special issue on polarization. J Chem Theor Comput, 2007, 12: 1877–1877 Nakano T, Kaminuma T, Sato T, Akiyama Y, Uebayasi M, Kitaura K. Fragment molecular orbital method: Application to polypeptides. Chem Phys Lett, 2000, 318: 614–618

25

26

27 28 29 30 31

32 33

34 35

36

Gao JL. Toward a molecular orbital derived empirical potential for liquid simulation. J Phys Chem B, 1997, 101: 657–663 Fedorov DG, Kitaura K, Li H, Jensen JH, Gordon MS. The polarizable continuum model (PCM) interfaced with the fragment molecular orbital method (FMO). J Comput Chem, 2006, 27: 976–985 Zhang DW, Xiang Y, Zhang JZH. New advance in computational chemistry: Full quantum mechanical ab initio computation of streptavidin-biotin interaction energy. J Phys Chem B, 2003, 107: 12039– 12041 Chen XH, Zhang YK, Zhang JZH. An efficient approach for ab initio energy calculation of biopolymers. J Chem Phys, 2005, 122: 184105 He X, Zhang JZH. The generalized molecular fractionation with conjugate caps/molecular mechanics method for direct calculation of protein energy. J Chem Phys, 2006, 124: 184703 Deev V. Collins MA. Approximate ab initio energies by systematic molecular fragmentation. J Chem Phys, 2005, 122: 154102 Li SH, Li W, Fang T. An efficient fragment-based approach for predicting the ground-state energies and structures of large molecules. J Am Chem Soc, 2005, 127: 7215–7226 Bettens RPA, Lee AM. A new algorithm for molecular fragmentation in quantum chemical calculations. J Phys Chem A, 2006, 110: 8777– 8785 Wang B, Merz KM. A fast QM/MM (Quantum Mechanical/Molecular Mechanical) approach to calculate nuclear magnetic resonance chemical shifts for macromolecules. J Chem Theor Comput, 2006, 2: 209– 215 Xie WS, Gao JL. Design of a next generation force field: The X-POL potential. J Chem Theor Comput, 2007, 3: 1890–1900 Xie WS, Song LC, Truhlar DG, Gao JL. Incorporation of a QM/MM buffer zone in a variational double self-consist field method. J Phys Chem B, 2008, 112: 14124–14131 Gao JL, Cembran A, Mo YR. Generalized X-POL theory and charge delocalization states. J Chem Theor Comput, 2010, 6: 2402–2410 Duan LL, Mei Y, Zhang DW, Zhang QG, Zhang JZH. Folding of a helix at room temperature is critically aided by electrostatic polarization of intraprotein hydrogen bonds. J Am Chem Soc, 2010, 132: 11159–11164 Ji CG, Mei Y, Zhang JZH. Developing polarized protein-specific charges for protein dynamics: MD free energy calculation of pKa shifts for Asp26/Asp20 in thioredoxin. Biophys J, 2008, 95: 1080– 1088 Ji CG, Zhang JZH. Protein polarization is critical to stabilizing AF-2 and helix-2' domains in ligand binding to PPAR. J Am Chem Soc, 2008, 130: 17129–17133 Duan LL, Mei Y, Zhang QG, Zhang JZH. Intra-protein hydrogen bonding is dynamically stabilized by electronic polarization. J Chem Phys, 2009, 130: 115102 Tong Y, Ji CG, Mei Y, Zhang JZH. Simulation of NMR data reveals that protein's local structures are stabilized by electronic polarization. J Am Chem Soc, 2009, 131: 8636–8641 Ji CG, Zhang JZH. NMR scaling coupling constant reveals that intraprotein hydrogen bonds are dynamically stabilized by electronic polarization. J Phys Chem B, 2009, 113: 13898–13900 Tong Y, Mei Y, Li YL, Ji CG, Zhang JZH. Electrostatic polarization makes a substantial contribution to free energy of avidin-biotin binding. J Am Chem Soc, 2010, 132: 5137–5142 Marqusee S, Baldwin RL. Helix stabilization by glu…lys+ salt bridges in short peptides of de novo design. Proc Natl Acad Sci USA, 1987, 84: 8898–8902 Marqusee S, Robbins VH, Baldwin RL. Unusually stable helix formation in short alanine-based peptides. Proc Natl Acad Sci USA, 1989, 86: 5286–5290 Baldwin RL. In search of the energetic role of peptide hydrogen bonds. J Biol Chem, 2003, 278: 17581–17588 Ghosh T, Garde S, Garcia AE. Role of backbone hydration and salt-bridge formation in stability of α-helix in solution. Biophys J, 2003, 85: 3187–3193 Chowdhury S, Zhang W, Wu C, Xiong G, Duan Y. Breaking nonnative hydrophobic clusters is the rate-limiting step in the folding of

Mei Y, et al.

37

38

39

40 41 42 43

44

45

46

47

48 49 50 51

52

53

54

Sci China Chem

an alanine-based peptide. Biopolymers, 2003, 68: 63–75 Wang WZ, Lin T, Sun YC. Examination of the folding of a short alanine-based helical peptide with salt bridges using molecular dynamics simulation. J Phys Chem B, 2007, 111: 3508–3514 Jang S, Kim E, Pak Y. Direct folding simulation of α-helices and -hairpins based on a single all-atom force field with an implicit solvation model. Proteins, 2007, 66: 53–60 Xiong K, Asciutto EK, Madura JD, Asher SA, Salt dependence of an -helical peptide folding energy landscapes. Biochemistry, 2009, 48: 10818–10826 Dzubiella J. Salt-specific stability and denaturation of a short saltbridge-forming -helix. J Am Chem Soc, 2008, 130: 14000–14007 Dzubiella J. Salt-specific stability of short and charged alanine-based -helices. J Phys Chem B, 2009, 113: 16689–16694 Yan ZQ, Wang J, Wang W. Folding and dimerization of the ionic peptide EAK 16-IV. Proteins, 2008, 72: 150–162 Zou DW, Tie ZX, Lu CM, Qin M, Lu XM, Wang M, Wang W, Chen P. Effects of hydrophobicity and anions on self-assembly of the peptide EMK16-II. Biopolymers, 2010, 93: 318–329 Hawkins GD, Cramer CJ, Truhlar DG. Pairwise solute descreening of solute charges from a dielectric medium. Chem Phys Lett, 1995, 246: 122–129 Hawkins GD, Cramer CJ, Truhlar DG. Parametrized models of aqueous free energies of solvation based on pairwise descreening of solute atomic charges from a dielectric medium. J Phys Chem, 1996, 100:19824–19839 Tsui V, Case DA. Molecular dynamics simulations of nucleic acids with a generalized born solvation model. J Am Chem Soc, 2000, 122: 2489–2498 Tsui V, Case DA. Theory and applications of the generalized born solvation model in macromolecular simulations. Biopolymers, 2000, 56: 275–291 Zhang DW, Zhang JZH. Molecular fractionation with conjugate caps for full quantum mechanical calculation of protein-molecule interaction energy. J Chem Phys, 2003, 119: 3599–3605 Mei Y, Zhang DW, Zhang JZH. New method for direct linear-scaling calculation of electron density of proteins. J Phys Chem A, 2005, 109: 2–5 Gao AM, Zhang DW, Zhang JZH, Zhang YK. An efficient linear scaling method for ab initio calculation of electron density of proteins. Chem Phys Lett, 2004, 394: 293–297 Bayly CI, Cieplak P, Cornell WD, Kollman PA. A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: The RESP model. J Phys Chem, 1993, 97: 10269– 10280 Cieplak P, Cornell WD, Bayly C, Kollman PA. Application of the multimolecule and multiconformational RESP methodology to biopolymers: Charge derivation for DNA, RNA, and proteins. J Comput Chem, 1995, 16: 1357–1377 Cornell WD, Cieplak P, Bayly CI, Kollman PA. Application of RESP charges to calculate conformational energies, hydrogen bond energies, and free energies of solvation. J Am Chem Soc, 1993, 115: 9620– 9631 Case DA, Darden TA, Cheatham III TE, Simmerling CL, Wang J, Duke RE, Luo R, Crowley R, Walker RC, Zhang W, Merz KM, Wang B, Hayik S, Roitberg A, Seabra G, Kolossvary I, Wong KF, Paesani F, Vanicek J, Wu X, Brozell SR, Steinbrecher T, Gohlke H, Yang L, Tan C, Mongan J, Hornak V, Cui G, Seetin MG, Sagui C,

December (2011) Vol.54 No.12

55

56

57

58 59

60

61

62 63

64

65 66 67

68

1981

Babin V, Kollman PA. Amber 10. 2008 Rocchia W, Sridharan S, Nicholls A, Alexov E, Chiabrera A, Honig B. 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, 2002, 23: 128–137 Cornell WD, Cieplak P, Bayly CI, Gould IR, Merz KM, Ferguson DM, Spellmeyer DC, Fox T, Caldwell JW, Kollman PA. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J Am Chem Soc, 1995, 117: 5179–5197 Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, Montgomery Jr. JA, Vreven T, Kudin KN, Burant JC, Millam JM, Iyengar SS, Tomasi J, Barone V, Mennucci B, Cossi M, Scalmani G, Rega N, Petersson GA, Nakatsuji H, Hada M, Ehara M, Toyota K, Fukuda R, Hasegawa J, Ishida M, Nakajima T, Honda Y, Kitao O, Nakai H, Klene M, Li X, Knox JE, Hratchian HP, Cross JB, Bakken V, Adamo C, Jaramillo J, Gomperts R, Stratmann RE, Yazyev O, Austin AJ, Cammi R, Pomelli C, Ochterski JW, Ayala PY, Morokuma K, Voth GA, Salvador P, Dannenberg JJ, Zakrzewski VG, Dapprich S, Daniels AD, Strain MC, Farkas O, Malick DK, Rabuck AD, Raghavachari K, Foresman JB, Ortiz JV, Cui Q, Baboul AG, Clifford S, Cioslowski J, Stefanov BB, Liu G, Liashenko A, Piskorz P, Komaromi I, Martin RL, Fox DJ, Keith T, Al-Laham MA, Peng CY, Nanayakkara A, Challacombe M, Gill PMW, Johnson B, Chen W, Wong MW, Gonzalez C, Pople JA. Gaussian 03, Revision E.01. 2004 McDonald IK, Thornton JM. Satisfying hydrogen-bonding potential in proteins. J Mol Bio, 1994, 238: 777–793 Mongan J, Case DA, McCammon JA. Constant pH molecular dynamics in generalized born implicit solvent. J Comput Chem, 2004, 25: 2038–2048 Onufriev A, Bashford D, Case DA. Exploring protein native states and large-scale conformational changes with a modified generalized born model. Proteins, 2004, 55: 383–394 Feig M, Onufriev A, Lee MS, Im W, Case DA, Brooks CL. Performance comparison of generalized Born and Poisson methods in the calculation of electrostatic solvation energies for protein structures. J Comput Chem, 2004, 25: 265–284 Lwin TZ, Zhou RH, Luo R. Is Poisson-Boltzmann theory insufficient for protein folding simulations? J Chem Phys, 2006, 124: 034902 Ryckaert JP, Ciccotti G, Berendsen HJC. Numerical integration of the Cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes. J Comput Phys, 1977, 23: 327 Srinivasan J, Trevathan MW, Beroza P, Case DA. Application of a pairwise generalized born model to proteins and nucleic acids: Inclusion of salt effects. Theor Chem Acc, 1999, 101: 426–434 Leach AR. Molecular Modelling: Principles and Applications. 2nd ed. New York: Prentice Hall, 2001 Berendsen HJC, Grigera JR, Straatsma TP. The missing term in effective pair potentials. J Phys Chem, 1987, 91: 6269–6271 Horn HW, Swope WC, Pitera JW, Madura JD, Dick TJ, Hura GL, Head-Gordon T. Development of an improved four-site water model for biomolecular simulations: Tip4p-EW. J Chem Phys, 2004, 120: 9665–9678 Shao JY, Tanner SW, Thompson N, Cheatham III TE. Clustering molecular dynamics trajectories. 1. Characterizing the performance of different clustering algorithms. J Chem Theor Comput, 2007, 3: 2312–2334