Article pubs.acs.org/JPCB

Hydration Behavior at the Ice-Binding Surface of the Tenebrio molitor Antifreeze Protein Uday Sankar Midya and Sanjoy Bandyopadhyay* Molecular Modeling Laboratory, Department of Chemistry, Indian Institute of Technology, Kharagpur 721302, India ABSTRACT: Molecular dynamics (MD) simulations have been carried out at two different temperatures (300 and 220 K) to study the conformational rigidity of the hyperactive Tenebrio molitor antifreeze protein (TmAFP) in aqueous medium and the structural arrangements of water molecules hydrating its surface. It is found that irrespective of the temperature the ice-binding surface (IBS) of the protein is relatively more rigid than its nonice-binding surface (NIBS). The presence of a set of regularly arranged internally bound water molecules is found to play an important role in maintaining the flat rigid nature of the IBS. Importantly, the calculations reveal that the strategically located hydroxyl oxygens of the threonine (Thr) residues in the IBS influence the arrangements of five sets of ordered waters around it on two parallel planes that closely resemble the basal plane of ice. As a result, these waters can register well with the ice basal plane, thereby allowing the IBS to preferentially bind at the ice interface and inhibit its growth. This provides a possible molecular reason behind the ice-binding activity of TmAFP at the basal plane of ice.

1. INTRODUCTION The organisms living in subzero temperatures usually contain a particular class of proteins, known as antifreeze proteins (AFPs).1,2 The presence of AFPs allows such organisms to survive by avoiding cell damage due to freezing of the body fluid. AFPs can noncolligatively lower the freezing point of a solution without actually changing its melting point. The difference between the two temperatures termed as “thermal hysteresis” depends on the concentration of the AFP.3 In addition to lowering the freezing point, AFPs often induce faceted ice crystal morphology with unchanged shape and size of the crystal within the hysteresis range.4 The structural characteristics and activities of a large variety of AFPs found in insects,5−12 plants,13,14 polar fishes,15,16 bacteria,17,18 etc., have been extensively studied. These proteins are broadly classified into two groups based on their activities, namely, moderately active and hyperactive AFPs.19 Generally, fish AFPs are categorized as moderately active, while the AFPs in other species belong to the hyperactive class. Due to structural diversity, there is no single consensus theory that can explain AFP activity. In an early study, Raymond and DeVries20 proposed that despite such diverse structural features AFPs in general function through a common adsorption−inhibition mechanism. According to this mechanism, it is believed that AFPs first adsorb irreversibly at the growing ice surfaces and then allow the ice to grow in between the adsorbed proteins. This leads to an increase of curvature on the ice surface, thereby reducing the local freezing temperature (Kelvin effect21). Recently, Braslavsky and co-workers22,23 showed from fluoroscence microscopy and microfluidic experiments that moderately active type III and hyperactive AFPs from insect © 2014 American Chemical Society

Tenebrio molitor indeed bind irreversibly to ice surfaces. However, such an irreversible adsorption−inhibition mechanism could not explain the binding of moderately active AFP from fish winter flounder to ice surfaces as obtained from NMR studies.24 In fact, this study provides evidence of a reversible binding mechanism. Despite significant efforts, a molecular level understanding of protein−ice interaction and the relative molecular arrangements between them that govern the mechanism of this recognition process are not clear. As a result, this is considered to be one of the most complex recognition problems in biology.25 In an important early work, Knight et al.26 carried out ice etching experiments to study the binding of the AFP present in fish winter flounder to ice. The α-helical AFP present in this species has a repetitive sequence of 11 amino acid units ended with threonine (Thr) along the axis of the helix. They found that the spacing between these Thr residues (16.5 Å) matches with the interoxygen distance (16.7 Å) on the {202̅1} pyramidal planes of ice along the ⟨011̅2⟩ directions. On the basis of this, it was proposed that the binding between the AFP and ice occurs through hydrogen bonds between the protein side chains and ice. Wen and Laursen27 showed that unidirectional binding of winter flounder AFP to ice along ⟨1̅102⟩ on {202̅1} pyramidal planes of ice allows the formation of protein−water hydrogen bonds between Thr and other polar residues present in the repeat segments of the protein and water. Molecular modeling and simulation studies28,29 on this α -helical AFP have also Received: December 21, 2013 Revised: March 6, 2014 Published: April 11, 2014 4743

dx.doi.org/10.1021/jp412528b | J. Phys. Chem. B 2014, 118, 4743−4752

The Journal of Physical Chemistry B

Article

Although the structural aspects of AFPs and their hydration properties have received immense attention in the recent past, several aspects of the problem are still elusive. In particular, interactions between AFPs and water and their preferential binding at the ice−water interface have not been understood properly. As a result, a microscopic knowledge of the antifreeze activities of these proteins is yet to emerge. In this work, we report results obtained from atomistic MD simulations of aqueous solutions of the Tenebrio molitor antifreeze protein (TmAFP) at two different temperatures (300 and 220 K). Efforts have been made to explore the differential conformational rigidity of the ice-binding and nonice-binding surfaces (IBS and NIBS) of the protein and the structure and ordering patterns of water molecules hydrating the two interfaces. The rest of the article has been organized as follows. In Section 2, we provide a brief description of the system setups and the simulation methods employed. The results obtained from our analyses are presented and discussed in the following section (Section 3). Finally, the important findings from the study and the conclusions reached therefrom are highlighted in Section 4.

