1. Structural Biology and Molecular Biophysics
  2. Computational and Systems Biology
Download icon

The biological function of an insect antifreeze protein simulated by molecular dynamics

  1. Michael J Kuiper  Is a corresponding author
  2. Craig J Morton
  3. Sneha E Abraham
  4. Angus Gray-Weale
  1. The University of Melbourne, Australia
  2. St Vincent's Institute of Medical Research, Australia
Research Article
  • Cited 35
  • Views 4,960
  • Annotations
Cite this article as: eLife 2015;4:e05142 doi: 10.7554/eLife.05142

Abstract

Antifreeze proteins (AFPs) protect certain cold-adapted organisms from freezing to death by selectively adsorbing to internal ice crystals and inhibiting ice propagation. The molecular details of AFP adsorption-inhibition is uncertain but is proposed to involve the Gibbs–Thomson effect. Here we show by using unbiased molecular dynamics simulations a protein structure-function mechanism for the spruce budworm Choristoneura fumiferana AFP, including stereo-specific binding and consequential melting and freezing inhibition. The protein binds indirectly to the prism ice face through a linear array of ordered water molecules that are structurally distinct from the ice. Mutation of the ice binding surface disrupts water-ordering and abolishes activity. The adsorption is virtually irreversible, and we confirm the ice growth inhibition is consistent with the Gibbs–Thomson law.

https://doi.org/10.7554/eLife.05142.001

eLife digest

Water expands as it freezes. If this happens to the water inside plants and animals, the resulting ice crystals can rupture cells. To prevent this, many plants and animals that live in cold climates have evolved ‘antifreeze proteins’. When a small particle of ice first starts to form, the antifreeze proteins bind to it and prevent the water around it freezing, hence preventing the growth of an ice crystal.

There are many different types of antifreeze protein, and some are more active than others. For example, some insects including the spruce budworm are exposed to extremely cold temperatures—sometimes below −30°C—and these insects have antifreeze proteins that are highly active.

It is not fully understood how different antifreeze proteins interact with ice and prevent the growth of ice crystals. This is largely because, as yet, there are no experimental techniques that make it possible to see how antifreeze proteins and water molecules arrange themselves at the surface of a growing particle of ice. Instead, scientists have developed computer simulations to investigate this process. While many of these studies have provided valuable information, the computational methods used have only recently become powerful enough to analyze how the antifreeze proteins approach the surface of the ice particle.

Kuiper et al. carried out simulations involving a highly active antifreeze protein from the spruce budworm. The results of these simulations revealed that this antifreeze protein does not bind directly to ice; instead, water molecules at the surface of the protein act as a bridge between the protein and the ice. These water molecules are highly ordered and though they have similarities with how water is structured in the ice, they are distinct from the ice lattice itself. Furthermore, this arrangement appears to be important for allowing the spruce budworm antifreeze protein to interact with the ice.

This study provides detailed insights as to how a highly active antifreeze protein helps to prevent ice crystals forming. In the future, the computational simulations used here may be extended to study the dynamics of other antifreeze proteins, and also how crystals of other materials form.

https://doi.org/10.7554/eLife.05142.002

Introduction

Many species of fish, insects, plants and micro-organisms living in cold environments produce antifreeze proteins (AFPs) whose main function is to target and modify the growth of regular ice (Raymond and DeVries, 1977; Davies and Hew, 1990). By binding to ice, AFPs are thought to both lower the freezing point and raise the melting point of ice through the Gibbs–Thomson effect whereby adsorbed AFPs cause increased micro-curvature of the ice front (Yeh and Feeney, 1996). The resulting difference between the melting and freezing points, known as the thermal hysteresis gap (TH), can range from a few tenths to over 6°, depending on the specific AFP (Lin et al., 2011).

Regular hexagonal ice (Ih) has three equivalent a-axes (a1a3) which are parallel to three respective prism planes and perpendicular to the c-axis. The ice basal plane lies perpendicular to the c-axis. Remarkably AFPs have evolved independently numerous times to bind stereo-selectively to ice in various orientations, resulting in a structurally diverse protein family (Hew et al., 1985; Ko et al., 2003; Pentelute et al., 2008; Hakim et al., 2013; Middleton et al., 2013). The potency of different AFPs reflects their adsorption characteristics and environmental adaptation. Polar fish, for example, do not normally experience temperatures much lower than the freezing point of sea water (−1.9°C), and have AFPs that provide protection to just a few tenths of a degree below this. Fish AFPs often bind oblique to ice prism planes, typically forming hexagonal bipyramid crystals at temperatures within the TH gap which rapidly grow as elongated spicules along the c-axis below the TH gap (Davies and Hew, 1990). In contrast insects which are exposed to much lower terrestrial temperatures, (often below −30°C), require AFPs with higher activity, possibly by selecting more effective prism and basal plane binding (Lin et al., 2011). Despite a wide selection of available structural and experimental data, a complete molecular description of how an AFP binds to its respective ice planes and inhibits crystal growth remains elusive.

The difficulty in characterizing the molecular mechanisms of AFPs arises largely due to their unusual relationship with water, which in this case serves as both the protein's solvent and, as ice, its target (Jia and Davies, 2002). Slowly freezing solutions of dilute AFP will incorporate AFPs into the growing ice phase at roughly their original concentration, while non-AFPs and solutes are largely excluded (Marshall et al., 2004). Many structural studies have shown a correlation between the periodicity of the ice-binding surface elements of AFPs and their respective ice-binding planes, though it has not been entirely clear whether AFPs bind directly to ice or through some ordered water intermediate (Sharp, 2011). Observation of the molecular arrangements of AFPs and surrounding water as they bind and modify the surface curvature of ice is currently technically unattainable. Given these experimental limitations, computer simulations have played a crucial role in expanding our understanding of the AFP phenomena.

Early computational simulations of winter-flounder AFP (wfAFP) by Wen and Laursen (1992); Jorgensen et al. (1993); Madura et al. (1994) firmly established the correlation of the distance between regularly spaced polar threonine residues and that of water molecules of the ice lattice of the {2021} pyramidal ice plane in the [112] direction by searching for energetically stable conformations. This was in agreement with previous experimental work of Knight et al which showed the preferential binding orientation of wfAFP on single crystal ice hemisphere (Knight et al., 1991) and also the later work of Sicheri and Yang who determined the wfAFP crystal structure (Sicheri and Yang, 1995). Cheng and Merz further expanded this work using molecular dynamics simulations of solvated wfAFP in the bound conformation to the {2021} pyramidal plane to propose key hydrogen bond interactions between polar residues and the ice surface, suggesting that hydrogen bonds were the primary driving force of ice adsorption (Cheng and Merz, 1997). However, this model was to be short lived as mutational work by the groups of Harding and Laursen (Zhang and Laursen, 1998; Haymet et al., 1999) raised questions about the nature of the wfAFP-ice interaction; conservative substitution of threonine residues for serine abolished activity, while threonine to valine substitution, which removed the proposed threonine hydroxyl–ice interactions surprisingly retained AFP activity.

Spurred by the paradoxical mutant wfAFP results both Dalal and Sönnichsen (2000) and Jorov et al. (2004) employed Monte Carlo simulations of the wfAFP on the pyramidal plane, both groups finding that hydrophobic interactions contributed significantly to the adsorption to ice. Further related work by Wierzbicki et al. (2007) proposed that wfAFP doesn't bind directly to ice but rather accumulates at the ice-water interface with preference for the hydrophobic over the hydrophilic face facing the ice. Additional computational studies of different AFPs have followed, including type II sea-raven AFP (Wierzbicki et al., 1997), type III eel pout AFP (Madura et al., 1996) and an insect AFPs from Tenebrio molitor (Liu et al., 2005). Nada and Furukawa (Nada and Furukawa, 2011) were amongst the first to provide a model of the spruce budworm AFP (sbwAFP) at the prism ice-water interface. Using rigid models of sbwAFP they showed that two pre-determined binding conformations were able to affect local ice growth kinetics, inferring that reduced ice growth and induced curvature was indicative freezing point depression consistent of the Gibbs–Thomson effect.

