# Supporting Information for “On the Engulfment of Antifreeze Proteins by Ice”

Aniket U. Thosar  Yusheng Cai  Sean M. Marks  Zachariah Vicars  Jeongmoon Choi
Akash Pallath  Amish J. Patel<EMAIL_ADDRESS>Department of
Chemical and Biomolecular Engineering, University of Pennsylvania

##  1 Molecular Simulations of Antifreeze Proteins (AFPs)

###  1.1 Simulation parameters

All molecular dynamics simulations were performed using GROMACS version 2016.3
1. The equations of motion were integrated using the leap-frog algorithm with
a time step of 2 fs, and periodic boundary conditions were used in all
dimensions. To simulate systems at a particular temperature, the stochastic
velocity-rescale thermostat 2 was employed with a time constant of 0.5 ps,
whereas system pressure was maintained using the anisotropic Parrinello-Rahman
barostat 3, 4 with a time constant of 4 ps; fluctuations in the simulation box
dimensions were allowed only in the $z$-direction. Water molecules were
represented explicitly using the TIP4P/Ice model 5 and the AFPs were modeled
using the Amber99SB force field 6. The Lorentz-Berthelot mixing rules were
used to obtain the cross-parameters for Lennard-Jones interactions. Both the
Lennard-Jones and the short-range electrostatic interactions were truncated at
1 nm, and long-range electrostatic interactions were treated using the
particle mesh Ewald (PME) method 7. Water hydrogens were constrained using the
SETTLE algorithm 8, whereas bonds involving protein hydrogens were constrained
using the LINCS algorithm 9, 10. To restrain the position and orientation of
the AFP, the positions of all protein heavy atoms within 0.5 nm of its center
of mass were restrained using a relatively soft spring constant of 1000
kJ/mol/nm2. By restraining only a small fraction of heavy atoms (roughly 20
atoms or 1 % of the total) near the protein’s center of mass, we are able to
localize the AFP without restraining any of its surface atoms, i.e., all
protein atoms on the IBS and the NBS can fluctuate freely.

Figure 1:  (A) Simulation snapshot of sbwAFP (pink) bound to the primary prism
plane of ice. The water oxygens shown in blue and purple are position
restrained, resulting in liquid-like and ice-like boundary slabs,
respectively, which ensure that ice grows in the $+z$-direction. The
unrestrained (mobile) ice and liquid water molecules are shown in cyan and
red, respectively. (B) Simulation snapshot of sbwAFP engulfed by ice, which
serves as the initial configuration for biased simulations. The snapshots were
visualized using the Visual Molecular Dynamics (VMD) package 11.

###  1.2 Simulation set-up

We studied the following four AFPs: (i) sbwAFP from the spruce budworm
(Choristoneura Fumiferana, Isoform 501; PDB ID: 1M8N), (ii) TmAFP from
Tenebrio molitor (PDB ID: 1EZG), (iii) opAFP from ocean pout (Zoarces
americanus, PDB ID: 1HG7), and (iv) LpAFP from Lolium perenne (PDB ID: 3ULT)
12, 13, 14, 15. The protein systems were prepared for simulation using the
‘pdb2gmx’ utility of the GROMACS molecular dynamics package 1. The majority of
our calculations were performed using sbwAFP with the simulation box
dimensions, $L_{x}$ and $L_{y}$, chosen to be roughly 6 nm (Figure 1A) to
ensure that the protein heavy-atoms were at least 1 nm away from the box
edges. Simulation boxes with $L_{x}$ and $L_{y}$ ranging from 6 nm to 9 nm
were used to characterize the variation of interfacial curvature with AFP
separation (Figures 3B and 4). To set up our simulations, the GenIce package
16 was first used to generate a 2 nm thick slab of ice (Figure 1A, cyan) with
its primary prism ($10\bar{1}0$) plane normal to the $z$-direction. The
Packmol package 17 was then used to judiciously place the AFP of interest in
the simulation box. In particular, the AFP was centered in the $x$ and $y$
directions, and along the $z$-direction, the AFP was placed less than 0.5 nm
from the ice slab with its IBS facing the primary prism plane; the azimuthal
orientation of the AFP was chosen to ensure that the chemical groups on the
IBS (e.g., threonine hydroxyls) were aligned with the hydroxyls of the ice
slab. The protein was then solvated in a roughly 7 nm thick slab of water
(Figure 1A, red) using the ‘solvate’ utility of GROMACS; sodium or chloride
ions were added, if necessary, to maintain the charge-neutrality of the
system. Lastly, two 1 nm thick position-restrained slabs of liquid water
(Figure 1A, blue) and ice (Figure 1A, purple) were introduced at the bottom of
the simulation box to break translational symmetry and ensure that ice grows
in the +$z$ direction. To prevent intrusion of the mobile water molecules
(Figure 1A, cyan and red) into these boundary slabs, the spring constants for
the restraining potentials were chosen to be 40,000 kJ/mol/nm2 and 4,000
kJ/mol/nm2 for the liquid water and ice slabs, respectively. The solvated AFP
was then position-restrained, as described above, and energy-minimized using
the steepest-descent algorithm, followed by a 0.5 ns long NVT simulation (at
$T=298$ K), and a 2 ns long NPT simulation (at $T=298$ K and $P=1$ bar). To
facilitate binding of the AFP to the ice-water interface, the position
restraints on the AFP were relaxed, and a 20 ns long simulation was run at
$T=298$ K and $P=1$ bar in the presence of a biasing potential that prevents
the ice slab from melting, as described in sec. 2B below. The final
configuration is shown in Figure 1A.

##  2 Characterizing the Free Energetics of AFP Engulfment

To prepare the initial configurations for our simulations as well as to
subsequently characterize the free energetics of AFP engulfment, we bias the
number of ice-like mobile waters in our simulations.

###  2.1 Biasing the approximate number of ice-like waters.

To bias the number of ice-like waters using molecular dynamics simulations,
the order parameter being biased must not only be capable of discriminating
between liquid water and ice, but also be a continuous and differentiable
function of particle positions 18. To this end, we define the order parameter,
$M$, corresponding to the approximate number of ice-like molecules within a
observation volume $v$ in our system as follows 19:

$M\equiv\sum_{i=1}^{N_{\rm
w}}\tilde{h}_{v}(i)\cdot\tilde{h}_{\bar{q}}(\bar{q}_{6}(i)),$ (1)

where the sum is performed over the $N_{\rm w}$ mobile waters in the
simulation box; the coarse-grained function $\tilde{h}_{v}(i)$ indicates
whether a water molecule $i$ is in an observation volume of interest, $v$; the
indicator function, $\tilde{h}_{\bar{q}}(\bar{q}_{6}(i))$, is used to classify
whether a water molecule, $i$, is ice-like, and its argument, $\bar{q}_{6}(i)$
is closely related to the local order parameter introduced by Lechner and
Dellago 20. In particular, $\bar{q}_{6}(i)$ is defined as:

$\bar{q}_{6}(i)\equiv\sqrt{\frac{4\pi}{13}\sum_{m=-6}^{6}|\bar{q}_{6m}(i)|^{2}}\text{,}$
(2)

