Structures, dynamics, and hydrogen-bond interactions of antifreeze proteins in TIP4P/Ice water and their dependence on force fields

Tenebrio molitor antifreeze protein (TmAFP) was simulated with growing ice-water interfaces at a realistic melting temperature using TIP4P/Ice water model. To test compatibility of protein force fields (FFs) with TIP4P/Ice water, CHARMM, AMBER, and OPLS FFs were applied. CHARMM and AMBER FFs predict more β-sheet structure and lower diffusivity of TmAFP at the ice-water interface than does OPLS FF, indicating that β-sheet structure is important for the TmAFP-interface binding and antifreeze activity. In particular, CHARMM FF more clearly distinguishes the strengths of hydrogen bonds in the ice-binding and non-ice-binding sites of TmAFP than do other FFs, in agreement with experiments, implying that CHARMM FF can be a reasonable choice to simulate proteins with TIP4P/Ice water. Simulations of mutated TmAFPs show that for the same density of Thr residues, continuous arrangement of Thr with the distance of 0.4~0.6 nm induces the higher extent of antifreeze activity than does intermittent arrangement of Thr with larger distances. These findings suggest the choice of CHARMM FF for AFP-TIP4P/Ice simulations and help explain the relationship between Thr-residue arrangement and antifreeze activity.


Introduction
Antifreeze proteins (AFPs), which consist of polypeptides with various sizes and structures, are found in arctic or antarctic organisms such as bacteria [1], fungi [2], plants [3], insects [4], and fish [5][6][7] that can survive at temperatures even lower than -30˚C. AFPs have shown great potential for industrial applications such as cryopreservation [8], food processing [9], and hydrate inhibition [10], since they bind to specific planes of the growing ice crystal and inhibit ice growth, showing a noncolligative property [11][12][13], which can overcome the limitation of synthetic polymer-based antifreezes that typically require high concentrations. To increase the antifreeze efficiency of AFP, the structure of AFP and its interaction with specific ice plane need to be understood, which has motivated many experimental and theoretical studies.
The DeVries group pioneered experimental and theoretical studies on the interactions between AFPs and ice-water interfaces. They found that AFPs adsorb to ice surfaces and inhibit ice growth through a noncolligative process, which was explained by a mechanism of a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 the melting-point depression on the curved ice-water interface, called the "Gibbs-Thomson" (or Kelvin) effect [7,14,15]. Graether observed that the mutation of the ice-binding site of AFP significantly reduces antifreeze activity [4], and Davies and coworkers showed the effects of AFP structure and ice-plane type on the AFP-ice binding [16][17][18], which were interpreted by hydrogen-bond interactions. Meister et al. observed different hydration layers on the icebinding sites of different AFPs [19,20], and Olijve et al. showed that the efficiency of AFP cannot be determined by any single factor, but by a combination of ice-plane type, adsorption rate of AFP, surface coverage, and kinetics of ice nucleation [21]. Liu et al. showed that ice nucleation is depressed by the non-ice-binding face (bulky hydrophobic and charged groups) of AFP but promoted by the ice-binding face (hydroxyl groups) of AFP [22]. While these experiments have shown the dependence on the AFP structure and ice-plane type and proposed the mechanism of the Gibbs-Thomson effect, the results are not always easy to interpret at the level of individual molecules or specific hydrogen-bond interactions between AFPs and water molecules because of the limited resolution of experimental techniques.
To support and complement experimental observations, molecular dynamics (MD) simulations have been performed. McDonald et al. simulated a winter flounder AFP (wfAFP) with pyramidal planes, showing the adsorption of wfAFP to ice surface because of the interactions between Threonine (Thr) residues and oxygen atoms in the ice lattice [23]. Cheng and Merz also found that aspartic acid (Asp), asparagine (Asn), and Thr residues of wfAFP form hydrogen bonds with ice surfaces, to an extent dependent on the distance and ordering of those residues [24]. Dalal et al. observed that the extent of hydrogen bond does not significantly differ in the ice-water interface and bulk water, implying that the protein-ice binding is mainly due to van der Waals forces [25]. Wierzbicki et al. calculated free energies and found that wfAFPs do not directly bind to the ice surface but instead interact with the ice-water interface in the semisolidsemiliquid phase [26]. To more realistically simulate the AFP-ice binding, Nada and Furukawa simulated wfAFP with a growing ice-water interface composed of six-site water molecules at the actual melting point of water and showed that wfAFPs become partially surrounded by the grown ice, leading to a decrease in the velocity of ice growth [27], which supports the Gibbs-Thomson effect. They also observed the dependence of the wfAFP-ice binding on the ice plane, in agreement with experiments [11]. Recently, Kuiper et al. simulated a spruce budworm AFP (sbwAFP) with prism planes, showing that sbwAFP binds to ordered water molecules near the prism face, which can be disrupted by replacing Thr with Leu [28]. Chakraborty and Jana found that water molecules near the ice surface are anchored to hydroxyl groups of Thr [29,30]. These simulations have captured hydrogen-bond interactions between AFPs and ordered water molecules near the curved ice-water interface [31,32], although the dependence of antifreeze activity on Thr arrangement, which might be also influenced by different force fields (FFs) in simulation studies, has not been systematically studied through computation.
As a further step toward understanding the effect of Thr arrangement on antifreeze activity, we here report all-atom MD simulations of Tenebrio molitor AFP (TmAFP) with ice-growing surfaces at a realistic melting temperature using TIP4P/Ice water model [33]. We first test compatibility of TIP4P/Ice water with different protein FFs such as CHARMM [34,35], AMR-BER [36,37], and OPLS [38,39]. Secondary structures, hydrogen-bond lifetimes, and antifreeze activities of TmAFPs are compared to suggest the FF that is most compatible with TIP4P/Ice water. Then, some of Thr residues of TmAFP are mutated and simulated, showing the effect of Thr-residue arrangement on antifreeze activity, which is rationalized by considering hydrogen-bond interactions between Thr residues and ordered water at the ice-water interface. We will show that these results help explain experimental observations regarding the AFP-ice binding strength and antifreeze activity, to an extent that depends on the arrangement of Thr residue.

