Comparison of popular force fields for molecular modeling of proteins applied to ice binding of the tenebrio molitor antifreeze protein

Ice Binding Proteins (IBPs) are a class of proteins that affect the melting and freezing temperatures of ice. While substantial research has been conducted to understand the properties of IBPs, experimental and molecular dynamics simulations have not satisfactorily explained the mechanism of the irreversible binding of protein to ice. Parallel rows of threonine residues play an important role in the activity of hyperactive proteins. Past experimental and computational studies show a peculiar stability of the protein structure and rigidity of the two rows of threonine. The hydroxyl groups on the two threonine rows are known to closely match the dimensions of hexagonal ice, but what about this geometry allows a protein to bind irreversibly to ice? To investigate this mechanism, we use molecular dynamics simulations to examine the Tenebrio molitor antifreeze protein (TmAFP), which is generally representative of the family of hyperactive antifreeze proteins. In particular, when designing simulations of these proteins, force field parameters that govern the molecular system must be examined closely to ensure an accurate physical representation. We survey four common protein system force fields, two each from the CHARMM and AMBER families. Trials using AMBER force fields (ff12SB and ff14SB) as well as one newer CHARMM force field (CHARMM36) are not able to accurately mimic the ice binding face, and simulations with these force fields fail to dock with ice. An older CHARMM force field, CHARMM22, is able to reproduce threonine rigidity. We measured hydroxyl distances between neighboring threonine residues and the threonine dihedral angles with each force field and demonstrate that the experimentally observed rigidity is mimicked only in the Tenebrio molitor antifreeze protein (TmAFP), which is generally representative of the family of hyperactive antifreeze proteins. In particular, when designing simulations of these proteins, force field parameters that govern the molecular system must be examined closely to ensure an accurate physical representation. We survey four common protein system force fields, two each from the CHARMM and AMBER families. Trials using AMBER force fields (ff12SB and ff14SB) as well as one newer CHARMM force field (CHARMM36) are not able to accurately mimic the ice binding face, and simulations with these force fields fail to dock with ice. An older CHARMM force field, CHARMM22, is able to reproduce threonine rigidity. We measured hydroxyl distances field and demonstrate that the experimentally observed rigidity is mimicked only in the simulations using CHARMM22. Further, a simulation with CHARMM22 shows successful binding of the TmAFP with ice, while an identical simulation with AMBER ff12SB does not bind. We propose that the threonine rigidity is essential to ice binding, and that the hardness of the CHARMM22 force field offers sufficient rigidity. Future simulations must therefore use a force field able to reproduce this rigidity, as in CHARMM22.