supported the hydrogen bond hypothesis for its binding along the ⟨1̅102⟩ directions. On the other hand, Knight and coworkers30 showed that on mutating the Thr residues on the icebinding surface (IBS) of the winter flounder AFP by serines (Ser) that preserve the γ-hydroxyl groups the protein looses its residual activity. However, the protein was found to retain some degree of activity on similar mutations by valine (Var) residues that preserve the γ-methyl groups. In another study with a type III AFP, mutations of the alanines (Ala) on the IBS by Thr residues were found to reduce the activity of the protein.31 Similar effects of mutations on the activity of AFPs have also been shown recently by Braslavsky and co-workers.32 Thus, in contrast to earlier studies, these experiments showed that formation of hydrogen bonds may not be a necessary requirement for antifreeze activity. Additionally, it is also not clear why the IBS residues of an AFP would prefer to form hydrogen bonds with ice water instead of liquid water which is present in large excess. MD simulations carried out by Madura and co-workers33 in water and at an ice−water interface also did not reveal any noticeable increase in hydrogen bonds involving the protein at the interface. Interestingly, using a random network model for water,34 Gallagher and Sharp35 showed that the polar groups present at the IBS of type III antifreeze proteins exhibit nonpolar-like hydration characteristics. It is further proposed that the flat IBS of the AFPs prefers ice binding through surface complementarity via van der Waals interaction between the protein and ice.36 However, with the ice−water interface being diffused in nature,37 such complementary interaction may not be that strong at the interface. In contrast, Nutt and Smith38 recently showed from their simulation study that the hyperactive AFP found in spruce budworm first prefers to order its solvation waters into quasiice-like structure during the recognition process. This is followed by merging of the protein solvation shell with the ice−water interface and subsequent freezing of the merged zone and binding of the protein on the ice surface. Simulated water density distribution revealed that the hydration pattern on IBS resembles the basal and primary prism planes of ice, which agrees well with experimental findings.19 Simulation studies by Smolin and Daggett39 also showed similar ice-like water structure on the IBS of type III AFPs. Recent solid-state NMR experiments40,41 support the simulation data. These experiments showed that on freezing the solution containing type III AFP the interfacial waters around the IBS of the protein turn to ice, whereas those around the NIBS (nonicebinding surface) remain in liquid form. Most of the hyperactive AFPs have a regularly arranged two-dimensional array of Thr residues on the flat IBS.6,8,12,18 In addition, some of the AFPs also contain water molecules located in between the ordered Thr residues in their cystalline forms.6,12,18 The separations between these Thr residues and that between the crystal waters are found to match closely with the separations between the oxygens at the ice surface. These regularly spaced waters are believed to play important roles in maintaining the rigidity of the protein structure and its activity in solution. Interestingly, X-ray crystallography and fluorescent imaging studies12,18,42 showed nonpreferential adsorption of these AFPs on both the prism and the basal planes of ice. On the other hand, icegrowth studies19,42,43 below the nonequilibrium freezing point showed that hyperactive AFPs allow the growth of ice on the prism plane, thereby indicating strong preference of these AFPs to bind with the basal plane at the ice−water interface.

2. SIMULATION DETAILS We have carried out two simulations, designated as S1 and S2, with the protein TmAFP in aqueous medium at temperatures 300 and 220 K, respectively. The initial coordinates of TmAFP were taken from its crystal structure as reported by Davies and co-workers6 and available in the Protein Data Bank (PDB code: 1EZG). TmAFP is a pseudorectangular-shaped right-handed parallel β-helix protein with 84 amino acids. It contains seven loops with a series of 12 residue repeats with sequence TCTxSxxCxxAx. The IBS of the protein contains the TCT (threonine-cysteine-threonine) motifs of the loops (except the N-terminal loop) and is on one side forming a flat β-sheet. We denote the remaining part of the protein as the NIBS. Among the crystal waters, there are five that are internally bound and located regularly between the six C-terminal loops close to the conserved alanine (Ala) residues and six others that are externally bound to the IBS and located regularly between the Thr residue side chains of the β-strands. The protein crystal structure was reported with an iodinated tyrosine (Tyr) residue at the 71st position and two end residues missing at the Cterminus of the sequence. We removed the iodine and added the missing residues, and the initial structure of the protein is shown in Figure 1. The simulations were performed using the NAMD package.44 The protein with its crystal waters was first minimized using the conjugate gradient energy minimization method as implemented in the NAMD code.44 The minimized protein structure was used to set up the two simulation systems (S1 and S2) in aqueous medium at two temperatures. The systems were prepared by inserting the protein with its crystal waters in two cubic cells containing equilibrated water molecules at respective temperatures. To avoid any unfavorable contact between the protein and water molecules, the insertion process in each case was carried out by carefully removing those water molecules that were found within 2 Å from the protein atoms. Two sodium ions (Na+) were then added in each system to neutralize the overall charge. Finally, the systems S1 and S2 contained the 84-residue protein dissolved in cubic cells with initial edge lengths of 65.49 and 65.60 Å containing 8895 and 8901 water molecules with 2 Na+ counterions, respectively. Once again, the two systems in the presence of solvent were minimized following the same protocol as mentioned above. 4744

dx.doi.org/10.1021/jp412528b | J. Phys. Chem. B 2014, 118, 4743−4752

The Journal of Physical Chemistry B

Article

equilibrated for 10 ns each under NPT ensemble conditions. The temperatures of the systems were controlled by the Langevin dynamics method with a friction constant of 1 ps−1, while the pressures were controlled by the Nosé−Hoover Langevin piston method.45 In both cases, the volumes of the simulation cells were allowed to fluctuate isotropically during this time period. At this stage, the cell volumes attained steady values with edge lengths of 65.37 and 65.35 Å for systems S1 and S2, respectively. We then fixed the individual cell volumes, and the simulation conditions were changed from that of NPT to NVT ensemble (constant volume and temperature). The equilibration of each of the two systems was then continued under NVT conditions for another 10 ns. This was followed by long NVT production runs of 130 ns duration for each of the two systems. The simulations were carried out with an integration time step of 1 fs, while the trajectories were stored with a time resolution of 500 fs for subsequent analyses. All bonds involving the hydrogen atoms were constrained by the SHAKE algorithm.46 The periodic boundary conditions and the minimum image convention47 were employed to calculate the short-range Lennard-Jones interactions with a spherical cutoff distance of 12 Å and a switch distance of 10 Å. The longrange electrostatic interactions were calculated using the particle-mesh Ewald (PME) method. 48 The all-atom CHARMM22 force field and potential parameters for proteins with CMAP corrections49,50 and the TIP4P model51 for water were employed in the calculations. It has been proved recently that though the CHARMM force fields were originally developed using the TIP3P model51 the TIP4P model that provides better structural properties of water in the bulk state

Figure 1. Crystal structure of the TmAFP protein6 (drawn as cartoons) with two missing residues added at the C-terminus. Two different orientations of the protein are shown, where the β-strands forming the IBS are drawn in blue, while the NIBS comprised of the remaining parts of the protein are drawn in green. Five internally and six externally bound waters at the IBS are shown as red spheres, while the side chains of the two sets of threonine (T) residues in the IBS are drawn as sticks. The primary amino acid sequence of the protein (in one-letter code) and the N- and C-terminal residues (Q(1) and H(84)) are marked for convenience.