where

$\bar{q}_{6m}(i)\equiv\frac{1}{1+\tilde{N}_{\mathrm{nn}}(i)}\cdot\Bigg{[}q_{6m}(i)+\sum_{k=1}^{N_{w}}\tilde{h}_{v_{i}}(k)q_{6m}(k)\Bigg{]}.$
(3)

Here, $\tilde{h}_{v_{i}}(k)$ indicates whether water $k$ is a nearest neighbor
of water $i$, $\tilde{N}_{\mathrm{nn}}(i)\equiv\sum_{k\neq i}^{N_{\rm
w}}\tilde{h}_{v_{i}}(k)$ is the coarse-grained number of nearest neighbors of
water $i$, and

$q_{6m}(i)\equiv\frac{1}{\tilde{N}_{{\rm{nn}}}(i)}\sum_{k=1}^{N_{\rm
w}}\tilde{h}_{v_{i}}(k)Y_{6m}(\mathbf{r}_{ik}),$ (4)

where $Y_{6m}(\mathbf{r}_{ik})$ is the spherical harmonic function, and its
argument, $\mathbf{r}_{ik}$, is the vector from the oxygen of water $i$ to
that of water $k$. To ensure that the indicator functions, $\tilde{h}_{v}(i)$,
$\tilde{h}_{v_{i}}(k)$, and $\bar{h}_{\bar{q}}$ are continuous and
differentiable functions of particle positions, they are obtained by
integrating Gaussian functions of width $\sigma$ that are truncated and
shifted at $2\sigma$ following refs. 21, 19. The indicator functions are
plotted in Figure 2. The observation volume, $v$, is chosen to be box-spanning
in the $x$ and $y$ directions, and spans from $z_{\rm min}=1.5$ nm to $z_{\rm
max}=8.5$ nm along the $z$ direction, so that $\tilde{h}_{v}(i)$ increases
continuously from 0 to 1 (over a distance of $z_{\rm min/max}\pm 2\sigma$ with
$\sigma=0.01$ nm) as water molecule $i$ approaches the boundaries of $v$, as
shown in Figure 2A; $z_{i}$ denotes the $z$-co-ordinate of water $i$. Figure
2B shows the indicator function $\tilde{h}_{v_{i}}(k)$ used to determine the
nearest neighbors of water $i$. As the distance, $r_{ik}$, between water $i$
and its putative neighbor $k$ increases above $r_{\rm max}=0.35$ nm, the value
of $\tilde{h}_{v_{i}}(k)$ decreases smoothly from 1 to 0 (over $r_{\rm max}\pm
2\sigma$ with $\sigma=0.01$ nm). Figure 2C shows the indicator function
$\tilde{h}_{\bar{q}}$, which determines whether the local orientational order
of water $i$ resembles that in bulk ice. As $\bar{q}_{6}(i)$ crosses the value
of $q_{\rm min}=0.352$, which occurs with equal probability in bulk water and
ice at 273 K, the value of $\tilde{h}_{\bar{q}}$ increases smoothly from 0 to
1 (over $q_{\rm min}\pm 2\sigma$ with $\sigma=0.01$). The approximate number
of ice-like waters, $M$, was biased, and the corresponding forces were
distributed across the water molecules in the system using an in-house plug-in
that was developed to interface with the PLUMED package 22, 23, 24.

Figure 2:  Functional forms of the indicator functions used to calculate the
number of ice-like molecules, $M$, in an observation volume, $v$, of interest.
(A) The indicator function $\tilde{h}_{v}(i)$ determines whether or not water
molecule $i$ is inside $v$. (B) The indicator function $\tilde{h}_{v_{i}}(k)$
determines whether water $k$ is a nearest neighbor of water $i$. (C) The
indicator function $\tilde{h}_{\bar{q}}(\bar{q}_{6}(i))$ determines whether
the local orientational order of water $i$ resembles that of ice.

###  2.2 Biased Simulations

We performed biased simulations to prepare the initial configurations for
subsequent calculations as well as to characterize the free energetics of AFP
engulfment. In particular, biased simulations were used to obtain initial
configurations with the AFP bound to ice (Figure 1A) and the AFP engulfed by
ice (Figure 1B). To obtain the former, a linear biasing potential, $\phi M$,
was applied (with $\phi=-0.548$ kJ/mol) to prevent the melting (or growth) of
the ice slab (Figure 1A, cyan). To engulf the AFP, we initialized the system
in the bound state configuration (Figure 1A), and ran a 20 ns long biased
simulation using a linear potential, $\phi M$ (with $\phi=-1.1$ kJ/mol), which
resulted in the engulfed configuration shown in Figure 1B. Subsequent biased
simulations for characterizing the free energetics of engulfment were
initialized using the engulfed configuration (Figure 1B), and employed
harmonic potentials, $\frac{\kappa}{2}(M-M^{*})^{2}$, with $M^{*}$ being
varied systematically to sample the pertinent range of $M$-values. To mitigate
the formation of defects, $M^{*}$ is decreased linearly to its desired value
over 3 ns. All biased simulation are run for 15 ns, which includes 3 ns for
varying $M^{*}$ to its final value, the subsequent 5 ns for equilibration, and
the final 7 ns for production.

###  2.3 Estimating the true number of ice molecules.

Although $\bar{q}_{6}(i)$ and similar order parameters are excellent at
discriminating between liquid water and ice, a small fraction of water
molecules in the liquid phase will nevertheless have high
$\bar{q}_{6}(i)$-values, and be incorrectly classified as ice-like, giving
rise to false positives. Moreover, classifying interfacial water molecules,
i.e., those near solid surfaces, as ice-like and liquid-like is more
challenging than classifying waters in the bulk because the interfacial waters
are neighbor-deficient, i.e., they have fewer than the four nearest neighbors
found in the bulk phases. Thus, the order parameter, $M$, defined in Equation
1 above, is unable to identify ice-like molecules in the vicinity of the AFP,
resulting in false negatives. To more accurately estimate the true number of
ice molecules, $\lambda$, we first used the CHILL+ algorithm 25 to identify
the ice-like waters (hexagonal, cubic and interfacial ice) in our system and
classified the remainder as liquid-like waters. As shown in Figure 3A,
although CHILL+ is excellent at identifying molecules that belong to ice, it
nevertheless misidentifies a few molecules that are clearly in bulk water. To
filter out the misidentified molecules, we located the instantaneous ice-water
interface following ref. 26, and water molecules that were determined to be
above the interface using the Moller-Trumbore algorithm 27 were reclassified
as being liquid-like. As shown in Figure 3B, this procedure was able to
successfully filter out the misidentified ice-like molecules seen in Figure
3A. Finally, to identify ice-like interfacial waters, i.e., those near the
AFP, we used the following procedure: if an AFP hydration water (i.e., water
within 0.6 nm of a protein heavy atom), which has been classified as liquid-
like, is found to be within 0.55 nm of 8 or more ice-like waters, it is
reclassified as being an ice-like water; this procedure is then repeated
recursively until no further reclassification occurs. As shown in Figure 3C,
this protocol enables identification of the ice-like waters near the AFP
(yellow).