As highlighted by an excellent review by Nada and Furukawa. (Nada and Furukawa, 2012), in real world systems, AFPs bind to an ice-water interface rather than ice crystal planes alone, thus earlier AFP simulations are unable to provide complete mechanistic details of the AFP-ice interaction. Even though a number of AFP studies have simulated ice-water interactions (McDonald et al., 1995; Dalal et al., 2001; Wierzbicki et al., 2007; Nada and Furukawa, 2008), Nada and Furukawa also point out the need for simulations to target a growing ice-water interface. Ideally, simulations would allow unconstrained AFPs to migrate from the water phase to interact freely the growing ice-water interface without introducing prior binding orientation bias. The computational requirements to do so are significant, and only recently became relatively accessible for researchers.

The highly active AFP from the spruce budworm, Choristoneura fumiferana, (sbwAFP) provides an excellent model for antifreeze simulations and has been used previously in computational studies (Nutt and Smith, 2008; Nada and Furukawa, 2011). SbwAFP has a β-helix structure with a triangular cross-section forming three parallel β-sheet faces and a hydrophobic core as shown in Figure 1. The ice binding face has a regular array of repeating threonine-x-threonine (T-X-T) motifs whereby the spacing between threonine residues match the repeat spacing parallel to the c-axis of ice of approximately 7.4 Å and motifs on adjacent β-strands closely match the a-axis repeating dimensions of the prism face of ice of 4.5 Å. The second row of threonine residues of the ice-binding surface is particularly sensitive to mutation, with replacement of the threonines drastically reducing sbwAFP's potency as an AFP (Graether et al., 2000). In addition to inducing a non-colligative freezing point depression of ice by up to 6°C, sbwAFP also demonstrates a significant melting inhibition (super-heating) of ice of up to 0.44°C (Celik et al., 2010). Here we employ multiple and extended molecular dynamics simulations of a potent sbwAFP isoform 501 (Leinala et al., 2002) with a novel ice/water interface model designed to create continual ice growth conditions, providing new insight into molecular specificity and the resultant Gibbs–Thomson induced thermal hysteresis.

Cartoon representations of the spruce budworm AFP (sbwAFP) 501 isoform.

This shows the triangular cross section of parallel β-helix structure. The ice binding surface is drawn as sticks.

https://doi.org/10.7554/eLife.05142.003

Results and discussion

Building a model of the ice-water interface

The TIP4P water model (Jorgensen et al., 1983) was used for simulations due to its ability to give realistic phase transformation behavior and superior freezing characteristics compared with the default TIP3P model, its support by the NAMD software (Phillips et al., 2005), and its compatibility with biological CHARMM parameters (Glass et al., 2010). Although the melting point (Tm) of the TIP4P water model is known to be 230.5 K, (García Fernández et al., 2006) it is highly unlikely such a water model will spontaneously freeze or remain frozen in simulation, particularly in the few degrees temperature range below Tm where AFPs exhibit their biological activity, due to the equilibrium radius R of the critical ice nucleus governed by the Gibbs–Thomson effect. This effect is a well known phenomenon whereby the chemical potential across an interface varies with curvature. For positive interfacial energies smaller crystals, with their higher curvature, are only able to be in equilibrium with their melt at lower temperatures than larger crystals. AFPs are thought to exploit this property by effectively subdividing a large ice surface of relatively low curvature into smaller surfaces with greater curvature thus giving rise to lower freezing temperatures, as shown in Figure 2. The critical radius at a given temperature Tr can be approximated with the following expression adapted from Pereyra et al. (2011):

(1) RMwσLρiAgT0(T0Tr),

where Mw is the molecular weight of water, σ is the surface tension of the ice-water interface, L is the latent heat of fusion, ρi is the density of ice, T0 is the melting temperature at atmospheric pressure, Ag is a constant of 2 or 1 for spherical or cylindrical symmetry conditions respectively. This relationship implies large ice crystals are required to maintain growth. For example, assuming spherical geometry, for a single degree of supercooling the critical radius of an ice embryo is about 518 Å containing approximately 18 million water molecules, which is currently computationally prohibitive to model (see ‘Materials and methods’ for calculation). To simulate AFP-ice interactions more efficiently a specially arranged model of an ice-water interface was developed to provide a continual ice growth mechanism while mimicking an infinite ice plane by taking advantage of periodic boundary conditions (PBCs). The ice/water system employed here consisted of 17,434 water molecules in a 51.5 × 126.6 × 83.9 Å box with PBCs, of which 1692 water molecules were harmonically constrained to an ideal prism ice lattice. This served as the seed crystal while an additional 1049 water molecules below the seed were constrained in a 5 Å thick disordered layer to act as an ice barrier preventing downward ice growth, leaving a small wedge of mobile water between constrained layers as shown in Figure 3. A perspective view of the model is shown in Figure 3—figure supplement 1. The seed ice crystal included an additional small prism step to initiate step layer growth and the entire seed was angled such that the primary prism face was approximately 4.3° with respect to the long axis. The angle of the seed is such that step layer growth across the prism face reappears, due to PBCs, one layer higher to where it started, thus providing a continual ice growth mechanism. Unimpeded, this mechanism will freeze the area above the seed crystal at any temperature below Tm. This arrangement mimics step growth as might be produced by screw dislocations, as originally proposed by Frank (1949) and later observed on ice surfaces by Ketcham and Hobbs (1968).

Schematic figure of antifreeze protein (AFP) ice inhibition via the Gibbs–Thomson effect.

(A) A mixture of AFP and non-AFP proteins ahead of a growing ice front. (B) AFPs selectively bind to the ice, while non-AFPs are excluded. AFP adsorption to the ice surface subdivides into smaller growth fronts which increases local curvature. If the induced curvature of the new ice fronts is less than the critical radius (R) of an ice embryo at a given temperature then the ice remains in equilibrium with the melt.

https://doi.org/10.7554/eLife.05142.004
Figure 3 with 1 supplement see all
Section view of the sbwAFP simulation setup.

The seed ice crystal is shown in dark blue, constrained disordered water coloured cyan and unrestrained water is coloured red. The sbwAFP is centrally positioned 20 Å above the seed ice. The prism face of the seed ice crystals is angled approximately 4.3° with respect to the long axis of the simulation box, such that the periodic image of the seed ice is consistent with itself at the edges.

https://doi.org/10.7554/eLife.05142.005

As the simulation proceeds, water molecules join the ice phase through hydrogen bonding and can be easily identified by their tetrahedral arrangement, which appears hexagonal when viewed directly down the c-axis of ice. Within this ice/water box the unconstrained sbwAFP model was initially placed 20 Å above the seed, with its ice binding face directed towards the seed approximately aligned in the expected binding orientation. In our preliminary work the seed ice crystal was arranged with the prism face parallel with the axis of the simulation boundaries. Crystals from these arrangements did not grow well, extending the ice front very slowly if at all. We realized that this was because one had to effectively wait for a two-dimensional nucleation event to occur on the prism plane and reach a critical size before the next layer of ice would grow.

Characterization of the sbwAFP binding mode