The temperatures of the two systems were then gradually increased to their target values within short MD runs of 100 ps each under isothermal−isobaric ensemble (NPT) conditions at a constant pressure of 1 atm. The systems were then

Figure 2. Snapshots of a few representative configurations of the TmAFP protein as obtained at 10, 50, 80, 110, and 150 ns from simulations S1 at 300 K (a−e) and S2 at 220 K (f−j). The images are shown in two different orientations, and the protein coloring scheme is the same as that in Figure 1. 4745

dx.doi.org/10.1021/jp412528b | J. Phys. Chem. B 2014, 118, 4743−4752

The Journal of Physical Chemistry B

Article

could be used along with the CHARMM parameters without loss of generality in protein−water interactions.52

3. RESULTS AND DISCUSSION 3.1. Structural Features of TmAFP. It is necessary to explore the conformational flexibility of a protein to understand its interaction with the surrounding solvent at a microscopic level. For that, we first display snapshots of a few configurations of the protein as obtained from the two simulations in Figure 2. A comparison between these configurations with the crystal structure (Figure 1) shows close resemblance between the simulated solution structure of the protein and that of its crystalline form. This is particularly true for the IBS of the protein. A close examination of the trajectories reveals minor differences between the crystal structure and the simulated protein configurations around the terminal loop segments at the nonice-binding surface (NIBS) at room temperature. This indicates relatively higher flexibility of the protein at its two termini at 300 K. Note that the flexibility at the terminal segments of the AFP at room temperature is not surprising as the protein has evolved to exhibit its function under subzero temperatures. However, such behavior is consistent with NMR experiments, where it is found that the N-terminal loop of the protein in aqueous solution is not well-defined at room temperature.53 The conformational flexibilities of TmAFP at the two temperatures have been probed further in a more quantitative manner by measuring the root-mean-square deviations (RMSDs) of the simulated protein configurations with respect to the crystal structure. RMSD calculations can provide first hand information on local conformational jumps often exhibited by complex biomolecules. Here, the calculations are carried out by considering the non-hydrogen atoms for the whole protein (except the two terminal residues) as well as separately for only the residues forming the IBS and the NIBS of the protein. Time evolutions of the RMSD data are shown in Figure 3. Significant rigidity of the protein in aqueous environment with minor deviations from the crystalline form is evident from the figure. Expectedly, the deviations are relatively higher at 300 K (S1) as compared to that at 220 K (S2). Interestingly, at both temperatures the IBS is found to be relatively more rigid than the NIBS. This is an important observation as one would expect that for efficient ice-binding activity the IBS should avoid large fluctuations and maintain its flat rigid nature. The estimated average RMSD values as obtained from the two trajectories are listed in Table 1. Note that irrespective of the temperature the average RMSD of the IBS of the protein is more than two times smaller than that of the NIBS. In addition, the effect of temperature on the average RMSD value is found to be relatively less significant for the IBS, thereby indicating its inherent rigid nature. It is known that in addition to eight disulfide bonds the presence of five internally bound water molecules plays an important role in controlling the rigidity of TmAFP.6 Our calculations reveal that these five water molecules remain in their locations during the entire simulation periods in the two systems. To explore the importance of these internally bound waters, we have carried out an additional simulation of the protein without those waters for about 60 ns duration at 300 K following the same protocols as described before. Calculation of water distribution inside the protein during the time scale of the simulation reveals nondetectable population of waters around these specific internal locations within the protein. This shows

Figure 3. Time evolutions of the RMSDs based on the non-hydrogen atoms of the whole TmAFP protein and its two surfaces (IBS and NIBS) with respect to the crystal structure as obtained from simulations (a) S1 at 300 K and (b) S2 at 220 K. The corresponding data obtained from the simulation of the protein at 300 K without the internally bound water molecules are shown in the inset.

Table 1. Average RMSD Values (in Å) Based on the NonHydrogen Atoms of the Whole TmAFP Protein and That of Its Two Surfaces (IBS and NIBS) as Obtained From Simulations S1 at 300 K and S2 at 220 Ka simulation

whole

IBS

NIBS

S1 S2

1.14 (0.11) 0.90 (0.08)

0.60 (0.09) 0.48 (0.06)

1.24 (0.12) 0.98 (0.10)

a

The values in the parentheses are the corresponding standard deviations.

that the penetration of waters inside the protein is minimum. The time evolutions of the RMSDs of the protein without these specific bound waters with respect to the experimental structure are shown in the inset of Figure 3a. It can be seen that although the overall flexibility of the protein increases significantly in the absence of the internally bound waters the IBS still maintains its relatively higher rigidity as compared to the NIBS. The average RMSD values of the IBS and the NIBS of the protein in the absence of the bound waters are found to be 1.91 (±0.10) Å and 4.1 (±0.20) Å, respectively. Therefore, it is clear that these internally bound water molecules strategically located inside the β-helical core of the protein play an important role in reducing the protein flexibility, which in turn is important for its activity. In Figure 4 we show the variation of per residue solventaccessible surface area (SASA) of the protein as a function of time as obtained from simulations S1 and S2 at the two temperatures. Again, the calculations are carried out separately for the two surfaces (IBS and NIBS) and shown in the figure. Note that due to differential number of residues at the two surfaces the results are obtained by averaging over the number of residues in each case. To compute SASA, spheres are first drawn around the protein atoms with their radii increased by the solvent probe radius of 1.4 Å.54 The SASA values are then computed from the nonoverlapping areas of such larger spheres 4746

dx.doi.org/10.1021/jp412528b | J. Phys. Chem. B 2014, 118, 4743−4752

The Journal of Physical Chemistry B

Article

last 130 ns of equilibrated NVT trajectories as obtained from simulations S1 and S2. Structural arrangements of water molecules around the IBS and NIBS of the protein have been probed by calculating the pairwise correlation function, s(r), between waters and the Cα atoms of the two surfaces. The calculations are similar to that of the radial distribution function (rdf), except that the minimum distance condition is used to identify the tagged water with respect to a Cα atom of one of the two surfaces. In addition, the s(r) values are normalized by using the actual volume available to water around a surface (IBS or NIBS) within a shell at a distance r with thickness δr instead of by the entire volume of the corresponding spherical shell, as usually done in standard rdf calculations. This allows one to avoid possible misinterpretation of hydration patterns obtained from standard rdf calculations around heterogeneously rough protein surfaces. We have used a Monte Carlo integration method using random number insertions as proposed by Astley et al.58 to calculate the available volume of water around the IBS and NIBS of the protein. The results as obtained from our calculations at the two temperatures are shown in Figure 5(a,b). Differential

