A first-principles examination of the ice–cellulose interface: towards bioinspired antifreeze design

Abstract

Examination of binding of cellulose to ice using ab initio modeling reveals that new C–O bonds are formed on the basal ice surfaces, where some of the O atoms are exposed at the surface due to missing H bonds. Further analysis suggests that the cellulose unit binds in such a way as to form a tetrahedral arrangement at the ice surface, evidenced by a geometric measure of tetrahedrality. This hypothesis is further validated for both primary and secondary prismatic planes. This leads us to conclude that in the case of cellulose molecules, binding at ice is dependent on preserving its tetrahedral bonding arrangement. Our findings suggest that the idea of tetrahedrality is very widely applicable to coordination ranging from water to ice-binding proteins, highlighting a design criterion for novel ice-binding/antifreeze proteins/materials.

---

Design, System, Application

Antifreeze proteins (AFPs), found in cold-adapted organisms such as arctic fish and bacteria, have the ability to bind to ice and inhibit its growth. These proteins can serve as a source of bio-inspiration for developing new antifreeze materials aimed at preventing infrastructure damage caused by freeze–thaw cycles. This is of paramount importance since the AFPs derived from these organisms are not suited for large-scale applications due to their limited availability, synthesis challenges, and high production costs. A promising alternative lies in the use of cellulosic biopolymers, which are abundant, environment friendly, and possess water retention similar to AFPs. Ab initio modeling offers a pathway to explore this bioinspired design strategy by providing a molecular-level understanding of the interactions of cellulose with ice, and more importantly whether ice-binding is possible. In this work, we provide a first of its kind first-principles probe into the nature of the interaction of cellulose with the ice surface by considering various basal and prismatic planes. We identify a key geometrical descriptor that can provide a qualitative and quantitative picture of the cellulose–ice binding leading to a framework that could be utilized to design sustainable antifreeze materials.

1 Introduction

The freezing of water and melting of ice are among the most common natural phenomena known to us, playing a crucial role in shaping landscapes, influencing soil stability¹ and infrastructure durability.² Freeze–thaw cycles, where water infiltrates soils and structures, freezes, expands, and then melts, exert significant mechanical stresses. These repeated cycles can lead to soil destabilization,³ increased erosion,⁴ and degradation of roads, bridges, and buildings. In colder regions, extensive damage to infrastructure occurs annually due to freeze–thaw-induced cracking and foundation instability. The associated repair and maintenance costs exceed billions of dollars each year and are projected to increase as temperature fluctuations become more extreme.^5,6 One key approach to mitigating the detrimental effects of freeze–thaw cycles is the development of advanced materials designed to control ice formation and reduce structural damage. Novel antifreeze⁷ and ice-nucleating⁸ materials can help regulate phase changes, preventing excessive ice expansion in soils and construction materials. Designing such materials is crucial for enhancing the longevity of infrastructure and maintaining the stability of built environments in regions affected by severe freeze–thaw conditions.

In nature, several organisms such as bacteria⁹ and Arctic fish¹⁰ survive extreme cold by preventing ice formation even under sub-zero conditions. They achieve this by producing protein molecules that bind to specific planes of ice crystals, thereby inhibiting their growth. Understanding these interactions has been a key focus in the design of novel antifreeze materials.^11,12 Building on this foundation, our work aims to uncover the principles behind designing novel antifreeze biopolymers that mimic this natural behavior. Typically, antifreeze glycoproteins (AFGPs)¹³ consist of sugar derivatives such as polysachharide units connected by glycosidic bonds. In general, a mechanistic understanding of antifreeze proteins (AFPs) remains far from being completely understood, although various models have been put forth to explain how the ice-growth is inhibited.¹⁴

The widespread use of AFGPs is limited by challenges in their synthesis and issues related to scalability.¹³ To overcome these limitations, the design of AFGP analogs must incorporate both the use of naturally abundant materials and environmentally friendly processes. Sustainability—guided by the principles of the circular economy¹⁵—is a key consideration in this approach, emphasizing the use of materials that can be safely recycled back into the environment.¹⁵ One possible candidate, cellulose, (C₆H₁₀O₅)_n, containing the β (1–4) glycosidic bond between their glucopyranose rings, has recently shown to possess ice recrystallization inhibition¹⁶ as well as ice nucleation activity.¹⁷ As a potential antifreeze material design strategy promoting sustainability, we therefore focus on cellulose, which is the most abundant biopolymer on the planet, can be recovered from agricultural waste, and is a great prototype for the circular economy¹⁸ for various applications ranging from biological¹⁹ to electronics.²⁰

In this work, we focus on the interaction of cellulose with ice to gain a fundamental understanding of how cellulosic materials could potentially be designed to act as ice nucleating/antifreeze agents. Ice, the solid phase of water, is one of the most widely studied materials, and even to this date, an understanding of its structure is far from complete, with new ice phases being reported by researchers regularly.²² However, hexagonal ice (I_h, sometimes also denoted as I_h), remains the most widely observed and studied. Therefore, we focus on investigating the binding of cellulose with hexagonal ice.