To characterize the binding mode of sbwAFP, 32 independent, unconstrained (apart from the waters constrained in the ice seed and ice-barrier) simulations of at least 150 nanoseconds (ns) to as long as 250 ns duration were carried out at 228 K. Definitive ice binding was observed 26 times with each binding event occurring with equivalent orientation with respect to the ice lattice, though adsorbing at a number of different step heights. The remaining 6 simulations either diffused away from the growing ice front (3 times), or had partial engagement at alternate angles (3 times). Of the partial engagements, the AFP still had significant rotational and translational movement with respect to the ice lattice so was not considered truly bound. Each of the definitive ice binding events incorporated a linear array of 6 water molecules bridging between the prism face and a row of conserved threonine residues (Thr 7, 23, 39, 54, 69, 84 and 101) as shown in Figure 4. These water molecules had the same periodicity of the prism face, spaced approximately 4.5 Å apart. Interestingly they were distinct from the bulk ice crystal, bridging between two prism face water molecules parallel along the c-axis which would normally be bridged by two water molecules in regular ice as shown in Figure 5. The sbwAFP-bound water molecules had fully satisfied tetrahedral hydrogen bonding arrangements; two hydrogen bonds to the prism face ice water molecules and two hydrogen bonds to the hydroxyl groups of adjacent threonine residues of the sbwAFP ice binding face.

The ice-binding surface of sbwAFP.

Ordered water molecules are shown in yellow hydrogen bonded between adjacent threonine residues of the top row.

https://doi.org/10.7554/eLife.05142.007
The ice binding arrangement of spruce budworm AFP in relation to ice.

(A) c-axis view of ice (cyan spheres) and the arrangement of the sbwAFP ordered water (yellow spheres). (B) a-axis view of the sbwAFP binding orientation showing the relative arrangement of the ordered water to ice. The ice lattice repeats approximately 4.5 Å and 7.4 Å along the a axis and c axis respectively.

https://doi.org/10.7554/eLife.05142.008

To further investigate the function of the ordered waters of the ice binding surface 16 independent simulations of a sbwAFP mutant (with 4 Thr to Leu substitutions at positions 7, 21, 39, 69) known to disrupt activity (Graether et al., 2000), were performed, each of 250 ns length. No ice binding events were observed with all mutant AFPs diffusing ahead of the growing ice surface in distinct contrast to the wild type protein. Figure 6 shows the relative height coordinates of both wild type and mutant sbwAFP throughout the docking simulations clearly demonstrating the former AFP adsorbing and fixing to the ice surface while the latter, unable to bind, diffuses ahead of the growing ice front. The Thr to Leu substitutions disrupt the binding of the ordered water to sbwAFP by removing the stabilizing hydroxyl side chain hydrogen bond interactions of the threonine residues. These results indicate ordered waters are crucial in initiating sbwAFP recognition and binding to ice and are suggestive of longer range protein-water interactions consistent with recent experimental findings by Meinster et al, of a similar beetle AFP from Dendroides canadensis, measured with terahertz spectroscopy (Meister et al., 2013). Representative simulations of the binding and non-binding events of the respective wild type and mutant sbwAFP are shown in supplementary Videos 1, 2.

Simulations of ice docking of wildtype and inactive mutant sbwAFP against a growing ice surface.

(A) Representative side view of final docking orientation of wild type sbwAFP. (B) Overlay of 16 independent wild type sbwAFP docking simulations mapping the centre of mass z-axis coordinate of the AFP over time, an unchanging z-coordinate shows ice binding. (C) Representative side view of mutant sbwAFP in simulation. (D) Overlay of 16 independent mutant sbwAFP docking simulations showing all AFPs diffusing ahead of the ice front and no ice binding events. The mutant form of the sbwAFP consisted of 4 Thr to Leu mutations at positions 7, 21, 39 and 69 which disrupts the ice binding surface.

https://doi.org/10.7554/eLife.05142.009
Video 1
Representative simulation of spruce budworm AFP (sbwAFP) docking to ice.

Side view of unconstrained sbwAFP docking to the prism face of a growing ice surface at 228 K. (250 nanosecond total simulation.) The seed ice crystal is in dark blue.

https://doi.org/10.7554/eLife.05142.010
Video 2
Representative simulation of mutant sbwAFP.

Unconstrained sbwAFP mutant simulated at 228 K for 250 nanoseconds. Ice binding ability has been lost by mutation of 4 threonine residues (7, 21, 39 and 69) of the ice binding surface to leucine.

https://doi.org/10.7554/eLife.05142.011

Adsorption inhibition

Three contiguous ∼1 μs MD simulations of sbwAFP adsorbed to an active ice front at temperatures of 225 K, 230 K and 232 K were run to test and simulate the adsorption-inhibition mechanism by observing subsequent ice growth as shown in Figure 7. Representative cross-section snapshots of the simulation at the three temperatures is shown in Figure 7—figure supplement 1. Though the Gibbs–Thomson phenomena is widely accepted as the basis of AFP activity, to the best of our knowledge it has not been observed directly at the ice-water interface. In the first simulation the sbwAFP adsorbed from the solvent phase to the advancing ice front at 225 K (5 K below the model system Tm). The sbwAFP adsorbed to the ice front within 40 ns and within 150 ns the ice front displayed marked convex curvature creating cylindrical ice fronts as shown in Figure 8. By 250 ns the ice front had reached equilibrium with an averaged cylindrical radius of approximately 50 Å. The extent of the inhibited ice front fluctuated, advancing and retreating in the range of 5–8 Å showing the reversible arrangement of water molecules to the crystal ice lattice, (this is an important indication of reaching equilibration and adequate sampling given the reduced water diffusion at these low temperatures). The ice fraction (defined as a proportion of total ice-like water molecules above the ice seed crystal), was measured over the 280 to 1000 ns time frame and found to equilibrate at 0.46 ± 0.01 (SD). In contrast AFP-free simulations froze almost completely within 150 ns. The vacillating nature of the ice front in simulation makes quantitative measurement of the curvature difficult, however we find good qualitative agreement between theoretical Gibbs–Thomson curvature and the averaged ice front as shown in Figure 8 where the theoretical radius is overlaid on to the averaged simulation ice cross-section. This is good indication that the simulation ice curvature is due to the expected Gibbs–Thomson phenomena.

Figure 7 with 1 supplement see all
Freezing and melting inhibition of ice by adsorbed sbwAFP.

Ice fraction of 3 contiguous ∼1 microsecond simulation temperatures 225 K (approx. 5 K below melting point Tm), 230 K (approx. at Tm) and 232 K (approx. 2 K above Tm). Linear regression of equilibrated states shows no appreciable growth or melting (red lines). A simulation containing no AFP (blue line on left) shows rapid freezing of the entire model.

https://doi.org/10.7554/eLife.05142.012
Side and end views of steady-state equilibrated ice fronts inhibited by adsorbed sbwAFP.

Cyan spheres represent the averaged ice front (defined as water occupation of greater than 85% averaged over approximately 1 microsecond simulation for each temperature). Dark blue spheres represent the constrained ice seed. (A) At 225 K the ice front shows curvature around the adsorbed sbwAFP. The middle shaded circle segment has a radius of 46 Å equal to the theoretical Gibbs–Thomson cylindrical curvature of TIP4P water at 5 K supercooling. (B) At 230 K the ice retreats, but is stabilized directly beneath the bound sbwAFP. (C) At 232 K, (approximately 2°C above Tm), the bound sbwAFP remains in place stabilizing an ice protrusion, pinning back further ice retreat. The side profile views are made from two adjacent simulation periodic images.

https://doi.org/10.7554/eLife.05142.014

When the temperature was raised to 230 K, (the approximate system melting point Tm) the simulation showed a significant reduction of the ice fraction to an equilibrated value of 0.21 ± 0.02, (calculated over the 1100 to 2130 ns time frame) shown in Figure 8B. The sbwAFP remained bound to ice, stabilizing an ice protrusion directly beneath the AFP approximately 7 Å high relative to the main ice front.