Figure 4. Time evolutions of per residue solvent accessible surface areas (SASAs) averaged over the whole TmAFP protein and its two surfaces (IBS and NIBS) as obtained from simulations (a) S1 at 300 K and (b) S2 at 220 K.

using spherical coordinates. Near-uniform variations of SASA with minor fluctuations in all cases can be seen from the figure. This is consistent with the overall rigid nature of the protein in aqueous medium as discussed before. Importantly, irrespective of the temperature the results demonstrate noticeably lower SASA for the IBS of the protein. The lower per residue SASA for the IBS suggests its relatively flat nature with reduced curvature and less exposed residues to solvent. The average SASA values as obtained from the two simulated trajectories are listed in Table 2. The data show that the change of temperature Table 2. Average per Residue SASA Values (in Å2) of the Whole TmAFP Protein and That of Its Two Surfaces (IBS and NIBS) as Obtained From Simulations S1 at 300 K and S2 at 220 Ka simulation

whole

IBS

NIBS

S1 S2

50.94 (0.92) 50.31 (0.69)

38.29 (1.41) 37.86 (1.18)

54.39 (1.17) 53.71 (0.76)

Figure 5. Surface pair correlation function, s(r), of water molecules as a function of distance from the Cα atoms of the residues of the two surfaces (IBS and NIBS) of the TmAFP protein as obtained from simulations (a) S1 at 300 K and (b) S2 at 220 K. The Voronoi volume distribution per water, P(V), averaged over the water molecules that are present in the first layers around the two surfaces of the protein as obtained from simulations (c) S1 at 300 K and (d) S2 at 220 K.

a

The values in the parentheses are the corresponding standard deviations.

has only a minor influence on the average SASA values for IBS and NIBS. This correlates well with the temperatureindependent structural rigidity of the protein, as already discussed. 3.2. Structure and Ordering of Surface Water. Protein structure is intrinsically heterogeneous at the microscopic level. Due to such heterogeneity, protein−water interactions can be nonuniform, thereby leading to differential hydration characteristics at the surface.55−57 Such nonuniformity in hydration patterns around the active and nonactive sites of a protein can often be correlated with its binding activity.55 In this section we explore such differential hydration behavior, if any, at the surface of TmAFP by comparing the distribution patterns and structural arrangements of water molecules at two of its surfaces (IBS and NIBS) at two different temperatures. The results presented in this section are carried out by averaging over the

solvent distributions around the two surfaces of the protein can be easily seen. We observe noticeable structuring of waters with several well-defined solvation shells and distance extending up to about 9 Å from the IBS. In contrast, no significant water structuring beyond the first hydration shell (within 5 Å) has been observed around the NIBS. Such nonuniform arrangements of water molecules near the two surfaces exist even at 300 K; however, the difference is more at the low temperature (220 K). This demonstrates that the flat IBS of the protein with reduced conformational oscillations has a strong influence on the nearby water molecules and tends to align them in an ordered manner over a longer distance from the surface. Similar 4747

dx.doi.org/10.1021/jp412528b | J. Phys. Chem. B 2014, 118, 4743−4752

The Journal of Physical Chemistry B

Article

long-range structural ordering of water has also been noticed from earlier MD simulations around the IBS of winter flounder AFP.59 Interestingly, it may be noted that such a hydration pattern around the rigid flat IBS of the AFPs is similar to extended layered solvent structures at solid−water interfaces.60−63 Our calculations reveal that layering of waters around the IBS is an inherent feature of TmAFP as evident from the results obtained at 300 K. This is an important observation as such long-range ordered arrangement of water molecules is expected to act as a precursor for the ice-binding activity of the protein. It is worth mentioning here that Havenith and coworkers64,65 have shown from recent terahertz spectroscopy measurements that the organized water layer around the IBS of an AFP can extend even further (up to 20 Å or so) than that observed in the present study. We have also calculated the Voronoi volume47,66 distributions per water molecule by averaging over those waters that are present in the first layers around the two surfaces of the protein at two temperatures, as shown in Figure 5(c,d). Voronoi volume is defined as the volume of a polyhedron formed by the closest points with respect to an atom (water oxygen in this case) present at the center of the polyhedron. In the calculation, a surface water molecule is selected as that which has one of the faces of the polyhedron in common with a non-hydrogen atom of the protein. It is apparent from the results that a water molecule present at the IBS occupies to some extent larger volume than a corresponding water present at the NIBS. This is true at both the temperatures and indicates lower local density of water at the IBS. The calculated average Voronoi volumes per water at the IBS and NIBS for system S1 (300 K) are found to be 30.64 and 30.06 Å3, while those for system S2 (220 K) are 29.78 and 29.29 Å3, respectively. Note that though the differences are marginal they are meaningful in the context of a small density difference between ice and water. Relatively lower water density at the IBS indicates the preference of these water molecules to reorient themselves in an ice-like fashion around it. So far, we have shown how the structural arrangements of water molecules are influenced in a nonuniform manner at the two surfaces of the protein. It is demonstrated that the water molecules near the IBS exhibit layered structuring with apparent lower density than that around the NIBS. We now explore how the regular tetrahedral ordering of water in the bulk state is modified around the two surfaces at two temperatures. For that, we have calculated the tetrahedral order parameter, qt, defined as67−69 qt = 1 −

3 8

3

Figure 6. Distributions of the tetrahedral order parameter, P(qt), for the water molecules that are present within 5 Å from the nonhydrogen atoms of the two surfaces (IBS and NIBS) of the TmAFP protein as obtained from simulations (a) S1 at 300 K and (b) S2 at 220 K. The corresponding distributions for pure TIP4P water at identical temperatures are included for comparison.

