Langmuir 2004, 20, 5353-5357
5353
Molecular Dynamics Study of the Interactions of Ice Inhibitors on the Ice {001} Surface Surya Devarakonda,† James M. B. Evans,‡,§ Alfred Y. Lee,‡ and Allan S. Myerson*,‡ Department of Chemical Engineering, Polytechnic University, 6 Metrotech Center, Brooklyn, New York 11201, and Department of Chemical Engineering, Illinois Institute of Technology, 3301 South Dearborn Street, Chicago, Illinois 60616 Received March 14, 2003. In Final Form: March 24, 2004 The interactions of antifreeze protein (AFP) type I, antifreeze glycoproteins, polyvinyl pyrrolidone (PVP), and various amino acids with ice are investigated using Cerius2, a molecular modelling tool. Binding energies of these additives to a major ice crystal face {001} are computed. Binding energy comparison of threonine molecules (by themselves) and as threonine residues within AFP type I demonstrate their role in improving AFP’s binding ability to the ice crystal face. The shifts in onset points of ice crystallization with AFP type I, PVP, and amino acids are measured using differential scanning calorimetry. These values when correlated with their respective binding energies reveal a direct proportionality and demonstrate AFP’s effectiveness in inhibiting growth and nucleation of ice, over amino acids.
I. Introduction In the past two decades or so, significant progress has been made in characterizing several classes of proteins capable of ice interaction. These include two functionally distinct and opposite classes of proteins: antifreeze proteins (AFPs), which inhibit ice crystal growth, and ice nucleation proteins,1 which promote ice growth. To date, three main classes of AFPs, types I, II, and III,2-4 have been identified with the unique capability of inhibiting ice crystal growth via an “adsorption inhibition mechanism”,2 where the AFPs adsorb preferentially with specific planes of ice crystals, as opposed to binding to all surfaces. It is this preferential adsorption which changes the growth rate of the various surface forms and, hence, the morphology of the ice crystal. Previous works1,5 reports that a variety of ice habits and growth rates in the presence/ absence of AFPs have been observed. In recent years, this potential to reduce the growth and the damage of freezing has attracted increasing interest because the ability to control the rate and growth of ice crystals has a variety of industrial applications. A review of the current literature showed that the majority of the investigations have been carried out on the HPLC6 winter flounder isoform (an AFP type I protein).6-8 HPLC6 has a helical structure that is 37 amino acids long, and consists of three complete 11 amino acid * To whom correspondence should be addressed. Tel.: +1 312 567 3163. Fax: +1 312 567 5205. E-mail:
[email protected]. † Polytechnic University. Current address: Dr. Reddy’s Labs, Bulk Actives Unit II, Bollaram, Hyderabad 502 325, India. ‡ Illinois Institute of Technology. § Current address: Glaxo Wellcome Manufacturing Pte, Ltd., 1 Pioneer Sector 1, Singapore 628413, Singapore. (1) Abe, K.; Watabe, S.; Emori, Y.; Watanabe, M.; Arai, S. FEBS Lett. 1989, 258, 2, 297-300. (2) Raymond, J. A.; DeVries, A. L. Proc. Natl. Acad. Sci. U.S.A. 1977, 74, 2589-2593. (3) Gronwald, W.; Loewen, M. C.; Lix, A. J.; Daugulis, So¨nnichsen, F. D.; Davies, P. L.; Sykes, B. D. Biochemistry 1998, 37, 4712-4721. (4) Jia, Z.; Deluca, Ci; Chao, H.; Davies, P. L. Nature 1996, 384, 6606, 285-289. (5) Scholander, P. F.; Maggert, J. E. Cryobiology 1971, 8, 371-374. (6) Chao, H.; Houston, M. E.; Hodges, R. D.; Lay, C. M.; Syker, B. D.; Loewen, M. C.; Davies, P. L.; So¨nnichsen, F. D. Biochemistry 1997, 36, 14652-14660.
repeats of Thr-X2-Asx-X7 (where X is generally alanine) and additional N- and C-capping sequences. A number of molecular dynamics (MD) and Monte Carlo simulations have been undertaken to identify the nature of the interactions and the key ice planes. These initial studies9-13 proposed that the binding of AFP type I to an ice surface occurs because of the hydrogen bonding between the threonine residue and the oxygen atoms and the repelling of the water by the hydrophobic alanine residue. More recent studies6,14,15 postulate a reduced role for hydrogen bonding in the binding of AFPs to an ice surface and that the system will be dominated by van der Waals interactions. Experimental support for this was obtained by Chao et al.6 and Haymet et al.16 who showed that the binding residues of HPLC6 could be mutated with serine, valine, alanine, and glycine substituted for the threonine, and these mutated systems showed varying degrees of ice inhibition. The strongest evidence for non-hydrogen bonding interactions controlling the binding was obtained with the substitution of serine (LSAAN), which should have a similar hydrogen bonding potential and valine (LVAAN) for threonine.6 In LSAAN, there was a catastrophic loss in activity, yet in LVAAN, there was only a minor loss of activity. The change in activity was found to be associated not with a hydroxyl group and its associated hydrogen bonding but with the presence or absence of a γ-methyl group. The key ice planes that an AFP type I binds to have been identified14 as the [2021] (7) Zhang, W.; Laursen, R. A. J. Biol. Chem. 1998, 273, 52, 3480634812. (8) Wen, D.; Laursen, R. A. J. Biol. Chem. 1992, 267, 14102-14108. (9) Chou, C. K. J. Mol. Biol. 1992, 223, 509-517. (10) Jorgensen, H.; Mori, M.; Matsui, H.; Kaaoka, M.; Yanagi, H.; Yabusaki, Y.; Kikuzono, Y. Protein Eng. 1993, 6, 19. (11) McDonald, S. M.; Shawn, M.; Brady, J. W.; Clancy, P. Biopolymers 1993, 33, 1481-1504. (12) Cheng, A.; Merz, K. M. Biophys. J. 1997, 73 (6), 2851-2873. (13) Hong, J. S.; Jung, D. H.; Jhon, M. S. Mol. Simul. 1998, 20, 303. (14) Dalal, P.; So¨nnichsen, F. D. J. Chem. Inf. Comput. Sci. 2000, 40, 1276-1284. (15) Dalal, P.; Knickelbein, J.; Haymet, A. D. J.; So¨nnichsen, F. D.; Madura, J. D. PhysChemComm 2001, 1, 32-37. (16) Haymet, A. D. J.; Ward, L. G.; Harding, M. M. J. Am. Chem. Soc. 1999, 121, 941-948.
10.1021/la0344377 CCC: $27.50 © 2004 American Chemical Society Published on Web 05/20/2004
5354
Langmuir, Vol. 20, No. 13, 2004
Devarakonda et al.
Table 1. Energies (Given in kcal/mol) and Thermal Hysteresis of the Additives with the Ice Crystal Surface {001} additive
conformational energy
binding energy
∆Tf (°C)
valine tyrosine threonine proline lysine arginine serine aspartate glutamate tryptophan histidine HPLC6 (1000 ppm) PVP24
31.2 30.1 32.6 29.3 35.2 34.9 54.0 28.8 39.4 79.3 26.8 116.0 102.2 (Mw 11K)
11.9 3.2 -8.2 -17.7 -22.5 -24.0 -28.4 -32.9 -42.7 -49.1 -58.2 -203.1 -179.8
-0.9 -0.7 -0.5 0.2 1.2 1.6 1.8 1.9 2.8 4.2 4.9 5.2 6.0
plane, the prism [1010] plane, the basal [0001] plane, and a secondary basal [1120] plane. The purpose of this work was to investigate the potential of various L-amino acids (Table 1) to act as growth inhibitors of the {001} ice surface via their binding energies, which were calculated via MD simulations. The {001} surface was of interest because Scholander and Maggert5 and Raymond and DeVries2 showed that in the presence of AFPs ice will grow as a fine needle and that the inhibition occurs because the AFP adsorbs preferentially onto the prism face {001}. In addition to this MD, simulations were also undertaken on HPLC6 and polyvinyl pyrrolidone (PVP). Differential scanning calorimetry (DSC) was used to measure thermal hysteresis, which is a measure of the degree of crystal growth inhibition, and was used to enable a realistic comparison with the MD simulations. II. Materials and Apparatus The amino acids and PVP were purchased from Sigma Chemical Co. and used at a 5000 ppm concentration in double-distilled water. AFP type I was purchased from A/F Protein, Inc. A Perkin-Elmer DSC-7 was used to determine the thermal hysteresis. Cerius2 molecular modelling software17 was used to optimize the molecular structures and predict their binding energies. The procedures used in this study will be presented in two sections: first, the modelling procedure used for calculating the binding energies and, second, measurement of the thermal hysteresis using DSC. III. Modelling Section III.1. Binding Energy Calculation. The molecular modelling procedure employed Cerius2 (ref 17) to compute binding energies is presented in the following section. The force field Drieding 2.2118 was used for all the molecular mechanics and dynamics simulations. The procedure consists of three general steps. III.2. Building the Additive Molecules. The chemical structures of the amino acid molecules were either obtained from the Cambridge Structural Database19 or built using the sketcher module available in Cerius2.17 The amino acids were then geometrically optimized, and their partial atomic charges were calculated via MOPAC20 (17) Cerius2; Molecular Simulations, Inc.: San Diego, CA, 1999. (18) Mayo, S. L.; Olafson, B. D.; Goddard, W. A., III. J. Phys. Chem. 1990, 94, 8897-8909. (19) Cambridge Structural Database; Cambridge Crystallographic Data Center: Cambridge, 1999. (20) Stewart, J. J. P. MOPAC. J. Comp.-Aided Mol. Des. 1990, 4, 1-45.
Figure 1. HPLC6 isoform (37 amino acid monomeric R-helical protein with three 11 amino acid repeats: D TASDAAAAAAL TAANAKAAAEL TAANAAAAAAA TAR).
using a modified neglect of diatomic overlap (MNDO) Hamiltonian approximation.21 Molecular mechanics and dynamics were then carried out on each of the structures (thus, ensuring that the minimum energy configuration had been achieved). The minimum energy configuration was then saved, and the conformational energy was recorded for further use. PVP was built using the threedimensional sketcher module; the first step was to build a single unit of PVP (110 g/mol), and this was then used to simulate the polymer of a desired molecular weight. Given the high molecular weight of PVP (360 000 g/mol) and the fact that this would be computationally intensive to simulate, a smaller PVP of molecular weight 11 000 g/mol was used in the molecular simulations. The HPLC6 isoform, 1WFA,22 was obtained from the protein databank (Figure 1) and loaded directly into Cerius2 (ref 17), where it was geometrically optimized and the minimum energy conformation was obtained via molecular mechanics and dynamics. III.3. Building the Ice Unit Cell and Morphology Prediction. The crystal structure of ice (form Ih) was obtained from the inorganic section of the Cambridge Structural Database.19 Attachment energy calculations were used to identify the morphologically important faces {001}, {010}, and {0-21}. In this study, only the {001} surface was of interest. III.4. Binding Energy Calculation. The {001} face was selected and cleaved out of the crystal to a depth of 5 Å, then the surface was extended to 8 × 8, and a nonperiodic superstructure was created. The additive molecule was then placed onto the cleaved crystal surface and orientated in such a way as to increase the probability of hydrogen bonding or have atoms of opposing polarity in close proximity. The structure was then minimized for energy using molecular mechanics. Molecular mechanics used a conjugate gradient23 for energy minimization. The energy minimization was terminated once the residual (21) Dewar, M. J. S.; Thiel, W. J. J. Am. Chem. Soc. 1977, 99, 48994907. (22) Sicheri, F.; Yang, D. S. Ice binding structure and mechanism of an antifreeze protein from winter flounder. Nature, 1995, 375, 427431. (23) Fletcher, R.; Reeves, C. M. Comput. J. 1964, 7, 149-154.
Interactions of Ice Inhibitors on the Ice {001} Surface
Langmuir, Vol. 20, No. 13, 2004 5355
root-mean-square force on the structure falls below 0.1 kcal/(mol Å). MDs were performed for a minimum of 20 ps to find a “local” minimum for the given surface size, position, and orientation of the additive. Running the simulation for this length of time normally ensures that the structure has equilibrated. MDs simulation applies an external force on the molecule and solves the Newtonian equations of motion to compute the new atomic positions. During this process, the atoms can move from a higher energy state to a lower one. The algorithm used here was constant NVT (number, volume, and temperature) with the temperature set at 273 K. Ewald summation was used to enable the long-range electrostatic interactions to be incorporated into the binding energy calculation. The procedure was repeated with different additive positions and orientations until a constant total energy was found, thus, ensuring a “global minimum” (term loosely used here) position for the additive on the surface, within the limits of trial and error. A more detailed description of the procedure and definition of the molecular mechanics (MM) and MDs used in this study can be found in Devarakonda.24 The total energy includes all nonbond and H-bond (if any) interactions between the surface and the additive molecule. The binding energy is the difference between the total energy and the conformational energy of the additive molecule alone:
Ebind ) Etotal - Eadditive
(1)
The larger the absolute value of the binding energy, the higher the affinity of the additive molecule toward the crystal face. Thus, the impurities or additives, which have higher binding energies, are more effective in adsorbing to that particular face. By adsorbing selectively to the face, the additive can inhibit or promote the growth of a certain face. In addition to the binding energy, MD simulations enable the nature of the specific binding interactions to be evaluated, that is, hydrogen bonding, van der Waals, and electrostatic interactions. III.5. Assumptions and Limitations. Only one additive molecule per surface is used. The idea here is to compare the relative numbers of the binding energies but not the absolute values. The {001} ice surface in these simulations was fixed, and while it is recognized that a majority of MD simulations use flexible surfaces, previous works14,26 has shown that it is feasible to carry out simulations using a fixed or partially fixed surface. While the crystal surface is fixed, the additive was free to rotate and move along the surface. This should allow the additive molecule to orient itself on the surface and optimize the potential interactions, while MM and MD simulations are carried out. All calculations are carried out in the vacuum. IV. Experimental Section IV.1. Thermal Hysteresis. The thermal hysteresis of the various additives was found using DSC. A solution of known concentration (up to 5000 ppm) of additive in double-distilled water is prepared. About 20 µL of the solution is introduced into a sample pan using a micropipet and sealed using a press to prevent any leaks. The sample is then placed into a sample holder and cooled, at the rate of 3 °C/min, to a few degrees below the onset of crystallization. The onset temperature is noted and the process is repeated at least twice more for statistical accuracy. The cooling rate employed here is a tradeoff between the (24) Devarakonda, S. Ph.D. Thesis, Polytechnic University, Brooklyn, New York, 1997. (25) Harding, M. M.; Ward, L. G.; Haymet, A. D. J. Eur. J. Biochem. 1999, 264, 653-665. (26) Storr, M. T.; Taylor, P. C.; Monfort, J.-P.; Rodger, P. M. J. Am. Chem. Soc. 2004, 126 (5), 1569-1576.
sensitivity of the instrument and the accuracy of the measurement. At lower cooling rates, the sensitivity of the instrument is sacrificed, while at higher cooling rates, the accuracy of the measurement is lost. This procedure was repeated with all of the amino acids, PVP (Mw 360 000 g/mol), and the HPLC6 isoform at various concentrations in water. The difference between the onset point of ice crystallization with (T′) and without (T) an additive is the thermal hysteresis:
∆Tf ) T - T′
(2)
A positive value is desired, and the higher the value the more effective the additive becomes at inhibiting the crystal growth.
V. Results and Discussion The binding energies of the amino acids, PVP, and HPLC6 with the {001} surface are presented in Table 1. The ice surface used in the MD simulations was set at 8 × 8 (the use of a standard size for all simulations enables a more realistic comparison of the binding energies than the traditional approach where the surface size changes as a function of the additive size). The results in Table 1 show that, in general, the degree of crystal growth inhibition increases (∆Tf) as the affinity between the ice surface and the additive increases. The correlation between the amino acid binding energies and the thermal hysteresis is considered in more detail in Figure 2, which shows that there is a strong relationship between the degree of inhibition and the binding energy (R ) 96.8%). This is an important observation because it supports ones intuition and ties modelling to experiments. However, it should be noted that not all of the amino acids have the desired effect on lowering the freezing point of water. Valine, tyrosine, and threonine show a negative shift, that is, a negative ∆Tf. A negative ∆Tf suggests growth promotion.25 The results in Table 1 and Figure 2 also show that both HPLC6 isoform and PVP inhibit the crystallization of ice and that PVP is a potential alternative to AFPs in controlling the morphology, nucleation, and crystal size of ice. There was some initial concern with the fact that the PVP thermal hysteresis was significantly larger than the HPLC6 isoform value (Table 1); an extended thermal hysteresis series for PVP (Table 2), showed that the difference was purely a function of concentration and that at the same concentrations, the HPLC6 isoform has a larger thermal hysteresis than PVP. While the results in Table 1 provide an accurate comparison between the various amino acids, they are limited because they only compare single molecule additives against the AFP which is thought to have four amino acid residue binding sites. To perform a more realistic comparison between the AFP and threonine, (which is the usual binding residue for the AFP), four threonines were individually bound to an 8 × 8 surface and the binding energy was calculated (Table 3). The results in Table 3 show that increasing the number of threonine molecules that are individually bound to the ice surface does result in an increase in the binding energy. The HPLC6 isoform still has a significantly higher binding energy of -203.138 kcal/mol. This shows that the AFP has a much stronger affinity to bind with the ice surface (via the threonine residuals) than the threonine molecules by themselves. The difference in the binding energy of the HPLC6 isoform and the four threonine molecules primarily arises from the γ-methyl functional group of the threonine and the i + 4 and i + 8 alanine residues that are in close
5356
Langmuir, Vol. 20, No. 13, 2004
Devarakonda et al.
Figure 2. Thermal hysteresis as a function of the binding energy (kcal/mol). Table 2. Thermal Hysteresis of PVP as a Function of Concentration PVP concentration (ppm)
∆Tf (°C)
PVP concentration (ppm)
∆Tf (°C)
500 1000
2.7 3.3
2000 5000
5.5 6.0
Table 3. Breakdown of the Interactions between the Additives and the {001} Ice Surface (kcal/mol) additive threonine (1) threonine (4) HPLC6
binding energy van der Waals electrostatic hydrogen -8.2 -38.8 -203.1
6.6 34.5 238.0
15.6 58.6 -711.5
-0.5 -2.1 -403.9
proximity and contact with the ice surface.27,28 These residues are found on the hydrophobic surface of the HPLC6 helix, which has been shown to orientate itself toward the ice surface. A secondary cause of the increase in the binding energy is thought to be the fact that AFPs are shaped in such a way that a significant proportion of their surface will come into contact with the ice surface. The final stage of the analysis was to consider the specific nature of the nonbonded interactions between the additives and the ice surface. The results in Table 4 tend to suggest that, for the majority of additives, the binding to the ice surface is driven by electrostatic interactions. In terms of the amino acid series, the dominance of electrostatic interactions is hardly surprising because under normal conditions, most amino acids carry a charge (zwitterions). This is even true of some amino acid residues, for example, histidine and lysine, which have been incorporated into protein structures, in addition to this in a normal protein: The amino end carries a positive charge (-NH2+), the carboxyl end carries a negative charge (-CO2-), and ice has a large dipole moment. While the charges of the amino and carboxyl ends are important, it is thought that the majority of the ice/amino interactions are controlled/determined by the functionality of the side chain (see Table 4) because the chemistry of this group (27) Davies, P. L.; Baardsnes, J.; Kuiper, M. J.; Walker, V. K. Philos. Trans. 2002, 357 (1423), 927-935. (28) Baardsnes, J.; Kondejewski, L. S.; Hodges, R. S.; Chao, H.; Kay, C.; Davies, P. L. FEBS Lett. 1999, 463 (1-2), 87-91.
is known to have a significant impact on host/impurity interactions.29 For example, Table 4 shows that in general a larger electrostatic interaction value can be expected for systems where the side chains carry a formal electric charge. These are the acidic aspartate and glutamate and basic arginine, lysine, and histidine amino acids. The remainder of the interaction energy essentially comes from van der Waals interactions. There are three exceptions to this, valine, tyrosine, and threonine, and it is interesting to note that of the three systems where van der Waals interactions dominate the system, two of these systems (valine and tyrosine) are seen to promote the growth of ice, as defined by a negative thermal hysteresis (Table 1). As with the inhibiting systems, the promotion of ice growth is thought to primarily arise from the functionality of the side chains which for valine and tyrosine are hydrophobic. Thus, these two amino acids are more likely to repel water/ice than other amino acids which interact with and disrupt the growth of the ice surface. In terms of hydrogen bonding, the results in Table 4 suggest that, apart from the HPLC6 isoform, this interaction does not play a significant role in inhibiting the growth of the {001} ice surface and, even here, the system still appears to be governed by electrostatic interactions. The larger hydrogen bonding component of the HPLC6 isoform is thought to arise from several factors: first, the relative size of the AFP in comparison to the amino acids and, second, the fact that the isoform tends to orientate itself such that a significant portion of the molecule is in contact with the ice surface. VI. Conclusions Binding energy results show that the HPLC6 isoform of the AFP type I has the most affinity for the ice and that PVP is a potential alternative to the AFPs for controlling the morphology, nucleation, and crystal size of ice. The binding energies of the amino acids can be correlated to the thermal hysteresis, thus, tying together modelling (29) Rastall, R. A. Protein Structure and Function. http://www.fst.rdg.ac.uk/courses/fs916/lect1/lect1.htm (accessed Feb 2004).
Interactions of Ice Inhibitors on the Ice {001} Surface
Langmuir, Vol. 20, No. 13, 2004 5357
Table 4. Breakdown of the Nonbond Interactions between the Additives and the {001} Ice Surface (All Values Are Given in kcal/mol) additive
binding
van der Waals
electrostatic
hydrogen
side chain
valine tyrosine threonine proline lysine arginine serine aspartate glutamate tryptophan histidine HPLC6
11.9 3.2 -8.2 -17.7 -22.5 -24.0 -28.4 -32.9 -42.6 -49.1 -58.2 -203.1
12.9 12.6 6.6 6.4 9.4 13.9 4.4 6.8 5.8 6.7 6.4 238.0
1.0 2.6 15.6 -17.6 -14.8 -33.9 -0.7 -32.1 -29.7 -21.9 -22.6 -711.5
0.0 0.0 -0.5 -0.1 × 10-1 -0.9 × 10-4 0.0 -0.1 -0.9 -7.1 × 10-2 0.0 0.0 -403.8
aliphatic aromatic hydroxyl secondary amino basic basic hydroxyl acidic acidic aromatic basic N/A
and experimental results. These results also demonstrate the effectiveness of the additive inhibition capacity via binding. The procedure and results presented in this article demonstrate a simple yet effective way to compare modelling and experimental results, while providing some insight into how additives (especially AFPs) interact with the ice crystal surface. Although progress has been made in identifying the molecular interactions between the
additives and the ice surface, it is clear that the antifreeze mechanism is not fully understood, as we report for the first time qualitative values of the electrostatic interactions obtained from MD simulations. From these values, it is apparent that electrostatic interactions play a significant role in ice inhibition. LA0344377