Previous ab initio modeling studies have focused on water adsorption on bulk cellulose substrates,²³ but to the best of our knowledge the behavior of cellulose adsorption on ice/water has not been given much attention. Studies of pure ice surfaces²⁴ as well as metal and inorganic molecules such as NH₃ (ref. 25) adsorption on the ice surface have been carried out and also point to the importance of proton order/disorder at the ice surface as investigated by Buch et al.²⁶ While the amphiphilic nature of cellulose has been discussed in the literature for some time,²⁷ Alqus et al.²⁸ conducted a thorough investigation of the amphiphilicity of cellulose and its binding to water in the presence of graphene using molecular dynamics simulations. They found that the –CH– abundant face (100) remained relatively intact with their interface with graphene, while a significant rearrangement was observed for the –OH dominant face (010), such that the reoriented cellulose again presented –CH– groups to the graphene sheet, strongly outlining the partly hydrophobic nature of cellulose and providing an indirect evidence of its amphiphilicity. Recently, Bruel et al.²⁹ have experimentally validated the amphiphilic nature of cellulose to be connected to its anisotropy in the case of cellulose nanocrystals (CNCs). Nanocellulose have already shown to be active ice recrystallization inhibition (IRI) agents¹⁶ with potential applications as cryoprotectants in food (e.g. ice cream³⁰) industry, as well as tissue engineering.³¹

While there is enough evidence about the potential of cellulose to inhibit ice-growth, the exact mechanism of its antifreeze behavior is very poorly understood. Nanocrystalline cellulose containing beta glycosidic units have –OH groups attached to the pyranose ring similar to commonly used and well-studied antifreeze agents, e.g. polyvinyl alcohol (PVA). PVA is known to recognize and bind to the prismatic ice planes inhibiting ice recrystallization.^32,33 As a first step to explore the antifreeze behavior of cellulose, we present a first study of cellulose binding to ice-surfaces. Through this work, we also aim to develop a preliminary framework to explore the binding characteristics of antifreeze/ice-nucleating agents on ice surfaces towards building a model that would include hydration in subsequent work.

Cellulose exists in nature in its most stable form Iβ, which is composed of H-bonded chains of D-glucopyranose units. In cellulose Iβ, there are two distinct chains, each of them forming a sheet of parallel chains (Fig. 1a), that stack together to give the 3D monoclinic crystal structure, space group P112₁. The other allotrope, Iα consists of one unique chain that stacks to form the 3D triclinic crystal structure.

Fig. 1 (a) Relaxed cellulose Iβ unit cell, and (b) the H-bond network present in the cellulose sheet and the cellulose unit (right) extracted and the atomic labeling shown as in Wada et al.²¹ C, H, and O atoms are shown in brown, pink, and red colors as presented in the legend in the top right.

In studies on ice, the question of proton order/disorder has been often addressed in several modeling studies in the last two decades.^24,26 Pan et al.²⁴ elucidated that the distribution of the energetics of the various proton/disordered ice surfaces spanned more than an order of magnitude greater than that of the bulk ice, due to the sensitivity of the charge–charge interactions of the dangling –OH bonds, which was also argued by Fletcher.³⁹ In this work, we extend the question to the role of adsorbents, and ask: are the antifreeze properties related to the dangling –OH bonds on the ice surface? Are van der Waals (vdW) interactions the most important in deciding the nature of binding of such antifreeze agents? We find that besides vdW bonding due to new H bonds between cellulose and ice, new non-covalent ionic C–O bonds are equally important in deciding the strong hydrophilic binding between cellulose and ice surfaces. We also note that excess protons on the ice surface impede this bond formation, and also lead to rejection of protons, and in some cases water molecules, from the ice surface.

We focus on whether the concentration of protons at the surface impacts the binding of cellulose to the ice basal surface. For this investigation, we include the Fletcher surface,³⁹ which is the striped surface (Fig. 2a) that is the most stable until a critical temperature (30 K for basal planes), when the disordered surface becomes more stable. In addition, we also consider the fully hydroxylated basal surface where all the H bonds of the terminating O atoms are satisfied, to compare it to the pristine basal surface, with O atoms missing H bonds in the out-of-plane direction. Finally, we also discuss the binding of cellulose to prismatic planes, and compare it to the observed behavior for the basal plane. We generalize our findings to argue that cellulose can serve as an antifreeze agent, and being a naturally occurring polymer, can find critical applications in sustainable solutions to thwart the infrastructure damage posed by irregular freeze–thaw cycles across the Arctic region.

Fig. 2 Arrangement of H bonds in the terminating bilayer present on the basal surface of (a) fletcher phase, (b) pristine ice surface, as well as with (c) cellulose adsorbed on the pristine surface, and (d) the fully hydroxylated ice surface used in this work.

It has been established that the microscopic structural properties of water are connected to its second shell.⁴⁰ Russo et al.⁴¹ found that the local translation order of the second shell can be captured by considering the directional intermolecular H bonds, suggesting the tetrahedral ordering of water to be central to this structural characterization. More recently Klatt et al.⁴² demonstrated that the local density fluctuations in liquid water and its higher order moments, all detect tetrahedrality, which essentially is a many-body property. Given the role of tetrahedrality in describing water properties and anomalies, we apply this idea to explore the binding of cellulose to ice surfaces. We show that tetrahedrality can not only detect this binding process but can also be used to quantify the binding energies for ice surfaces.

2 Computational methods