higher qt values (Figure 6b). Importantly, the calculation reveals similar distribution patterns around the two surfaces of the protein that closely resemble that of pure bulk water at a particular temperature. However, the height of the distribution curve for the IBS water molecules is found to be slightly more at 220 K. This shows a small but noticeably increased fraction of tetrahedrally arranged waters near the IBS, which approaches that of water in the pure bulk state at lower temperature. The results presented above indicate low-density ice-like water arrangements near the IBS of the protein. To obtain a better understanding of the spatial arrangements of water molecules, we have computed their two-dimensional probability distributions around the protein. For such calculations, first the coordinates of all the atoms along the trajectories are transformed in such a manner so that the helical axis of the protein remains aligned along the z axis, and the IBS remains parallel to the xz plane. Next, to probe the influence of the central four β-strands that take part in ice-binding activity on water distribution, the water molecules that are within −12.5 Å ≤ z ≤ 12.5 Å are mapped on the plane perpendicular to the z axis (i.e., xy plane). The distributions at the two temperatures are calculated with a grid size of 0.125 Å × 0.125 Å, and the results are shown in Figure 7. Irrespective of the temperature of the system, the distribution is characterized by five ordered high intensity peaks appearing at the IBS of the protein (marked as a to e). These peaks originate from sets of water molecules aligned parallel to the z axis near the IBS. The spacings between the peaks a and b and between c and e are found to be ∼7.5 Å, which is close to the lattice separations in the basal (7.83 Å) and the primary prism (7.35 Å) planes of ice. This is an interesting observation that demonstrates the tendency of these sets of water molecules to form planes at the IBS that are analogous to the primary prism or the basal plane of ice even at room temperature. To explore such behavior further, in Figure 8 we plot the three-dimensional water density isosurface (with local densities of 0.08, 0.10 water Å−3 for systems S1 and S2, respectively) around the central four β-strands of the IBS. For

4

2 ⎛ 1⎞ ⎜cos ψ + ⎟ jk ⎝ 3⎠ k=j+1

∑ ∑ j=1

(1)

where ψjk is the angle between the bond vectors rij and rik, where j and k are the four atoms that are nearest neighbors to the i-th water. The calculations are carried out by assuming that the non-hydrogen atoms of the protein can act as neighboring sites for the first layer of water molecules around the protein. We have calculated the order parameter distribution, P(qt), for those water molecules that are present within 5 Å from the protein non-hydrogen atoms. The results are shown in Figure 6. For comparison, the corresponding distributions for pure bulk water as obtained from separate MD simulations of TIP4P water at identical temperatures are shown in the figure. As expected, on lowering the temperature tetrahedral ordering of waters increases noticeably. This is reflected in narrower distributions with maximum peak positions shifted toward 4748

dx.doi.org/10.1021/jp412528b | J. Phys. Chem. B 2014, 118, 4743−4752

The Journal of Physical Chemistry B

Article

between the Thr hydroxyl oxygens on the IBS are close to the lattice separations in the primary prism and basal planes of ice. However, the relative arrangement between these two sets of waters and the Thr hydroxyl oxygens, as shown schematically in Figure 9, closely resemble arrangement of waters in the basal

Figure 7. Two-dimensional probability distributions of water molecules around the TmAFP protein on the xy plane (the protein helical axis is aligned along the z axis) as obtained from simulations S1 at 300 K and S2 at 220 K. Note the presence of five ordered high intensity peaks (marked as a to e) appearing at the IBS of the protein at both the temperatures. Figure 9. Two-dimensional schematic representation of relative arrangements of water oxygens in the basal (left) and primary prism (right) planes of ice. The oxygens shown in blue form a plane that lies above the plane formed by the oxygens shown in gray. How the arrangement of the Thr oxygens (blue spheres labeled with “T”) at the IBS and the sets a and b waters (blue spheres labeled with “a” and “b”, respectively) can preferentially map with the ice basal plane oxygens is indicated in the figure.

plane of ice. It is also revealed from Figure 8 that the water molecules that belong to the sets c, d, and e form another plane that is above and parallel to the IBS (xz plane), with sets c and e waters located above the two sets of Thr hydroxyl oxygens and set d waters located nearly above the set b waters. The relative arrangement and separations between these bound waters (sets c, d, and e) also match well with that corresponding to the basal plane of ice. These are important observations that demonstrate that irrespective of the temperature TmAFP exhibits two parallel layers of ice basal plane-like hydration pattern at the IBS. This allows the IBS hydration waters to align appropriately and register with the quasi-ice-like interfacial waters on the ice basal plane. The protein thus recognizes the ice surface in this manner and binds with it, thereby preventing further growth of ice. Davies and co-workers19,42 have recently used ice crystal morphology studies and demonstrated that hyperactive AFPs (including TmAFP) prevent ice growth at the basal plane. Further, they concluded from fluorescent microscopy analysis that the hyperactive AFPs in general have affinities to bind at the basal plane of ice, thereby protecting ice growth in the direction perpendicular to such plane.42 Therefore, our simulation results provide direct support to the binding affinity of TmAFP as predicted from these experiments. To the best of our knowledge, we have been able to provide for the first time a microscopic explanation of the mechanism behind the preferential binding of TmAFP, an important class of hyperactive AFP, at the basal plane of ice. Importantly, such ice-like hydration pattern of the protein seems to exist even at room temperature, as evident from the present study.

Figure 8. Three-dimensional water density isosurface plots (shown in green) around the IBS of the TmAFP protein as obtained from simulations S1 at 300 K and S2 at 220 K. The helical axis of the protein is oriented perpendicular to the plane in the top panel, while in the bottom panel it is oriented parallel to the plane. Four central threonine-cysteine-threonine (TCT) motifs of the IBS are drawn as sticks, and the locations of the five sets of periodically arranged water molecules with high densities aligned parallel to the IBS around the threonine hydroxyl oxygens are marked (a to e). As a reference, the positions of the crystal waters that are externally bound at the IBS and located in between the threonine residue side chains of the central four β-strands are shown as red spheres.