In the final contiguous simulation the temperature was raised to 232 K, (approximately 2 K above the simulation melting point) and showed further reduction in the ice fraction which equilibrated to an average of 0.15 ± 0.01 (calculated over 2190 to 3160 ns time frame). The sbwAFP remained attached to the ice protrusion approximately 8.5 Å above the main ice plane which developed slight concave curvature with respect to the adsorbed protein (Figure 8C). A continuation of the simulation at 235 K resulted in the detachment of the sbwAFP from the ice surface within 10 ns allowing further melting of the ice to a residual fraction of 0.08 ± 0.01. Steady-state representations of the equilibrated ice fronts for each temperature were generated averaging water density over each simulation segment, by mapping water occupancy of greater than 85% using the ‘volmap’ module of the visualization software VMD (Humphrey et al., 1996). Extended simulations showing the induced curvature of the ice surface are shown in Videos 3, 4, 5. These simulations clearly demonstrate the Gibbs–Thomson effect showing how adsorbed protein induces ice-surface curvature and subsequent melting and freezing inhibition of ice. An important aspect of AFP activity that remains to be resolved involves the question of the strength of binding of AFP to ice. Previous estimations of the binding energy of AFPs to ice have varied greatly in the literature. Initial calculations from Cheng and Merz (1997) put the wfAFP binding energy of −157 kcal/mol while Dalal and Sönnichsen (2000) put the figure at −67 kcal/mol. Later estimates by Jorov et al. (2004) and Wierzbicki et al. (2007) both put the figure for wfAFP much lower at below −5 kcal/mol. Such low values of binding, however, imply an equilibrium with appreciable amounts of exchange of the AFP from the ice surface and its surroundings. Contrary to this, experimental evidence suggests that the binding energy must be relatively high and virtually irreversible. Microfluidic experiments by Celik et al. (2013) show that ice growth remains inhibited despite the exchange of the solution surrounding AFP-inhibited ice crystals with AFP-free buffer, suggesting there is little or no exchange of AFPs from the ice surface once bound. This finding is further strengthened by recent field observations of Cziko et al. where it was found Antarctic notethenoid fishes internally accumulated ice does not melt as expected over summer warming periods (Cziko et al., 2014). Our own observations in this study suggest relatively strong binding of a magnitude similar to that proposed by Dalal and Sönnichsen (2000), with the sbwAFP never detaching once bound to ice, (apart from deliberate melting events), even over the course of the extended millisecond range simulations at above melting temperatures. We intend to address the magnitude of the binding energy in our future work.

Video 3
Ice-bound sbwAFP simulated for 1000 nanoseconds at 232 K, showing stabilization of the ice front by the presence of the bound sbwAFP.
https://doi.org/10.7554/eLife.05142.015
Video 4
Simulation of sbwAFP binding to the prism ice front at 225 K.

The adsorption of the sbwAFP to the prism face causes curvature to the ice consistent with the Gibbs–Thomson effect.

https://doi.org/10.7554/eLife.05142.016
Video 5
Continuation of Video 4 showing the ice curvature with periodic boundary conditions applied.

The second half of the video highlights the positions of the bound water molecules and shows how they migrate to the ice binding surface from the beginning of the simulation.

https://doi.org/10.7554/eLife.05142.017

AFP activity and ordered water

Our simulations have shown that sbwAFP activity depends on binding and arranging water in a similar periodicity to ice which in turn bind to the prism face of ice, and that the ordered water plays a crucial role in determining ice plane specificity and integrating into the ice lattice. This conclusion is supported by previous work by Nutt and Smith (2008) whose molecular dynamics simulations of sbwAFP in solution indicated water was more ice-like on its binding face at low temperatures and speculated that this preconfigured water was important for initial recognition and binding events. Recent work by Midya and Bandyopadhyay (2014) simulating insect AFP from Tenebrio molitor with a similar ice binding surface to sbwAFP in water have also found equivalent arrangements of bound water molecules emphasizing the importance of the hydroxyl groups of the ice-binding threonine residues. Contrary to this, as mentioned previously, it is noteworthy that the putative ice-binding threonine residues of the alpha helical wfAFP can be substituted to structurally similar but hydrophobic valine residue with little effect on activity (Haymet et al., 1999). This shows wfAFP ice binding does not depend on the threonine hydroxyl group, even suggesting a hydrophobic mode of interaction. Though tempting to make comparisons with sbwAFP threonine residues, we should be reminded that these proteins are completely different structures with significantly different ice-binding properties. It will be important to characterize both the wfAFP/ice interaction as well as analogous sbwAFP threonine to valine mutations before valid interpretations can be made.

Many of the solved AFP crystal structures have been found to contain ordered waters which have been postulated to integrate directly into the ice lattice (Liou et al., 2000). In a recently determined structure for a highly active bacterial AFP from Marinomonas primiryensis (mpAFP), bound waters exhibited an extensive ‘clathrate-like’ network that matched ice lattice dimensions (Garnham et al., 2011). Interestingly, unlike most AFPs, mpAFP appears to bind to all ice plane orientations, rather than specific ones. This may result from extended ice-like ordering of water from the ice binding surfaces interacting with ice in non-specific orientations, perhaps in a process similar to ice sintering (Kuroiwa, 1961).

As noted in a recent review by Nada and Furukawa (2012), previous AFP simulations have focused on modeling prearranged AFPs at the ice-water interface at equilibrium and have not yet tackled free AFP interacting with an actively growing ice interface which would reveal more mechanisms of the molecular process and be less biased to initial positioning. Our freely diffusing MD simulations of AFP spontaneously adsorbing and inhibiting an active ice surface model have demonstrated that simulations can be used to provide detailed insight into molecular processes like crystal inhibition and modification.

We have presented simulations of the sbwAFP spontaneously interacting with a growing ice surface and demonstrated key aspects of the AFP mechanisms; namely the process of molecular surface recognition and irreversible adsorption inhibition of the ice surface. We have shown that the Gibbs–Thomson formula agrees closely with the simulated inhibition of the growing ice front, and demonstrated irreversibility of binding by free energy analysis. This is one of the first computer simulations to provide a complete in silico demonstration of a protein's biological function.

Materials and methods

Molecular modelling

Molecular simulations were performed using NAMD2.9 (Wierzbicki et al., 1997) with CHARM22 (MacKerell et al., 1998) forcefield employing a TIP4P water model supported on both IBM BlueGene/P and BlueGene/Q architecture. The sbwAFP model was based on the pdb structure 1M8N (Leinala et al., 2002). Simulations were run with PBCs using the NPT ensemble at various temperatures (225, 228, 230, 232 and 235 K) and 1 bar pressure employing Langevin dynamics. The PBCs were constant in the XY dimensions. Long-range Coulomb forces were computed with the Particle Mesh Ewald method with a grid spacing of 1 Å. 2 fs timesteps were used with non-bonded interactions calculated every 2 fs and full electrostatics every 4 fs while hydrogens were constrained with the SHAKE algorithm. The cut-off distance was 12 Å with a switching distance of 10 Å and a pair-list distance of 14 Å. Pressure was controlled to 1 atmosphere using the Nosé-Hoover Langevin piston method employing a piston period of 100 fs and a piston decay of 50 fs. Trajectory frames were captured every 100 picoseconds. Ice fractions were determined by measuring the ratio of mobile to immobile water molecules. Water molecules were determined to be immobile, (i.e., ice), if the average oxygen position of water had moved less than 0.8 Å over three successive trajectory frames (i.e., 200 ps). The initial pdb model of the sbwAFP at the ice/water interface is included in the supplementary data (Supplementary file 1).

Calculating critical size of spherical ice embryo at 1 K supercooling

Using Equation 1 for the case of spherical geometry (Ag = 2),

RMwσLρi2T0(T0Tr),