All density functional theory (DFT) calculations were carried out using a plane-wave basis set within the Projector-Augmented Wave (PAW) approach^43,44 with the exchange-correlation treated at the generalized gradient approximation level of Perdew–Becke–Ernzehof (GGA-PBE)⁴⁵ as implemented in the Vienna Ab initio Simulation Package (VASP 5.4.4),^46–49 except for the initial optimization of cellulose Iβ structure. The geometry optimization of cellulose Iβ taken from Nishiyama et al.³⁴ was first carried out using the Quantum Espresso 6.8 (ref. 50 and 51) DFT package using norm-conserving pseudopotentials⁵² and GGA-PBE exchange-correlation with a kinetic energy cut-off of 90 Ry for the plane-wave expansion and a Γ-centered Monkhorst–Pack⁵³ K-grid of 8 × 8 × 8 until the total energy was converged to within 0.01 meV per atom and the forces on the atoms were within 0.02 eV Å⁻¹. The obtained structure was subsequently relaxed using VASP with a kinetic energy cut-off of 900 eV, and a Γ-centered Monkhorst–Pack K-grid of 8 × 8 × 8, with a convergence criteria of total energy change for the self-consistency (SCF) cycle as 10⁻⁵ eV per supercell and the norm of the forces on each atom being less than 0.01 eV Å⁻¹. All calculations in this work used the D3 correction with Becke–Johnson⁵⁴ damping to incorporate the dispersion interactions that are important for ice and cellulose.

For the hexagonal ice (I_h), we use a 16 water molecule orthorhombic supercell as presented by Morrison et al.³⁷ We choose this supercell to be able to include the inherent proton disorder in bulk ice I_h using a relatively larger cell that has a zero net dipole moment. This structure also presents the terminating H-atoms in the ‘up’ position, that has been experimentally detected to be the case for the basal plane of ice I_h using vibrational spectroscopy,⁵⁵ and most recently by direct visualization of the basal (0001) surface of ice I_h using atomic-resolution imaging.⁵⁶ The kinetic energy cut-off of the plane-wave expansion was 900 eV. For the relaxation of the supercell geometry, a Γ-centered Monkhorst–Pack grid of 8 × 8 × 8 k-points was used. The optimization was continued until the total energy change of the supercell for the SCF cycle was less than 10⁻⁵ eV and the norm of the forces on each atom was less than 0.01 eV Å⁻¹. The ice–cellulose interfaces were optimized at constant volume with the cellulose unit on the relaxed ice surfaces having a vacuum length of 25 Å normal to the interface, with an energy cut-off of 900 eV with the convergence criteria of total energy change of SCF cycle at 10⁻⁵ eV and the norm of the forces on each atom less than 0.01 eV Å⁻¹ at the Γ K-point for the evaluation of the low-energy configurations for all cases. The various basal and prismatic surfaces were constructed by a 2 × 2 repetition of the optimized ice I_h supercell in the in-plane direction and containing 3–4 bilayers (section 3.1.1) in the out-of-plane direction to the surface with vacuum. We used the fletcher surface from Engel et al.⁵⁷ to create the ice interface with cellulose.

2.1 Cellulose model

In this work, we represent cellulose by a two D-glucopyranose unit (C₆H₁₀O₅)₂ structure containing the β-glycosidic linkage extracted from a cellulose chain (Fig. 1b) of the relaxed cellulose Iβ unit cell. These two D-glucopyranose rings are the smallest repeating units of cellulose and combine to form longer chains, as well as sheets to assemble into fibrils, and even crystalline phases. Most importantly, during treatment of cellulose, such as pyrolysis, it is widely accepted in the literature that the breakdown of cellulose can proceed with the breakage of either the C₁–O or the C₄–O bonds^58,59 along with a more thermodynamically favored concerted intermediate mechanism. In fact, it was shown that this activation of the glycosidic bond is probable even at low-temperatures using a catalyzed mechanism involving nearby –OH groups.⁶⁰ This ‘active cellulose’ unit that can be periodically repeated to create chains, sheets, and crystalline cellulose (Iα and Iβ), is missing a usual O atom at the C₁ site (Fig. 1b, right). Recently, it was also pointed out that the partially charged C and O atoms in cellulose might be responsible for its several properties.⁶¹ Keeping all such ideas into consideration in the present work, we explore the implications of such a reactive cellulose unit binding on ice surfaces.

2.2 Binding energy

The ice–cellulose interfaces were generated by placing the cellulose unit on the basal, the primary prismatic, and the secondary prismatic surfaces of ice I_h, for different translations of cellulose (schematic of the supercell shown in Fig. 2c). For all the interfacial calculations, the supercell was held at a constant volume with a vacuum length of 25 Å, and all the atoms were allowed to relax. The binding energy of cellulose on the ice-surface can be calculated using the following relationship:

Eⁱ_hetero = E_ice + E_cellulose + Eⁱ_bind (1)

where Eⁱ_hetero, E_ice, and E_cellulose are the total energies of an interfacial configuration i (with vacuum), hexagonal ice (with vacuum), and our cellulose unit (C₆H₁₀O₅)₂ (in vacuum) respectively, and Eⁱ_bind is the binding energy of cellulose on the ice surface for the particular configuration i.

Similarly, we can compute the binding energy of cellulose for a different configuration j as,

E^j_hetero = E_ice + E_cellulose + E^j_bind (2)

It can be seen from eqn (1) and (2) that the relative binding energy between two positions i and j can be expressed as