visual clarity the isosurface plots are displayed for water molecules that are present within −8.0 Å ≤ x ≤ 8.0 Å, −18.0 Å ≤ y ≤ 2.0 Å, and −12.5 Å ≤ z ≤ 12.5 Å with a grid dimension of 0.25 Å × 0.25 Å × 0.25 Å. As a reference, the positions of the crystal waters that are externally bound at the IBS and located in between the threonine residue side chains of the central four β-strands are also shown in the figure. It is evident from Figures 7 and 8 that the set a and set b water molecules are located in between the β-strands parallel to the xz plane (i.e., the IBS). Note that the set b water molecules that are separated by ∼4.6 Å and remain in plane with the Thr hydroxyl oxygens correspond to the crystal water positions.6 Thus, the oxygens of set a and set b waters and the Thr hydroxyl oxygens together form a plane at the IBS. Note that the average distance between the intrastrand Thr hydroxyl oxygens is around 7.44 Å, while that between the interstrand Thr hydroxyl oxygens at equivalent positions is around 4.64 Å.6 Therefore, the separations between the set a and set b waters and that

4. CONCLUSIONS In this work, we have carried out atomistic MD simulations of the antifreeze protein TmAFP in aqueous medium at two different temperatures, 300 and 220 K. Conformational rigidity 4749

dx.doi.org/10.1021/jp412528b | J. Phys. Chem. B 2014, 118, 4743−4752

The Journal of Physical Chemistry B

Article

(4) Harrison, K.; Hallett, J.; Burcham, T. S.; Feeney, R. E.; Kerr, W. L.; Yeh, Y. Ice Growth in Supercooled Solutions of Antifreeze Glycoprotein. Nature (London) 1987, 328, 241−243. (5) Graham, L. A.; Liou, Y. C.; Walker, V. K.; Davies, P. L. Hyperactive Antifreeze Protein from Beetles. Nature 1997, 388, 727− 728. (6) Liou, Y. C.; Tocilj, A.; Davies, P. L.; Jia, Z. Mimicry of Ice Structure by Surface Hydroxyls and Water of a β -Helix Antifreeze Protein. Nature 2000, 406, 322−324. (7) Tyshenko, M. G.; Doucet, D.; Davies, P. L.; Walker, V. K. The Antifreeze Potential of the Spruce Budworm Thermal Hysteresis Protein. Nat. Biotechnol. 1997, 15, 887−890. (8) Graether, S. P.; Kuiper, M. J.; Gagné, S. M.; Walker, V. K.; Jia, Z.; Sykes, B. D.; Davies, P. L. β -Helix Structure and Ice-Binding Properties of a Hyperactive Antifreeze Protein from an Insect. Nature 2000, 406, 325−328. (9) Graham, L. A.; Davies, P. L. Glycine-Rich Antifreeze Proteins from Snow Fleas. Science 2005, 310, 461−461. (10) Pentelute, B. L.; Gates, Z. P.; Tereshko, V.; Dashnau, J. L.; Vanderkooi, J. M.; Kossiakoff, A. A.; Kent, S. B. H. X-ray Structure of Snow Flea Antifreeze Protein Determined by Racemic Crystallization of Synthetic Protein Enantiomers. J. Am. Chem. Soc. 2008, 130, 9695− 9701. (11) Kristiansen, E.; Ramløv, H.; Højrup, P.; Pedersen, S. A.; Hagen, L.; Zachariassen, K. E. Structural Characteristics of a Novel Antifreeze Protein from the Longhorn Beetle Rhagium Inquisitor. Insect Biochem. Mol. Biol. 2011, 41, 109−117. (12) Hakim, A.; Nguyen, J. E.; Basu, K.; Zhu, D. F.; Thakral, D.; Davies, P. L.; Isaacs, F. J.; Modis, Y.; Meng, W. Crystal Structure of an Insect Antifreeze Protein and Its Implications for Ice Binding. J. Biol. Chem. 2013, 288, 12295−12304. (13) Sidebottom, C.; Buckley, S.; Pudney, P.; Twigg, S.; Jarman, C.; Holt, C.; Telford, J.; McArthur, A.; Worrall, D.; Hubbard, R.; Lillford, P. Heat-Stable Antifreeze Protein from Grass. Nature 2000, 406, 256− 256. (14) Worrall, D.; Elias, L.; Ashford, D.; Smallwood, M.; Sidebottom, C.; Lillford, P.; Telford, J.; Holt, C.; Bowles, D. A Carrot LeucineRich−Repeat Protein That Inhibits Ice Recrystallization. Science 1998, 282, 115−117. (15) DeVries, A. L.; Wohlschlag, D. E. Freezing Resistance in Some Antarctic Fishes. Science 1969, 163, 1073−1075. (16) Sicheri, F.; Yang, D. S. C. Ice-Binding Structure and Mechanism of an Antifreeze Protein from Winter Flounder. Nature 1995, 375, 427−431. (17) Gilbert, J. A.; Hill, P. J.; Dodd, C. E. R.; Laybourn-Parry, J. Demonstration of Antifreeze Protein Activity in Antarctic Lake Bacteria. Microbiology 2004, 150, 171−180. (18) Garnham, C. P.; Campbell, R. L.; Devies, P. L. Anchored Clathrate Waters Bind Antifreeze Proteins to Ice. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 7363−7367. (19) Scotter, A. J.; Marshall, C. B.; Graham, L. A.; Gilbert, J. A.; Garnham, C. P.; Davies, P. L. The Basis for Hyperactivity of Antifreeze Proteins. Cryobiology 2006, 53, 229−239. (20) Raymond, J. A.; DeVries, A. L. Adsorption Inhibition as a Mechanism of Freezing Resistance in Polar Fishes. Proc. Natl. Acad. Sci. U.S.A. 1977, 74, 2589−2593. (21) Wilson, P. W. Explaining Thermal Hysteresis by The Kelvin Effect. Cryo-Lett. 1993, 14, 31−36. (22) Pertaya, N.; Marshall, C. B.; DiPrinzio, C. L.; Wilen, L.; Thomson, E. S.; Wettlaufer, J. S.; Davies, P. L.; Braslavsky, I. Fluorescence Microscopy Evidence for Quasi-Permanent Attachment of Antifreeze Proteins to Ice Surfaces. Biophys. J. 2007, 92, 3663−3673. (23) Celik, Y.; Drori, R.; Pertaya-Braun, N.; Altan, A.; Barton, T.; Bar-Dolev, M.; Groisman, A.; Davies, P. L.; Braslavsky, I. Microfluidic Experiments Reveal That Antifreeze Proteins Bound to Ice Crystals Suffice to Prevent Their Growth. Proc. Natl. Acad. Sci. U.S.A. 2013, 110, 1309−1314. (24) Ba, Y.; Wongskhaluang, J.; Li, J. Reversible Binding of the HPLC6 Isoform of Type I Antifreeze Proteins to Ice Surfaces and the