where R is the radius, Mw the molecular weight of water (0.018 kg/mol), L the latent heat of melting of water ice (6.02 × 103 J/mol) (Vega et al., 2005), σ is the ice-water surface energy 29.1 × 10−3 J/m2 (Handel et al., 2008), ρi is the density of ice (917 kg/m3), T0 the normal melting temperature (273.15 K), and Tr is the melting temperature depressed by curvature (272.15 K). These values give R = 518 Å. A sphere of this radius has a volume of 5.82 × 108 Å3 and would contain 1.8 × 107 water molecules.

Using parameters derived from the TIP4P water model (Vega et al., 2005; Handel et al., 2008), the corresponding values are T0 of 230.5 K, σ of 23 × 10−3 J/m2, ρi of 944 kg/m3 and L of 4.4 × 103 J/mol. Using the cylindrical version of this Gibbs–Thomson equation (Ag = 1) and a freezing point depression of 5 K, we find the critical cylindrical radius to be 46 Å.

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
    Biochemistry of fish antifreeze proteins
    1. PL Davies
    2. CL Hew
    (1990)
    FASEB Journal 4:2460–2468.
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
    Direct calculation of solid-liquid interfacial free energy for molecular systems: TIP4P ice-water interface
    1. R Handel
    2. RL Davidchack
    3. J Anwar
    4. A Bruhko
    (2008)
    Physical Review Letters 100:036104.
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
    A peek at ice binding by antifreeze proteins
    1. KA Sharp
    (2011)
    Proceedings of the National Academy of Sciences of USA 108:7281–7282.
    https://doi.org/10.1073/pnas.1104618108
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53

Decision letter

  1. Gerhard Hummer
    Reviewing Editor; The Max Planck Institute of Biophysics, Germany

eLife posts the editorial decision letter and author response on a selection of the published articles (subject to the approval of the authors). An edited version of the letter sent to the authors after peer review is shown, indicating the substantive concerns or comments; minor concerns are not usually shown. Reviewers have the opportunity to discuss the decision before the letter is sent (see review process). Similarly, the author response typically shows only responses to the major concerns raised by the reviewers.

Thank you for sending your work entitled “The biological function of an insect antifreeze protein simulated by molecular dynamics” for consideration at eLife. Your article has been favorably evaluated by John Kuriyan (Senior editor), a Reviewing editor, and two reviewers.

The Reviewing editor and the reviewers discussed their comments before we reached this decision, and the Reviewing editor has assembled the following comments to help you prepare a revised submission.

This paper explores the molecular mechanism of an antifreeze protein with the help of molecular dynamics simulations. The simulations show that spruce budworm antifreeze protein (sbwAFP) stalls ice growth at the ice prismatic plane in a manner consistent with the Gibbs-Thomson effect. While the study is clearly interesting, a number of serious criticisms have been raised, in particular concerning novelty (in light of earlier simulation studies of ice-growth inhibition by the same protein), the choice of the water model (with a TIP4P having a very low freezing point), the assumed ice growth mode, and the calculation of binding free energies. The points listed below would have to be convincingly addressed for the paper to be publishable in eLife.

Major comments:

1) Novelty. The present study is closely related to the earlier work of Furukawa and Nada. In Nada, 2001, the same protein was studied, and very similar conclusions were reached. In particular, the ice melting point depression due to the Gibbs-Thomson effect has already been reported by Nada and Furukawa (Nada, 2008 and Nada, 2001). Specific questions to address are: What distinguishes the present study from these earlier simulations? Is the claim that “to the best of our knowledge [the Gibbs-Thomson phenomena] has not been observed directly” justified? What are the differences to these earlier studies? Does the broader temperature range, the consideration of mutant proteins, and the longer simulations lead to substantially different conclusions or a deeper understanding? It has to be clear that the present paper constitutes a major scientific advance.

2) Water model. The simulations are performed with the TIP4P model. However, TIP4P gives a only a poor rendering of the water phase diagram. In particular, the ice melting point is far too low. As a result, the simulations were carried out at ∼230 K. However, the diffusion of water and protein at that temperature is very slow, resulting in poor sampling. There are several models in which the reproduction of the phase diagram is much better than in the TIP4P model, e.g., the six-site model and TIP5P/ice model. Could the authors comment about the dependence of the simulation results on the water potential model? It would be reassuring if simulations with a different water model gave similar results.

3) Ice growth mode. The simulations are set up assuming a step growth mode on the ice prismatic plane. Is this realistic? At 1 atm pressure, ice crystals grown from water do not have prismatic plane facets. This means that prismatic planes should not grow by a step-growth mode. Prismatic plane facets appear on grown ice crystal shapes only when pressure is extremely high (M. Maruyama, J Cryst. Growth 275, 2005, 598). The formation of hexagonal column ice crystals in the presence of sbwAFP suggests that sbwAFP binds to the prismatic plane and inhibit ice growth on the plane, but it does not necessarily mean that the prismatic plane growth occurs by a step growth mode. Would the results be affected if a different, more typical growth mechanism were considered?

4) Estimate of binding free energy and enthalpy. The calculated binding free energies are questionable because of several issues: (i) The variance in the calculated work values is very large (several kT), which makes the application of perturbative methods problematic; (ii) Entropic effects (e.g., free rotation of the molecule in the dissociated state) are probably not sampled correctly; (iii) It is not clear how the work is calculated; (iv) It is not clear if the free energy contains any contributions from the harmonic restraint that would have to be subtracted for a proper binding free energy (see, e.g., Hummer, Szabo, Acc. Chem. Res. 2005); (v) The enthalpy estimate in terms of hydrogen bond numbers does not seem to include hydrogen bonds of free water with ice. To get a more controlled estimate of the binding free energy, it would be important to include also the reverse process of binding and use (variants of) the Crooks fluctuation theorem. In its present form, the binding free energy calculations seem to be more confusing than helpful. They may be more suitable for a separate technical publication, in which details can be explored more carefully.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled “The biological function of an insect antifreeze protein simulated by molecular dynamics” for further consideration at eLife. Your revised article has been evaluated by John Kuriyan (Senior editor) and a member of the Board of Reviewing Editors, and the two reviewers who read your original submission. The manuscript has been improved but there are major remaining issues that need to be addressed before acceptance, as outlined below in the summary statement and review prepared by the Reviewing Editor.

It should be possible for you to address the issues raised by this round of review by editing your manuscript in response to the issues raised, without the need for any further calculations. Please pay particular attention to the fact that the Reviewing Editor asks for the elimination of the data and discussion concerning the calculation of ice-binding free energies. Our proceeding further with a revised manuscript is contingent on this being done, because the editor and the reviewers have identified gaps in this analysis.

Summary statement from the Reviewing Editor:

Molecular dynamics simulations are used to study the trapping of an advancing ice-water interface by antifreeze proteins. The simulations are insightful and interesting, and shed light on an important but poorly understood process at the interface of biology and physics. In their rebuttal and revision, the authors have made a good effort to address the concerns raised in the previous reviews. However, there are important issues that have not been adequately addressed.

To make the paper acceptable for publication in eLife will require (1) a considerably more nuanced discussion of preceding work, and (2) the elimination of the data and discussion of the calculation of ice-binding free energies. Detailed reasons are given below. In addition, (3) the discussion of the semi-quantitative agreement with the Gibbs-Thomson formalism, as given in equation 1 should better reflect evident uncertainties in this comparison. Also, (4) a clearer justification of the choice of ice-growth model should be included in the paper, not just in the rebuttal. In the following, I discuss these four points in detail.

First, I find the discussion of earlier studies inadequate. I am particularly troubled by statements such as “previous AFP simulations have fallen short of complete, convincing molecular descriptions of an AFP-ice interaction” and “recent advances in computational science now make complete simulations feasible and we report them here.” What do the authors mean by their simulations being “complete”? Is there nothing to be done in this field after this study? This is certainly not the case, with many questions remaining open, e.g., concerning the water model with its freezing point of ∼240 K, the quantification of the phenomenological theory in equation 1, the mode of ice growth, or the calculation of a binding free energy. A fairer discussion of earlier results and a more modest presentation of the new results are in place. Such a discussion should explicitly mention the work of Nada and Furukawa on the use of MD simulations to study ice growth kinetics around AFPs, as influenced by the Gibbs-Thomson effect.