Methods
All simulations and analyses were performed using the GROMACS5. , and then the mutated TmAFPs were generated by replacing some of Thr residues with Leu using the Swiss-Pdb viewer (Table 1) [44]. For water and ice, TIP4P/Ice model [33] was used because this model reproduces more realistic melting temperature of~260 K than do other water models such as TIP3P (146 K) [45], TIP4P (232 K) [45], TIP4P/Ew (245 K) [46], SPC (190 K) [47], and SPC/E (215K) [48]. Note that TIP5P [45] and six-site [49] water models also show realistic melting temperatures, but they are not considered in this work because they are computationally very expensive. For the ice-crystal plane, Drori et al.'s recent experiments showed that TmAFP adsorb to both prism and basal planes, but the adsorption rate is much faster at the prism plane than at the basal plane [50]. Thus, the prism plane was used in this work, and its coordinate was obtained from Kuiper et al.'s work [28]. To test the ability of different FFs to predict the structure of TmAFPs and their interactions with TIP4P/Ice water, TmAFP was modeled with different FFs such as CHARMM [34,35], AMBER [36,37], and OPLS FFs [38, 39].
A single TmAFP was placed above the ice-water interface with a z-directional distance of 1.8 nm between the protein center and the ice-water interface, and 2 counterions (Na + ) were added to establish electro-neutrality in a periodic box of size 5.1 × 12.6 × 8.6 nm 3 (Fig 1). A  T26L T38L T50L T62L  3   TmAFP-m2  CHARMM  T16L T38L T40L T62L T64L  3   TmAFP-m3  CHARMM  T16L T26L T28L T50L T52L T62L T64L  3 https://doi.org/10.1371/journal.pone.0198887.t001 temperature of 260 K and a pressure of 1 bar were maintained by using a velocity-rescale thermostat [51] and Parrinello-Rahman barostat [52] in the NPT ensemble. A real space cutoff of 1.2 nm was applied for Lennard-Jones and electrostatic forces with the inclusion of particle mesh Ewald for long-range electrostatics [53]. The LINCS algorithm was used to constrain the bond lengths [54,55]. To obtain more samples, three or five simulations were performed for each system for 120 ns with a time step of 2 fs on computer facilities supported by the National Institute of Supercomputing and Networking/Korea Institute of Science and Technology Information with supercomputing resources including technical support (KSC-2017-C3-061).
The last-30 ns trajectories were used for analyses.