Figure 3:  Estimating the true number of ice molecules, $\lambda$. (A) Ice-
like waters identified using the CHILL+ algorithm 25 are shown (cyan) for a
partially engulfed AFP (pink, front view), whereas liquid waters are hidden
for clarity. A few molecules in bulk water are misidentified as ice-like. (B)
Ice molecules determined using Moller-Trumbore algorithm 27 to be above the
instantaneous ice-water interface 26 are filtered out. (C) Waters near the AFP
(top view) that were reclassified to being ice-like are shown here in yellow.

###  2.4 Reweighting

Although the approximate number of ice-like waters, $M$, is biased using a
harmonic potential, $\frac{\kappa}{2}(M-M^{*})^{2}$ in our simulations, the
correlation between $M$ and the true number of ice molecules, $\lambda$,
enables us to indirectly sample a wide range of $\lambda$-values across our
biased simulations 18. We then reweight the biased simulations to obtain the
free energetics of engulfment, $G(\lambda;\Delta T)$, and the corresponding
resistance to engulfment, $g_{\lambda}\equiv dG/d\lambda$. For biased
distributions with overlap in their sampled $M$ (and $\lambda$) values,
$G(\lambda;\Delta T)$ and $g_{\lambda}$ were obtained using standard weighted
histogram reweighting methods 28, 29, 18. To estimate $g_{\lambda}$ more
efficiently for different separations, $L$, and AFP types, we additionally
made use of the sparse sampling approach 30, which does not require overlap
between adjacent biased distributions. In particular, we first obtained
$g_{M}\equiv dG/dM$ at select values of $M=\langle M\rangle_{\kappa,M^{*}}$
using the relationship $g_{M}\big{|}_{M=\langle
M\rangle_{\kappa,M^{*}}}=\kappa(M^{*}-\langle M\rangle_{\kappa,M^{*}})$ 30, as
shown in Figure 4A. We then averaged $\lambda$ over our biased simulations to
obtain estimates of $\langle\lambda\rangle_{\kappa,M^{*}}$, and found a strong
linear correlation between $\langle M\rangle_{\kappa,M^{*}}$ and
$\langle\lambda\rangle_{\kappa,M^{*}}$, as shown in Figure 4B. Accordingly, we
posited that $M$ and $\lambda$ were linearly related to one another with the
proportionality constant, $M_{\lambda}$, being well-approximated by the slope
of the linear fit to $\langle M\rangle_{\kappa,M^{*}}$ vs
$\langle\lambda\rangle_{\kappa,M^{*}}$ (Figure 4B). Finally, we estimated
$g_{\lambda}$ for $\lambda$-values corresponding to
$\lambda=\langle\lambda\rangle_{\kappa,M^{*}}$ as:

$g_{\lambda}\big{|}_{\lambda=\langle\lambda\rangle_{\kappa,M^{*}}}=g_{M}\big{|}_{M=\langle
M\rangle_{\kappa,M^{*}}}\cdot M_{\lambda}.$ (5)

In Figure 4D, we compare our estimates of $g_{\lambda}$ obtained using sparse
sampling (Equation 5) against those obtained using WHAM 28, 29, and find
excellent agreement between them. The errors in our sparse sampling estimates
of $g_{\lambda}$ were obtained using error propagation with contributions from
the standard errors in $\langle M\rangle_{\kappa,M^{*}}$ as well as errors in
$M_{\lambda}$ (obtained using block bootstrapping).

Figure 4:  Estimating $g_{\lambda}$ for sbwAFP with $L\approx 6$ nm and
$T=298$ K. (A) At select values of $M=\langle M\rangle_{\kappa,M^{*}}$,
$g_{M}$ is shown as a function of $M$. (B) We find $\langle
M\rangle_{\kappa,M^{*}}$ to be linearly correlated with
$\langle\lambda\rangle_{\kappa,M^{*}}$ with the slope of the best linear fit
being 0.92. (C) Estimates of $g_{\lambda}$ obtained using sparse sampling
(symbols) are in good agreement with those obtained using WHAM (line).

##  3 Non-equilibrium Simulations

To interrogate whether the supercooling at which sbwAFP (with $L\approx 6$ nm)
is engulfed spontaneously by ice corresponds to the supercooling, $\Delta
T^{*}=18$ K, at which engulfment barriers in $G(\lambda;\Delta T)$ vanish, we
performed unbiased non-equilibrium simulations wherein an sbwAFP molecule
bound to the ice-water interface is quenched to different temperatures. The
system is initialized with the AFP bound to ice, as shown in Figure 1A. Then,
at time $t=0$, the system is supercooled to either $\Delta T=20$ K or $\Delta
T=16$ K, and the corresponding non-equilibrium simulations are run for 240 ns;
the final configurations of these simulations are shown in Figures 5A and 5B,
respectively, and the time-dependence of the number of ice molecules during
these simulations is shown in Figure 5C. These results highlight that sbwAFP
arrests ice growth for $\Delta T<\Delta T^{*}$, but is engulfed by ice for
$\Delta T>\Delta T^{*}$.

Figure 5:  (A) When sbwAFP bound to the ice-water interface is quenched to a
supercooling of $\Delta T=20$ K, it is engulfed by ice, but (B) when it is
quenched to a supercooling of $\Delta T=16$ K, it arrests ice growth. (C) The
number of ice molecules, $\lambda$, is shown as a function of time for non-
equilibrium simulations with a supercooling of $\Delta T=20$ K (blue) and
$\Delta T=16$ K (red).

##  4 Macroscopic Interfacial Thermodynamics

According to macroscopic interfacial thermodynamics, the free energetics of a
system, which has an ice-water interface profile, $h(x,y)$, contains $\lambda$
ice-like waters, and is at supercooling, $\Delta T$, is given by:

$G_{\rm th}([h],\lambda;\Delta T)=-\Delta\mu\rho V[h]+\gamma A[h]+\Delta\gamma
A_{\rm nbs}+\mathcal{L}_{\lambda}(\lambda-\rho V[h]),$ (6)

where $V[h]=\iint h\,dx\,dy$ is the volume of ice in the system and
$A[h]=\iint\sqrt{1+h_{x}^{2}+h_{y}^{2}}\,dx\,dy$ is the ice-water interfacial
area. As discussed in the main text, $\Delta\mu$ is the chemical potential
difference between liquid water and ice, $\rho$ is the density of ice,
$\gamma$ is the ice-water interfacial tension, $\Delta\gamma$ quantifies NBS
ice-phobicity, $A_{\mathrm{nbs}}$ is the surface area of the NBS in contact
with ice; and $\mathcal{L}_{\lambda}$ is a Lagrange multiplier that constrains
$\rho V$ to $\lambda$. Because we are interested in studying AFP engulfment,
we assume that the AFP is bound to ice for all interfacial profiles, $h(x,y)$,
and $\lambda$-values under consideration. above. Moreover, ice-phobicity can
vary across different regions of the NBS, so $\Delta\gamma$ represents an
average over the area of the NBS that is in contact with ice. Finally, because
$A_{\rm nbs}$ does not depend on the shape of the interfacial profile, but
only on the location of its boundaries, the functional derivative of $G_{\rm
th}$ with respect to $h(x,y)$ is:

$\frac{\delta G}{\delta h}=-\rho(\mathcal{L}_{\lambda}+\Delta\mu)\frac{\delta
V}{\delta h}+\gamma\frac{\delta A}{\delta h}$ (7)

Therefore, the mean-field interface, $\bar{h}(x,y)$, must obey:

$\frac{\delta G}{\delta h}\bigg{|}_{\bar{h}}=0\implies\frac{\delta A/\delta
h}{\delta V/\delta
h}\bigg{|}_{\bar{h}}=\frac{\rho(\mathcal{L}_{\lambda}+\Delta\mu)}{\gamma}.$
(8)

By taking the functional derivatives of $V[h]$ and $A[h]$, we get $\delta
V/\delta h=1$ and:

$\frac{\delta A}{\delta
h}=-\frac{[1+h_{x}^{2}]h_{yy}-2h_{x}h_{y}h_{xy}+[1+h_{y}^{2}]h_{xx}}{[1+h_{x}^{2}+h_{y}^{2}]^{3/2}}=2\kappa_{\lambda}(x,y),$
(9)

where $\kappa_{\lambda}(x,y)$ is the negative of the mean interfacial
curvature. Combining Equations 8 and 9, we obtain:

$\kappa_{\lambda}(x,y)\bigg{|}_{\bar{h}}=\frac{\rho}{2\gamma}\cdot(\mathcal{L}_{\lambda}+\Delta\mu).$
(10)

Because the right hand side of this equation is independent of $x$ and $y$,
the mean-field ice-water interface must be a constant mean curvature surface.

To illustrate the physical significance of the Lagrange multiplier,
$\mathcal{L}_{\lambda}$, we recognize that the mean-field free energy of the
system is given by $\bar{G}(\lambda;\Delta T)=G_{\rm
th}([\bar{h}],\lambda;\Delta T)$, and the corresponding resistance to
engulfment, $g_{\lambda}\equiv d\bar{G}(\lambda;\Delta T)/d\lambda$ is:

$g_{\lambda}=\frac{d}{d\lambda}G_{\rm th}([\bar{h}],\lambda;\Delta
T)=\int\int\cancelto{0}{\frac{\delta G_{\rm th}}{\delta
h}\bigg{|}_{\bar{h}}}\cdot\frac{d\bar{h}}{d\lambda}\cdot
dxdy+\frac{\partial}{\partial\lambda}G_{\rm th}([\bar{h}],\lambda;\Delta
T)=\mathcal{L}_{\lambda}+\frac{\partial\mathcal{L}_{\lambda}}{\partial\lambda}\cdot\cancelto{0}{(\lambda-\rho
V[\bar{h}])}=\mathcal{L}_{\lambda}.$ (11)

Finally, combining Equations 10 and 11, we obtain an expression for the
magnitude of the interfacial curvature:

$\kappa_{\lambda}=\frac{\rho}{2\gamma}\cdot(g_{\lambda}+\Delta\mu).$ (12)

##  5 Characterizing the shape of the ice-water interface

In this section, we describe the methods used to characterize the shape of the
ice-water interface, observed in our biased simulations, and to obtain the
corresponding mean and normalized Gaussian curvature fields. We illustrate the
techniques using partially engulfed sbwAFP molecules that are separated by
roughly 8 nm.

###  5.1 Locating the interface of ice with water and/or the AFP

For every configuration in our biased simulations, we first identified the
ice-like molecules using the procedure described in section 2C. Following ref.
26, we then obtained the coarse-grained density of ice-like molecules at every
position, $\mathbf{r}$, in our simulation box position as:

$\rho(\mathbf{r})=\sum_{i}\phi(|\mathbf{r}-\mathbf{r_{i}}|),$

where $\mathbf{r_{i}}$ corresponds to the position of the $i^{\rm th}$ ice-
like water and the coarse graining function $\phi(r)$ is defined as:

$\phi(r)=(2\pi\xi^{2})^{-3/2}\exp(-r^{2}/2\xi^{2}),$

with the coarse-graining length, $\xi=0.24$ nm. The instantaneous ice density
field, $\rho(\mathbf{r})$, was then averaged over the simulation trajectory,
to obtain the averaged ice density field, $\bar{\rho}(\mathbf{r})$, and the
Marching cubes algorithm was used to obtain the interface of ice with liquid
water and/or the AFP. In particular, the interface was defined as the
$\bar{\rho}(\mathbf{r})=s$ iso-surface, with $s$ chosen to be 15 nm-3 (i.e.,
roughly half the density of bulk ice). As shown in Figure 6A, the interface
(blue) separates ice (cyan) from liquid water (red) and the AFP (pink).

Figure 6:  (A) The (average) interface (blue) separating ice (cyan) from
liquid water (red) and the AFP (pink) is superimposed with a simulation
snapshot. (B) The mean curvature map of the smoothed interface.

###  5.2 Smoothing the interface and estimating interfacial curvature

Because curvature is sensitive to rough interfacial features, we first
smoothed the interface using the implicit fairing algorithm, which retains the
geometric characteristics of the interface (e.g., curvature) 31. In
particular, 10 iterations of the implicit fairing algorithm were executed with
a step size of 0.1 nm2. We then used polynomial fitting to compute the mean
curvature at every point on the smoothed interface 32. The mean curvature map
thus obtained is shown in Figure 6B; its central region corresponds to the
AFP-ice interface, whereas its peripheral regions correspond to the ice-water
interface. To obtain an estimate of the (constant) mean curvature associated
with the ice-water interface, $\kappa_{\lambda}^{\rm int}$, we averaged the
interfacial curvature over the regions farthest from the AFP (dashed lines
shown in Figure 6B).

###  5.3 Identifying the three-phase contact line

Although the above procedure provides a faithful representation of the ice-
water interface far from the AFP, the shape of the ice-water interface near
the AFP is obfuscated by the proximity of the AFP-ice interface. To address
this challenge, we first identify the three-phase contact line,
$\partial\mathcal{I}$. To this end, we averaged the number of ice molecules
within 1 nm of every AFP surface heavy atom, and normalized it by the
corresponding quantity for the fully engulfed AFP. A simulation snapshot of
the partially engulfed sbwAFP is shown in Figure 7A with every protein heavy
atom $i$ colored according to the normalized number of ice molecules,
$\phi_{i}$, in its vicinity; atom $i$ was deemed to be a contact line atom if
$0.5\leq\phi_{i}\leq 0.6$. The $z$-coordinates of the contact line atoms were
tightly clustered, so the atoms were fitted to a contact plane that is normal
to the $z$-direction, as shown in Figure 7B. The intersection between the
contact plane and the protein surface (Figure 7C) then provided an initial
approximation of the contact line (Figure 7D, purple). To obtain the protein
surface mesh, an iso-surface of the average coarse-grained density field of
all mobile waters, i.e., those belonging to both liquid water and ice, was
used, as described in section 5A. We note that an iso-surface of the coarse-
grained density field of protein heavy atoms can also be used to obtain the
protein mesh, and it results in a nearly identical estimate of the contact
line. Because the contact line resides on the ice-water interface (and not on
the protein), we extended our initial estimate of the contact line outward by
0.32 nm (corresponding to the Lennard Jones $\sigma$ of the water oxygen) to
obtain our final estimate of the contact line, $\partial\mathcal{I}$ (Figure
7D, red).