Second, there are major issues with the calculation of the protein-ice binding free energy. Quoting Reviewer #2, in revising their manuscript, the authors have substantially increased the details of their simulation methodology and have now included a much-expanded discussion on the calculation of the binding free energy. Unfortunately these details now render these calculations and the conclusions derivable from them somewhat questionable. These calculations could be presented separately in a more technical journal to do justice to the many challenges involved in such demanding calculations.

The instantaneous switching of the control parameters (either the spring constant or harmonic minima) is expected to yield reliable free energy estimates only when such changes result in dissipative work values on the order of kT (where the exponential of the work is well approximated by perturbation theory). The accumulated work values in Figure 9 place this study far outside of this regime, and it is likely the poor sampling of exponentially rare work values that lead to systematic errors and ultimately the disagreement between various protocol estimates of the binding energy discussed in the appendix. Contrary to the authors' claims in their revision, simulating both forward and backward protocols to converge the binding energy would be hugely beneficial precisely because of the irreversibility of binding that the authors allude to. Only in instances of microscopically reversible dynamics do the two procedures yield redundant information. The authors may, e.g., consider Dellago and Hummer, 2013, for guidelines concerning the application of these equations (Dellago, C. and Hummer, G., “Computing Equilibrium Free Energies Using Non-Equilibrium Molecular Dynamics”, Entropy 16.1 (2013): 41-61).

The discussion of the binding free energy in terms of the number of broken hydrogen bonds is also too vague to be useful. As pointed out by Reviewer #2, the number listed for the hydrogen bond enthalpy seems to be judiciously picked to yield agreement with the calculated binding energy (the reference pointed to in the manuscript lists a slightly different value) and moreover hydrogen bonding energies are known to be rapidly varying functions of bonding geometry. The authors' recent correction of this estimate to reflect an activation energy rather than the binding energy presupposes that the binding is barrierless for the equivalence of these two quantities to be thermodynamically meaningful and that it is necessary to break hydrogen bonds along the reaction coordinate for this estimate to be kinetically meaningful. From the data provided and cited, it is not clear that either of those assumptions is valid.

I want to add to the comments of Reviewer #2 that even the definition of the binding free energy is unclear. Binding equilibria from a bulk phase to a surface have dimensional contributions that require a choice of standard state. Furthermore, does the calculated free energy, fully or partially, include the entropic contribution from a substantial increase in rotational freedom? Probably not, since the time to sample full rotation in the dissociated state has been too small. On the technical side, I wonder how the work was calculated (i.e., which formula was used and implemented), in particular for the cases in which the spring constant was changed with time.

Third, Reviewer #2 points out that the near quantitative agreement of the simulation results with the radius of curvature calculated from Equation 1 seems fortuitous. Is there a reason to believe that the surface tension or enthalpy of fusion for the water model used are the same as in experiment? Can experimental parameters taken from a website (http://www.its.caltech.edu/∼atomic/snowcrystals/ice/ice.htm) be used to explain the observations of the simulations? Such a validation procedure seems questionable. There are also technical issues with the application of Equation 1 to a molecular simulation system. As discussed previously (Findenegg et al., “Freezing and melting of water confined in silica nanopores”, ChemPhysChem 9.18 (2008): 2651-2659), how exactly the interface and the radius of curvature are calculated is not unique. The implications on the calculated value of the radius would need to be discussed for this quantitative comparison to be meaningful.

Finally, Reviewer #1 asks for further clarification in the paper concerning the ice-growth mode. In the rebuttal, the authors mention: “Our preliminary work using ‘perfect’ crystals did not grow well, extending the ice front very slowly if at all. We realised this was because one had to effectively wait for a 2 dimensional nucleation event on the prism plane to occur and reach a critical size before the next layer of ice would grow.” The authors should include this explanation in the paper to justify their assumption of a step-growth mode on the prismatic plane.

https://doi.org/10.7554/eLife.05142.019

Author response

1) Novelty. The present study is closely related to the earlier work of Furukawa and Nada. In Nada, 2001, the same protein was studied, and very similar conclusions were reached. In particular, the ice melting point depression due to the Gibbs-Thomson effect has already been reported by Nada and Furukawa (Nada, 2008 and Nada, 2001). Specific questions to address are: What distinguishes the present study from these earlier simulations?

Our studies are the first to show reproducible ice binding of an AFP as it migrates freely from the aqueous phase to the ice surface and the first to show total ice inhibition consistent with the Gibbs-Thomson mechanism while also revealing surprising details of the water arrangements at the interface layer. It is also one of the first, by including mutant control proteins, to demonstrate how crucial the local water binding arrangement is for AFP activity.

Nada and Furukawa performed superficially similar simulations in 2011, however their study has serious shortcomings. Firstly, their two AFPs initial conformations was pre-determined using a docking routine on an idealized ice prism surface without liquid water. As their simulations were also quite short, (only 2 simulations of 6 nanoseconds), the pre-positioning of AFP at the ice interface compounds the biases by limiting the sampling of interfacial water.

Secondly, Nada and Furukawa have not clearly demonstrated total ice inhibition, but only a slowing of the ice growth rate (shown in their Figure 5). Their ice growth never plateaued and their simulations are too short to allow proper sampling or protein diffusion at the ice interface. It was perhaps premature on their part to attribute the curvature to the Gibbs-Thomson effect; indeed, one could argue that over their simulation time-scale a similar sized non-AFP protein may have caused an equivalent ice-distortion, but this was not tested with a control.

Nada and Furukawa based their conclusions on limited sampling of 2 simulations of only 6 nanoseconds length. Our study has the fortune to be based on sampling of a long 3000 nanosecond simulation plus 32 wildtype binding simulations (each of 150 to 250 nanoseconds in length) and 16 simulations of the mutant form (each of 250 nanoseconds in length) totalling over 12 microseconds for the binding/inhibition part of our experiments alone.

We have modified our text to emphasise these advances over earlier studies (Introduction):

“Impeded by computational constraints, previous AFP simulations have fallen short of complete, convincing molecular description of an AFP-ice interaction. […] Recent advances in computational science now make complete simulations feasible and we report them here.”

Is the claim that “to the best of our knowledge [the Gibbs-Thomson phenomena] has not been observed directly” justified?

Upon wider reading, we did find references that refer to observing the Gibbs-Thomson effect in nanomaterials research dealing with semiconductors. Consequently, strictly speaking the original wording would be invalid. We further consulted with an expert in the field, Dr. Charlie Knight, who agreed that in context of AFPs, resolving the ice/water interface at the scale of a few tens of nanometers is nearly impossible, even with atomic force microscopy.

We have thus modified the subject text to be more specific (subsection titled “Adsorption inhibition”):

“Though the Gibbs-Thomson effect is widely accepted as the basis of AFP activity, resolving these nanoscale features at the ice/water interface has so far not been attained experimentally.”

2) Water model. The simulations are performed with the TIP4P model. However, TIP4P gives a only a poor rendering of the water phase diagram. In particular, the ice melting point is far too low. As a result, the simulations were carried out at ∼230 K. However, the diffusion of water and protein at that temperature is very slow, resulting in poor sampling. There are several models in which the reproduction of the phase diagram is much better than in the TIP4P model, e.g., the six-site model and TIP5P/ice model. Could the authors comment about the dependence of the simulation results on the water potential model? It would be reassuring if simulations with a different water model gave similar results.