Eⁱ_hetero − E^j_hetero = Eⁱ_bind − E^j_bind (3)

In other words, the total energy difference between two configurations of cellulose on an ice surface is equal to their binding energy difference.

3 Results and discussion

There has been considerable attention devoted to the structure of ice (0001) surfaces for over two centuries. Since our focus is to understand the binding of cellulose to different ice planes, we focus on the relative binding of cellulose on the basal and prismatic planes for our ice model in this work. However, to ensure we take into account the role of unterminated-O at the surface (denoted as dangling-O by Buch²⁶), or non H-bond forming H; i.e., a hydroxylated surface (denoted as d-H similarly by Buch²⁶), we include the case of the basal (0001) surface when all the terminating O have 4 bonds (by adding a H atom above the O atom), as shown in Fig. 2d. This compared to the case, when some of the terminating O lack one H bond (partially hydroxylated, Fig. 2b) can elucidate the role of surface protons in the binding of absorbents as was determined in the cases of NH₃ (ref. 62) and CH₃OH.⁶³

3.1 Bulk ice and cellulose

The optimized lattice parameters of cellulose are given in Table 1, while those of the ice Ih model chosen in this work are provided in Table 2. Glycosidic torsional angles, Φ, Ψ, χ for the cellulose chains are the dihedral angles as described in Table 1. The PBE-D3 vdW correction is shown to improve the predictions of lattice parameters as compared to D2 and D corrections.

Table 1 Lattice parameters of cellulose Iβ

	PBE-D (Li^23 2011)	PBE-D2 (Bučko^35 2011)	PBE-D2 (Kubicki^36 2013)	PBE-D3 (this work)	Experiment^34
a Φ is defined as O5′–C1′–O1–C4. b Ψ as C1′–O1–C4–C5. c χ as O6–C6′–C5′–O5′. d χ′ as O6–C6′–C5′–C4′.
a (intersheet)(Å)	7.85	7.57	7.55	7.82	7.78
b (interchain)(Å)	8.18	8.14	8.14	8.16	8.20
c (intrachain)(Å)	10.47	10.39	10.40	10.41	10.38
γ(°)	96.6	96.5	96.4	96.6	96.5
Φ(°)^a	95, 95	93, 94	93, 94	94, 94	99, 89
Ψ(°)^b	144, 143	143, 145	143, 145	143, 144	142, 147
χ(°)^c	164, 165	168, 165	—	167, 165	170, 158
χ′(°)^d	78, 77	74, 77	—	75, 76	70, 83

---

Table 2 Optimized lattice parameters of ice I_h orthorhombic cell taken from ref. 37 based on the Bernal–Fowler rules³⁸

	PBE-D3
a (Å)	8.70
b (Å)	7.53
c (Å)	7.11

---

3.1.1 Ice surface energy. The surface energy for basal and prismatic surfaces of hexagonal ice is computed using the following relation:where E^slab_n and E^bulk are the total energies of the slab (with n bilayers), and bulk ice (per atom) respectively, and N^slab_n is the total number of atoms in the slab (with n bilayers). A is the surface area of the slab.

From Fig. 3, it can be seen that the surface energy is rapidly converged within 2–3 bilayers for all ice surfaces, and increasing the number of bilayers does not provide an improvement in the accuracy of the surface energy. The ice basal (0001) surface energy is around 14 meV Å⁻², within the predicted range of 13–15 meV Å⁻² by other PBE results in the literature.²⁴ The slight oscillation seen in the case of the prismatic surfaces around 16–17 meV Å⁻² is due to the different configurations of the H atoms on the two surfaces for different number of bilayers in the ice slab. We note that this is an extremely small variation. Therefore, we use 3 bilayers of basal and primary prismatic planes to construct the adsorbate for the cellulose unit in the present work. In the case of the secondary prismatic surface, the odd and the even number of bilayers lead to a switching of orientations of the surface H atoms such that in the first case (odd number of bilayers) some O atoms do not have out-of-plane H atoms, while in the second case (even number of bilayers) all O atoms have some out-of-plane H atoms. Therefore, for the secondary prismatic ice slab, we choose 4 bilayers to ensure both the surfaces in the supercell have similarly oriented H atoms in both in-plane and out-of-plane directions.

Fig. 3 Surface energy variation for different number of bilayers in the model.

3.2 Ice–cellulose interface

3.2.1 Basal planes. We compare the relative binding energy of cellulose for the (0001) surface using eqn (3). To examine the dependence of this binding energy on the surface H distribution, we consider two additional cases, shown in Fig. 2. First, we consider the fletcher surface³⁹ (Fig. 2a), where the O and the H atoms form alternating rows. Theoretically, the fletcher surface is the most stable basal surface until at ∼32 K, when the ordered arrangement is lost. Second, we consider another case, where we fully hydroxylate the terminating basal surface by adding additional H atoms for the missing O atoms on the pristine basal surface as a result of the cleavage of the plane. We create contour maps of the binding energy change for different configurations of cellulose on the ice surfaces (Fig. 4a) and plot this for all cases in Fig. 4b.

Fig. 4 (a) Schematic of various configurations spanned by cellulose translations on the (0001) plane; relative binding energy plot for b) Fletcher, c) pristine, and d) H excess (0001) ice surfaces; and e) scatter plot of all relative binding energies of cellulose on the (0001) basal surface.