Effect of different force fields on the interaction between TmAFP and TIP4P/Ice water
Since the TIP4P/Ice water model reproduces the realistic melting temperature of water while still being computationally less expensive than TIP5P or six-site water, this water model has been often used for simulations of AFPs, although the compatibility of TIP4P/Ice water with different protein FFs has not yet been studied. To test the ability of different FFs to predict the strength of hydrogen bonds between AFPs and TIP4P/Ice water, simulations of TmAFP with TIP4P/Ice water were performed using different FFs such as CHARMM, AMBER, and OPLS. TmAFP and the ice-water interface, which is more significantly observed in CHARMM and AMBER FFs than in OPLS FF. This implies that CHARMM and AMBER FFs can more accurately reproduce the binding of TmAFP with the ice interface and its effect on the suppression of ice growth than does OPLS FF.
To quantify this effect of TmAFP, the z-directional movement of TmAFP was analyzed by calculating the z-directional distance between the center of mass (COM) of TmAFP and its initial COM as a function of time. In Fig 3, although there are some exceptions, simulations with CHARMM and AMBER FFs show that most TmAFPs migrate toward the disordered-water region within a distance of 0.5 nm from its initial position and then do not move much after 90 ns, while simulations with OPLS FF show that some TmAFPs (systems 1 and 2) drastically move until the end, indicating that those TmAFPs do not bind to the ice-water interface and thus keep pushed toward the disordered-water region. Recall from Fig 2 that CHARMM and AMBER FFs show the order water adjacent to the ice-binding site of TmAFP, but OPLS FF does not. These results imply that CHARMM and AMBER FFs can more accurately reproduce the binding of TmAFP with the ice interface and its effect on the suppression of ice growth than does OPLS FF.
To confirm these different extents of the TmAFP-ice binding, diffusion coefficients of TmAFPs were calculated by obtaining the slopes of the mean-square displacements (MSD) of the COM of TmAFP for the initial (0~30 ns) and final (90~120 ns) simulation times. Note that the MSD usually needs to be corrected for finite size effects using Yeh and Hummer's analytic formula for cubic boxes [57], but for systems in non-cubic boxes, the analytic equation of Yeh and Hummer is not necessarily applicable. Also, the size effects should not significantly affect the comparison of diffusivities within the same box geometry, and thus the finite size effect is not corrected here. In Fig 4, diffusivities significantly decrease at the end of simulations with CHARMM and AMBER FFs, apparently because TmAFPs bind to the ice-water interface and inhibit ice growth, as observed in Figs 2 and 3. For simulations with OPLS FF, the diffusivity does not decrease much, showing that TmAFPs keep migrating toward the disordered-water region. This, combined with Figs 2 and 3, indicates that CHARMM and AMBER FFs can more effectively capture the antifreeze activity of TmAFP than does OPLS FF.
To understand this dependence of the antifreeze effect on FFs, secondary structures of TmAFP were calculated using the DSSP program [58]. Table 2 shows that CHARMM and AMBER FFs predict more β-sheet and less random-coil structures than does OPLS FF. Note that the percentage of β-structure is lower than the experimental value of 40% [59], since solvent conditions differ in experiments and our simulations, which precludes a quantitative comparison of secondary structures. Despite this, CHARMM and AMBER FFs predict more β-sheet structure than does OPLS FF and thus more favorably compares with experimental results. These, combined with Figs 2-4, indicate that TmAFP must retain β-sheet structure to suppress ice growth, which can be captured more effectively by CHARMM and AMBER FFs than by OPLS FF.

Strength of hydrogen-bond interactions: comparison of different force fields
Experimental results have shown that the ice-binding site of AFP consists of many Thr residues and thus forms more hydrogen bonds with the surrounding water than does the non-ice-binding site of AFP [22], implying that the strengths of those hydrogen bonds may also differ. To compare the strengths of hydrogen bonds on the ice-binding and non-ice-binding sites and test the ability of different FFs to predict this, we calculated average lifetimes of hydrogen bonds between Thr residues and water molecules, which is the inverse of the rate constant for hydrogen-bond kinetics (breaking) described by the autocorrelation function of the hydrogen-bonding existence functions (either 0 or 1) averaged over all hydrogen bonds [60,61]. Here, we