The TIP4P water model does indeed have a well-known discrepancy with regards to its freezing temperature. We would have preferred to use a more sophisticated water model such as the TIP5P/ice, however the only water models that were supported by the NAMD software we used were TIP3P and TIP4P. Other MD programs such as Amber can support more water models, but at significantly slower performance. NAMD was the only software we had that would scale efficiently on our BlueGene clusters to get the required timescales. Also the TIP4P model has been used in a number of notable ice simulation papers, namely Mochizuki et al., (Nature 498.7454 (2013): 350-354) and Matsumoto et al., (Nature 416.6879 (2002): 409-413).

Also, we do agree the diffusion of water is slow at ∼230K but with our longer simulations we have seen the fluctuations in the ice front; observing the front grow and retreat. The fact that we show the molecules are able to reversibly arrange themselves into a crystal structure shows sufficient sampling the relevant regions of phase space. It is certainly a good point to be wary of, so we have added an additional note in the text (in the subsection “Adsorption inhibition”):

“The extent of the inhibited ice front fluctuated, advancing and retreating in the range of 5 to 8 Å showing the reversible arrangement of water molecules to the crystal ice lattice, (this is an important indication of reaching equilibration and adequate sampling given the reduced water diffusion at these low temperatures).”

Overall, we are confident that the TIP4P model can be used to give good representative behaviour of AFP activity, as long as sufficient sampling is done to compensate for slower water diffusion.

3) Ice growth mode. The simulations are set up assuming a step growth mode on the ice prismatic plane. Is this realistic?

Ice crystals are typically riddled with imperfections as shown by various vacuum etching studies, which serve as growth mechanisms such as screw dislocations (Sinha, N. K. (1977). Dislocations in ice as revealed by etching. Philosophical Magazine, 36(6), 1385-1404).

Our preliminary work using “perfect” crystals did not grow well, extending the ice front very slowly if at all. We realised this was because one had to effectively wait for a 2 dimensional nucleation event on the prism plane to occur and reach a critical size before the next layer of ice would grow. The tilt of the seed crystal and resultant step growth ensures a constant ice growth mechanism with similar effect to a screw dislocation. We believe this to be not only a realistic model of ice growth, but immensely practical to introduce as it greatly reduces the computational time necessary to achieve a certain rate of growth.

At 1 atm pressure, ice crystals grown from water do not have prismatic plane facets. This means that prismatic planes should not grow by a step-growth mode. Prismatic plane facets appear on grown ice crystal shapes only when pressure is extremely high (M. Maruyama, J Cryst. Growth 275, 2005, 598).

Our model ice face is only centred on the prism ice face for practical purposes. It would be incorrect to take this as representative of a prism plane of macroscopic ice crystal.

The formation of hexagonal column ice crystals in the presence of sbwAFP suggests that sbwAFP binds to the prismatic plane and inhibit ice growth on the plane, but it does not necessarily mean that the prismatic plane growth occurs by a step growth mode. Would the results be affected if a different, more typical growth mechanism were considered?

The step growth model we have developed here is a convenient way to test the AFP activity in the presence of a growing ice front. With a much larger system, (and much more computational overhead!), one could create other crystal growth initiators conditions, such as a screw dislocation for similar layer-by-layer growth mechanism. We would expect similar results, as steps originating from screw dislocations are still essentially step growth. It would certainly be interesting to test AFP activity around such a dislocation in simulation, but that would currently require yet another massive increase in computational capability.

4) Estimate of binding free energy and enthalpy. The calculated binding free energies are questionable because of several issues: (i) The variance in the calculated work values is very large (several kT), which makes the application of perturbative methods problematic.

If the error bounds seem large, it is because we have been most careful and conservative in our estimation of the free energy. The estimates are as derived and analysed by Gore et al.. We are aware of no superior method of error estimation. We do not see that larger error and bias estimates raise the possibility of sources of error not captured by the estimates used, though we have carefully examined the literature on the use of the Jarzynski theorem.

Further to these error estimates, we have examined the effect of varying the speed of the protein, and checked that calculation of ΔG for the motion of the protein through pure water is appropriately close to 0. We used two independent sets of trajectories and identified their strengths and weaknesses. Our paper depends only on the demonstration that the binding is irreversible. We submit that this is clearly shown by our results. We have now also expanded and clarified the discussion of the Jarzynski method.

(ii) Entropic effects (e.g., free rotation of the molecule in the dissociated state) are probably not sampled correctly.

In the dissociated state, the protein is surrounded by water on all sides, to a distance of at least the correlation length in pure water. This state on average is symmetrical, and so sampling more orientations would not alter the thermodynamics of this state. More broadly, and not only for the purposes of the Jarzynski calculation, equilibration is important. As discussed in the main paper, we see stabilisation of the ice front and protein in each structure examined.

(iii) It is not clear how the work is calculated.

We have greatly expanded and clarified the report of this calculation. There are too many changes for these to be set out in full here, but for example the relevant Methods section now opens with an outline of its components, and how they fit together.

Other salient improvements are in the Methods section, where the role of the error and bias estimates is discussed, as well as the use of tests of the bias method.

(iv) It is not clear if the free energy contains any contributions from the harmonic restraint that would have to be subtracted for a proper binding free energy (see, e.g., Hummer, Szabo, Acc. Chem. Res. 2005).

Fortunately the free energy of a harmonic oscillator is exactly and easily calculable, and the force constants we use correspond to free energies smaller than the random errors.

The work by Hummer and Szabo pertains to the application of Jarzynski's theorem to single molecule force spectroscopy, rather than to a simulation. The important correction that they consider concerns the measurement of force as a function of time, when the work W in Jarzynski's theorem is the integral of force with respect to distance. This is a subtle issue, but one that we circumvent in a simulation because the exact value of the work W is available for each trajectory.

(v) The enthalpy estimate in terms of hydrogen bond numbers does not seem to include hydrogen bonds of free water with ice.

Estimating the hydrogen bonds of free water with ice is problematic as it is difficult to define which molecule is water and which is ice. We have clarified our original estimate and argument to consider the initial step of removing the protein from the ice as requiring simultaneously breaking the 12 hydrogen bonds we observe made from the ordered water to the ice and relating this to the activation free energy for the removal of the protein from the ice.

To get a more controlled estimate of the binding free energy, it would be important to include also the reverse process of binding and use (variants of) the Crooks fluctuation theorem.

The individual trajectories are done deterministically and so doing more trajectories is the same as studying also the reverse process. To use the Crooks theorem on top of Jarzynski would be redundant.

In its present form, the binding free energy calculations seem to be more confusing than helpful. They may be more suitable for a separate technical publication, in which details can be explored more carefully.

We have made improvements on the clarity of the binding free energy calculations as detailed above.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

To make the paper acceptable for publication in eLife will require (1) a considerably more nuanced discussion of preceding work, and (2) the elimination of the data and discussion of the calculation of ice-binding free energies. Detailed reasons are given below. In addition, (3) the discussion of the semi-quantitative agreement with the Gibbs-Thomson formalism, as given in equation 1 should better reflect evident uncertainties in this comparison. Also, (4) a clearer justification of the choice of ice-growth model should be included in the paper, not just in the rebuttal. In the following, I discuss these four points in detail.

First, I find the discussion of earlier studies inadequate. I am particularly troubled by statements such as “previous AFP simulations have fallen short of complete, convincing molecular descriptions of an AFP-ice interaction” and “recent advances in computational science now make complete simulations feasible and we report them here.” What do the authors mean by their simulations being “complete”? Is there nothing to be done in this field after this study? This is certainly not the case, with many questions remaining open, e.g., concerning the water model with its freezing point of ∼240 K, the quantification of the phenomenological theory in equation 1, the mode of ice growth, or the calculation of a binding free energy. A fairer discussion of earlier results and a more modest presentation of the new results are in place. Such a discussion should explicitly mention the work of Nada and Furukawa on the use of MD simulations to study ice growth kinetics around AFPs, as influenced by the Gibbs-Thomson effect.