The binding energy change for the three basal surface configurations are shown in Fig. 4c–e, where the darker regions suggest greater binding between the cellulose and the ice surface. The configurations with the highest relative binding energies for various translations for each of the ice surfaces are annotated by numbers. The range of the binding energy difference is ∼3 eV for the fletcher (c) and the pristine (d) basal surfaces, whereas it spans about 20 eV for the fully hydroxylated (e) surface as depicted in Fig. 4b. We find that the structures with the highest relative binding energy (structures 1–6) all show the formation of a C₁(cellulose)–O (ice) bond at the interface (Fig. 5a), resulting in a configuration that is ∼3 eV lower in energy than the other relaxed configurations of cellulose for different translations on the ice surface. Furthermore, we find that this C–O bond formation facilitates the formation of a tetrahedral geometry at the cellulose–ice interface as portrayed in Fig. 5b. As water and ice both have tetrahedral geometry, we postulate that the affinity of cellulose to the ice surface is associated with its ability to retain the tetrahedral coordination of ice by forming this bond. Electronic density of states (DOS) of one of the preferred binding configurations (structure 5) shows that the cellulose unit dopes electrons uniformly in the valence band region, including right at the valence band edge (Fig. S1).

Fig. 5 (a) Binding of cellulose unit on the basal surface (structure 5) before and after the energy minimization presenting the change of the cellulose shape as it binds to the surface, (b) top view of (a) showing the tetrahedral geometry formed as a result of this binding.

Upon further examination of configurations of the highest relative binding energy for the latter case (structures 7 and 8), we find that the dissociation of the surface layers/cellulose in the form of H atoms and molecules (Fig. S2) takes place. It is important to note that Buch²⁶ had also pointed out that an all proton terminated ice I_h surface should be unlikely based on their molecular dynamics simulations, and our findings suggest the same to hold even in the case of a cellulose molecule in the process of binding to it.

3.2.2 Role of hydroxylation. Lee et al.⁶⁴ studied the binding of βD-glucose and cellobiose on the hydroxylated kaolonite (Si₂Al₂O₅(OH)₄) (0001) surface and found that the binding is primarily caused by the formation of new H-bonds during the attachment of the adsorbed cellulose on the surface. They found that the hydroxylated surface promoted greater binding than the siloxane surface due to the formation of a greater number of H bonds with cellulose.

This explains the binding of cellulose unit to the fletcher and the pristine (partially hydroxylated) basal surface, where new H-bonds are seen. Specifically, we note the formation of 4–5 H bonds in each of their most stable configurations, where 3 H bonds are formed with the O atoms from cellulose, while 1–2 H-bonds form with the O atoms present on the ice surface (Table S1). While we observe the formation of new H-bonds between cellulose and ice surface in all our basal surfaces, it is noted that the hydroxylated surface do not lead to the stable cellulose–ice bound configurations, as they prohibit the formation of tetrahedral coordination between the adsorbed cellulose and the ice surface, as observed for the fletcher and the pristine surfaces, where the presence of O atoms (lacking a H bond above) at the surface enables the formation of this geometry (Fig. 5).

On the contrary, in the case of fully hydroxylated ice Ih (0001), we observe dissociation of this surface which suggests that it is thermodynamically more favorable for cellulose to bind to the partially hydroxylated ice surfaces than to the fully hydroxylated one where no unbonded O atoms exist. Thus it can be concluded that the binding of cellulose is dependent on both a) the ability to form new H-bonds and b) the formation of tetrahedral coordination via the C₁(cellulose)–O(ice) bond when possible (Fig. 5). It should be pointed out that the same binding mechanisms have been demonstrated to be active in the case of several AFPs.⁶⁵

3.2.3 Prismatic planes. We also extend our investigation to the primary and the secondary prismatic planes as shown in Fig. 6. Fig. 6a and b show the top view of the primary and secondary prism planes with the inset showing the side views of their corresponding terminating layers. The relative binding energy following the same approach as for the basal planes are depicted in Fig. 6c for the primary prismatic (101̄0) and Fig. 6d for the secondary prismatic (1̄21̄0) planes. For both prismatic planes, we again observe that the configurations associated with the largest relative binding of cellulose show the formation of a new C–O bond at the interface. Furthermore, as in the case of the basal planes, we note the formation of 5–6 H bonds (Table S1). The range of the binding energy difference is of the same order as in the case of basal surfaces, suggesting that cellulose might possess a relatively equal degree of affinity with the prismatic ice surfaces.

Fig. 6 a) Primary prismatic ice surface with inset showing the side view of the terminating layer with grey (out) and yellow (in) O atoms; b) secondary prismatic ice surface with inset showing the side view of the terminating layer with grey O atoms; relative binding energy plots for cellulose binding on c) primary prismatic; and d) secondary prismatic ice surfaces.