Figure 7:  (A) Simulation snapshot denoting ice molecules (cyan) and AFP heavy
atoms colored according to their $\phi_{i}$-values; the atoms colored in white
($0.5\leq\phi_{i}\leq 0.6$) approximate the location of the three-phase
contact line where the ice-water interface is pinned to the AFP. (B) The
positions of the atoms that fall on the three-phase contact line (blue) are
fitted to a plane (red). (C) The intersection of the contact plane (red) with
the surface mesh of protein (blue). (D) The intersection of contact plane with
protein surface mesh (purple) is extended outwards by the size of a water
molecule to obtain the three-phase contact line (red).  Figure 8:  Identifying
the ice-water interface near the contact line. (A) A periodic flat mesh
excluding the region enclosed by the contact line is generated using
constrained Delaunay triangulation 33. (B) The mesh is pinned at the contact
line and evolved into a constant mean curvature surface using the Lagrange
process 34; the curvature map of the converged mesh confirms that it is indeed
a constant mean curvature surface. (C) The evolved surface (red) and the
smoothed interface (blue, Figure 6A) agree well with one another. The root
mean square difference (RMSD) between the surfaces is less than 0.1 nm, which
is roughly 3% of the range of interface heights, $\Delta h$.

###  5.4 Determining the ice-water interface near the contact line

To determine the shape of the ice-water interface near the contact line,
$\partial\mathcal{I}$, we used constrained Delaunay triangulation 33 to first
construct a flat 2D mesh, which excludes the region enclosed by the contact
line, and is periodic in the $x$ and $y$ directions (Figure 8A). We then
employed a Lagrange process 34 to evolve the flat mesh into a constant mean
curvature surface using $\kappa_{\lambda}^{\rm int}$ determined above (sec.
5B). In particular, keeping the mesh pinned at the contact line, the positions
of the rest of the mesh vertices were evolved using:

$\frac{d\bar{r}_{i}}{dt}=-(\kappa_{i}-\kappa_{\lambda}^{\rm
int})\leavevmode\nobreak\ \hat{n}_{i},$ (13)

where $\bar{r}_{i}$ is the position of vertex $i$, $\kappa_{i}$ is the
magnitude of the mean curvature at that position, and $\hat{n}_{i}$ is the
corresponding unit vector normal to the interface. Upon convergence (i.e.,
$d\bar{r}_{i}/dt=0$), the generated interface satisfies
$\kappa_{i}=\kappa_{\lambda}^{\rm int}$ for every mesh vertex $i$. The
converged mean curvature field, thus obtained, is shown in Figure 8B. As shown
in Figure 8C, the corresponding ice-water interface (red) agrees well with the
smoothed interface (blue, Figure 6A). We also computed the root mean square
difference between the two surfaces (RMSD), and found it to be less than 0.1
nm or roughly 3 % of the range of heights spanned by the ice-water interface,
$\Delta h$ (Figure 8C).

Using the above procedure to characterize the shape of the ice-water
interface, both near and far from the AFP, we were able to obtain the
interface profiles and curvature fields for a wide range of $\lambda$-values
as sbwAFP molecules, separated by roughly 6 nm, undergo engulfment; see Figure
9.

Figure 9:  For sbwAFP molecules separated by roughly 6 nm, as $\lambda$ is
increased, the evolution of the average ice-water interface profiles,
$\bar{h}_{\lambda}(x,y)$, is shown in the top row, the corresponding mean
curvature fields, $\kappa_{\lambda}(x,y)$, are shown in the middle row, and
the corresponding Gaussian curvatures, $\mathcal{K}_{\lambda}$, normalized by
the square of the mean curvature, are shown in the bottom row. In each case,
$\kappa_{\lambda}(x,y)$ is independent of $(x,y)$ (middle row), highlighting
that the ice-water interfaces are constant mean curvature surfaces. However,
the interfaces have complex shapes and display regions of both positive and
negative Gaussian curvature (bottom row).

##  6 Variation of $g_{\rm m}$ with $\Delta T$

According to Equation 2 of the main text: $g_{\rm
m}=-\Delta\mu+(2\gamma/\rho)\cdot\kappa_{\rm m}$. Upon taking the derivative
of this equation with respect to the supercooling, $\Delta T$, we obtain:

$\frac{dg_{\rm m}}{d\Delta T}=-\frac{d\Delta\mu}{d\Delta T}+\frac{d}{d\Delta
T}\bigg{(}\frac{2\gamma}{\rho}\bigg{)}\kappa_{\rm m}\approx\frac{\Delta
h}{T_{\rm m}}-2\kappa_{\rm m}\frac{d}{dT}\bigg{(}\frac{\gamma}{\rho}\bigg{)},$
(14)

where we use the fact that $\Delta\mu\approx-(\Delta h/T_{\rm m})\Delta T$ and
that $d\Delta T=-dT$. Because $\Delta h=6.01$ kJ/mol and $T_{\rm m}=273.15$ K,
the first term on the right hand side is: $\Delta h/T_{\rm m}=0.022$ kJ/mol-K.
Similarly, for our water model, i.e., TIP4P/Ice 5, $\Delta h=5.6$ kJ/mol and
$T_{\rm m}=271$ K, so that $\Delta h/T_{\rm m}=0.02$ kJ/mol-K. To estimate the
second term on the right hand side, we note that the maximum curvature
obtained for sbwAFP separated by $L=6$ nm was roughly $\kappa_{\rm m}=0.22$
nm-1 (Figure 2C), and that $\kappa_{\rm m}$ only decreases for larger $L$
(varying as $L^{-2}$). Moreover, $d(\gamma/\rho)/dT=-0.0049$ nm-kJ/mol-K 35
(and -0.00533 nm-kJ/mol-K for our water model 36), so that the second term
must be smaller than 0.0022 kJ/mol-K (for $L\geq 6$ nm). Thus, the second term
on the right hand side of Equation 14 is roughly an order of magnitude smaller
than the first, suggesting that $(2\gamma/\rho)\kappa_{\rm m}$ is a weak
function of $T$ relative to $\Delta\mu$, and that the variation of $g_{\rm m}$
with $\Delta T$ stems primarily from the (linear) variation of $\Delta\mu$
with $\Delta T$.

These arguments provide a basis for understanding the linear relationship,
$g_{\rm m}(\Delta T)=-m_{g}\Delta T+c_{g}$, observed in Figure 1C, and suggest
that the fit parameters must be given by: $m_{g}=(\Delta h/T_{\rm m})$ and
$c_{g}=(2\gamma/\rho)\kappa_{\rm m}$. We find that
$m_{g}=0.02\leavevmode\nobreak\ {\rm kJ/mol-K}$ (Figure 1C) is indeed in
excellent agreement with $(\Delta h/T_{\rm m})=0.02\leavevmode\nobreak\ {\rm
kJ/mol-K}$ of the TIP4P/Ice water model 5. Similarly, $c_{g}=0.36$ kJ/mol
(Figure 1C) agrees well with $(2\gamma/\rho)\kappa_{\rm m}=0.324$ kJ/mol,
where $\gamma=3.7\times 10^{-5}$ kJ/m2 and $\rho=5\times 10^{4}$ mol/m3 for
TIP4P/Ice 5, 36, and $\kappa_{\rm m}=0.22$ nm-1 (Figure 2C).