We have greatly elaborated on prior work as suggested by the reviewers, including citing an additional 8 references and adding three paragraphs in relation to the historical context of computational simulation of AFPs, with particular mention of the work from Nada and Furukawa.

Second, there are major issues with the calculation of the protein-ice binding free energy. Quoting Reviewer #2, in revising their manuscript, the authors have substantially increased the details of their simulation methodology and have now included a much-expanded discussion on the calculation of the binding free energy. Unfortunately these details now render these calculations and the conclusions derivable from them somewhat questionable. These calculations could be presented separately in a more technical journal to do justice to the many challenges involved in such demanding calculations.

The instantaneous switching of the control parameters (either the spring constant or harmonic minima) is expected to yield reliable free energy estimates only when such changes result in dissipative work values on the order of kT (where the exponential of the work is well approximated by perturbation theory). The accumulated work values in Figure 9 place this study far outside of this regime, and it is likely the poor sampling of exponentially rare work values that lead to systematic errors and ultimately the disagreement between various protocol estimates of the binding energy discussed in the appendix. Contrary to the authors' claims in their revision, simulating both forward and backward protocols to converge the binding energy would be hugely beneficial precisely because of the irreversibility of binding that the authors allude to. Only in instances of microscopically reversible dynamics do the two procedures yield redundant information. The authors may, e.g., consider Dellago and Hummer, 2013, for guidelines concerning the application of these equations (Dellago, C. and Hummer, G., “Computing Equilibrium Free Energies Using Non-Equilibrium Molecular Dynamics”, Entropy 16.1 (2013): 41-61).

The discussion of the binding free energy in terms of the number of broken hydrogen bonds is also too vague to be useful. As pointed out by Reviewer #2, the number listed for the hydrogen bond enthalpy seems to be judiciously picked to yield agreement with the calculated binding energy (the reference pointed to in the manuscript lists a slightly different value) and moreover hydrogen bonding energies are known to be rapidly varying functions of bonding geometry. The authors' recent correction of this estimate to reflect an activation energy rather than the binding energy presupposes that the binding is barrierless for the equivalence of these two quantities to be thermodynamically meaningful and that it is necessary to break hydrogen bonds along the reaction coordinate for this estimate to be kinetically meaningful. From the data provided and cited, it is not clear that either of those assumptions is valid.

I want to add to the comments of Reviewer #2 that even the definition of the binding free energy is unclear. Binding equilibria from a bulk phase to a surface have dimensional contributions that require a choice of standard state. Furthermore, does the calculated free energy, fully or partially, include the entropic contribution from a substantial increase in rotational freedom? Probably not, since the time to sample full rotation in the dissociated state has been too small. On the technical side, I wonder how the work was calculated (i.e., which formula was used and implemented), in particular for the cases in which the spring constant was changed with time.

We have removed all results and discussions regarding the free-energy calculations. We agree that this aspect of the work is difficult and perhaps better suited in a more technical journal where we can explore technical and theoretical details fully.

We greatly appreciate the detailed critique of the free-energy methods we originally employed. Since we have removed the controversial components raised by the reviewers, specifically instantaneous switching of control parameters and binding energy in terms of broken hydrogen bonds from this manuscript, we shall address them and binding energy definitions in subsequent work with and not elaborate them here.

Third, Reviewer #2 points out that the near quantitative agreement of the simulation results with the radius of curvature calculated from Equation 1 seems fortuitous. Is there a reason to believe that the surface tension or enthalpy of fusion for the water model used are the same as in experiment? Can experimental parameters taken from a website (http://www.its.caltech.edu/∼atomic/snowcrystals/ice/ice.htm) be used to explain the observations of the simulations? Such a validation procedure seems questionable. There are also technical issues with the application of Equation 1 to a molecular simulation system. As discussed previously (Findenegg et al., “Freezing and melting of water confined in silica nanopores”, ChemPhysChem 9.18 (2008): 2651-2659), how exactly the interface and the radius of curvature are calculated is not unique. The implications on the calculated value of the radius would need to be discussed for this quantitative comparison to be meaningful.

We have recalculated the radius of curvature of ice using parameters derived from the TIP4P model from literature values of Vega et al., and Handel et al., leading to a smaller radius value of 46 Å compared to our original 59 Å figure. We agree the quantitative comparison is problematic and have dropped the direct reference in favour for a qualitative comparison and a brief mention of the difficulties in measuring curvature from a simulation. We have replaced the shaded 60 Å radius circular segment in Figure 8, with that of a 46 Å radius segment equal to that of the theoretical radius of curvature value derived from TIP4P model parameters.

Finally, Reviewer #1 asks for further clarification in the paper concerning the ice-growth mode. In the rebuttal, the authors mention: “Our preliminary work using ‘perfect’ crystals did not grow well, extending the ice front very slowly if at all. We realised this was because one had to effectively wait for a 2 dimensional nucleation event on the prism plane to occur and reach a critical size before the next layer of ice would grow.” The authors should include this explanation in the paper to justify their assumption of a step-growth mode on the prismatic plane.

We have further clarified the step growth mode of the ice model by mentioning our earlier experience of slow growth of ‘perfect’ crystals and referencing work of Frank and Hobbs in relation to step propagation as a mechanism crystal growth.

https://doi.org/10.7554/eLife.05142.020

Article and author information

Author details

  1. Michael J Kuiper

    Victorian Life Sciences Computation Initiative, The University of Melbourne, Carlton, Australia
    Contribution
    MJK, Conception and design, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article, Contributed unpublished essential data or reagents
    For correspondence
    mkuiper@unimelb.edu.au
    Competing interests
    The authors declare that no competing interests exist.
  2. Craig J Morton

    ACRF Rational Drug Discovery Centre, St Vincent's Institute of Medical Research, Fitzroy, Australia
    Contribution
    CJM, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  3. Sneha E Abraham

    School of Chemistry, The University of Melbourne, Melbourne, Australia
    Contribution
    SEA, Conception and design, Acquisition of data, Analysis and interpretation of data
    Competing interests
    The authors declare that no competing interests exist.
  4. Angus Gray-Weale

    School of Chemistry, The University of Melbourne, Melbourne, Australia
    Contribution
    AG-W, Conception and design, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.

Funding

Australian Research Council (ARC) (DP110103388)

  • Sneha E Abraham
  • Angus Gray-Weale

University of Melbourne (Victorian Life Sciences Computation Initiative grant number VR0064)

  • Michael J Kuiper

Victorian Government

  • Michael J Kuiper

Victorian Partnership for Advanced Computing

  • Michael J Kuiper

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This research was supported by a Victorian Life Sciences Computation Initiative (VLSCI) grant number VR0064 on its Peak Computing Facility at the University of Melbourne, an initiative of the Victorian Government, Australia, and by a grant from the Australian Research Council (DP110103388) and computational resources from the Victorian Partnership for Advanced Computing. We thank Dr Christina Hall for discussion and manuscript preparation.

Reviewing Editor

  1. Gerhard Hummer, The Max Planck Institute of Biophysics, Germany

Publication history

  1. Received: October 12, 2014
  2. Accepted: May 6, 2015
  3. Accepted Manuscript published: May 7, 2015 (version 1)
  4. Version of Record published: May 26, 2015 (version 2)

Copyright

© 2015, Kuiper et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 4,960
    Page views
  • 829
    Downloads
  • 35
    Citations

Article citation count generated by polling the highest count across the following sources: Scopus, Crossref, PubMed Central.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Biochemistry and Chemical Biology
    2. Structural Biology and Molecular Biophysics
    Ching-Ju Tsai et al.
    Research Article Updated
    1. Biochemistry and Chemical Biology
    2. Structural Biology and Molecular Biophysics
    Vivekanandan Shanmuganathan et al.
    Research Article Updated