However, a closer examination of the preferred binding configurations (structures 9–12) reveals that the formation of this new C–O bond is accompanied by a release of a water molecule (Fig. S3) that was not observed in the case of basal planes. One plausible reason for this observation could be a different mechanism of ice-binding where some water molecules are displaced to aid the process, as in the case for the AFP in the bacterium Marinomonas primoryensis (MpAFP) binding on the primary prism planes shown by Davies.⁶⁶ Another possible reason for this phenomenon may be the destabilization of the ice surface by the hydrophobic part of cellulose, and therefore attributed to the greater exhibition of the amphiphilic nature of cellulose on prismatic than basal planes. It is interesting to note that for proteins that bind to multiple planes, slight changes in their ice-binding characteristics can be consequential in ice-shaping as seen to be the case for MpAFP that binds to both basal and prismatic planes,⁶⁷ and ice morphologies resulting from different ice growth kinetics on the basal and the prismatic planes as was noted in the case of Fragilariopsis cylindrus.⁶⁸ Therefore, we believe cellulose binding to ice should lead to varying ice-growth patterns depending on the proportion of binding on the basal and the prismatic planes. Importantly, our calculations suggest that cellulose binding on both basal and prismatic planes are energetically favored, with the mechanism of ice-binding appearing to be more complex for the prismatic planes.

3.3 Degree of tetrahedrality

The tetrahedral arrangement of molecules in water and ice is the key to their properties. As demonstrated in ref. 42, the signatures of the tetrahedral arrangement of water can be used to understand and even predict water properties at higher temperatures.

We observe that a tetrahedral coordination is formed around the interfacial C–O bond at the interface between the C and H atoms of cellulose and O and H atoms of ice as shown in Fig. 5b. To quantify this binding characteristic, we compute the degree of tetrahedrality q, as discussed by Errington and Debenedetti⁶⁹ originally based on the parameter introduced by Chau and Hardwick.⁷⁰ The degree of tetrahedrality is expressed as

where the value of q corresponds to 0 for an ideal gas and 1 for a perfect tetrahedral coordination. In this work, we compute q for all the lowest energy configurations of cellulose–ice interfaces where a C–O bond is observed such that cos Ψ_jk corresponds to the angles around the C atom with its neighbors j and k. For all these configurations, q is almost 1 suggesting a near-tetrahedron geometry around the cellulose–ice interface. For comparison, we also relaxed 4, 6, and 8 water molecules on the (0001) surface and noticed the average q to be very close to 1. These are presented in Fig. 7. Interestingly, Nutt and Smith⁷¹ had noted that the hydration water near the ice-binding face of the antifreeze protein of the spruce budworm, Choristoneura fumiferana can rearrange to enhance the local tetrahedrality to facilitate the ice-binding.

Fig. 7 Computed values of the degree of tetrahedrality q for the configurations with the highest relative binding energy from Fig. 4 and 6.

In several AFPs found in cold-adapted organisms, anchored clathrate motifs were determined to form between the AFP and the ice-binding surfaces in the cases of the bacteria Marinomonas primoryensis,⁹ and the insect Tenebrio molitor.^72,73 Furthermore, it has been established that the hydration water at the ice binding surface of antifreeze protein has an enhanced tetrahedral geometry that has also been seen in the hydration water of hydrophobic groups of alcohols.⁷⁴ In the case of the antifreeze protein, however, the tetrahedral character is much more pronounced than in the case of alcohols and, is ice-like.⁷⁵ In the case of the AFP found in ocean pout, the tetrahedral water cluster at the hydrophobic region of the ice-binding surface was experimentally determined by joint neutron and X-ray diffraction,⁷⁶ further supporting the principal argument of tetrahedrality. From an energetic analysis, however, it is difficult to establish this mechanism as ab initio calculations employed in the present work are limited by the system size that they can explicitly model thus prohibiting us from including water molecules in our model. This is primarily due to the fact that any atomistic model that includes water molecules needs to carefully account for the statistics associated with the configuration of the water layer. We do aim to fill this gap in our next work where we are developing a model that accounts for the entropy of water molecules on top of the ice surface. In the view of these established results, we strongly suggest that the enhanced tetrahedral geometry observed at the cellulose–ice interfaces when the hydrophobic face of cellulose presents itself to the various ice planes, is an indirect evidence of the favored tetrahedral coordination of the hydration/interfacial water at the AFP–ice interface, especially prominent in the case of cellulose binding to prismatic planes leading to displacement of a water molecule. Therefore, our results are a first ab initio validation of this phenomenon.

Conclusions

In this work, cellulose binding to ice surfaces has been investigated in detail by considering the basal, primary prismatic, and secondary prismatic planes. The relative binding energies for each ice surface with cellulose suggests that cellulose likes to bind to the ice surface forming a tetrahedral geometry. For the case of the basal plane, we also assess the role of surface hydroxylation, and find that the cellulose–ice interface having partial H and partial O is preferred over the fully hydroxylated case. Our primary focus of the present work was to understand the initial binding of cellulose to ice at the molecular design level. We note that modeling the actual cellulose structure such as nanocellulose, and features such as the chirality of the larger structure, will require more coarse-grained length scales to be computationally accessible. We do plan, in future work, to develop such models, informed by the results of our present model, to understand these larger scale effects.

Our findings pave a way towards the design principles of cellulose based antifreeze materials inspired by arctic organisms, and suggest that cellulose binds to ice in order to preserve its tetrahedral bonding (ice-like) pattern. As cellulose also exhibits a dominant network of H-bonds, as in the case of water and ice, it explains the natural affinity between ice and cellulose to bind that extends the tetrahedral bonding in ice; even though cellulose is considered an amphiphilic material possessing both hydrophilic and hydrophobic groups. As hydrophobic components are commonly found in the AFPs/AFGPs that bind to ice, we believe this to be the case in cellulose and propose that the degree of tetrahedrality could be a critical design parameter for sustainable antifreeze materials.