Fig 3. Distances between the center of mass (COM) of TmAFP and its initial COM in the z-direction as a function of time.
Since five simulations were performed for each force field, five values are represented with different colors, where black, red, green, blue, and yellow lines respectively correspond to systems 1, 2, 3, 4, and 5 in Fig 2. https://doi.org/10.1371/journal.pone.0198887.g003 assume that a hydrogen-bonding interaction exists when the donor-acceptor distance is less than 0.35 nm and the angle of the hydrogen-donor-acceptor triplet is less than 30˚ [62]. Analyses were performed using different distance and angle criterions, including 0.3 nm-25˚and 0.35 nm-25˚, which showed that these stricter criteria produce similar qualitative trends. Thus, the 0.35 nm-30˚criterion was used for our hydrogen-bonding analysis of all the simulations. 11 Thr residues, which are regularly distributed along two lines on the bottom side of TmAFP, are considered to be an ice-binding site, while the remaining 6 residues are considered a non-icebinding site (Fig 1). In Fig 5, simulations with CHARMM FF show that average hydrogen-bond lifetimes are 1.54~2.63 and 0.37~1.83 ns, respectively, for the ice-binding and non-ice-binding sites of TmAFP, indicating the formation of more persistent hydrogen bonds on the ice-binding site than on the non-ice-binding site, in agreement with experiments [22]. The similar trend was also observed for other FFs, but hydrogen-bond lifetimes of the ice-binding and non-icebinding sites are more clearly distinguished in CHARMM FF than in other FFs. Also, recall that β-sheet structure is more predominantly observed in CHARMM FF than in other FFs (Table 2). These results indicate that when TIP4P/Ice water is used for AFP-ice simulations, CHARMM FF can be a reasonable choice to predict the experimentally observed β-sheet structure of AFPs and their hydrogen-bond interactions with the ice-water interface.

Effect of the arrangement of Thr residues on antifreeze activity
Experimentally, Graether et al. showed that the mutation (Thr to Leu) of the ice-binding site of AFP significantly reduces antifreeze activity [4], which was supported by Kuiper et al.'s simulations [28], indicating that Thr-residue arrangement may influence antifreeze activity. To confirm this, TmAFP was mutated by replacing some of Thr residues with Leu in various patterns and simulated with TIP4P/Ice water using CHARMM FF. Mutations of TmAFP are described in Fig 6 and Table 1, showing that the ice-binding site of TmAFP-m1 has Thr residues along only a single line, and the ice-binding sites of TmAFP-m2 and TmAFP-m3 have intermittent arrangement of Thr residues along two lines. In Fig 6, final snapshots of simulations show that TmAFP-m1 and TmAFP-m2 induce the formation of curved ice-water interfaces, while TmAFP-m3 does not. Fig 6 also shows that the water molecules adjacent to the ice-binding site of proteins are more ordered in TmAFP-m1 than in other mutated ones, indicating the stronger binding of TmAFP-m1 with the ice-water interface. To quantify the extent of antifreeze activity, diffusion coefficients of TmAFP and mutated TmAFPs were compared. Fig 7 shows that diffusivities of TmAFP and TmAFP-m1 become very close at the end of simulations, while diffusivity of TmAFP-m3 at the end of simulation is almost same as that of TmAFP at the beginning of simulation. This indicates that TmAFP-m1 strongly binds to the ice-water interface and thus does not diffuse much, similar to TmAFP, while TmAFP-m3 does not bind to the ice-water interface and thus keeps pushed toward the water region, leading to a significant decrease in antifreeze activity. For TmAFP-m2, diffusivity is higher than TmAFP-m1 and lower than TmAFP-m3. These indicate that the extent of the TmAFP-ice binding strength and antifreeze activity is higher in the order of TmAFP and TmAFP-m1 > TmAFP-m2 > TmAFP-m3, indicating that antifreeze activity can be reduced by replacing Thr with Leu, as observed in previous experiments [4] and simulations [28].