##  7 Estimating $\eta_{\rm m}$ for sbwAFP

In Equation 4 of the main text, we showed that the maximum interfacial
curvature, $\kappa_{\rm m}$, that can be sustained by the AFP is given by:

$\kappa_{\rm m}=\frac{\eta_{\rm m}p_{\rm m}}{2(L^{2}-a_{\rm
m})}=\frac{1}{2}\eta_{\rm m}p_{\rm m}\frac{L^{-2}}{(1-a_{\rm m}L^{-2})},$ (15)

where $\eta_{\rm m}$ is the effectiveness of the NBS at resisting engulfment
by ice, $a_{\rm m}$ and $p_{\rm m}$ are the area enclosed by the three phase
contact line and its perimeter, respectively, and $L$ is the geometric mean of
the simulation box dimensions, $L_{x}$ and $L_{y}$. To estimate $\eta_{\rm m}$
for the AFP of interest, first we systematically varied the AFP separation,
$L$ (by varying the dimensions of our simulation box), and uncovered how it
influences the maximum curvature, $\kappa_{\rm m}$, that can be sustained by
the NBS. Then, we obtained estimates of $p_{\rm m}$ and $a_{\rm m}$ using the
contact line (Figure 7D, red). Finally, we fit the $\kappa_{\rm m}$ vs $L$
relationship to Equation 15, and obtained $\eta_{\rm m}$ as the parameter,
which provided the best fit to the data.

###  7.1 Characterizing the variation of $\kappa_{\rm m}$ with $L$

To characterize how $\kappa_{\rm m}$ varies with $L$, we performed biased
simulations of AFP engulfment using different simulation boxes sizes. To
determine $\kappa_{\lambda}$ for a given simulation box size, $L$, we used
Equation 12, i.e., $\kappa_{\lambda}^{\rm
th}=\frac{\rho}{2\gamma}\cdot(g_{\lambda}+\Delta\mu)$. In particular, by
leveraging the fact that $\kappa_{\lambda}^{\rm th}$ is relatively insensitive
to temperature, and recognizing that relaxation times decrease precipitously
with $T$ (necessitating longer simulation times), we used biased simulations
at $T=298$ K to determine the resistance to engulfment, $g_{\lambda}$, using
the procedure outlined in section 2D. We obtained $\kappa_{\lambda}$ by
plugging our estimate of $g_{\lambda}$ in Equation 12 along with $\rho=5\times
10^{4}$ mol/m3, $\gamma=3.7\times 10^{-5}$ kJ/m2, and $\Delta\mu=-0.548$
kJ/mol for TIP4P/Ice at $T=298$ K 5, 36, 19. As shown in Figures 10B – E for
sbwAFP, the interfacial curvatures thus calculated, $\kappa_{\lambda}^{\rm
th}$, agree well with those estimated using interfacial shape,
$\kappa_{\lambda}^{\rm int}$ (using the procedure described in section 5) for
each of the four box sizes that we studied. The maximum curvature,
$\kappa_{\rm m}^{\rm th}$, was then plotted against $L^{-2}$ (Figure 10F).

###  7.2 Estimating $p_{\rm m}$ and $a_{\rm m}$

To obtain estimates of the perimeter, $p_{\rm m}$, of the three-phase contact
line corresponding to $\lambda=\lambda_{\rm m}$, and the area, $a_{\rm m}$,
that it encloses, we used the contact line, $\partial\mathcal{I}$, determined
in section 5C (Figure 7D, red). The perimeter, $p_{\rm m}$ was computed by
adding up the lengths of the line segments that make up the contact line,
i.e., $p_{\rm m}=\int_{\partial\mathcal{I}}ds$, whereas the area, $a_{\rm m}$,
was calculated using Green’s theorem, i.e., $a_{\rm
m}=\iint_{\notin\mathcal{I}}dxdy=0.5\int_{\partial\mathcal{I}}(ydx+xdy)$.
Estimates of $p_{\rm m}$ and $a_{\rm m}$ were obtained for each of the four
simulation box sizes used to study sbwAFP, and as seen in Table 1, they agree
well with one another. To minimize the uncertainty in their estimates, values
of $p_{\rm m}$ and $a_{\rm m}$ averaged over the different simulation box
sizes were used to fit the $\kappa_{\rm m}$ vs $L$ data using Equation 15
(Figure 10F).

Table 1: Estimates of $p_{\rm m}$ and $a_{\rm m}$ obtained for the different
simulation box dimensions used to study sbwAFP.

$L_{x}\times L_{y}$ (nm2) $p_{\rm m}$ (nm) $a_{\rm m}$ (nm2) 5.89 x 6.33 13.27
11.8 6.62 x 7.23 12.89 11.23 7.35 x 8.13 13.5 12.03 8.09 x 9.04 12.63 10.71

The best fit to the data was obtained for $\eta_{\rm m}=0.95\pm 0.05$
suggesting that sbwAFP is optimal at resisting engulfment by ice. To estimate
the error in $\eta_{\rm m}$, we randomly chose a $\kappa_{\rm m}$-value for
every AFP separation, $L$, drawing from a Gaussian distribution with its mean
and standard deviation corresponding to $\kappa_{\rm m}^{\rm th}$ and its
error, respectively. The set of $\kappa_{\rm m}$-choices was then fitted to
Equation 15 to obtain an $\eta_{\rm m}$-estimate. This procedure was repeated
1000 times and the standard deviation of the $\eta_{\rm m}$-estimates was
reported as the error in $\eta_{\rm m}$.

Figure 10:  Characterizing the effectiveness of wild-type sbwAFP at resisting
engulfment by ice. (A) Snapshot of sbwAFP (pink) bound to ice (cyan); liquid
water is hidden for clarity. (B-E) The variation of interfacial curvature,
$\kappa_{\lambda}$, with the number of ice-like waters, $\lambda$, is shown
for four different simulation box sizes: $5.88\times 6.33$ nm2 (B),
$6.62\times 7.23$ nm2 (C), $7.35\times 8.13$ nm2 (D), and $8.09\times 9.04$
nm2 (E). Estimates of $\kappa_{\lambda}^{\rm th}$, obtained from theory
(Equation 12, blue) agree well with those obtained from the shape of the ice-
water interface (red). (F) The maximum interfacial curvature, $\kappa_{\rm
m}^{\rm th}$, for each simulation box is plotted against $L^{-2}$ (symbols),
and the data is fit to Equation 15 (line) with $\eta_{\rm m}$ as a fit
parameter ($a_{\rm m}$ and $p_{\rm m}$ were determined from the shape of the
AFP). The simulation data agrees well with Equation 15 with the best fit being
provided by $\eta_{\rm m}=0.95$.

##  8 Data for other AFPs