Author contributions

Credit: Aakash Kumar: conceptualization, methodology, software, validation, formal analysis, investigation, visualization, writing – original draft, writing – review and editing; Dilip Gersappe: conceptualization, supervision, project administration, funding acquisition, resources, writing – review and editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI).

Supplementary information: density of states of the ice surface without and with cellulose adsorbed, relaxed ice–cellulose structures showing dissociation of molecules on the basal ice plane, ice–cellulose structures for prismatic plane, H-bond analysis. See DOI: https://doi.org/10.1039/d5me00137d.

Acknowledgements

The authors thank Prof. Benjamin Hsiao from Stony Brook University for helpful discussions. All DFT calculations were performed on the Seawulf cluster of the Institute for Advanced Computational Science (IACS) at Stony Brook University. This material is based upon work supported by the Broad Agency Announcement Program and the Cold Regions Research and Engineering Laboratory (ERDC-CRREL) under Contract No. W913E522C0001. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the Broad Agency Announcement Program and ERDC-CRREL.

References

1. E. Wang, R. M. Cruse, X. Chen and A. Daigh, Can. J. Soil Sci., 2012, 92, 529–536.
2. R. Wang, Z. Hu, Y. Li, K. Wang and H. Zhang, Constr. Build. Mater., 2022, 321, 126371.
3. E. C. Rooney, V. L. Bailey, K. F. Patel, A. R. Possinger, A. C. Gallo, M. Bergmann, M. SanClements and R. A. Lybrand, J. Geophys. Res.:Biogeosci., 2022, 127, e2022JG006889.
4. M. Deprez, T. De Kock, G. De Schutter and V. Cnudde, Earth-Sci. Rev., 2020, 203, 103143.
5. D. A. Streletskiy, S. Clemens, J.-P. Lanckman and N. I. Shiklomanov, Environ. Res. Lett., 2023, 18, 015006.
6. E. Manos, C. Witharana and A. K. Liljedahl, Commun. Earth Environ., 2025, 6, 221.
7. Y. Gao, C. Han and J. Wang, Mol. Syst. Des. Eng., 2025, 692–721.
8. S. Bogler and N. Borduas-Dedekind, Atmos. Chem. Phys., 2020, 20, 14509–14522.
9. C. P. Garnham, R. L. Campbell and P. L. Davies, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 7363–7367.
10. L. Chen, A. L. DeVries and C.-H. C. Cheng, Proc. Natl. Acad. Sci. U. S. A., 1997, 94, 3817–3822.
11. I. Voets, Soft Matter, 2017, 13, 4808–4823.
12. L. L. C. Olijve, K. Meister, A. L. DeVries, J. G. Duman, S. Guo, H. J. Bakker and I. K. Voets, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 3740–3745.
13. A. K. Balcerzak, C. J. Capicciotti, J. G. Briard and R. N. Ben, RSC Adv., 2014, 4, 42682–42696.
14. M. Schauperl, M. Podewitz, T. S. Ortner, F. Waibl, A. Thoeny, T. Loerting and K. R. Liedl, Sci. Rep., 2017, 7, 11901.
15. P. Morseletto, Resour., Conserv. Recycl., 2020, 153, 104553.
16. T. Li, Y. Zhao, Q. Zhong and T. Wu, Biomacromolecules, 2019, 20, 1667–1674.
17. N. Hiranuma, K. Adachi, D. M. Bell, F. Belosi, H. Beydoun, B. Bhaduri, H. Bingemer, C. Budke, H.-C. Clemen and F. Conen, et al. , Atmos. Chem. Phys., 2019, 19, 4823–4849.
18. K. I. Johnson, G. Ilacas, R. Das, H.-Y. Chang, P. R. Sharma, C. O. Dimkpa and B. S. Hsiao, Sustainability Sci. Technol., 2024, 1, 014001.
19. T. Abitbol, A. Rivkin, Y. Cao, Y. Nevo, E. Abraham, T. Ben-Shalom, S. Lapidot and O. Shoseyov, Curr. Opin. Biotechnol., 2016, 39, 76–88.
20. M. Giese and M. Spengler, Mol. Syst. Des. Eng., 2019, 4, 29–48.
21. M. Wada, Y. Nishiyama, H. Chanzy, T. Forsyth and P. Langan, Powder Diffr., 2008, 23, 92–95.
22. T. C. Hansen, Nat. Commun., 2021, 12, 3161.
23. Y. Li, M. Lin and J. W. Davenport, J. Phys. Chem. C, 2011, 115, 11533–11539.
24. D. Pan, L.-M. Liu, G. A. Tribello, B. Slater, A. Michaelides and E. Wang, J. Phys.: Condens. Matter, 2010, 22, 074209.
25. H. Ogasawara, N. Horimoto and M. Kawai, J. Chem. Phys., 2000, 112, 8229–8232.
26. V. Buch, H. Groenzin, I. Li, M. J. Shultz and E. Tosatti, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 5969–5974.
27. B. Medronho, A. Romano, M. G. Miguel, L. Stigsson and B. Lindman, Cellulose, 2012, 19, 581–587.
28. R. Alqus, S. J. Eichhorn and R. A. Bryce, Biomacromolecules, 2015, 16, 1771–1783.
29. C. Bruel, J. R. Tavares, P. J. Carreau and M.-C. Heuzey, Carbohydr. Polym., 2019, 205, 184–191.
30. M. Li, C. R. Luckett and T. Wu, Biomacromolecules, 2021, 23, 497–504.
31. A. K. Tamo, J. Mater. Chem. B, 2024, 7692–7759.
32. C. Budke and T. Koop, ChemPhysChem, 2006, 7, 2601–2606.
33. P. M. Naullage, L. Lupi and V. Molinero, J. Phys. Chem. C, 2017, 121, 26949–26957.
34. Y. Nishiyama, P. Langan and H. Chanzy, J. Am. Chem. Soc., 2002, 124, 9074–9082.
35. T. Bučko, D. Tunega, J. G. Ángyán and J. Hafner, J. Phys. Chem. A, 2011, 115, 10097–10105.
36. J. D. Kubicki, M. N.-A. Mohamed and H. D. Watts, Cellulose, 2013, 20, 9–23.
37. I. Morrison, J.-C. Li, S. Jenkins, S. S. Xantheas and M. C. Payne, J. Phys. Chem. B, 1997, 101, 6146–6150.
38. J. D. Bernal and R. H. Fowler, J. Chem. Phys., 1933, 1, 515–548.
39. N. H. Fletcher, Philos. Mag., 1962, 7, 255–269.
40. Z. Yan, S. V. Buldyrev, P. Kumar, N. Giovambattista, P. G. Debenedetti and H. E. Stanley, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2007, 76, 051201.
41. J. Russo and H. Tanaka, Nat. Commun., 2014, 5, 3556.
42. M. A. Klatt, J. Kim, T. E. Gartner, III and S. Torquato, J. Chem. Phys., 2024, 160, 204502.
43. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
44. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
45. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
46. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
47. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
48. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
49. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
50. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
51. P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. D. Jr, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. dela Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu and S. Baroni, J. Phys.: Condens. Matter, 2017, 29, 465901.
52. D. Hamann, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 085117.
53. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Solid State, 1976, 13, 5188.
54. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
55. Y. Nojima, Y. Shioya, H. Torii and S. Yamaguchi, Chem. Commun., 2020, 56, 4563–4566.
56. J. Hong, Y. Tian, T. Liang, X. Liu, Y. Song, D. Guan, Z. Yan, J. Guo, B. Tang and D. Cao, et al. , Nature, 2024, 1–6.
57. E. A. Engel, B. Monserrat and R. J. Needs, J. Chem. Phys., 2016, 145, 044703.
58. X. Zhang, W. Yang and C. Dong, J. Anal. Appl. Pyrolysis, 2013, 104, 19–27.
59. X. Yang, Z. Fu, D. Han, Y. Zhao, R. Li and Y. Wu, Renewable Energy, 2020, 147, 1120–1130.
60. V. Maliekkal, S. Maduskar, D. J. Saxon, M. Nasiri, T. M. Reineke, M. Neurock and P. Dauenhauer, ACS Catal., 2018, 9, 1943–1955.
61. M. C. Jarvis, Cellulose, 2023, 30, 667–687.
62. D. H. Lee and H. Kang, J. Phys. Chem. C, 2023, 127, 2885–2893.
63. G. McLaurin, K. Shin, S. Alavi and J. A. Ripmeester, Angew. Chem., Int. Ed., 2014, 53, 10429–10433.
64. S. G. Lee, J. I. Choi, W. Koh and S. S. Jang, Appl. Clay Sci., 2013, 71, 73–81.
65. A. Hudait, Y. Qiu, N. Odendahl and V. Molinero, J. Am. Chem. Soc., 2019, 141, 7887–7898.
66. P. L. Davies, Trends Biochem. Sci., 2014, 39, 548–555.
67. M. Bar Dolev, R. Bernheim, S. Guo, P. L. Davies and I. Braslavsky, J. R. Soc. Interface, 2016, 13, 20160210.
68. M. Bayer-Giraldi, G. Sazaki, K. Nagashima, S. Kipfstuhl, D. A. Vorontsov and Y. Furukawa, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 7479–7484.
69. J. R. Errington and P. G. Debenedetti, Nature, 2001, 409, 318–321.
70. P.-L. Chau and A. Hardwick, Mol. Phys., 1998, 93, 511–518.
71. D. R. Nutt and J. C. Smith, J. Am. Chem. Soc., 2008, 130, 13066–13073.
72. A. Hudait, N. Odendahl, Y. Qiu, F. Paesani and V. Molinero, J. Am. Chem. Soc., 2018, 140, 4905–4912.
73. A. Hudait, D. R. Moberg, Y. Qiu, N. Odendahl, F. Paesani and V. Molinero, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 8266–8271.
74. J. G. Davis, K. P. Gierszal, P. Wang and D. Ben-Amotz, Nature, 2012, 491, 582–585.
75. K. Meister, S. Strazdaite, A. L. DeVries, S. Lotze, L. L. Olijve, I. K. Voets and H. J. Bakker, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 17732–17736.
76. E. I. Howard, M. P. Blakeley, M. Haertlein, I. P. Haertlein, A. Mitschler, S. J. Fisher, A. C. Siah, A. G. Salvay, A. Popov and C. M. Dieckmann, et al. , J. Mol. Recognit., 2011, 24, 724–732.

---