The relationship between Thr position and the TmAFP-ice binding
As discussed above, antifreeze activity of TmAFP is influenced by the arrangement of Thr residues. Since Thr residues of the ice-binding site form strong hydrogen bonds with the ice-water interface and induce antifreeze activity, the strength of the AFP-ice binding may depend on the arrangement of Thr residues, although their relationship has not been systematically studied.
To resolve this, we calculated minimum distances between Thr residues. In Fig 8, TmAFP and TmAFP-m1 show that distances between Thr residues are less than 0.9 nm, while TmAFP-m2 and TmAFP-m3 show larger distances. Fig 9 plots the distributions of minimum distances between Thr residues for all simulated systems, showing that minimum distances for TmAFP and TmAFP-m1 are mostly distributed in the range of 0.4~0.6 nm, while those for TmAFP-m2 and TmAFP-m3 are in the broad range of 0.5~2 nm. Recall from Fig 7 that the extent of antifreeze activity is higher in the order of TmAFP % TmAFP-m1 > TmAFP-m2 > TmAFP-m3, indicating that antifreeze activity increases with decreasing the distance between Thr residues. Also, note that the distance between ordered water molecules of ice surface is~0.45 nm, which is in the minimum-distance range of 0.4~0.6 nm for TmAFP and TmAFP-m1, implying that Thr residues need to be densely spaced for the formation of strong hydrogen bonds. To check this, we calculated average lifetimes of hydrogen bonds between water molecules and Thr residues on the ice-binding sites of mutated TmAFPs. Fig 10 shows [28]. In particular, our simulations show that although the same number of Thr residues are distributed on the ice-binding sites of TmAFP-m1 and TmAFP-m2, the extent of antifreeze activity is higher for TmAFP-m1 than for TmAFP-m2. This indicates that the densely distributed arrangement of Thr with the distance of 0.4~0.6 nm more effectively induces antifreeze activity than does the intermittent arrangement of Thr with larger distances, implying the dependence on the distance between Thr residues. These simulations show the effect of different FFs on the secondary structure and antifreeze activity of TmAFP, and the strength of the TmAFP-ice surface binding, suggesting a choice of CHARMM FF especially for simulations of AFP and TIP4P/Ice water at a realistic melting temperature. In particular, simulations show the higher extent of antifreeze activity when Thr residues are more densely distributed with the distance of 0.4~0.6 nm, indicating the effect of Thr-residue arrangement on hydrogen-bond interactions and antifreeze activity. Note that different types of AFPs have specific binding affinities with various ice-plane types, so our results should be considered for specific AFP and ice plane with caution. To more precisely determine the antifreeze activity of different TmAFPs, it would obviously be interesting to calculate the free energy of binding at the AFP-ice interface, which we hope to report on elsewhere.

Conclusions
We performed MD simulations of TmAFP with growing ice-water interfaces using TIP4P/Ice water that can reproduce a realistic melting temperature while still being computationally less expensive than TIP5P or six-site water. Since the interactions between proteins and TIP4P/Ice water have not been verified, we first tested compatibility of TIP4P/Ice water with different FFs such as CHARMM, AMBER, and OPLS. TmAFP binds to the ice-water interface, induces curvature of ice surface, and suppresses ice growth, although their extents differ. CHARMM and AMBER FFs predict more β-sheet structure and lower diffusivity of TmAFP at the icewater interface than does OPLS FF, indicating the stronger TmAFP-ice surface binding and antifreeze activity. In particular, analyses of hydrogen-bond lifetime show that CHARMM FF more clearly distinguishes the strengths of hydrogen bonds in the ice-binding and non-icebinding sites of TmAFP than do other FFs, which favorably compares with experiments, implying that CHARMM FF can be reasonably chosen to reproduce experimental observations on the secondary structure of AFP and its interaction with TIP4P/Ice water.
To understand the relationship between Thr-residue arrangement and antifreeze activity of AFP, mutated TmAFPs were generated by replacing some of Thr residues with Leu and simulated with TIP4P/Ice water using CHARMM FF. Simulations of mutated TmAFP show that when densities of Thr are same, densely distributed arrangement of the distance of 0.4~0.6 nm more effectively induces antifreeze activity than does intermittent arrangement of Thr with larger distances. This indicates that distances between Thr residues must be small enough to maintain ordered structure of ice-surface water that can form persistent hydrogen bonds with TmAFP. These findings suggest that CHARMM FF can be reasonably chosen for simulations of AFP with TIP4P/Ice water at a realistic melting temperature, as well as help explain experimental observations regarding the difference of the binding strength in the ice-binding and non-ice-binding sites of TmAFP, and the dependence of antifreeze activity on the distance between Thr residues and hydrogen-bond interactions with ice surfaces.