Research Article  Open Access
Effects of Water Models on Binding Affinity: Evidence from AllAtom Simulation of Binding of Tamiflu to A/H5N1 Neuraminidase
Abstract
The influence of water models SPC, SPC/E, TIP3P, and TIP4P on ligand binding affinity is examined by calculating the binding free energy of oseltamivir carboxylate (Tamiflu) to the wild type of glycoprotein neuraminidase from the pandemic A/H5N1 virus. is estimated by the Molecular MechanicPoisson Boltzmann Surface Area method and allatom simulations with different combinations of these aqueous models and four force fields AMBER99SB, CHARMM27, GROMOS96 43a1, and OPLSAA/L. It is shown that there is no correlation between the binding free energy and the water density in the binding pocket in CHARMM. However, for three remaining force fields decays with increase of water density. SPC/E provides the lowest binding free energy for any force field, while the water effect is the most pronounced in CHARMM. In agreement with the popular GROMACS recommendation, the binding score obtained by combinations of AMBERTIP3P, OPLSTIP4P, and GROMOSSPC is the most relevant to the experiments. For wildtype neuraminidase we have found that SPC is more suitable for CHARMM than TIP3P recommended by GROMACS for studying ligand binding. However, our study for three of its mutants reveals that TIP3P is presumably the best choice for CHARMM.
1. Introduction
The determination of binding affinity is a central problem in computeraided drug design which is a useful tool to search for potential leads for various diseases. The accuracy of estimation of the binding free energy of ligand to receptor depends on computational methods and modeling of receptorligand complex. The docking method is usually used for locating binding sites and virtual screening of potential drug candidates from large data bases. This approach suffers from low accuracy and its results usually have to be refined by more sophisticated methods based on the molecular dynamics (MD) simulation. In many cases MD methods can reproduce reliable results on binding free energy having acceptable correlation with experimental data [1ā5]. Among them Molecular MechanicPoisson Boltzmann Surface Area (MMPBSA) [1], thermodynamic integration (TI) [2], linear interaction energy (LIE) [3], linear response approximation (LRA) [6, 7], free energy perturbation (FEP) [4], and steered MD [8, 9] methods are widely used. Each method should be considered carefully to compromise between CPU time efficiency and accuracy level.
In modeling of biosystems in aqueous environment, it is important to develop appropriate force fields and water models. Force fields, which are given in the form of empirical potential energy functions, have been developed by different groups. Today OPLS, CHARMM, AMBER, and GROMOS are the most popular force fields for allatom simulation of biomolecules. To describe aqueous environment, one can use various models such as SPC [10], SPC/E [11], TIP3P [12], and TIP4P [13]. The adjustment of parameters of these models is based on their ability to reproduce the enthalpy of vaporization and the density of water. SPC/E is especially accurate for capturing experimental properties of water such as the diffusion coefficient and dielectric constant. Because SPC, SPC/E, TIP3P, and TIP4P are relatively simple and able to provide reasonable results, they are often employed for simulation of peptides [14, 15] and proteins [16].
Previous studies [17, 18] have revealed that different force fields provide different estimations for . Recently, the role of water molecules in the binding process has been considered [19ā21]. GROMACS manual (http://www.gromacs.org/Support/Online_Manual) suggests that for allatom simulation of biomolecules water model TIP3P is suitable for AMBER and CHARMM, while TIP4P and SPC are more appropriate, respectively, for OPLS and GROMOS force fields. However, what water model is the best fit for a given force field in computation of remains largely unknown.
To shed light on this problem, in this paper we study the impact of combinations of four main water models SPC [10], SPC/E [11], TIP3P [12], and TIP4P [13] with AMBER99SB [24], CHARMM27 [25], OPLSAA [26], and GROMOS96 43a1 [27] force fields on the binding affinity of Tamiflu to the wild type (WT) of A/H5N1 neuraminidase (NA). We have chosen the NATamiflu complex because A/H5N1 virus causes great damage to live poultry markets [28], especially being recognized as human transmitted virus [29]. More importantly, the binding free energy of Tamiflu to NA has been experimentally determined [23] and this gives us the opportunity to compare theoretical estimates with the experimental ones.
Using MMPBSA method we have shown that combinations AMBERTIP3P, OPLSTIP4P, and GROMOSSPC are the best choice for estimation of of Tamiflu. This result is in agreement with the GROMACS recommendation, which is followed from force field development [24ā27]. Contrary to the GROMACS suggestion, it is shown that SPC is more suitable for CHARMM than TIP3P but this conclusion is valid for the wild type of NA. Our study of three mutants H274Y, N294S, and Y252H reveals that TIP3P is presumably the best choice for CHARMM as suggested by GROMACS.
It is found that obtained by the OPLS force field is much less sensitive to water models compared to other force fields. The difference in water models seems to have the drastic effect in CHARMM modulating both the receptorligand interaction energy and hydrogen bond network in binding area. For all studied force fields, SPC/E is worse than other aqueous models overestimating .
2. Materials and Methods
2.1. Crystal Structure of A/H5N1 NA and Parametrization of Tamiflu
The initial structures of A/H5N1 WT and mutants H274Y and N294S were obtained from Protein Data Bank with PDB ID 2HU4, 3CL0, and 3CL2, respectively [23]. Y252H was derived by the corresponding mutation in WT structure using the mutagenesis module, integrated in PyMOL package [30]. For Tamiflu we use the oseltamivir carboxylate type. Its charges and atom types, used for MD simulation, are described in detail in our previous work [18]. Namely, for unitedatom GROMOS96 43a1 force field, charges and atom types of oseltamivir were fully parametrized by Dundee PRODRG2.5 Server (Beta) [31]. For the remaining allatom force fields, atomic partial charges for Tamiflu were derived by ESP charge. To obtain optimal geometry for electrostatic potential calculations, its structure is first optimized with the help of Gaussian98 [32] using the B3LYP/631G* level of theory. Fitting charges to the electrostatic potential was subsequently done by the RESP method. Atom types for Tamiflu were derived from different modules to get along with each force field. ACPYPE [33] and MKTOP [34] were adjusted to provide suitable atom types in OPLSAA/L [26], while for AMBER99SB and CHARMM27 [25], atom types were named by ACPYPE and SwissParam [35], respectively.
2.2. Water Models
Water, known as an indispensable solvent in almost chemical and biological reactions, has been built in different ways to obtain reasonable models for computational study [10ā13]. A water molecule is characterized by its geometrical parameters such as bond lengths and angles which could be kept rigid or flexible during simulation. Each model is parametrized with atomic partial charges of oxygen and hydrogen and assigned with dispersion and repulsion forces approximated by LennardJones potential [22]. Water models are categorized by the number of points used to shape them, by rigid or flexible structures, and by integration or not with polarization effects [22]. There are 46 distinct models [36] that have from 3 to 6 sites. However, only 3 and 4site models (Figure 1) are often used in simulations of biological systems. For 3site model like SPC [10] and TIP3P [12], a water molecule is constructed by one oxygen atom and two hydrogen atoms. Each atom is assigned with atomic partial charge, but only oxygen is allowed to have the LennardJones interaction with other atoms. The van der Waals (vdW) interaction among hydrogen atoms was not parametrized yet. Threesite models are known as rigid and have the experimental geometry of water, except SPC which has the ideal tetrahedral angle of water as 109.47, but not 104.5Ā°. SPC/E [11], an updated version of SPC model, adds an average polarization correction to the potential energy function, resulting in the better density and diffusion constant than SPC model.
(a)
(b)
The foursite model TIP4P [13] is a rigid planar foursite interaction potential for water (Figure 1), having a similar geometry to the Bernal and Fowler model [37]. Here the negative charge is shifted from the oxygen atom to a point 0.15āĆ along the bisector between hydrogen atoms. In this paper, we just limit our study to only four frequently used water models SPC, SPC/E, TIP3P, and TIP4P (Figure 1). Their geometrical and physical characteristics are described in Table 1. Here , , and are the partial charges of hydrogen, oxygen, and lone pair, respectively. and are the HāOāH and lone pairOāH angles, while and are the well depth and vdW radius, respectively.

2.3. Molecular Dynamic Simulations
Complex of NATamiflu is placed in a triclinic box of around 12000 water molecules with 1ānm distance between the solute and box (a typical snapshot is shown in Figure S1 of supplementary material available online at http://dx.doi.org/10.1155/2014/536084). The receptor and ligand have 3832 and 5749 atoms in the united atom and allatom models, respectively. Periodic boundary condition is imposed in three directions. We use 1.4ānm and 1.0ānm cutoff for vdW and electrostatic interactions, respectively. Long range electrostatic interaction was computed by the particlemesh Ewald summation method [38]. Equations of motion were integrated using a leapfrog algorithm [39] with a time step 1āfs. The nonbonded interaction pairlist was updated every 10āfs with the cutoff of 1ānm. All systems were neutralized by adding counterions and then minimized to remove the local strain in protein upon addition of all hydrogen atoms and to remove bad vdW contacts with water. By using the conjugate gradient method for every 50 steps of steepest descent, minimization is converged when maximum force becomes smaller than 0.01ākJ/mol/nm. Then, all bonds of protein were restrained, leaving remaining parts to relax for 100āps to obtain evenly distributed system. The temperature was gradually heated to 300āK during 100āps with 5ākcal/mol harmonic restraints in all systems. The equilibration was next performed, coupled with temperature and pressure. Constant temperature 300āK was enforced using Berendsen algorithm [40] under 50āps NVT simulation with a damping coefficient of 0.1āps. ParrinelloRahman pressure coupling [41] was used in 100āps NPT run for 1āatm with the damping coefficient of 0.5āps. Final NPT simulations of 20āns were carried out for all force fields. Each force field is combined with four different water models SPC, SPC/E, TIP3P, and TIP4P, except GROMOS, which uses only SPC and SPC/E models. In total, we have 14 different models for the NATamiflu complex. All simulations have been carried out in the GROMACS suit with Gromacs4.5 package [42].
2.4. Binding Free Energy Calculation by MMPBSA
The details of MMPBSA method are given in our previous work [18]. Overall, in this method the binding free energy of ligand to receptor is defined as follows: where and are contributions from electrostatic and vdW interactions, respectively. and are nonpolar and polar solvation energies. The entropic contribution is estimated using the normal mode approximation [18]. From 20āns MD simulation output only snapshots collected in equilibrium are used to compute the binding free energy given by (1).
2.5. Measures Used in Data Analysis
The rootmeansquare deviation (RMSD) is employed to measure the deviation of receptor structures from the initial configuration. The hydrogen bond (HB) is assumed to be formed if the distance between proton donor (D) and proton acceptor (A) is less than 0.35ānm and the angle HDA is also less than 30Ā°.
2.6. Definition of Binding Site
The binding pocket is defined as a space surrounded by 50 amino acids as shown in Figure S2. Our definition is compatible with that of Cheng et al. [43]. The list of these amino acids is given in Table 2. Volume of the binding pocket is approximately estimated as volume of smallest convex hull which contains all of the fifty atoms [44ā46] (Figure 2). The number of water molecules inside the convex hull is considered as the number of water molecules in the binding site. The binding pocket volume and number of water molecules inside it are calculated by Matlab software [47]. The water density in binding site is a ratio of the number of water molecules to its volume.

3. Results and Discussion
In this section we present results obtained for WT NA if not otherwise stated.
3.1. Equilibration Time Scales Depend on Force Field and Water Model
The time dependence of RMSD of NA in complex with Tamiflu is shown in Figure 3 for different sets of force fields and water models. The equilibration of system is reached when RMSD gets saturation. In AMBER the equilibration time āns for SPC, while ā10āns is needed to equilibrate the system in other water models (Figure 3). SPC gives the largest RMSD in equilibrium state. In OPLS RMSD reaches saturation quite fast, after about 3āns for TIP3P and 5āns for remaining models (Figure 3). TIP3P provides a bit larger departure from the initial structure compared to remaining models.
(a)
(b)
(c)
(d)
For CHARMM āns for all sets (Figure 3). The effect of water on stability of the NATamiflu complex is at most pronounced in CHARMM, where TIP4P affects the stability to a greater extent than other models. In GROMOS one has different time scales for equilibration in SPC (āns) and SPC/E (āns), but in equilibrium the average values of RMSD are almost the same for both water models. Due to unitedatom nature GROMOS is the most unstable having average value of RMSD ānm against 0.13ānm of other force fields (Figure 3).
3.2. Estimation of Binding Free Energy by MMPBSA Method
3.2.1. Effect of Water Model on the ReceptorLigand Interaction Energy
The interaction energy between ligand and receptor is shown in Figure 4. The most pronounced dependence on water is observed for CHARMM as fluctuates at very different levels. Particularly, TIP3P and TIP4P models make the interaction energy highly unstable during the first 7āns, while it remains almost stable during the whole MD run for SPC and SPC/E (Figure 4). In equilibrium TIP3P and TIP4P give lower interaction energy than SPC and SPC/E. Average interaction energy , , , and ākcal/mol for SPC, SPC/E, TIP3P, and TP4P, respectively.
(a)
(b)
(c)
(d)
In GROMOS , obtained by using SPC, is higher than SPC/E. The effect of water modeling is also visible for AMBER, where SPC/E provides the lowest interaction energy in equilibrium. is almost the same in TIP3P and TIP4P (Figure 4). As in the RMSD case (Figure 3), the results, obtained by the OPLS force field, are not affected much by water models (Figure 4). The combination of GROMOS with SPC and SPC/E gives the highest receptorligand interaction energy among four force fields, while the lowest is obtained by CHARMMSPC and CHARMMSPC/E.
3.2.2. Electrostatic Interaction Dominates over vdW Interaction in All Models
The separate contributions of these two interactions are shown in Figures S3 and S4. Clearly, the electrostatic interaction is far superior than vdW in binding affinity of Tamiflu to NA. This observation was reported previously [18] for a few number of models, but the role of water has not been explored yet.
As follows from Tables 3ā6, the contribution of vdW interaction to the binding free energy is not sensitive to water models for all force fields except CHARMM where SPC makes markedly higher contribution compared to other models. The effect of environment on the electrostatic interaction is weak in OPLS (Table 4) leaving almost equal in 4 water models. For AMBER (Table 3) SPC/E gives the lowest estimation for , while the drastic water effect is observed in CHARMM (Table 5) and GROMOS (Table 6). In the latter case two models yield the difference in of about 30ākcal/mol, but SPC/E and TIP4P in CHARMM provide even larger difference of ā61.5kcal/mol.




3.2.3. Binding Free Energy Depends on Water Models
Apolar solvation energy and entropy contributions are not sensitive to force fields and water models (Tables 3ā6). is about 4.5ākcal/mol, while ā is 13ā15ākcal/mol for all models. The dependence of on water mainly comes from competition between the electrostatic energy and polar solvation energy . If they compensate each other as in the case of AMBER and GROMOS, then the absolute value of is small (Tables 3 and 6). For OPLS is far below the absolute value of leading to large . This result suggests that the charge parametrization of OPLS is not suitable for studying binding affinity of oseltamivir to NA and its mutations [18]. Since overestimation of was obtained by the MMPBSA method, it remains unclear whether other methods would change this conclusion. The similar noncompensation effect is observed in CHARMMSPC/E set (Table 5) where is also far below the experimental result (āā40.85ākcal/mol).
SPC/E generates the most negative values for both electrostatic and vdW interactions compared with other 3 models (Figure 4). Therefore, this model provides the highest binding affinity in all studied force fields (Tables 3ā6). This observation agrees with the previous study of Hu and Jiang [16] that the Coulomb interaction between water and lysozyme is more negative in SPC/E than in SPC and TIP3P since SPC/E has weaker selfdiffusion than others, but closer to the experiment. In GROMOS force fields SPC and SPC/E give and ākcal/mol, respectively (Table 6). Clearly, SPC result is closer to the experiment [23].
Averaging the binding free energy over water models, one has , , , and for AMBER, OPLS, CHARMM, and GROMOS, respectively. Thus, the strongest effect of water modeling is observed in CHARMM as the departure from the average value is about 8.6ākcal/mol, while the weakest influence is seen in OPLS with deviation of ā2.3ākcal/mol.
3.3. Recommendation for the Best Sets of Force Field and Water Model
To recommend the best combination one has to rely on the experimental results. The experiments of Collinās group [23] have shown that the binding free energy of Oseltaminir to A/H5N1 NA kcal/mol. Clearly, in AMBER TIP3P is the best fit to the experiments giving kcal/mol (Table 3). Thus, in accord with the GROMACS recommendation, AMBERTIP3P is the best choice for studying ligand binding affinity. The agreement with the GROMACSās suggestion has been also obtained for GROMOSSPC and OPLSTIP4P (Tables 6 and 4) having closer to the experiments than other sets. It should be noted that OPLSTIP4P is marginally better than OPLSTIP3P because their difference in is less than 1ākcal/mol. So OPLSTIP3P may not be a bad choice for estimation of binding affinity.
For CHARMM the closest to experiment result (kcal/mol) falls into SPC model (Table 5). TIP3P is ranked second having kcal/mol which is far from the experimental estimate. Thus, based on the results obtained for WT NA, one may recommend to use CHARMMSPC instead of CHARMMTIP3P suggested by the GROMACS. Since this conclusion is obtained for one system, to ascertain that SPC is the best choice for CHARMM we have computed for three more systems including mutants Y252H, N294S, and H274Y which have been studied experimentally [23]. We summarize the main results in Table 7 providing details of calculations for different water models in Supporting Information (SI) (Figure S5 and Tables S1āS4). The experiments show the ranking for binding affinity as Y252H WT N294S H274Y. This ranking is correctly captured by TIP3P and TIP4P but not SPC as well as SPC/E (Table 7). Comparing absolute values of with the experiments one can see that SPC and SPC/E are the best for WT and Y252H, while TIP3P is the best for both N294S and H274Y. Taken together, in accord with GROMACSās suggestion, TIP3P is most suitable for CHARMM.

3.4. Hydrogen Bond Network at the Binding Site
From previous MMPBSA results, the hydrogen bonding, which mostly contributes to the electrostatic energy, plays the key role in the interaction between Tamiflu and A/H5N1 NA [18, 48]. However, the role of water has not been explored yet. Figure S6 shows the time dependence of HBs obtained by different force fields and water models. The HB number not only levels significantly among force fields but also depends on aqueous environments.
3.4.1. Amber Force Field
Typical HB networks of four sets with AMBER are shown in Figure 5, where one has 7, 7, 6, and 7 HBs for SPC, SPC/E, TIP3P, and TIP4P, respectively. For all water models oseltamivir has the strong hydrogen bonding with residues E119, D151, R292, and R371 (lower panel of Figure 5). Within SPC/E the Hbonding with R152 is weaker than other models which have the population exceeding 80%. The strong interaction with E277 is observed only for this water model. Thus, in terms of individual contributions of ligand atoms SPC/E differs from other models. However, in equilibrium the average numbers of HBs are almost equal in all aqueous environments having , , , and 7.0 for SPC, SPC/E, TIP3P, and TIP4P, respectively (Figure S6).
3.4.2. OPLS Force Filed
For OPLS four aqueous models show nearly the same HB networks (Figure S6) (in equilibrium , , , and 7.2 for SPC, SPC/E, TIP3P, and TIP4P, resp.). This is not surprising because they also have little effect on the binding free energy as discussed in the previous section (Table 4). HB patterns are quite similar among various water models in AMBER and OPLS force fields (Figure 5 and Figure S7) implying that the geometry of ligand and area around the binding pocket does not depend much on water models. A slight difference is in population at residues R152 and E277 for two force fields.
3.4.3. CHARMM Force Filed
The situation becomes very different in the case of CHARMM where water has the strong effect on the HB network (Figure S8). The average number of HB in equilibrium varies between 5.5 for SPC and 4.4 for TIP4P (Figure S6). Contrary to AMBER and OPLS, only R371 remains the key residue for 4 aqueous models having the population more than 75%. SPC gives also strong Hbonding with E119 and D151 (population %), while D151, R152, and R292 are Hbonded with Tamiflu for the most simulation time with SPC/E (Figure S8). The ligand forms HB with R118 and E277 if one uses TIP3P but not other models (Figure S8). TIP4P shows the modest HB population with E119, R152, and Y347, while together with R371 residue R292 is conserved in this water model. The diversity of HB networks in CHARMM presumably causes strong variation of the binding free energy among 4 aqueous models (Table 5).
3.4.4. GROMOS Force Filed
As follows from Figure S6, due to the unitedatom approximation used for this force field the number of HBs is much lower ( and 1.6 for SPC and SPC/E) than other force fields. Consequently, HB networks are very poor (Figure S9). Hbonding in SPC/E is stronger than SPC leading to its higher binding affinity (Table 6). Residue R152 has the substantial population in this water model. For SPC Hbonding is weak for all residues from the binding pocket.
3.5. Effect of Hydration on Binding Affinity
3.5.1. AMBER Force Field
Volume of the binding pocket, estimated by approximate polyhedron (Materials and Methods), fluctuates during the course of MD simulation (Figure 6) depending on types of hydration. It gets saturated in equilibrium and the average volume is shown in Table S5. The largest volume is obtained for SPC/E, while TIP4P provides the smallest volume. The time fluctuation of the number of water molecules inside the pocket (Figure 7) indicates that water molecules keep going out and coming back (see Movie 1). The weak dependence on water models is observed for AMBER force field because in equilibrium the binding pocket contains 42ā46 water molecules (Table S5). SPC/E widens the binding site volume to a greater extent compared to other models providing the largest number of water molecules.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
Figure S10 shows the time dependence of the water density in the binding space for all cases. As expected, (Table S5) is lower than the standard density of 1ākg/L of water surrounding protein. It is well known that water weakens Hbonding leading to lower binding affinity than in vacuum. If this is true, then SPC/E model, for example, would provide the lowest binding affinity having the highest . However, this is not the case as this model provides the highest binding affinity (Table 4). In general, one has the strong correlation (correlation level ) between and (Figure 8) that the higher water density is, the higher binding affinity is. Since this correlation is at odds with the role of water in weakening Hbonding, one expects that HBs alone do not govern ligand binding affinity.
(a)
(b)
(c)
(d)
Using parameters of water models (Table 1), one can show that is not correlated with either the dielectric constant or dipole moment. Thus, one can not work out a unique factor that controls the binding affinity of Tamiflu to NA in AMBER. This is also true for other force fields.
3.5.2. OPLS Force Field
As in the AMBER case, binding site volume (Figure 6), number of water molecules (Figure 7), and water density (Figure S10) do not show much variations among water models. In equilibrium the volume fluctuates around 5000 for all water models. SPC shows the highest water density, while the lowest value of is given by TIP4P (Table S6). The latter model also has the smallest binding site. SPC/E and TIP3P have the same water density (Table S6) but different binding free energies (Table 4). For four water models there is a modest correlation () between and (Figure 8). Again, similar to the AMBER case, this correlation cannot explain the binding affinity through the influence of water on HB network.
3.5.3. CHARMM Force Field
The situation becomes entirely different in the case of CHARMM where the binding pocket volume (Figure 6) and number of water molecules inside it (Figure 7) are substantially higher than other force fields. TIP4P, which has the lowest volume in AMBER and OPLS, widens the pocket to the largest extent in CHARMM (Tables S5āS7). The number of water molecules in this model is about twofold larger than SPC/E. The packing of TIP4P water inside the binding site is also much denser (Figure S10 and Table S7) than other models having kg/L. Nevertheless, the corresponding binding free energy remains higher than SPC/E (Table 5). Overall, there is no correlation between and (Figure 8) in CHARMM.
3.5.4. GROMOS Force Field
In GROMOS the binding pocket volume (Figure 6), the number of water molecules (Figure 7), and water density (Figure S10) are lower than other force fields (Table S8). This presumably comes from unitedatom approximation. There is a pronounced difference in the binding free energies obtained by SPC and SPC/E due to different water densities. Thus, similar to AMBER and OPLS, decreases with .
4. Conclusions
Our previous study [18] on binding affinity of Tamiflu to variants of influenza A/H5N1 neuraminidase has revealed that AMBER99SB is the best among popular force fields as it reproduces with highest correlation and closest range with the experiments. In this paper, water models are tested with four force fields to examine their effects on the Tamiflu binding affinity toward the WT of A/H5N1 NA. As a result, TIP3P rather than SPC, SPC/E, and TIP4P goes along with AMBER99SB better than any combination of water models and other force fields.
Within one force field, in agreement with the GROMACS recommendation, combinations AMBER99SBTIP3P, OPLSTIP4P, CHARMMTIP3P, and GROMOSSPC are more suitable for simulation of ligand binding. Although this result has been obtained for NATamiflu complex and its validity for other systems is a subject for further investigation, the choice of these combinations is recommended. One has to bear in mind that the above combinations may not work in some cases. For example, CHARMMSPC is a better choice for WT NA than CHARMMTIP3P. In the case of Y252H, SCP/E works better than other water models if one utilizes CHARMM (Table 7). We have demonstrated that within one force field the binding free energy greatly varies for different combination of force fields and water models. For WT NA SPC/E always provides the lowest binding free energy among all water models regardless of force fields.
Contrary to the remaining 3 force fields, estimated by OPLSAA/L force field does not vary much among water models (Table 5). Due to strong electrostatic interaction, their values are too low compared with other force fields and experiments. It remains unclear if this is an artifact of MMPBSA approximation. Apolar solvation and entropy contributions are not affected by either force fields or water models (Tables 3ā6), but other terms are sensitive to them. The HB network between Tamiflu and NA changes little upon water models in OPLS and AMBER force fields, while it does strongly in GROMOS. CHARMM is a medium case. The pronounced influence of aqueous models on water density inside binding pocket has been observed in CHARMM force field.
Conflict of Interests
The authors declare that there is not conflict of interests in writing this paper.
Acknowledgments
The useful discussion with P. L. Chau is highly appreciated. The work was supported by Narodowe Centrum Nauki in Poland (Grant no. 2011/01/B/NZ1/01622) and Department of Science and Technology at Ho Chi Minh city, Vietnam.
Supplementary Materials
Additional figures including the initial structure for MD simulation, the binding site, the time dependence of the number of hydrogen bonds between Oseltamivir and receptor, typical snapshots of hydrogen bond networks, and time dependence of water density inside the binding pocket are presented. The tables on water density inside the binding site are provided.
References
 P. A. Kollman, I. Massova, C. Reyes et al., āCalculating structures and free energies of complex molecules: combining molecular mechanics and continuum models,ā Accounts of Chemical Research, vol. 33, no. 12, pp. 889ā897, 2000. View at: Publisher Site  Google Scholar
 J. G. Kirkwood, āStatistical mechanics of fluid mixtures,ā The Journal of Chemical Physics, vol. 3, no. 5, pp. 300ā313, 1935. View at: Google Scholar
 J. Aqvist, C. Medina, and J.E. Samuelsson, āA new method for predicting binding affinity in computeraided drug design,ā Protein Engineering, vol. 7, no. 3, pp. 385ā391, 1994. View at: Google Scholar
 R. W. Zwanzig, āHightemperature equation of state by a perturbation method. I. Nonpolar gases,ā The Journal of Chemical Physics, vol. 22, no. 12, pp. 1420ā1426, 1954. View at: Publisher Site  Google Scholar
 M. Karplus and J. A. McCammon, āMolecular dynamics simulations of biomolecules,ā Nature Structural Biology, vol. 9, no. 9, pp. 646ā652, 2002. View at: Publisher Site  Google Scholar
 F. S. Lee, Z.T. Chu, M. B. Bolger, and A. Warshel, āCalculations of antibodyantigen interactions: microscopic and semimicroscopic evaluation of the free energies of binding of phosphorylcholine analogs to McPC603,ā Protein Engineering, vol. 5, no. 3, pp. 215ā228, 1992. View at: Google Scholar
 U. Bren, J. Lah, M. Bren, V. Martnek, and J. Florin, āDNA duplex stability: the role of preorganized electrostatics,ā Journal of Physical Chemistry B, vol. 114, no. 8, pp. 2876ā2885, 2010. View at: Publisher Site  Google Scholar
 H. Grubmüller, B. Heymann, and P. Tavan, āLigand binding: molecular mechanics calculation of the streptavidinbiotin rupture force,ā Science, vol. 271, no. 5251, pp. 997ā999, 1996. View at: Google Scholar
 B. K. Mai, M. H. Viet, and M. S. Li J, āTop leads for swine influenza A/H1N1 virus revealed by steered molecular dynamics approach,ā Journal of Chemical Information and Modeling, vol. 50, no. 12, pp. 2236ā2247, 2010. View at: Publisher Site  Google Scholar
 J. H. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, and J. Hermans, āInteraction models for water in relation to protein hydration,ā in Intermolecular Forces, B. Pullmann, Ed., pp. 331ā342, 1981. View at: Google Scholar
 H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, āThe missing term in effective pair potentials,ā Journal of Physical Chemistry, vol. 91, no. 24, pp. 6269ā6271, 1987. View at: Google Scholar
 W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, āComparison of simple potential functions for simulating liquid water,ā The Journal of Chemical Physics, vol. 79, no. 2, pp. 926ā935, 1983. View at: Google Scholar
 W. L. Jorgensen and J. D. Madura, āTemperature and size dependence for Monte Carlo simulations of TIP4P water,ā Molecular Physics, vol. 56, no. 6, pp. 1381ā1392, 1985. View at: Publisher Site  Google Scholar
 B. Hess and N. F. A. van der Vegt, āHydration thermodynamic properties of amino acid analogues: a systematic comparison of biomolecular force fields and water models,ā Journal of Physical Chemistry B, vol. 110, no. 35, pp. 17616ā17626, 2006. View at: Publisher Site  Google Scholar
 P. Florova, P. Sklenovsky, P. Banas, and M. Otyepka, āExplicit water models affect the specific solvation and dynamics of unfolded peptides while the conformational behavior and flexibility of folded peptides remain intact,ā Journal of Chemical Theory and Computation, vol. 6, no. 11, pp. 3569ā3579, 2010. View at: Publisher Site  Google Scholar
 Z. Hu and J. Jiang, āAssessment of biomolecular force fields for molecular dynamics simulations in a protein crystal,ā Journal of Computational Chemistry, vol. 31, no. 2, pp. 371ā380, 2009. View at: Publisher Site  Google Scholar
 M. Almlof, B. O. Brandsdal, and J. Aqvist, āBinding affinity prediction with different force fields: examination of the linear interaction energy method,ā Journal of Computational Chemistry, vol. 25, no. 10, pp. 1242ā1254, 2004. View at: Publisher Site  Google Scholar
 T. T. Nguyen, B. K. Mai, and M. S. Li, āStudy of tamiflu sensitivity to variants of A/H5N1 virus using different force fields,ā Journal of Chemical Information and Modeling, vol. 51, no. 9, pp. 2266ā2276, 2011. View at: Publisher Site  Google Scholar
 P. Cozzini, M. Fornabaio, A. Marabotti, D. J. Abraham, G. E. Kellogs, and A. Mozzarelli, āFree energy of ligand binding to protein: evaluation of the contribution of water molecules by computational methods,ā Current Medicinal Chemistry, vol. 11, no. 23, pp. 3093ā3118, 2004. View at: Google Scholar
 J. Michel, J. TiradoRives, and W. L. Jorgensen, āPrediction of the water content in protein binding sites,ā Journal of Physical Chemistry B, vol. 113, no. 40, pp. 13337ā13346, 2009. View at: Publisher Site  Google Scholar
 L. Wang, B. J. Berne, and R. A. Friesner, āLigand binding to proteinbinding pockets with wet and dry regions,ā Proceedings of the National Academy of Sciences of the United States of America, vol. 108, no. 4, pp. 1326ā1330, 2011. View at: Publisher Site  Google Scholar
 M. Chaplin, āWater structure and science,ā 2000, http://www.lsbu.ac.uk/water/. View at: Google Scholar
 P. J. Collins, L. F. Haire, Y. P. Lin et al., āCrystal structures of oseltamivirresistant influenza virus neuraminidase mutants,ā Nature, vol. 453, no. 7199, pp. 1258ā1261, 2008. View at: Publisher Site  Google Scholar
 V. Hornak, R. Abel, A. Okur, B. Strockbine, A. Roitberg, and C. Simmerling, āComparison of multiple amber force fields and development of improved protein backbone parameters,ā Proteins, vol. 65, no. 3, pp. 712ā725, 2006. View at: Publisher Site  Google Scholar
 B. R. Brooks, C. L. Brooks, A. D. Mackerell et al., āCHARMM: the biomolecular simulation program,ā Journal of Computational Chemistry, vol. 30, no. 10, pp. 1545ā1614, 2009. View at: Publisher Site  Google Scholar
 W. L. Jorgensen and J. TiradoRives, āThe OPLS potential functions for proteins. Energy minimizations for crystals of cyclic peptides and crambin,ā Journal of the American Chemical Society, vol. 110, no. 6, pp. 1657ā1666, 1988. View at: Google Scholar
 W. F. van Gunsteren, S. R. Billeter, A. A. Eising et al., Biomolecular Simulation: The GROMOS96 Manual and User Guide, Vdf Hochschulverlag AG an der ETH Zurich, Zurich, Switzerland, 1996.
 R. G. Webster and E. A. Govorkova, āH5N1 influenza—continuing evolution and spread,ā The New England Journal of Medicine, vol. 355, no. 21, pp. 2174ā2177, 2006. View at: Publisher Site  Google Scholar
 The Writing Committee of the World Health Organization (WHO) Consultation on Human Influenza A/H5, āAvian influenza A (H5N1) infection in humans,ā The New England Journal of Medicine, vol. 353, no. 13, pp. 1374ā1385, 2005. View at: Publisher Site  Google Scholar
 Schrödinger, PyMOL: The PyMOL Molecular Graphics System, Version 1.3, Schrödinger, 2010.
 D. M. F. van Aalten, R. Bywater, J. B. C. Findlay, M. Hendlich, R. W. W. Hooft, and G. Vriend, āPRODRG, a program for generating molecular topologies and unique molecular descriptors from coordinates of small molecules,ā Journal of ComputerAided Molecular Design, vol. 10, no. 3, pp. 255ā262, 1996. View at: Publisher Site  Google Scholar
 M. J. Frisch, G. W. Trucks, H. B. Schlegel et al., Gaussian 03, Revision C.02, Gaussian, Wallingford, Conn, USA, 2004.
 A. W. S. D. Silva, W. F. Vranken, and E. D. Laue, āACPYPE—AnteChamber PYthon Parser interfacE,ā submitted. View at: Google Scholar
 A. S. T. R. Andre, A. C. H. Bruno, and B. A. Ricardo, āMKTOP: a program for automatic construction of molecular topologies,ā Journal of the Brazilian Chemical Society, vol. 19, no. 7, pp. 1433ā1435, 2008. View at: Publisher Site  Google Scholar
 V. Zoete, M. A. Cuendet, A. Grosdidier, and O. Michielin, āSwissParam: a fast force field generation tool for small organic molecules,ā Journal of Computational Chemistry, vol. 32, no. 11, pp. 2359ā2368, 2011. View at: Publisher Site  Google Scholar
 B. Guillot, āA reappraisal of what we have learnt during three decades of computer simulations on water,ā Journal of Molecular Liquids, vol. 101, no. 13, pp. 219ā260, 2002. View at: Publisher Site  Google Scholar
 J. D. Bernal and R. H. Fowler, āA theory of water and ionic solution, with particular reference to hydrogen and hydroxyl ions,ā The Journal of Chemical Physics, vol. 1, no. 8, pp. 515ā548, 1933. View at: Google Scholar
 T. Darden, D. York, and L. Pedersen, āParticle mesh Ewald: an N·log(N) method for Ewald sums in large systems,ā The Journal of Chemical Physics, vol. 98, no. 12, pp. 10089ā10092, 1993. View at: Google Scholar
 R. W. Hockney, S. P. Goel, and J. W. Eastwood, āQuiet highresolution computer models of a plasma,ā Journal of Computational Physics, vol. 14, no. 2, pp. 148ā158, 1974. View at: Google Scholar
 H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. Dinola, and J. R. Haak, āMolecular dynamics with coupling to an external bath,ā The Journal of Chemical Physics, vol. 81, no. 8, pp. 3684ā3690, 1984. View at: Google Scholar
 M. Parrinello and A. Rahman, āPolymorphic transitions in single crystals: a new molecular dynamics method,ā Journal of Applied Physics, vol. 52, no. 12, pp. 7182ā7190, 1981. View at: Publisher Site  Google Scholar
 B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, āGRGMACS 4: algorithms for highly efficient, loadbalanced, and scalable molecular simulation,ā Journal of Chemical Theory and Computation, vol. 4, no. 3, pp. 435ā447, 2008. View at: Publisher Site  Google Scholar
 L. S. Cheng, R. E. Amaro, D. Xu, W. W. Li, P. W. Arzberger, and J. A. McCammon, āEnsemblebased virtual screening reveals potential novel antiviral compounds for avian influenza neuraminidase,ā Journal of Medicinal Chemistry, vol. 51, no. 13, pp. 3878ā3894, 2008. View at: Publisher Site  Google Scholar
 C. B. Barber, D. P. Dobkin, and H. Huhdanpaa, āThe quickhull algorithm for convex hulls,ā ACM Transactions on Mathematical Software, vol. 22, no. 4, pp. 469ā483, 1996. View at: Google Scholar
 K. L. Clarkson, K. Menlhorn, and R. Seidel, āFour results on randomized incremental constructions,ā Computational Geometry, vol. 3, no. 4, pp. 185ā212, 1993. View at: Publisher Site  Google Scholar
 P. L. Chau, āWater movement during Ligand Unbinding from receptor site,ā Biophysical Journal, vol. 87, no. 1, pp. 121ā128, 2004. View at: Publisher Site  Google Scholar
 MATLAB Version 7.0.1 (R2007a), The MathWorks Inc., Natick, Mass, USA, 2007.
 O. Aruksakunwong, M. Malaisree, P. Decha et al., āOn the lower susceptibility of oseltamivir to influenza neuraminidase subtype N1 than those in N2 and N9,ā Biophysical Journal, vol. 92, no. 3, pp. 798ā807, 2007. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2014 Trang Truc Nguyen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.