of the ice-binding and nonice-binding surfaces (IBS and NIBS) of the protein and the structural arrangements of water molecules around those have been studied in detail. Efforts have been made to understand the microscopic origin of the ice-binding activity of the protein at its IBS. The calculations reveal that the overall structure of TmAFP is rigid in aqueous medium with features similar to that of its crystalline form. A closer examination of the simulated results shows that the IBS is relatively more rigid than the NIBS. This has important consequences as it is expected that for icebinding activity the IBS should avoid large fluctuations and maintain its flat rigid nature. It is demonstrated that the presence of five internally bound water molecules near the IBS is crucial in maintaining its rigidity, which in turn is important for the proteinʼs activity. Our investigations further reveal noticeably higher water structuring with several well-defined solvation shells over a longer distance and relatively lower local density around the IBS of the protein. It indicates preferential alignment of IBS water molecules toward an ice-like-ordered low-density arrangement. Besides, such layering of waters around the IBS is found to be an inherent feature of the protein as evident from similar structural preference of the surrounding solvent even at room temperature. Two-dimensional water distribution analysis and density isosurface calculations further show the existence of five sets of periodically arranged ordered water molecules aligned parallel to the IBS. Importantly, it is demonstrated that the strategically located hydroxyl oxygens of the Thr residues influence the arrangements of these sets of ordered waters on two parallel planes at the IBS that closely resemble the basal plane of ice. In addition, the IBS of the protein seems to retain such ice-like hydration pattern even at room temperature. We believe that the unique hydration behavior of TmAFP as observed in this study provides a molecular level understanding of the origin of its ice-binding activity at the basal plane of ice. It would be interesting to explore further how the nonuniform rigidity of the IBS and NIBS of TmAFP affects the dynamics of the surrounding water molecules. We are currently investigating these aspects in our laboratory.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This study was supported in part by grants from the Department of Science and Technology (DST) (SR/S1/PC23/2007), Government of India. Part of the work was carried out using the computational facilities created under DST-FIST programme (SR/FST/CSII-011/2005) and DST-IYC award. U.S.M. thanks the Council for Scientific and Industrial Research (CSIR), New Delhi, for providing a scholarship.

■

REFERENCES

(1) Duman, J. G. Antifreeze and Ice Nucleator Proteins in Terrestrial Arthropods. Annu. Rev. Physiol. 2001, 63, 327−357. (2) Jia, Z.; Davies, P. L. Antifreeze Proteins: An Unusual ReceptorLigand Interaction. Trends Biochem. Sci. 2002, 27, 101−106. (3) Yeh, Y.; Feeney, R. E. Antifreeze Proteins: Structures and Mechanisms of Function. Chem. Rev. 1996, 96, 601−617. 4750

dx.doi.org/10.1021/jp412528b | J. Phys. Chem. B 2014, 118, 4743−4752

The Journal of Physical Chemistry B

Article

(47) Allen, M. P.; Tidlesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, 1987. (48) Darden, T.; York, D.; Pedersen, L. J. Particle Mesh Ewald: An N·log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (49) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E., III; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (50) MacKerell, A. D., Jr.; Feig, M.; Brooks, C. L., III. Improved Treatment of the Protein Backbone in Empirical Force Fields. J. Am. Chem. Soc. 2004, 126, 698−699. (51) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (52) Nutt, D. R.; Smith, J. C. Molecular Dynamics Simulations of Proteins: Can the Explicit Water Model Be Varied? J. Chem. Theory Comput. 2007, 3, 1550−1560. (53) Daley, M. E.; Spyracopoulos, L.; Jia, Z.; Davies, P. L.; Sykes, B. D. Structure and Dynamics of a β -Helical Antifreeze Protein. Biochemistry 2002, 41, 5515−5525. (54) Le Grand, S. M.; Marz, K. M., Jr. Rapid Approximation to Molecular Surface Area Via the Use of Boolean Logic and Look-up Tables. J. Comput. Chem. 1993, 14, 349−352. (55) Bandyopadhyay, S.; Chakraborty, S.; Balasubramanian, S.; Pal, S.; Bagchi, B. Atomistic Simulation Study of the Coupled Motion of Amino Acid Residues and Water Molecules Around Protein HP-36: Fluctuations at and Around the Active Sites. J. Phys. Chem. B 2004, 108, 12608−12616. (56) Sinha, S. K.; Bandyopadhyay, S. Differential Flexibility of the Secondary Structures of Lysozyme and the Structure and Ordering of Surrounding Water Molecules. J. Chem. Phys. 2011, 134, 115101. (57) Pal, S.; Bandyopadhyay, S. Importance of Protein Conformational Motions and Electrostatic Anchoring Sites on the Dynamics and Hydrogen Bond Properties of Hydration Water. Langmuir 2013, 29, 1162−1173. (58) Astley, T.; Birch, G. G.; Drew, M. G. B.; Rodger, P. M.; Wilden, G. R. H. Effect of Available Volumes on Radial Distribution Functions. J. Comput. Chem. 1998, 19, 363−367. (59) Cui, J.; Battle, K.; Wierzbicki, A.; Madura, J. D. Investigations of Structure and Dynamics of Water Solvation of the Type I Antifreeze Protein. Int. J. Quantum Chem. 2008, 109, 73−80. (60) Argyris, D.; Tummala, N. R.; Striolo, A.; Cole, D. R. Molecular Structure and Dynamics in Thin Water Films at the Silica and Graphite Surfaces. J. Phys. Chem. C 2008, 112, 13587−13599. (61) Argyris, D.; Cole, D. R.; Striolo, A. Hydration Structure on Crystalline Silica Substrates. Langmuir 2009, 25, 8025−8035. (62) Godawat, R.; Jamadagni, S. N.; Garde, S. Characterizing Hydrophobicity of Interfaces by Using Cavity Formation, Solute Binding, and Water Correlations. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 15119−15124. (63) Tan, O. Z.; Tsai, K. H.; Wu, M. C. H.; Kuo, J. L. Structural and Dynamic Properties of Water on the GaN Polar Surface. J. Phys. Chem. C 2011, 115, 22444−22450. (64) Ebbinghaus, S.; Meister, K.; Prigozhin, M. B.; DeVries, A. L.; Havenith, M.; Dzubiella, J.; Gruebele, M. Functional Importance of Short-Range Binding and Long-Range Solvent Interactions in Helical Antifreeze Peptides. Biophys. J. 2012, 103, L20−L22. (65) Meister, K.; Ebbinghaus, S.; Xu, Y.; Duman, J. G.; DeVries, A.; Gruebele, M.; Leitner, D. M.; Havenith, M. Long-Range Protein Water Dynamics in Hyperactive Insect Antifreeze Proteins. Proc. Natl. Acad. Sci. U.S.A. 2013, 110, 1617−1622. (66) Finney, J. L. A Procedure for the Construction of Voronoi Polyhedra. J. Comput. Phys. 1979, 32, 137−143.