The analysis described in section 7 was repeated for TmAFP, opAFP and LpAFP,
and the corresponding results are shown in Figures 11, 12 and 13.

Figure 11:  Characterizing the effectiveness of TmAFP at resisting engulfment
by ice. (A) Snapshot of TmAFP (pink) bound to ice (cyan); liquid water is
hidden for clarity. (B-E) The variation of interfacial curvature,
$\kappa_{\lambda}$, with the number of ice-like waters, $\lambda$, is shown
for four different simulation box sizes: $5.88\times 6.33$ nm2 (B),
$6.62\times 7.23$ nm2 (C), $7.35\times 8.13$ nm2 (D), and $8.09\times 9.04$
nm2 (E). Estimates of $\kappa_{\lambda}^{\rm th}$, obtained from theory
(Equation 12, blue) agree well with those obtained from the shape of the ice-
water interface (red). (F) The maximum interfacial curvature, $\kappa_{\rm
m}^{\rm th}$, for each simulation box is plotted against $L^{-2}$ (symbols),
and the data is fit to Equation 15 (line) with $\eta_{\rm m}$ as a fit
parameter ($a_{\rm m}$ and $p_{\rm m}$ were determined from the shape of the
AFP). The simulation data agrees well with Equation 15 with the best fit being
provided by $\eta_{\rm m}=0.88$.  Figure 12:  Characterizing the effectiveness
of opAFP at resisting engulfment by ice. (A) Snapshot of opAFP (pink) bound to
ice (cyan); liquid water is hidden for clarity. (B-E) The variation of
interfacial curvature, $\kappa_{\lambda}$, with the number of ice-like waters,
$\lambda$, is shown for four different simulation box sizes: $5.88\times 6.33$
nm2 (B), $6.62\times 7.23$ nm2 (C), $7.35\times 8.13$ nm2 (D), and $8.09\times
9.04$ nm2 (E). Estimates of $\kappa_{\lambda}^{\rm th}$, obtained from theory
(Equation 12, blue) agree well with those obtained from the shape of the ice-
water interface (red). (F) The maximum interfacial curvature, $\kappa_{\rm
m}^{\rm th}$, for each simulation box is plotted against $L^{-2}$ (symbols),
and the data is fit to Equation 15 (line) with $\eta_{\rm m}$ as a fit
parameter ($a_{\rm m}$ and $p_{\rm m}$ were determined from the shape of the
AFP). The simulation data agrees well with Equation 15 with the best fit being
provided by $\eta_{\rm m}=0.94$.  Figure 13:  Characterizing the effectiveness
of LpAFP at resisting engulfment by ice. (A) Snapshot of LpAFP (pink) bound to
ice (cyan); liquid water is hidden for clarity. (B-E) The variation of
interfacial curvature, $\kappa_{\lambda}$, with the number of ice-like waters,
$\lambda$, is shown for four different simulation box sizes: $5.88\times 6.33$
nm2 (B), $6.62\times 7.23$ nm2 (C), $7.35\times 8.13$ nm2 (D), and $8.09\times
9.04$ nm2 (E). Estimates of $\kappa_{\lambda}^{\rm th}$, obtained from theory
(Equation 12, blue) agree well with those obtained from the shape of the ice-
water interface (red). (F) The maximum interfacial curvature, $\kappa_{\rm
m}^{\rm th}$, for each simulation box is plotted against $L^{-2}$ (symbols),
and the data is fit to Equation 15 (line) with $\eta_{\rm m}$ as a fit
parameter ($a_{\rm m}$ and $p_{\rm m}$ were determined from the shape of the
AFP). The simulation data agrees well with Equation 15 with the best fit being
provided by $\eta_{\rm m}=0.82$.

##  9 Studying an sbwAFP mutant to interrogate how charged residues on the NBS
influence $\eta_{\rm m}$

###  9.1 Preparation of mutant sbwAFP structure

The swapaa command in UCSF Chimera 37 was used to make the following mutations
from charged residues to alanines in sbwAFP: K19A, D28A, D34A, R40A, D45A,
K49A, K50A, D55A, D59A, K60A, R79A, R102A, and K112A. To assess the stability
of the mutated structure, an unbiased simulation with the hydrated sbwAFP
mutant is run at $T=298$ K for 50 ns, as shown in Figure 14. The root mean
square deviation (RMSD) of the heavy atoms was computed over the 50 ns
simulation trajectory (Figure 14B) and remained below 0.13 nm over the
duration of the simulation illustrating the stability of the mutated protein
structure.

Figure 14:  (A) The sbwAFP mutant retains its folded structure after a 50 ns
unbiased simulation. (B) The RMSD of protein backbone atoms with respect to
the structure at $t=0$ is roughly 0.1 nm, highlighting that the folded
structure of the mutated protein closely resembles that of the wild-type.

###  9.2 Estimating $\eta_{\rm m}$ for the mutant sbwAFP

The analysis described in section 7 was repeated for the sbwAFP mutant, and
the corresponding results are shown in Figure 15.

Figure 15:  Characterizing the effectiveness of the sbwAFP mutant at resisting
engulfment by ice. (A) Snapshot of sbwAFP mutant (pink) bound to ice (cyan);
the mutated residues are colored blue, and liquid water is hidden for clarity.
(B-E) The variation of interfacial curvature, $\kappa_{\lambda}$, with the
number of ice-like waters, $\lambda$, is shown for four different simulation
box sizes: $5.88\times 6.33$ nm2 (B), $6.62\times 7.23$ nm2 (C), $7.35\times
8.13$ nm2 (D), and $8.09\times 9.04$ nm2 (E). Estimates of
$\kappa_{\lambda}^{\rm th}$, obtained from theory (Equation 12, blue) agree
well with those obtained from the shape of the ice-water interface (red). (F)
The maximum interfacial curvature, $\kappa_{\rm m}^{\rm th}$, for each
simulation box is plotted against $L^{-2}$ (symbols), and the data is fit to
Equation 15 (line) with $\eta_{\rm m}$ as a fit parameter ($a_{\rm m}$ and
$p_{\rm m}$ were determined from the shape of the AFP). The simulation data
agrees well with Equation 15 with the best fit being provided by $\eta_{\rm
m}=0.8$.