COMPARISON OF POPULAR FORCE FIELDS FOR MOLECULAR MODELING OF PROTEINS APPLIED TO ICE BINDING OF THE TENEBRIO MOLITOR ANTIFREEZE PROTEIN by Nathan Odendahl A Senior Honors Thesis Submitted to the Faculty of The University of Utah In Partial Fulfillment of the Requirements for the Honors Degree in Bachelor of Arts In Chemistry Approved: ______________________________ Valeria Molinero Thesis Faculty Supervisor _____________________________ Cynthia Burrows Chair, Department of Chemistry _______________________________ Thomas Richmond Honors Faculty Advisor _____________________________ Sylvia D. Torti, PhD Dean, Honors College April 2016 ABSTRACT Ice Binding Proteins (IBPs) are a class of proteins that affect the melting and freezing temperatures of ice. While substantial research has been conducted to understand the properties of IBPs, experimental and molecular dynamics simulations have not satisfactorily explained the mechanism of the irreversible binding of protein to ice. Parallel rows of threonine residues play an important role in the activity of hyperactive proteins. Past experimental and computational studies show a peculiar stability of the protein structure and rigidity of the two rows of threonine. The hydroxyl groups on the two threonine rows are known to closely match the dimensions of hexagonal ice, but what about this geometry allows a protein to bind irreversibly to ice? To investigate this mechanism, we use molecular dynamics simulations to examine the Tenebrio molitor antifreeze protein (TmAFP), which is generally representative of the family of hyperactive antifreeze proteins. In particular, when designing simulations of these proteins, force field parameters that govern the molecular system must be examined closely to ensure an accurate physical representation. We survey four common protein system force fields, two each from the CHARMM and AMBER families. Trials using AMBER force fields (ff12SB and ff14SB) as well as one newer CHARMM force field (CHARMM36) are not able to accurately mimic the ice binding face, and simulations with these force fields fail to dock with ice. An older CHARMM force field, CHARMM22, is able to reproduce threonine rigidity. We measured hydroxyl distances between neighboring threonine residues and the threonine dihedral angles with each force field and demonstrate that the experimentally observed rigidity is mimicked only in the ii simulations using CHARMM22. Further, a simulation with CHARMM22 shows successful binding of the TmAFP with ice, while an identical simulation with AMBER ff12SB does not bind. We propose that the threonine rigidity is essential to ice binding, and that the hardness of the CHARMM22 force field offers sufficient rigidity. Future simulations must therefore use a force field able to reproduce this rigidity, as in CHARMM22. iii TABLE OF CONTENTS ABSTRACT ii INTRODUCTION 1 ICE BINDING PROTEINS 1 STRUCTURAL PATTERS IN ICE-BINDING PROTEINS 2 CONTROVERSY IN THE MECHANISM OF BINDING 5 MODELS AND METHODS 9 THE TENEBRIO MOLITOR ANTIFREEZE PROTEIN 9 COMPUTATIONAL METHODS 10 RESULTS AND DISCUSSION 13 CONCLUSIONS AND OUTLOOK 20 REFERENCES 22 iv 1 INTRODUCTION ICE BINDING PROTEINS A huge variety of novel and powerful protein functions exist, despite all proteins sharing only about two-dozen amino acid building blocks. Because proteins are too small to observe with a microscope, researchers are forced to use a variety of techniques to examine protein shape and function. Computational simulations can model a solvated protein on timescales of up to a few microseconds, and so offer a unique avenue of investigation. Ice Binding Proteins (IBPs) bind to the ice crystal lattice.9 Discovered about thirty years ago, IBPs are surprisingly good at manipulating the freezing and melting temperature of ice. Two main classes of IBP exist: antifreeze proteins (AFPs) and ice nucleating proteins (INPs). AFPs lower the freezing temperature of water by preventing ice from growing, while INPs encourage ice to form more readily.4 AFPs are also small, ~50x10 Å, while INPs can be up to 30 times larger.10 Antifreeze proteins have been shown to inhibit ice growth according to the GibbsThomson effect (or the Kelvin effect).17 The Gibbs-Thomson effect (or the Kelvin effect) observes that ice growth is slower when the ice front is curved. The mostly flat ice front becomes pocked where AFPs bind and prevent further ice growth. The “button mattress” result causes enough curvature between binding sites to inhibit ice growth by a few degrees.17 The Gibbs-Thomson effect caused by AFPs has been observed by ETM, where measured grooves were consistent with assumed AFP binding.12 Subsequent simulations and experiments have confirmed the accuracy of the Gibbs-Thomson effect and the related equations.6,17,26 2 The equilibrium melting point and temperature of freezing of water are not necessarily identical. At equilibrium, ice melts and freezes at 273 K. The difference between the melting point and inhibited freezing point is called the Thermal Hysteresis (TH). The macro-scale mechanism of antifreeze proteins is well understood, but the direct protein-ice interaction is still unclear. The structure of ice binding proteins is an important consideration in understanding their mechanism. STRUCTURAL PATTERNS IN ICE BINDING PROTEINS Since the first moderate AFP was discovered, a variety of AFPs structures have been characterized. The earliest AFPs found are produced in fish, and have since been classified as “moderate AFPS,” with Type I, Type II, and Type III subclasses. Characterized as small proteins with either globular or alpha helix structure, fish AFPs have been studied in depth, partially because they were discovered early, and partially because they tend to be small in size. The hyperactive AFP family, found mostly in insects, contrasts with moderate AFPs, found primarily in fish. Hyperactive AFPs show antifreeze activity at a lower concentration and are stronger inhibitors than moderate AFPs.4,48 The Gibbs-Thomson effect depresses the freezing point is depressed by up to 2 °C, for moderate AFPs, or over 5 °C, for hyperactive AFPs.14,48 The hyperactive class, while more efficient, also takes longer to effect a change, and requires more time to bind as temperatures drop.5,7 Although a variety of IBP structures have been observed, a few traits are common to them all. The ice-active surface area is called the ice-binding site (IBS) or ice-binding 3 face (IBF). All IBPs are highly organized, and nearly all have a very regular, repeated sequence along the ice-binding face.10,18,29 Moderate AFPs, which are less effective, have either a globular or alpha helical structure throughout, while hyperactive AFPs are composed of highly regular beta sheets.4,10 The organized beta sheets of the hyperactive family are also characteristic of ice nucleating proteins.11,29 By forming dimers or complexes, some IBPs increase the size of their ice binding face and therefore the favorability of binding. This phenomenon has been theorized,10 demonstrated,31 and observed in vivo.38 The ice binding faces of hyperactive AFPs are unusually flat and hydrophobic, and usually contain some hydroxyls that are able to interact with the solvent waters.4,10 The non-ice binding faces (nIBF) of AFPs show less organization and repetition, and have been observed to interrupt the order of water to discourage ice growth.22,26 Repetition allows the formation of a helical structure with regularly arranged residues to populate the ice binding face. Commonly, hyperactive AFPs contain two rows of organized threonine residues along the ice binding face that match the structure of ice and influence the dynamics of water.10,22,23 Threonine side chains along the IBS are exceptionally rigid and show lower fluctuations than non-IBS threonine.3 Waters along the threonine channels have been shown in experiments and simulations to slow when near the organized ice binding face.3,22,27 Figure 1 is taken directly from a 2014 review by Peter L. Davies and demonstrates the variety of known AFP structures.4 The proteins A-D are the moderate AFPs found in fish and characterized by alpha helical or globular structure. Protein E is a plant AFP. J and K are two microbial AFPs. Proteins F-I are insect AFPs, or hyperactive 4 proteins, characterized by beta sheets and a large, flat ice-binding site. Two other main AFP families have been observed. Antifreeze glycoproteins (AFGP) connect globular proteins to highly organized sugar chains that closely mimic ordered ice structure. The snow flea AFP (sfAFP, Figure 1 [I]) has a unique polyproline II Figure 1 – Antifreeze protein structures. (A) Maxi AFP; (B-D) Type I-III AFPs; (E) Plant AFP; (F-H) Hyperactive insect AFPs; (I) sfAFP; (J-K) microorganism AFPs. Figure from Davies.4 helical structure that is structurally different than other known AFPs, but which still shows high order matching the dimensions of ice.34,35 Ice nucleating proteins (INPs) perform the function opposite AFPs by encouraging ice nucleation at temperatures well above the homogeneous freezing temperature of water.4,20,24 The main structural difference between AFPs and INPs is the larger INP size.10 Kobashigawa et al. showed that an INP cut into smaller pieces exhibits antifreeze activity.16 Similarly, when a group of AFPs are immobilized onto a surface, Charpentier et al. demonstrated ice nucleation ability.2 These two experiments in conjunction show that beyond structural similarity, the two classes of IBP function according to similar principles: irreversible binding, the Gibbs-Thomson effect to produce antifreeze activity, antifreeze activity when proteins are free in solution, and heterogeneous ice nucleation when the protein is sufficiently large. 5 While IBPs are interesting from an academic perspective, there are important applications as well. INPs have been shown to play a role in atmospherics and precipitation.1,13 Pseudomonas syringae, perhaps the most well known producer of an INP, is a bacterium that causes significant frost damage on farm crops.13 Snomax® is a commercialized form of an INP that increases the yield of snow machines. Ice cream crystallization can be manipulated by AFPs, resulting in a smoother ice cream product. Understanding the mechanisms of ice growth and manipulation of ice inhibition offers many potential applications. A better understanding of fundamentals of ice is essential in cryogenics, in organ transplant transportation, and in producing frost resistant surfaces. A variety of organisms, ranging from fish and insects to bacteria, independently developed AFPs to reduce death by freezing.4 CONTROVERSY IN THE MECHANISM OF BINDING Despite 30 years of research, controversy still remains on the mechanism of ice binding. Researchers have put forth a number of theories to explain the surprisingly good binding ability of IBPs. Because of the close match the hydroxyls of threonine residues on the IBS make with the lattice structure of ice, the threonine residues were theorized to effectively replace waters within ice, forming direct hydrogen bonds with the ice waters.3,4,27,36 Mutation studies that replace essential threonine hydroxyls with the similarly shaped but non-polar valine on moderate AFPs show only slightly decreased ice binding activity, while mutation to the much differently shaped lysine on hyperactive AFPs destroys most ice binding activity.17,45 6 Irreversible AFP binding to ice has been shown experimentally.7,28 Irreversibility of the binding event remains an essential component in the mechanism of protein-ice binding. If proteins regularly unattached and reattached to the growing ice face, ice advancement would not be inhibited according to the Gibbs-Thomson effect because curvature would not be retained between attachment events. After researchers realized that the ice binding face of IBPs are mostly hydrophobic, an entropy driven mechanism was proposed. It was proposed that waters are constrained by the ice binding face and have reduced entropy.22,23,27,37 By replacing the constrained waters with the already constrained ice lattice, the system would experience a gain in entropy. This entropy advantage seems to play a role in IBPs’ ability to bind.3,4,8,14,17,23,26 A quantitative value of the entropy driven energetics has still not been calculated. Computational techniques offer a complimentary method of analysis to that of experimentation. A molecular dynamics (MD) simulation can offer detailed spatial information and statistics, but is limited by accuracy and time length. As molecular dynamics simulations grow in capability, simulations including protein, ice, and sufficient solvent play an increasing role in research projects. Simulations that investigate the structure and mechanism of ice binding are challenging to interpret. Some published results seem to show that the protein hydroxyls do not bind directly into the ice crystal lattice. A transition layer, called “ice-like water” by some researchers, may connect the protein to the ice surface.30,34 The ice-like water idea, however, has been criticized as an oversimplification, as not representative of all 7 proteins, or as flatly wrong.19 The water between the protein and ice may actually be less ice-like than the bulk water, with less tetrahedral character.8 The “anchored clathrate” mechanism is another explanation put forward by experts trying to explain the structure of water between the protein and ice.9 The proposal made by Garnham et al. suggests that direct hydrogen bonds are formed between the layer of water and threonine hydroxyls, but that a clathrate-like cage forms around the non-polar methyl groups. The hydrogen bonds to the hydroxyls serve as anchors, binding the protein to ice, while the cages around the methyl groups can connect with relatively little error into the ice lattice.4,9,10 The anchored clathrate model has only been demonstrated on one AFP so far. Perhaps a geometric argument is too complicated to classify or varies too widely between proteins to ever settle on a single theory. The energetics of entropy and the number of distinct hydrogen bonds may make a more compelling argument for the irreversibility of ice binding. Duboue-Dijon and Laage show (in liquid phase) that the hydrogen bonding along the ice-binding face of an antifreeze protein, CfAFP, is stronger than hydrogen bonding in bulk water due to steric confinement.8 The abnormally strong hydrogen bonding could explain why ice might prefer binding to the protein over bulk water. The energetics of this situation will doubtlessly be investigated in more depth. Ice crystal structure is mainly driven by tetrahedral interactions. Large-scale ice grows chiefly in a hexagonal prism, as shown in Figure 2, when at 1 bar and near 273 K. The two main face types are indicated by color. The two Figure 2 – Ice grows in a hexagonal prism. The basal face is indicated in blue and the prismatic faces in white. Ice growth is most favorable along the prismatic faces. 8 hexagonal basal faces are opposite one another, with the one visible basal face colored blue. The primary and secondary prismatic (or pyramidal) faces are indicated in white. Ice grows most favorably along the prismatic faces, with slower growth along the basal face. Numerous studies demonstrate that moderate and hyperactive AFPs show opposite binding preference, with superior binding to prismatic and basal faces, respectively. Initially determined by ice etching, the technique requires an ice crystal to be immersed into a solution of AFP and allowed to grow, then removed and left to sublime. Etching along the prismatic faces becomes apparent when moderate antifreeze proteins are used.46 Similar results were obtained using fluorescent labeling. Moderate and hyperactive proteins tagged with a fluorescent molecule were bound to an ice crystal, and the progress of growth imaged. A bipyramidal shape indicated growth along the basal plane only and inhibition along the prismatic faces, characteristic of moderate AFPs. Hyperactive proteins show fluorescence along all ice planes, but with a brighter signal at the basal faces, indicating strongest binding.4,6,7 Face selectivity is also predicted by the side chain spacing, which matches the ice lattice spacing of the prismatic planes, for moderate AFPs, or of the basal planes, for hyperactive AFPs. The hyperactive AFP side chains also form a less-ideal match with the prismatic ice faces, explaining the observed hyperactive binding along all faces of ice, but preferential binding to the basal face. The close interface and small size of the protein-ice connection complicates the investigation of IBP behavior. AFPs are very small, on the scale of 50x10 Å. Researchers have struggled to find effective probes of the ice-protein interface.4,22,25 9 Sum-frequency generation spectroscopy offers an experimental technique to investigate interface waters,19 but the method still does not effectively describe the binding mechanism. Mutation studies, where individual amino acids in the protein sequence are modified, help to probe the protein function. Combined with molecular dynamics simulations, key mutations have been shown to destroy ice-binding activity.17 Molecular dynamics simulations provide an invaluable tool for exploring the ice-protein dynamics, offering statistical answers at small cost.8,22,23 While experimental techniques continue to improve in detail and sophistication, computational techniques are recently powerful enough to offer useful information about complex and nano-scale systems. By utilizing molecular dynamics simulations, we hope to elucidate the mechanism of AFPs. We simulate the Tenebrio molitor antifreeze protein binding to ice. MODELS AND METHODS THE TENEBRIO MOLITOR ANTIFREEZE PROTEIN The Tenebrio molitor AFP (TmAFP) is a commonly studied and structurally representative hyperactive protein. Composed of beta sheets with two rows of parallel threonine residues, the TmAFP shows important similarities to various other hyperactive proteins (like the also common CfAFP). TmAFP also shows an exciting similarity to the Pseudomonas family of ice nucleating proteins.10 The parallel rows of threonine in TmAFP play an essential role in ice binding and have been investigated in the past.3,22,32 Water was observed to move more slowly within the channel created by the threonine side chains, which may indicate an important factor in the ice binding mechanism.22,23,32 10 The threonine residues organized into two rows play an essential role in the TmAFP. A threonine side chain makes a Y-shape with the hydroxyl and methyl groups, as shown in Figure 3. The residue attaches at its base to the protein backbone. A central carbon atom links the backbone with the two top branches: a nonpolar methyl group and a polar hydroxyl group. The distance between the neighboring Figure 3 – Threonine. Hydrogen (white), carbon (teal), oxygen (red), and nitrogen (blue). The nitrogen, center carbon, and rightmost carbon-oxygen pair are part of the protein backbone (-NH-CH-C=O-). The Yshape formed by the remaining large atoms (CH3-CH-OH) spins freely. threonine residues has been shown to match ice dimensions, and each threonine residue along the line is abnormally rigid.3 Figure 4 shows a wide view of a TmAFP monomer from the 1EZG crystal structure with only the ice binding face Figure 4 – TmAFP crystal structure (from 1EZG) showing the regular threonine residues on the ice binding face. Hydrogen atoms are not included in crystal structure. The rigidity of the threonine residues is essential to ice-binding activity, and must be reproduced in computational models. threonine residues explicitly displayed. COMPUTATIONAL METHODS Molecular dynamics (MD), the most common method for modeling proteins, requires a step-by-step force and velocity calculation of every atom at each femtosecond time interval. Many processors run in parallel to compute the forces acting on the tens of thousands of atoms contained in the box, and then evaluate the new position of each particle after the femtosecond time step. Multiple software packages are available that 11 carry out the necessary calculations. In this work, we use the AMBER14 software package.40 Powerful computers are needed to process the forces inside of the simulation box. MD simulations follow classical physics, which means that calculations ignore less impactful quantum mechanical effects. Atomistic simulations explicitly include interactions between every atom to preserve a greater level of accuracy, but do so at a high computational cost. Determining the interaction potentials between all pairs of atoms is no small task. By combining mathematical, numerical, and experimental models, researchers have made a number of “best fit” parameters to simulate physical systems. These parameters are referred to as “force fields” and contain information on the reaction of the most common elements when involved in different bonds. Force fields become increasingly complicated as more exact interactions are defined. The AMBER family of force fields are commonly used for simulating protein systems, including AMBER ff12SB (2012) and AMBER ff14SB (2014).43 CHARMM is a separate family of force fields for proteins which uses a slightly different premise, and includes CHARMM36 (2012)49 and CHARMM22 (1998)44. In AMBER and CHARMM, Van der Waals’ interactions are represented by a Leonard-Jones potential with electrostatic restraints and partial charges. A bond between a pair of atoms is harmonic springs according to the interaction potential. Angles and torsions are also controlled by simplified functions. We compare the performance of four force fields in modeling the TmAFP-waterice system: AMBER ff12SB43, AMBER ff14SB43, CHARMM3649, AND CHARMM2244. Because TmAFP is known to bind to ice, these force fields should 12 ideally be able to reproduce that binding behavior. All simulations were run simulations in the NPT ensemble at 1 bar pressure using under periodic box conditions. Temperatures ranged from 239 to 275 K using a Langevin thermostat. Atom motions were calculated with time steps of 2 fs. To preserve computation time, atom-atom interaction calculations were cutoff after 9 Å throughout each simulation. We visualized all resulting trajectories using Visual Molecular Dynamics (VMD).47 Our analysis of dihedral and hydroxyl distances used the in-built measurement features of VMD. We used the yellow mealworm (Tenebrio molitor) antifreeze protein structure [1EZG] found by Liou et al.39 as a starting point for each simulation. Since 1EZG is published as a dimer in the protein databank file (.pdb), we separated the first chain, added hydrogen atoms, and solvated the protein in TIP4P/2005 water.50 After minimization and equilibration, we concatenated TIP4P/2005 water as ice to the bottom of the system, which totaled 86,006 atoms. We ran all simulations using the AMBER pmemd package accelerated for GPU.40,41 Because the equilibrated protein was not in a favorable orientation relative to the ice lattice, we manually rotated the protein to align the ice binding face with the ice lattice and moved the protein to approximately 5 Å above the basal face, with the appropriate waters replaced. We then repeated minimization and equilibration of the system. We used the AMBER force fields ff14SB and ff12SB and the CHARMM force fields CHARMM36 and CHARMM22 within AMBER14 and CHAMBER (CHARMM enabled AMBER) respectively.40,42,43,49 We collected simulation information including dihedral angles and hydroxyl distances. We directly examined the parameter information 13 for CHARMM36 and ff14SB and determined those fields to carry the same flaw as ff12SB. Because of the observed flaws, we only ran full ice-water-protein simulations for ff12SB and CHARMM22. RESULTS AND DISCUSSION We investigated the capability of various force fields to mimic the ice binding behavior of the TmAFP and concluded that multiple newer force fields are too soft. Our initial simulations with the TmAFP did not bind to ice. The presented research explores and identifies the likely reason for poor or absent ice binding: threonine rigidity. Simulations run at 239 K and 245 K showed a few layers of ice growth over the 30 ns timescale. The protein, which we manually set close to the ideal binding configuration, was excluded from the ice lattice when parameterized with ff12SB, moving upwards away from the crystal face. After a few layers of ice growth, the protein became trapped within the ice lattice because of finite box size, but was not in an orientation representative of binding. Simulation with CHARMM22 confined the motions of threonine residues and led to successful binding. We know from a variety of experiments that the TmAFP binds readily to ice. Failure to bind to ice is unphysical and must derive from some flaw in the simulation parameters. The distance between neighboring threonine residues was measured as a probe of the stability of the ice binding face. Experiments have shown that the ice binding face quickly becomes rigid as temperature drops.3 14 Figure 5 shows the distance between two neighboring hydroxyls (residue index 51 and 63) when the protein is modeled in liquid water. The possible distances are confined to three general modes: -OH facing -OH (short distance, ~2.7 Å), parallel (medium distance, ~4.7 Å), and -OH facing -CH3 (long distance, ~7.2 Å). Dotted guidelines are included in Figure 5 to Figure 5 – Hydroxyl separation distance in Å over simulation time. THR51 and THR63 hydroxyls are representative of behavior over the ice-binding face. Variation is confined to three main modes, corresponding to shifting threonine orientations. Experiments show that only the parallel (4.5 Å) mode should manifest at low temperatures. emphasize the three modes. The parallel mode (~4.7 Å) matches ice geometry, while the facing modes (2.7 or 7.2 Å) do not. The trimodal behavior shown in Figure 5 is not conducive to ice binding because of the high variability and time spent in non-ice-like configuration. Because the threonine residues are close to one another, there is some communication between neighboring orientations, further disrupting the ice-like structure. Figure 6 shows a roughly parallel mode (~70° off of parallel), but the mode still does not mimic ice binding because of the anti-parallel configuration. The unphysical behavior of threonine lacking rigidity must be attributed to one or both of the two parameter sets that were used in the initial simulations. TIP4P/2005 is a widely used water model with generally good agreement with experimental values.50 The relatively large partial charges intrinsic to the model, however, could cause unpredictable interactions, effectively pushing or pulling the threonine residues out of alignment. The 15 AMBER force field ff12SB was the first parameter set applied to the protein interactions. As with all force fields used in molecular dynamics, the ff12SB estimates interactions to preserve computing power at the cost of some accuracy. To assess whether TIP4P/2005 the source of the unphysical IBS behavior, we simulated the gas phase (unsolvated) protein. The trimodal fluctuations were still observed in the gas phase, indicating the TIP4P/2005 is not the source of the unphysical fluctuations. We conclude, then, that the AMBER force field ff12SB contains an approximation that fails to reproduce the threonine rigidity that we theorize is essential to ice binding. A dihedral (or torsion) angle describes the orientation of four consecutive atoms within a molecule. Figure 6 shows two of three likely dihedrals in threonine, and Figure 7 shows the third orientation. Strictly, the dihedral angle describes the angle Figure 6 – Two threonine residues. These demonstrate two dihedral values (through the hydroxyl) of about 200° (right) and 80° (left). These two hydroxyls would be called parallel, corresponding to the 4.5 Å mode, but these could not bind to ice because of the opposite dihedrals. between two planes. Four consecutive atoms (we will call them A, B, C, and D) form three bonds (A-B, BC, and C-D), two angles (ABC, BCD), and one dihedral (ABCD). If we make two groups of three points, ABC and BCD, then each group is contained Figure 7 – Another threonine residue is displayed showing the third likely dihedral orientation, about 330°. in one unique plane. The unique planes formed by ABC and by BCD will necessarily 16 intersect. The smaller angle formed by the intersecting planes is defined as the dihedral angle. Force fields like ff12SB contain parameters that control dihedral motions. The threonine dihedral parameters in ff12SB are the main forces controlling the rotation of the threonine side chain. It seems likely, then, that the dihedrals do not Figure 8 – Histogram of dihedral angles of threonine residues along the ice binding face over the course of a simulation. This was produced with the ff12SB force field, but ff14SB would produce similar results. Three modes are seen centered on 80°, 190°, and 330°. None of these modes matches the 308° ideal seen in the protein crystal structure. appropriately model the rotation. Figure 8 shows the dihedral angles during the course of the simulation modeled with ff12SB. Crystal structure of TmAFP and experiments indicate that a dihedral of 308° is ideal for the rigid IBS. Three dihedrals, none matching the 308° ideal, dominate in the simulations with the ff12SB force field. A newer AMBER force field, ff14SB, was also tested. Similarly unphysical fluctuations were observed, indicating that the protein would not be able to bind to ice. CHARMM22,44 developed in 1998, is a very simplified force field that, nonetheless, accurately reproduces the proper threonine dihedral activity. Figure 11 compares the dihedral force parameters of CHARMM22 and ff12SB. Although CHARMM22 is relatively simple, the extra detail and softness in the newer force fields actually prevent the threonine residues from aligning properly. Figure 9 shows the dihedral values of the 17 CHARMM22 solvated protein. The dihedrals are confined almost entirely to 308°, in good agreement with the dihedral values from the xray crystallography obtained structure. Figure 10 shows the OH-OH distances, which again closely match the predicted values of the Figure 9 – Histogram of threonine dihedral angles along the ice binding face measured for the CHARMM22 simulation. The dihedrals are confined almost entirely around 308°, showing a strong match with the TmAFP crystal structure dihedral angles, and indicating that ice binding activity will occur. ice lattice. The deviations of the dihedral and OH-OH values are also greatly decreased in the CHARMM22 simulations, indicating higher rigidity. The dotted guidelines and the y-axis scale are preserved between Figures 5 and 10 to highlight the different OH-OH behavior between the two force fields. Only one mode can be observed in the CHARMM22 hydroxyl distances Figure 10 – A representative hydroxyl distance measurement from a solvated system when modeled with CHARMM22. The guidelines and y-axis match those shown in Figure 5 (for ff12SB). Less deviation, one mode, and high agreement with the TmAFP crystal structure all demonstrate that CHARMM22 preserves the needed rigidity. and dihedral. The two probes, OH-OH distance and dihedral of the threonine, both agree to show that CHARMM22 preserves threonine rigidity, while ff12SB does not. 18 CHARMM22 AMBER ff12SB Figure 11 – Summary of the dihedral forces and simulation results. The lines (not to scale with one another) show the dihedral force distribution along the 360° rotation of threonine, as in Figure 8. The bars are two normalized histograms of dihedral angles measured throughout two simulations, as in Figures 7 and 9. The CHARMM22 (red top line and red bars) remains localized around the key 308° of the crystal structure. The AMBER ff12SB (gray bottom line and gray bars) shows three modes around 80°, 190°, and 330°. Figure 11 offers a summary of the dihedral analysis. The two curves represent the dihedral force exhibited on the threonine side chain. Three roughly evenly spaced minima for each force field correspond to the three staggered conformations of threonine, which are the most energetically stable torsions. The curves on Figure 11 are not to scale with one another along the y-axis and are only offset for readability. The histogram of the CHARMM22 simulation matches closely with the CHARMM22 minimum found at the 308° staggered conformation. Importantly, the dihedral never oscillates between three minima offered by the CHARMM22 force field. This inability to jump between minima and therefore fluctuate in conformation we refer to as “hard” behavior. The dihedral from the AMBER ff12SB simulation conversely jumps between all three minima, exhibiting flexibility in what we call “soft” behavior. The softness of the 19 dihedral can therefore be linked directly to the fluctuating threonine residues and the protein’s inability to bind to ice. The simulations with CHARMM22 are able to bind to ice. Figure 12 shows the geometry of the TmAFP interfacing into the ice crystal basal plane. From a geometric analysis, TmAFP successfully binds to ice. The hydroxyl groups create hydrogen bonds with the ice lattice, while the ice-waters form half-cages around the hydrophobic methyl group. The threonine hydroxyls form a pentagon Top View (internal) or hexagon (external) that is able to interface directly into the ice lattice. From the binding geometry of the TmAFP, it appears that the threonine residues must stay planar and rigid in order to bind to ice. Side View Figure 12 – Top and side view of the central TmAFP threonine residue side chains binding to water. Atoms are represented as spheres, water oxygen (blue), threonine oxygen (red), and threonine methyl (green). Light blue cylinders represent hydrogen bonds, while the thin green cylinders represent covalent bonds within the threonine side chain. The protein (not pictured) connects above the threonine side chains, and the ice (partially pictured) connects below. The side view hexagonpentagon geometry transitions directly into hexagonal ice, indicating binding. These snapshots are taken from the simulation that used CHARMM22. The regularity of the ice crystal geometry must be matched by the threonine residues, as the residues show in the TmAFP crystal structure. Too much flexibility in the threonine 20 residues will break the peculiar ice-like spacing and prevent the protein from interfacing with ice. CONCLUSIONS AND OUTLOOK The results of this work highlight severe limitations of certain popular protein force fields in modeling ice-binding proteins. The flexibility of threonine residues found in AMBER ff12SB, AMBER ff14SB, and CHARMM36 arises from the equations governing dihedral motions and yields a fatal error in this particular protein system. Only the incidentally less specific parameterization in CHARMM22 results in an ice binding site that is not only rigid, but also causes the threonine hydroxyl distances to match those of the TmAFP crystal structure, allowing the protein to bind to ice in simulation. From the binding behavior of CHARMM22, we purport that rigidity of the threonine-lined channel is an essential attribute to the hyperactive antifreeze protein TmAFP. The present research naturally leads into further investigation of the mechanism of ice binding. More statistical analysis and investigation of the energetics of additional simulations of ice binding will help to probe and quantify the still debated mechanism of binding. Energy and entropy play a role in this mechanism, but more work is needed to elucidate the impact of these two contributions to binding. Mutation of the protein is an interesting tool to further investigate the antifreeze protein binding interaction. The methyl groups of the threonine, for example, seem to fill the ice lattice and modulate the geometry of hydrogen bonding in water. Replacement of 21 a threonine with a serine removes the crowding methyl but preserves the hydroxyl for continued hydrogen bonding. It would be interesting to examine how such a mutation could affect the ice binding ability. Alternatively, the methyl groups may be important because they help preserve the rigidity of the threonine-lined channel. Atomistic simulations still do not provide perfect accuracy, and all forces are still estimated. Many researchers utilize “coarse graining” to further simplify the calculations necessary to simulate a collection of particles. Often, coarse graining involves removing the least important atoms, like the extremely common but small hydrogen, in order to greatly reduce the computational cost. In order to model an ice binding protein system larger than the Tenebrio molitor antifreeze protein, such as the interesting family of ice nucleating proteins, the development of a coarse grained model will be necessary because of limitations in computing. Computational chemistry quickly becomes a balancing act. Using simple functions to govern pair interactions can increase processing power by orders of magnitude compared to Ab-Initio simulations, but also reduces precision. Since INPs are presently too large to simulate economically, useful nucleators, like the Pseudomonas borealis, are still confined to experimental investigation. The inaccurate parameters in many of the newest (and softest) force fields are an important consideration in the development of molecular models. 22 References [1] Amato, P., Joly, M., Schaupp, C., Attard, E., Möhler, O., Morris, C. E., Brunet, Y., and Delort, A.-M. (2015) Survival and ice nucleation activity of bacteria as aerosols in a cloud simulation chamber. Atmospheric Chemistry and Physics Atmos. Chem. Phys. 15, 6455–6465. [2] Charpentier, T. V., Neville, A., Millner, P., Hewson, R., and Morina, A. (2013) An Investigation of Freezing of Supercooled Water on Anti-Freeze Protein Modified Surfaces. Journal of Bionic Engineering 10, 139–147. [3] Daley, M. E., and Sykes, B. D. (2003) The role of side chain conformational flexibility in surface recognition by Tenebrio molitor antifreeze protein. Protein Sci. Protein Science 12, 1323–1331. [4] Davies, P. L. (2014) Ice-binding proteins: a remarkable diversity of structures for stopping and starting ice growth. Trends in Biochemical Sciences 39, 548–555. [5] Drori, R., Celik, Y., Davies, P. L., and Braslavsky, I. (2014) Ice-binding proteins that accumulate on different ice crystal planes produce distinct thermal hysteresis dynamics. Journal of The Royal Society Interface 11, 20140526–20140526. [6] Drori, R., Davies, P. L., and Braslavsky, I. (2015) Experimental correlation between thermal hysteresis activity and the distance between antifreeze proteins on an ice surface. RSC Adv. 5, 7848–7853. [7] Drori, R., Davies, P. L., and Braslavsky, I. (2015) When Are Antifreeze Proteins in Solution Essential for Ice Growth Inhibition? Langmuir 31, 5805–5811. 23 [8] Duboué-Dijon, E., and Laage, D. (2014) Comparative study of hydration shell dynamics around a hyperactive antifreeze protein and around ubiquitin. The Journal of Chemical Physics J. Chem. Phys. 141. [9] Garnham, C. P., Campbell, R. L., and Davies, P. L. (2011) Anchored clathrate waters bind antifreeze proteins to ice. Proceedings of the National Academy of Sciences 108, 7363–7367. [10] Garnham, C. P., Campbell, R. L., Walker, V. K., and Davies, P. L. (2011) Novel dimeric β-helical model of an ice nucleation protein with bridged active sites. BMC Structural Biology BMC Struct Biol 11, 36. [11] Graether, S. P., and Jia, Z. (2001) Modeling Pseudomonas syringae Ice- Nucleation Protein as aβ-Helical Protein. Biophysical Journal 80, 1169–1173. [12] Grandum, S., Yabe, A., Nakagomi, K., Tanaka, M., Takemura, F., Kobayashi, Y., and Frivik, P.-E. (1999) Analysis of ice crystal growth for a crystal surface containing adsorbed antifreeze proteins. Journal of Crystal Growth 205, 382–390. [13] Hartmann, S., Augustin, S., Clauss, T., Wex, H., Šantl-Temkiv, T., Voigtländer, J., Niedermeier, D., and Stratmann, F. (2013) Immersion freezing of ice nucleation active protein complexes. Atmospheric Chemistry and Physics Atmos. Chem. Phys. 13, 5751–5766. [14] Jia, Z., and Davies, P. L. (2002) Antifreeze proteins: an unusual receptor–ligand interaction. Trends in Biochemical Sciences 27, 101–106. [15] Kar, R. K., and Bhunia, A. (2015) Biophysical and biochemical aspects of antifreeze proteins: Using computational tools to extract atomistic information. Progress in Biophysics and Molecular Biology 119, 194–204. 24 [16] Kobashigawa, Y., Nishimiya, Y., Miura, K., Ohgiya, S., Miura, A., and Tsuda, S. (2005) A part of ice nucleation protein exhibits the ice-binding ability. FEBS Letters 579, 1493–1497. [17] Kuiper, M. J., Morton, C. J., Abraham, S. E., and Gray-Weale, A. (2015) The biological function of an insect antifreeze protein simulated by molecular dynamics. eLife 4. [18] Kumaki, Y., Kawano, K., Hikichi, K., Matsumoto, T., and Matsushima, N. (2008) A circular loop of the 16-residue repeating unit in ice nucleation protein. Biochemical and Biophysical Research Communications 371, 5–9. [19] Meister, K., Lotze, S., Olijve, L. L. C., Devries, A. L., Duman, J. G., Voets, I. K., and Bakker, H. J. (2015) Investigation of the Ice-Binding Site of an Insect Antifreeze Protein Using Sum-Frequency Generation Spectroscopy. J. Phys. Chem. Lett. The Journal of Physical Chemistry Letters 6, 1162–1167. [20] Li, Q., Yan, Q., Chen, J., He, Y., Wang, J., Zhang, H., Yu, Z., and Li, L. (2012) Molecular Characterization of an Ice Nucleation Protein Variant (InaQ) from Pseudomonas syringae and the Analysis of Its Transmembrane Transport Activity in Escherichia coli. Int. J. Biol. Sci. International Journal of Biological Sciences 8, 1097–1108. [21] Meister, K., Strazdaite, S., Devries, A. L., Lotze, S., Olijve, L. L. C., Voets, I. K., and Bakker, H. J. (2014) Observation of ice-like water layers at an aqueous protein surface. Proceedings of the National Academy of Sciences Proc Natl Acad Sci USA 111, 17732–17736. 25 [22] Midya, U. S., and Bandyopadhyay, S. (2014) Hydration Behavior at the Ice- Binding Surface of the Tenebrio molitor Antifreeze Protein. The Journal of Physical Chemistry B J. Phys. Chem. B 118, 4743–4752. [23] Modig, K., Qvist, J., Marshall, C. B., Davies, P. L., and Halle, B. (2010) High water mobility on the ice-binding surface of a hyperactive antifreeze protein. Phys. Chem. Chem. Phys. Physical Chemistry Chemical Physics 12, 10189. [24] Murakami, D., and Yasuoka, K. (2011) Molecular Dynamics Simulation of Quasi-Two-Dimensional Water Network on Ice Nucleation Protein. ASME/JSME 2011 8th Thermal Engineering Joint Conference. [25] Nada, H., and Furukawa, Y. (2012) Antifreeze proteins: computer simulation studies on the mechanism of ice growth inhibition. Polym J Polymer Journal 44, 690–698. [26] Nada, H., and Furukawa, Y. (2011) Growth inhibition at the ice prismatic plane induced by a spruce budworm antifreeze protein: a molecular dynamics simulation study. Phys. Chem. Chem. Phys. Physical Chemistry Chemical Physics 13, 19936. [27] Nutt, D. R., and Smith, J. C. (2008) Dual Function of the Hydration Layer around an Antifreeze Protein Revealed by Atomistic Molecular Dynamics Simulations. J. Am. Chem. Soc. Journal of the American Chemical Society 130, 13066–13073. [28] Pertaya, N., Marshall, C. B., Diprinzio, C. L., Wilen, L., Thomson, E. S., Wettlaufer, J., Davies, P. L., and Braslavsky, I. (2007) Fluorescence Microscopy Evidence for Quasi-Permanent Attachment of Antifreeze Proteins to Ice Surfaces. Biophysical Journal 92, 3663–3673. 26 [29] Schmid, D., Pridmore, D., Capitani, G., Battistutta, R., Neeser, J.-R., and Jann, A. (1997) Molecular organisation of the ice nucleation protein InaV from Pseudomonas syringae. FEBS Letters 414, 590–594. [30] Smolin, N., and Daggett, V. (2008) Formation of Ice-like Water Structure on the Surface of an Antifreeze Protein †. The Journal of Physical Chemistry B J. Phys. Chem. B 112, 6193–6202. [31] Stevens, C. A., Drori, R., Zalis, S., Braslavsky, I., and Davies, P. L. (2015) Dendrimer-Linked Antifreeze Proteins Have Superior Activity and Thermal Recovery. Bioconjugate Chem. Bioconjugate Chemistry 26, 1908–1915. [32] Sun, T., Gauthier, S., Campbell, R., and Davies, P. (2015) Revealing surface waters on an antifreeze protein by fusion protein crystallography. [33] Sun, T., Lin, F.-H., Campbell, R. L., Allingham, J. S., and Davies, P. L. (2014) An Antifreeze Protein Folds with an Interior Network of More Than 400 SemiClathrate Waters. Science 343, 795–798. [34] Todde, G., Whitman, C., Hovmöller, S., and Laaksonen, A. (2014) Induced Ice Melting by the Snow Flea Antifreeze Protein from Molecular Dynamics Simulations. The Journal of Physical Chemistry B J. Phys. Chem. B 118, 13527–13534. [35] Todde, G., Hovmöller, S., and Laaksonen, A. (2015) Influence of Antifreeze Proteins on the Ice/Water Interface. The Journal of Physical Chemistry B J. Phys. Chem. B 119, 3407–3413. [36] Yang, C., and Sharp, K. A. (2004) The mechanism of the type III antifreeze protein action: a computational study. Biophysical Chemistry 109, 137–148. 27 [37] Yang, Z., Zhou, Y., Liu, K., Cheng, Y., Liu, R., Chen, G., and Jia, Z. (2003) Computational Study on the Function of Water within a β-Helix Antifreeze Protein Dimer and in the Process of Ice-Protein Binding. Biophysical Journal 85, 2599–2605. [38] Lupi, L., Hudait, A., and Molinero, V. (2014) Heterogeneous Nucleation of Ice on Carbon Surfaces. J. Am. Chem. Soc. Journal of the American Chemical Society 136, 3156–3164. [39] Liou, Y., Tocilj, A., Davies, P. L., Jia, Z. (2000) Mimicry of ice structure by surface hydroxyls and water of a β-helix antifreeze protein. Letters to Nature 406, 322-324. [40] D.A. Case, J.T. Berryman, R.M. Betz, D.S. Cerutti, T.E. Cheatham, III, T.A. Darden, R.E. Duke, T.J. Giese, H. Gohlke, A.W. Goetz, N. Homeyer, S. Izadi, P. Janowski, J. Kaus, A. Kovalenko, T.S. Lee, S. LeGrand, P. Li, T. Luchko, R. Luo, B. Madej, K.M. Merz, G. Monard, P. Needham, H. Nguyen, H.T. Nguyen, I. Omelyan, A. Onufriev, D.R. Roe, A. Roitberg, R. Salomon-Ferrer, C.L. Simmerling, W. Smith, J. Swails, R.C. Walker, J. Wang, R.M. Wolf, X. Wu, D.M. York and P.A. Kollman (2015), AMBER 2015, University of California, San Francisco. [41] Scott Le Grand; Andreas W. Goetz; & Ross C. Walker. SPFP: Speed without compromise - a mixed precision model for GPU accelerated molecular dynamics simulations., Comp. Phys. Comm, (2013), 184, pp374-380, DOI: 10.1016/j.cpc.2012.09.022 [42] Crowley, M.; Williamson, M.; Walker, R. CHAMBER: Comprehensive support for CHARMM force fields within the AMBER software. Int. J. Quant. Chem., (2009),109, 3767–3772. 28 [43] Cornell, W.; Cieplak, P.; Bayly, C.; Gould, I.; Merz, K. Jr.; Ferguson, D.; Spellmeyer, D.; Fox, T.; Caldwell, J.; Kollman, P. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J. Am. Chem. Soc., (1995), 117, 5179–5197. [44] MacKerell Jr., A.; Bashford, D.; Bellott, M.; Dunbrack, R.; Evanseck, J.; Field, M.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D.; Prodhom, B.; Reiher, W.; Roux, B.; Schlenkrich, M.; Smith, J.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B, (1998), 102, 3586–3616. [45] Dalal, P., Knickelbein, J., Haymet, A. D. J., Sönnichsen, F. D., and Madura, J. D. (2001) Hydrogen bond analysis of Type 1 antifreeze protein in water and the ice/water interface. PhysChemComm 4, 32–36. [46] Knight, C., Cheng, C., and Devries, A. (1991) Adsorption of alpha-helical antifreeze peptides on specific ice crystal surface planes. Biophysical Journal 59, 409–418. [47] Humphrey, W., Dalke, A. and Schulten, K., "VMD - Visual Molecular Dynamics", J. Molec. Graphics, 1996, vol. 14, pp. 33-38. [48] Olijve, L. L. C., Meister, K., Devries, A. L., Duman, J. G., Guo, S., Bakker, H. J., and Voets, I. K. (2016) Blocking rapid ice crystal growth through nonbasal plane adsorption of antifreeze proteins. Proceedings of the National Academy of Sciences Proc Natl Acad Sci USA 113, 3740–3745. 29 [49] Best, R. B., Zhu, X., Shim, J., Lopes, P. E. M., Mittal, J., Feig, M., and Mackerell, A. D. (2012) Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the Backbone ϕ, ψ and Side-Chain χ 1 and χ 2 Dihedral Angles. J. Chem. Theory Comput. Journal of Chemical Theory and Computation 8, 3257–3273. [50] Abascal, J. L. F., and Vega, C. (2005) A general purpose model for the condensed phases of water: TIP4P/2005. The Journal of Chemical Physics J. Chem. Phys. 123, 234505. .