Antifreeze Mechanism Studied by Multiple Quantum Filtering-Spin Exchange NMR Experiment. J. Am. Chem. Soc. 2003, 125, 330−331. (25) Sharp, K. A. A Peek at Ice Binding by Antifreeze Proteins. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 7281−7282. (26) Knight, C. A.; Cheng, C. C.; DeVries, A. L. Adsorption of αHelical Antifreeze Peptides on Specific Ice Crystal Surface Planes. Biophys. J. 1991, 59, 409−418. (27) Wen, D.; Laursen, R. A. A Model for Binding of an Antifreeze Polypeptide to Ice. Biophys. J. 1992, 63, 1659−1662. (28) Madura, J. D.; Wierzbicki, A.; Harrington, J. P.; Maughon, R. H.; Raymond, J. A.; Sikes, C. S. Interactions of the D and L-Forms of Winter Flounder Antifreeze Peptide with the {201} Planes of Ice. J. Am. Chem. Soc. 1994, 116, 417−418. (29) Cheng, A.; Merz, K. M., Jr. Ice-Binding Mechanism of Winter Flounder Antifreeze Proteins. Biophys. J. 1997, 73, 2851−2873. (30) Haymet, A. D. J.; Ward, L. G.; Harding, M. M.; Knight, C. A. Valine Substituted Winter Flounder ‘Antifreezeʼ: Preservation of Ice Growth Hysteresis. FEBS Lett. 1998, 430, 301−306. (31) Deluca, C. L.; Davies, P. L.; Ye, Q.; Jia, Z. The Effect of Steric Mutations on the Structure of Type III Antifreeze Protein and its Interaction with Ice. J. Mol. Biol. 1998, 275, 515−525. (32) Bar, M.; Celik, Y.; Fass, D.; Braslavsky, I. Interactions of β -Helical Antifreeze Protein Mutants with Ice. Cryst. Growth Des. 2008, 8, 2954−2963. (33) Dalal, P.; Knickelbein, J.; Haymet, A. D. J.; Sönnichsen, F. D.; Madura, J. D. Hydrogen Bond Analysis of Type 1 Antifreeze Protein in Water and the Ice/Water Interface. PhysChemComm. 2001, 4, 32−36. (34) Henn, A. R.; Kauzmann, W. Equation of State of a Random Network, Continuum Model of Liquid Water. J. Phys. Chem. 1989, 93, 3770−3783. (35) Gallagher, K. R.; Sharp, K. A. Analysis of Thermal Hysteresis Protein Hydration Using the Random Network Model. Biophys. Chem. 2003, 105, 195−209. (36) Doxey, A. C.; Yaish, M. W.; Griffith, M.; McConkey, B. J. Ordered Surface Carbons Distinguish Antifreeze Proteins and Their Ice-Binding Regions. Nat. Biotechnol. 2006, 24, 852−855. (37) Karim, O. A.; Haymet, A. D. J. The Ice/Water Interface: A Molecular Dynamics Simulation Study. J. Chem. Phys. 1988, 89, 6889− 6896. (38) Nutt, D. R.; Smith, J. C. Dual Function of the Hydration Layer around an Antifreeze Protein Revealed by Atomistic Molecular Dynamics Simulations. J. Am. Chem. Soc. 2008, 130, 13066−13073. (39) Smolin, N.; Daggett, V. Formation of Ice-like Water Structure on the Surface of an Antifreeze Protein. J. Phys. Chem. B 2008, 112, 6193−6202. (40) Siemer, A. B.; McDermott, A. E. Solid-State NMR on a Type III Antifreeze Protein in the Presence of Ice. J. Am. Chem. Soc. 2008, 130, 17394−17399. (41) Siemer, A. B.; Huang, K. Y.; McDermott, A. E. Protein−Ice Interaction of an Antifreeze Protein Observed with Solid-State NMR. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 17580−17585. (42) Bar-Dolev, M.; Celik, Y.; Wettlaufer, J. S.; Davies, P. L.; Braslavsky, I. New Insights into Ice Growth and Melting Modifications by Antifreeze Proteins. J. R. Soc. Interface 2012, 9, 3249−3259. (43) Pertaya, N.; Marshall, C. B.; Celik, Y.; Davies, P. L.; Braslavsky, I. Direct Visualization of Spruce Budworm Antifreeze Protein Interacting with Ice Crystals: Basal Plane Affinity Confers Hyperactivity. Biophys. J. 2008, 95, 333−341. (44) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorsid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (45) Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. Constant Pressure Molecular Dynamics Simulation: The Langevin Piston Method. J. Chem. Phys. 1995, 103, 4613−4621. (46) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327−341. 4751

dx.doi.org/10.1021/jp412528b | J. Phys. Chem. B 2014, 118, 4743−4752

The Journal of Physical Chemistry B

Article

(67) Errington, J. R.; Debenedetti, P. G. Relationship Between Structural Order and the Anomalies of Liquid Water. Nature 2001, 409, 318−321. (68) Chau, P. L.; Hardwick, A. J. A New Order Parameter for Tetrahedral Configurations. Mol. Phys. 1998, 93, 511−518. (69) Sharma, R.; Chakraborty, S. N.; Chakravarty, C. Entropy, Diffusivity, and Structural Order in Liquids with Waterlike Anomalies. J. Chem. Phys. 2006, 125, 204501.

4752

dx.doi.org/10.1021/jp412528b | J. Phys. Chem. B 2014, 118, 4743−4752