A New Implicit Solvent Model for Brownian Dynamics Simulation: Solvent-Accessible Surface Area Dependent Effective Charge Model
Tadashi ANDO, Toshiyuki MEGURO and Ichiro YAMATO
To simulate a folding of a protein by molecular simulation technique, it is necessary to improve force fields that represent the protein-solvent system. The solvent environment has profound influences on the structures, thermodynamics, dynamics and functions of biological molecules . The proper representation of this solvation effect is one of the most important challenges for the simulation of biomolecule-water systems. Explicit solvent models have been extensively used in the simulations. In principle an explicit all-atom solvent model would include all solvation effects. However, the explicit treatment of thousands of solvent degrees of freedom requires several orders of magnitude more CPU time than the corresponding in vacuo simulation. This strongly inhibits the simulations of large-scale conformational transitions such as protein folding, which takes place in milliseconds or sub-milliseconds even for small proteins . In this respect, many implicit continuum solvent models have been developed [3 - 5].
The finite difference solution of the Poisson-Boltzmann (PB) equations based upon rigorous physics has been widely used for the calculation of the electrostatic contribution to the solvation energy as the most accurate continuum model . Until recently, use of the PB approach has been largely restricted to calculations involving static representations of molecular structure, but the recent development of simple methods to obtain solvation forces from the PB equation made it possible to use it in molecular dynamics (MD) simulations . However, obtaining accurate numerical solution of PB equation for a large system such as a protein is still too costly to permit useful long time dynamics of biological molecules to be routinely studied. Another simpler and faster implicit continuum model based on PB equation is the so-called generalized Born/solvent-accessible surface area (GB/SA) model introduced by Still and co-workers . This GB/SA model has been developed and applied for molecular dynamics simulations of proteins and nucleic acids [9 - 11]. In more simplified implicit solvent models, all contributions to the solvation free energy are assumed to be proportional to the solvent-accessible surface area [12 - 14] or Gaussian solvent exclusion volume .
Recently, we have developed an atomistic Brownian dynamics (BD) with multiple time step algorithm for the simulations that require long time calculations such as protein folding [16, 17]. BD simulation is one of the stochastic dynamics simulations in which the dynamical aspects of the solvent are treated as a dissipative random force in the equations of motion. Furthermore, by eliminating the less interesting high-frequency motions such as those arising from bond vibrations, one can use larger time steps in BD simulations [18 - 20]. So, BD simulation enables much longer simulations to be performed. In the BD algorithm, a protein was described by united-atom model with AMBER91 force field and solvation effect was introduced by distance-dependent dielectric and solvent-accessible surface area (DD/SA) model. This BD simulation of a peptide was 160 times faster than molecular dynamics simulation with explicit solvent and stable.
In this paper, we introduce a new simple implicit solvent model, effective charge (EC) model to improve the force field. In the EC model, the atomic charge is neutralized as a function of the solvent-accessible surface area of the atom. The GB/SA model was also applied for BD algorithm. In order to test the effectiveness of these implicit solvent models, BD simulations of a 28-mer bba folded peptide were carried out. A corresponding molecular dynamics simulation using explicit all-atom water model was also performed for comparison. In vacuo simulations are well known to suffer from several artifacts, such as smaller radius of gyration, shrinkage of the hydrophilic solvent-accessible surface area and smaller atomic fluctuations [1, 14]. We compared those structural and dynamics properties of the peptide in these simulations.
2 Theory and Methods
2. 1 Solvation Energy
In this study, four implicit solvation models were used: (1) distance-dependent dielectric (DD) model; (2) solvent-accessible surface area (SA) model; (3) effective charge (EC) model; (4) generalized Born/solvent-accessible surface area (GB/SA) model. We briefly summarize each implicit solvent model.
(1) DD model
DD model assumes the dielectric constant e = Crij where C is a constant value and rij is the distance between atoms i and j. This model is approximating the long-range electrostatic screening effects . In this study, C of 2 (i.e. e = 2rij) was used.
(2) SA model
In SA model, the solvation energy is given by:
where SAi(rN) is the solvent-accessible surface area of atom i and si is a solvation parameter depending on the atom type. This model is based on the assumptions that most of the solvation energy arises from the first water shell around the protein and that the energy of interaction of a solute with water can be considered as a sum of energies of atomic groups . The SAi(rN) was calculated using the approximate analytical expression of Hassel et al.  and the atomic solvation parameters determined by Wesson and Eisenberg ; s(C) = 12 cal/mol/A2, s(O, N) = -116 cal/mol/A2, s(S) = -18 cal/mol/A2, s(O-) = -175 cal/mol/A2 and s(N+) = -186 cal/mol/A2 were used.
(3) EC model
We introduce a new implicit solvent model, Effective Charge (EC) model. Water is a polar molecule with a dipole moment of 1.85 debye. Therefore, the point charge in water solution is surrounded by oriented water molecules and the point charge is shielded. To represent this shielding effect, atomic charge of atom i, qi, is neutralized as a function of solvent-accessible surface area of the atom, SAi(rN);
Here qi' is the effective charge of atom i, Si is the total solvent-accessible surface area of isolated atom i, aint is a shielding parameter against interior of the solute (wherein aint is set at unity) and aext is a shielding parameter for exterior water. Therefore, in this model effective atomic charge decreases to a factor of 1/aext linearly with the increase of solvent-accessible surface area. In this study, aext = 8 was used.
This EC model is similar to the solvent-accessibility modified Tanford-Kirkwood method developed by Matthew and Gurd . In their model, the electrostatic energy of pairwise interaction is calculated as the value of the Tanford-Kirkwood electrostatic interaction energy reduced by average solvent-accessibility (normalized by the corresponding accessible surface area of a model tripeptide) over atoms i and j.
(4) GB/SA model
In the GB/SA model , the total solvation energy (Gsolv) is given as the sum of a solvent-solvent cavity term (Gcav), a solute-solvent van der Waals term (GvdW) and a solute-solvent electrostatic polarization term (Gpol):
The first two terms are the non-electrostatic contributions and are approximately linearly related to their solvent-accessible surface area (SA). This is represented by the sum:
where gi is the atomic solvation parameter of atom i and SAi(rN) is the solvent-accessible surface area of atom i. This SAi(rN) was calculated using the same method as used in SA model. In this study, g(C(sp3), S) = 10 cal/mol/A2, g(C(sp2), C(sp)) = 7 cal/mol/A2 and g(N, O) = 0 cal/mol/A2 were used, in which the united atom model is used in SA calculation, thus g(H) = 0 cal/mol/A2 .
The third term, Gpol, is the electrostatic contribution and is expressed by the following form:
is the distance between atoms i and j with atomic charges qi and qj, respectively and e is the dielectric constant of the solvent (= 80). ai is the so-called effective Born radius of atom i. This is the generalized Born (GB) equation modified by Liu and Beveridge for improved estimates of solvation energies . The effective Born radii were calculated using the pairwise dielectric descreening procedure initially introduced by Hawkins et al.  and then developed by Liu and Beveridge  for AMBER91 united-atom force field. Derivatives of Eq. (5) were calculated by treating ai as a constant at every time step of BD simulations .
In this study, the three implicit solvent models described above were used. However, DD and EC models include only the electrostatic contribution to the free energy of solvation and neglect the non-electrostatic contribution. To consider both contributions, DD and EC models were combined with SA model, that is, DD/SA and DD/SA/EC models.
2. 2 BD Simulation
In order to test the effectiveness of the implicit solvent models described above, 5 ns multiple time step BD simulations of a designed bba fold 28-mer peptide pda8d  with the implicit solvent models as described in the previous paper  were performed. This peptide is one of the smallest proteins that contain two conventional secondary structures, b-strand and a-helix, and is a useful model peptide for developing new algorithms and parameters of molecular simulation; we previously reported its 5 ns MD simulation . For multiple time step algorithm, Dt of 5 fs and n of 8 (i.e. Dt = 40 fs) were used. The simulation temperature was 280 K. Cut-off method was not used. Coordinates and energies were recorded every 10 ps during the simulation. All calculations were performed on a 2.8 GHz Pentium4 processor based on Linux.
2. 3 MD Simulation
For the comparison, an MD simulation in explicit water molecules was also performed. The MD simulation was performed using the MD program AMBER 4.1  with united-atom force field , as described in the previous paper . The simulation temperature was the same as used for BD simulation.
3. 1 Dynamics of the Peptide Using Implicit Solvent Models
In this study, we tested four implicit solvent models, SA, DD/SA, DD/SA/EC and GB/SA, with the BD algorithm. First, the simulation results on structure and dynamics obtained by various implicit solvent simulations were compared.
Figure 1 shows Ca root mean square deviations (RMSD) of the peptide from the NMR structure as a function of the simulation time. A snapshot at 5 ns period of each trajectory is shown in Figure 2. By in vacuo simulation, RMSDs were around 4 A until 5 ns. On the other hand, other implicit solvent models appeared to give mostly stable trajectory, as judged by the smaller RMSD than those obtained from in vacuo simulation. However, while the DD/SA and DD/SA/EC models gave stable trajectories of the secondary and tertiary structures of the native peptide, the structures by the SA and the GB/SA models were collapsed slightly. Additionally, the DD/SA/EC model gave the smaller RMS deviations over 5 ns simulation than those obtained from DD/SA simulation.
Figure 1. Time evolutions of the Ca RMSD during simulations. The RMS deviations from the NMR structure versus time for 5 ns simulations using various solvent models are shown.
Figure 2. Ribbon representations of the snapshots at 5 ns in various solvent model simulations. On the left side, the initial structure of the simulations is also shown as a reference. The figures are generated with MOLMOL .
Figure 3 shows the radius of gyration (Rg) of the peptide as a function of time. As a reference, the average value of Rg calculated from the 32 structure models determined by NMR was shown. The explicit solvent simulation is, in general, most reliable at present. However, the Rg from the simulation were somewhat larger than those calculated from the NMR structures. The DD/SA and DD/SA/EC models also gave larger Rg than the value of the NMR structures during 5 ns simulation similarly to the explicit solvent model. On the other hand, the SA and GB/SA models gave somewhat smaller Rg values than those calculated from the NMR structures. Rg values obtained by in vacuo simulation were much smaller.
Figure 3. Time evolutions of the radius of gyration (Rg) during 5 ns simulations using various solvent models. The average value of Rg calculated from the 32 structure models determined by NMR is shown as a reference (dashed lines).
Figures 4, 5 show the hydrophilic solvent-accessible surface area (Sphi) and hydrophobic solvent-accessible surface area (Spho) of the peptide as a function of the simulation time, respectively. As a reference, the average values of Sphi and Spho calculated from the 32 structure models determined by NMR were also shown. By the SA and GB/SA models, Sphi values decreased drastically at early stage. Especially, in the GB/SA model, this decrease in Sphi was drastic and many salt bridges were observed. On the other hand, Sphi values generated by the DD/SA and DD/SA/EC models were somewhat larger than the value calculated from NMR structures and by explicit solvent model throughout 5 ns simulation. Spho values obtained by the SA, GB/SA models and in vacuo simulation were almost the same as the NMR structure value. By the DD/SA and DD/SA/EC models, Spho values were a little bit larger than those calculated by the NMR structures but were the closest to the values obtained by explicit solvent MD simulation.
Figure 4. Time evolutions of the total hydrophilic solvent-accessible surface area (Sphi) during 5 ns simulations using various solvent models. The average values of Sphi calculated from the 32 structure models determined by NMR are shown as a reference (dashed lines).
Figure 5. Time evolutions of the total hydrophobic solvent-accessible surface area (Spho) during 5 ns simulations using various solvent models. The average values of Spho calculated from the 32 structure models determined by NMR are shown as a reference (dashed lines).
Finally, we examined the Ca positional root mean square fluctuations (RMSF) of the peptide during simulation using various solvent models (Figure 6). The fluctuations observed in the simulation using the explicit solvent model were large, especially around the N and C termini of the peptide regions and the hairpin turn (residues 7 and 8), while those observed in the simulation using implicit solvent models and in vacuo were much smaller throughout the peptide.
Figure 6. The root means square fluctuations (RMSF) of Ca atom positions obtained from simulations using various solvent models. The RMS fluctuation values were calculated as the average values over 2.5 to 5 ns period of each simulation.
3. 2 Computational Cost
It is also important to examine the computational time required for each implicit solvent model. Table 1 lists the time required to run 1 ns simulation of pda8d with various solvent models on a 2.8 GHz Pentium4 processor in a personal computer and the relative times compared to that of multiple time step BD simulation with SA model. The BD simulation with the DD/SA/EC solvent models required about 7% more computational time than the DD/SA BD simulation, and 162 (= 2,057/13.7) times less than the MD simulation with explicit water molecules with cut-off radius of 9 A. The GB/SA model was 1.4 times slower than the DD/SA/EC BD simulation, but 115 (= 2,057/171.9) times faster than the explicit water model with cut-off radius of 9 A.
Table 1. Computation Time for 1 ns Dynamics of the 28-mer bba folded peptide, pda8d.
aAll calculations were performed using a Pentium4 2.8 GHz processor.
|Algorithm||Solvent model||Time(min)a||Relative time|
bThe simulations were performed using the BD with multiple time step algorithm developed here. All covalent bonds were constrained with LINCS algorithm. The short time step of 5 fs and n = 8 were used for multiple time step algorithm. Number of atoms was 304.
cThe simulation was performed using the MD program AMBER with united-atom force field. The peptide was solvated using a box extending at least 10 A in all directions. All covalent bonds were constrained with SHAKE algorithm. Cut-off radius was 9 A. The time step was 2 fs. Number of atoms was 7,681.
In this study we have evaluated the performance of the four implicit solvent models, generalized Born/solvent-accessible surface area (GB/SA), solvent-accessible surface area (SA), SA combined with distance-dependent dielectric (DD/SA) and effective charge model combined with DD/SA models, in combination with the atomistic Brownian dynamics algorithm based on AMBER united-atom force field. The computational times required for those implicit solvent models combined with BD algorithm were much smaller than that of MD simulation using explicit solvent water model (Table 1). Therefore, the BD simulations with those implicit solvent models permit long time simulations to be routinely studied.
The main artifacts often observed in in vacuo simulations, namely smaller radius of gyration (Rg) and shrinkage of the hydrophilic solvent-accessible surface area (Sphi) were considerably reduced when the DD/SA and DD/SA/EC solvent models using e = 2rij were applied for the BD simulation in order to include solvation effects into the simulation (Figures 3, 4). Judging from the smaller RMSD obtained from DD/SA/EC simulation throughout 5 ns than those calculated from simulations using other implicit solvent models, DD/SA/EC is the most effective for representing solvent effects among the four implicit solvent models evaluated in this study. Increase on Sphi in the simulation using DD/SA and DD/SA/EC models was different from the result of explicit solvent MD and the average value calculated from NMR structures. This disagreement may be attributed to the insufficient simulation time of MD and the initial model peptide structure determined by NMR in which structural models were calculated using nuclear Overhauser enhancement data in vacuo. The SA model was less effective for the reduction of the artifacts encountered in in vacuo simulations than that using DD and EC models. The BD simulation using the SA solvent model gave erroneous salt-bridges, smaller Rg (Figure 3) and shrinkage of Sphi (Figure 4). It does not seem that this solvent model can account for long-range screening effects.
The GB/SA simulation resulted in smaller Rg, decrease of Sphi especially due to many salt bridges (Figures 3, 4). The parameters used for the GB/SA model in this study were optimized using only one peptide containing all 20 amino acids whose structures were obtained from a short-time MD simulation (1 ns) and random generation . Recently, Zhu et al. pointed out the necessity of using structurally distinct several proteins and long-time trial simulations in the parametrization of implicit solvent models . Therefore, the insufficient parametrization might cause errors in calculation of solvation energy of the system in our GB/SA simulation. In addition, although this continuum approximation has the advantage of being fast to evaluate the solvation energy, the dielectric constant which is a macroscopic property of solvent cannot properly describe the atomistic nature of solvent at short distances . Zhu et al. noted the importance of maintaining the proper balance of non-bonded interactions, especially hydrogen-bond interactions, which are effective at short distances, in the parametrization of implicit solvent models . In their work, the treatment of the protons that can form hydrogen-bonding was modified. That is, the effective Born radius of any backbone hydrogen atom was taken to be the same as that of the corresponding nitrogen atom in the parametrization of a GB/SA solvent model and the simulations with this model were comparable to explicit solvent simulations . Very recently, Zhou also found that the Generalized Born model showed erroneous salt bridge effects between charged residues resulting in an over-weighting of a non-native structure of a peptide . He noted that even with the large successes of the GB continuum solvent model in many fields, such as pKa calculation, solvation free energy calculation and ligand-receptor bindings, one must be cautious with simulations involving large-scale conformational change. Furthermore, Sheinerman and Brooks proposed that water provides a lubrication mechanism during folding . Therefore, it is necessary to modify the estimation of solvation effects at short-range where the molecular property of the water strongly influences the protein folding. Owing to these reasons, the deviations of Sphi and Rg values and over-weighting salt bridges in the GB/SA model from the values calculated from NMR structures might have occurred. In this respect, our results in this study suggest that the DD/SA/EC model is the most preferable for long time simulation of peptides with BD.
Although some artifacts encountered in in vacuo simulations were reduced by using the DD/SA and DD/SA/EC solvent using e = 2rij model, the small atomic fluctuations observed in BD simulation using the implicit solvent models were different from those with explicit solvent model (Figure 6). The same was true for the SA and the GB/SA solvent models. On the other hand, Zhu et al. reported that RMS fluctuations obtained from stochastic dynamics (Langevin dynamics) simulations with a time step of 2 fs using a GB/SA and a DD/SA model based on GROMOS96 force field were not small and were almost in good agreement with those generated by the explicit solvent simulations . Therefore, this difference may be attributed to the BD algorithm using long time step, in which the less interesting high-frequency motions were omitted.
In EC model, although the one shielding parameter aext should be optimized, the model is very simple and requires less computational cost. From what has been discussed above, we can conclude that the DD/SA/EC model is more effective for representation of solvation effects than other implicit solvent models.
Proper representation of solvation effects in molecular simulations is a crucial point to describe the dynamics of biomolecules in aqueous systems, especially for the dynamics simulation using implicit solvent model such as Brownian dynamics. In this work, we developed a new implicit solvent model, effective charge (EC) model, to take into account solvation effects in the Brownian dynamics simulation. This EC model could reduce the artifacts typically encountered in in vacuo simulations considerably. Moreover, the EC model is simple and requires no additional computation time. Therefore, this new EC model is a good candidate for implementation in the molecular simulation that investigates the dynamic process required for long time.
6 Agreement for Using the Program
This Brownian dynamics program is freeware. Please contact us by e-mail when you want to use the program.
We thank Prof. Akihide Oguchi for helpful suggestions and discussions. This work was supported (I. Y.) by a Grant-in-Aid for Scientific Research on Priority Area (C) Genome Information Science from the Ministry of Education, Culture, Sports, Science and Technology of Japan.
[ 1] C. L. Brooks III, M. Karplus and B. M. Pettitt, Solvent influence on protein dynamics, Proteins: a theoretical perspective of dynamics, structure, and thermodynamics, ed. by I. Prigogine and S. A. Rice, John Wiley & Sons, Inc., New York (1988), pp. 137-174.
[ 2] A. R. Fersht, Structure and Mechanism in Protein Science, W. H. Freeman and Company, New York (1999).
[ 3] A. Warshel and A. Papazyan, Curr. Opin. Struct. Biol., 8, 211-217 (1998).
[ 4] B. Roux and T. Simonson, Biophys. Chem., 78, 1-20 (1999).
[ 5] C. J. Cramer and D. G. Truhlar, Chem. Rev., 99, 2161-2200 (1999).
[ 6] W. Rocchia, S. Sridharan, A. Nicholls, E. Alexov, A. Chiabrera and B. Honig, J. Comput. Chem., 23, 128-137 (2002).
[ 7] A. H. Elcock, M. J. Potter and J. A. McCammon, Application of Poisson-Boltzmann solvation forces to macromolecular simulations, Computer simulation of biomolecular system: theoretical and experimental applications, Vol. 3, ed. by W. F. van Gunsteren, P. K. Weiner and A. J. Wilkinson, Kluwer Academic Publishers, Netherlands (1997), pp 244-261.
[ 8] W. C. Still, A. Tempczyk, R. C. Hawley and T. Hendrickson, J. Am. Chem. Soc., 112, 6127-6129 (1990).
[ 9] B. N. Dominy and C. L. Brooks III, J. Phys. Chem. B, 103, 3765-3773 (1999).
 V. Tsui and D. A. Case, J. Am. Chem. Soc., 122, 2489-2498 (2000).
 J. Zhu, Y. Shi and H. Liu, J. Phys. Chem. B, 106, 4844-4853 (2002).
 J. Vila, R. L. Williams, M. Vasquez and H. A. Scheraga, Proteins, 10, 199-218 (1991).
 L. Wesson and D. Eisenberg, Protein Sci., 1, 227-235 (1992).
 F. Fraternali and W. F. van Gunsteren, J. Mol. Biol., 256, 939-948 (1996).
 T. Lazaridis and M. Karplus, Proteins, 35, 133-152 (1999).
 T. Ando, T. Meguro and I. Yamato, J. Comput. Chem. Jpn., 1, 115-122 (2002).
 T. Ando, T. Meguro and I. Yamato, Mol. Simulation, 29, 471-478 (2003).
 D. L. Ermak and J. A. McCammon, J. Chem. Phys., 69, 1352-1360 (1978).
 C. L. Brooks III, M. Karplus and B. M. Pettitt, Dynamical simulation methods, Proteins: a theoretical perspective of dynamics, structure, and thermodynamics, ed. by I. Prigogine and S. A. Rice, John Wiley & Sons, Inc., New York (1988), pp. 33-58.
 W. F. van Gunsteren, Molecular dynamics and stochastic dynamics simulation: a primer, Computer simulation of biomolecular system: theoretical and experimental applications, Vol. 2, ed. by W. F. van Gunsteren, P. K. Weiner and A. J. Wilkinson, ESCOM Science Publishers B. V., Leiden (1993), pp. 3-36.
 B. R. Gelin and M. Karplus, Biochemistry, 18, 1256-1268 (1979).
 D. Eisenberg and A. D. McLachlan, Nature, 319, 199-203 (1986).
 W. Hasel, T. F. Hendrickson and W. C. Still, Tetrahedron Comput. Meth., 1, 103-116 (1988).
 J. B. Matthew and F. R. N. Gurd, Methods Enzymol., 130, 413-436 (1986).
 D. Qui, P. S. Shenkin, F. P. Hollinger and W. C. Still, J. Phys. Chem. A, 101, 3005-3014 (1997).
 Y. Liu and D. L. Beveridge, Proteins, 46, 128-146 (2002).
 G. D. Hawkins, C. J. Cramer and D. G. Truhlar, Chem. Phys. Lett., 246, 122-129 (1995).
 B. I. Dahiyat, C. A. Sarisky and S. L. Mayo, J. Mol. Biol., 273, 789-796 (1997).
 W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, J. K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell, and P. A. Kollman, J. Am. Chem. Soc., 117, 5179-5197 (1995).
 S. J. Weiner, P. A. Kollman, D. A. Case, U. C. Singh, C. Ghio, G. Alagona, S. J. Profeta and P. Weiner, J. Am. Chem. Soc., 106, 765-784 (1984).
 P. E. Smith and W. F. van Gunsteren, Methods for the evaluation of long range electrostatic force in computer simulations of molecular systems, Computer simulation of biomolecular system: theoretical and experimental applications, Vol. 2, ed. by W. F. van Gunsteren, P. K. Weiner and A. J. Wilkinson, ESCOM Science Publishers B. V., Leiden (1993), pp. 182-212.
 R. Zhou, Proteins, 53, 148-161 (2003).
 F. B. Sheinerman and C. L. Brooks III, J. Mol. Biol., 278, 439-456 (1998).
 R. Koradi, M. Billeter and K. Wuthrich, J. Mol. Graphics, 14, 51-55 (1996).