## References

  * Berendsen et al. 1995 Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS: A message-passing parallel molecular dynamics implementation. _Comput. Phys. Commun._ 1995, _91_ , 43–56
  * Bussi et al. 2007 Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. _J. Chem. Phys._ 2007, _126_ , 014101
  * Nosé and Klein 1983 Nosé, S.; Klein, M. L. Constant pressure molecular dynamics for molecular systems. _Mol. Phys._ 1983, _50_ , 1055–1076
  * Parrinello and Rahman 1981 Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. _J. Appl. Phys._ 1981, _52_ , 7182–7190
  * Abascal et al. 2005 Abascal, J. L. F.; Sanz, E.; García Fernández, R.; Vega, C. A potential model for the study of ices and amorphous water: TIP4P/Ice. _J. Chem. Phys._ 2005, _122_ , 234511
  * Hornak et al. 2006 Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of multiple Amber force fields and development of improved protein backbone parameters. _Proteins_ 2006, _65_ , 712–725
  * Darden et al. 1993 Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N$\cdot$log(N) method for Ewald sums in large systems. _J. Chem. Phys._ 1993, _98_ , 10089–10092
  * Miyamoto and Kollman 1992 Miyamoto, S.; Kollman, P. A. Settle: An analytical version of the SHAKE and RATTLE algorithm for rigid water models. _J. Comput. Chem._ 1992, _13_ , 952–962
  * Hess et al. 1997 Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A linear constraint solver for molecular simulations. _J. Comput. Chem._ 1997, _18_ , 1463–1472
  * Hess 2008 Hess, B. P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation. _J. Chem. Theory Comput._ 2008, _4_ , 116–122
  * Humphrey et al. 1996 Humphrey, W.; Dalke, A.; Schulten, K. VMD – Visual Molecular Dynamics. _Journal of Molecular Graphics_ 1996, _14_ , 33–38
  * Leinala et al. 2002 Leinala, E. K.; Davies, P. L.; Doucet, D.; Tyshenko, M. G.; Walker, V. K.; Jia, Z. A $\beta$-Helical Antifreeze Protein Isoform with Increased Activity: STRUCTURAL AND FUNCTIONAL INSIGHTS*. _J. Biol. Chem._ 2002, _277_ , 33349–33352
  * Liou et al. 2000 Liou, Y. C.; Tocilj, A.; Davies, P. L.; Jia, Z. Mimicry of ice structure by surface hydroxyls and water of a beta-helix antifreeze protein. _Nature_ 2000, _406_ , 322–324
  * Antson et al. 2001 Antson, A. A.; Smith, D. J.; Roper, D. I.; Lewis, S.; Caves, L. S.; Verma, C. S.; Buckley, S. L.; Lillford, P. J.; Hubbard, R. E. Understanding the mechanism of ice binding by type III antifreeze proteins. _J. Mol. Biol._ 2001, _305_ , 875–889
  * Middleton et al. 2012 Middleton, A. J.; Marshall, C. B.; Faucher, F.; Bar-Dolev, M.; Braslavsky, I.; Campbell, R. L.; Walker, V. K.; Davies, P. L. Antifreeze protein from freeze-tolerant grass has a beta-roll fold with an irregularly structured ice-binding site. _J. Mol. Biol._ 2012, _416_ , 713–724
  * Matsumoto et al. 2018 Matsumoto, M.; Yagasaki, T.; Tanaka, H. GenIce: Hydrogen-Disordered Ice Generator. _J. Comput. Chem._ 2018, _39_ , 61–64
  * Martínez et al. 2009 Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: a package for building initial configurations for molecular dynamics simulations. _J. Comput. Chem._ 2009, _30_ , 2157–2164
  * Patel et al. 2011 Patel, A. J.; Varilly, P.; Chandler, D.; Garde, S. Quantifying density fluctuations in volumes of all shapes and sizes using indirect umbrella sampling. _J. Stat. Phys._ 2011, _145_ , 265–275
  * Marks et al. 2023 Marks, S.; Vicars, Z.; Thosar, A.; Patel, A. Characterizing Surface Ice-philicity Using Molecular Simulations and Enhanced Sampling. _J. Phys. Chem. B_ 2023, _127_ , 6125–6135
  * Lechner and Dellago 2008 Lechner, W.; Dellago, C. Accurate determination of crystal structures based on averaged local bond order parameters. _J. Chem. Phys._ 2008, _129_ , 114707
  * Patel et al. 2010 Patel, A. J.; Varilly, P.; Chandler, D. Fluctuations of water near extended hydrophobic and hydrophilic surfaces. _J. Phys. Chem. B_ 2010, _114_ , 1632–1637
  * Bonomi et al. 2009 Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. PLUMED: A portable plugin for free-energy calculations with molecular dynamics. _Comput. Phys. Commun._ 2009, _180_ , 1961–1972
  * Tribello et al. 2014 Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. _Comput. Phys. Commun._ 2014, _185_ , 604–613
  * PLUMED consortium 2019 PLUMED consortium, Promoting transparency and reproducibility in enhanced molecular simulations. _Nat. Methods_ 2019, _16_ , 670–673
  * Nguyen and Molinero 2015 Nguyen, A. H.; Molinero, V. Identification of Clathrate Hydrates, Hexagonal Ice, Cubic Ice, and Liquid Water in Simulations: the CHILL+ Algorithm. _J. Phys. Chem. B_ 2015, _119_ , 9369–9376
  * Willard and Chandler 2010 Willard, A. P.; Chandler, D. Instantaneous liquid interfaces. _J. Phys. Chem. B_ 2010, _114_ , 1954–1958
  * Möller and Trumbore 1997 Möller, T.; Trumbore, B. Fast, Minimum Storage Ray-Triangle Intersection. _Journal of Graphics Tools_ 1997, _2_ , 21–28
  * Shirts and Chodera 2008 Shirts, M. R.; Chodera, J. D. Statistically optimal analysis of samples from multiple equilibrium states. _J. Chem. Phys._ 2008, _129_ , 124105
  * Zhu and Hummer 2012 Zhu, F.; Hummer, G. Convergence and error estimation in free energy calculations using the weighted histogram analysis method. _J. Comput. Chem._ 2012, _33_ , 453–465
  * Xi et al. 2018 Xi, E.; Marks, S. M.; Fialoke, S.; Patel, A. J. Sparse sampling of water density fluctuations near liquid-vapor coexistence. _Mol. Simul._ 2018, _44_ , 1124–1135
  * Desbrun et al. 1999 Desbrun, M.; Meyer, M.; Schröder, P.; Barr, A. H. Implicit fairing of irregular meshes using diffusion and curvature flow. Proceedings of the 26th annual conference on Computer graphics and interactive techniques. USA, 1999; pp 317–324
  * Cazals and Pouget 2005 Cazals, F.; Pouget, M. Estimating differential quantities using polynomial fitting of osculating jets. _Comput. Aided Geom. Des._ 2005, _22_ , 121–146
  * The CGAL Project 2023 The CGAL Project, _CGAL User and Reference Manual_ , 5.5.2 ed.; CGAL Editorial Board, 2023
  * Lobaton and Salamon 2007 Lobaton, E. J.; Salamon, T. R. Computation of constant mean curvature surfaces: Application to the gas-liquid interface of a pressurized fluid on a superhydrophobic surface. _J. Colloid Interface Sci._ 2007, _314_ , 184–198
  * Ickes et al. 2015 Ickes, L.; Welti, A.; Hoose, C.; Lohmann, U. Classical nucleation theory of homogeneous freezing of water: thermodynamic and kinetic parameters. _Phys. Chem. Chem. Phys._ 2015, _17_ , 5514–5537
  * Espinosa et al. 2016 Espinosa, J. R.; Navarro, C.; Sanz, E.; Valeriani, C.; Vega, C. On the time required to freeze water. _J. Chem. Phys._ 2016, _145_ , 211922
  * Pettersen et al. 2004 Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. UCSF Chimera–a visualization system for exploratory research and analysis. _J. Comput. Chem._ 2004, _25_ , 1605–1612