Hydrogen bonding in glassy trehalose–water system: Insights from density functional theory and molecular dynamics simulations

Disaccharides such as trehalose are widely used as excipients in solid-state formulations of proteins and other biologics.^1–3 The protective mechanisms of disaccharides in a solid state include the replacement of water (the substance that can accelerate chemical reactions and physical degradation processes in biomolecules) and the formation of a glassy state, where all processes are expected to slow down dramatically. In practice, complete removal of water molecules from a carbohydrate matrix in processes like freeze-drying is impossible; there is always a certain amount of residual water.^4,5 The residual water can substantially influence the properties of the formulations, affecting their physical properties and stability.

Theoretical description of interactions between disaccharides and water at low hydration levels is a complicated task since, under these conditions, the matter becomes glassy with no rigid and ordered crystal structure. The enthalpy relaxation and entropic effects also become important. Standard theoretical approaches, which describe either molecules (in vacuum or in an environment) or solids with a regular structure, have to be adjusted for organic materials with properties intermediate between solid and liquid.

Experimental studies of the properties of residual water in amorphous formulations (and water in solid state carbohydrates generally) are challenging. First, the amounts of water can be very low, which limits the applications of certain techniques, especially keeping in mind the molecular similarities between the carbohydrate matrix and water. Second, the properties of the glassy state are dependent not only on the thermodynamic parameters during the experiment but also on the sample history. Whereas in cooling/heating experiments, the dependence on sample history can be described by existing models such as the Tool–Narayanaswamy–Moynihan (TNM) model,⁶ the dependence of glassy systems on sorption/desorption rates is considerably less understood. Therefore, experimental data on water properties in the glassy matrix can be strongly dependent on the experimental method design and previous sample history. Third, when experiments include the addition or removal of water during sorption or desorption, slow diffusion kinetics can influence experimental results.

An alternative way to address the problems described earlier is to use molecular simulations to model the hydrogen bonds of water with carbohydrates. Quantum calculations of trehalose conformations were described in the literature before,^14–16 as were trehalose–water interactions in the gas phase^17,18 and in solution,¹⁷ but interactions of trehalose with explicit water molecules in the solid state have not been reported yet. In this work, we perform density functional theory (DFT) quantum calculations of trehalose–water interactions in vacuum and in solid-state with periodic boundary conditions. Furthermore, the temperature dependence of these interactions and the glass transition in the system are modeled using molecular dynamics simulations. By gradually removing water molecules from the simulated system, we also consider water desorption from the glassy state. While the glass transition in the trehalose–water system was studied using molecular dynamics (MD) simulations before,^18,19 the issue of the enthalpy of water sorption by solid state trehalose has not, however, been explicitly addressed.

Quantum chemical calculations give a complete picture of strong covalent bonding and weak hydrogen bonding, with very high accuracy and with a minimal number of parameters. Such calculations can be performed for isolated molecules, molecules in liquid or solid-state environments, or with periodic boundary conditions. Since the target of the current study is the formation of weak hydrogen bonds, the accuracy of the method is essential. The results might be sensitive to the details of the calculations, such as the basis sets and the functionals. To justify the computational procedure, we studied the interaction between trehalose and water and tested it on a benchmark case—a water dimer.

This section briefly describes three different approaches that can be applied to the theoretical study of the interaction between trehalose and water: electronic structure calculation in vacuum, electronic structure calculation in solid state, and force field modeling.

In this study, we used two computational codes to describe the electronic structure of molecules in vacuum: OpenMX [employing basis sets in the form of a linear combination of localized pseudoatomic orbitals (LCPAOs)] and Molcas²⁰ [which utilizes basis sets obtained as a linear combination of atomic natural orbitals (ANOs)]. The former basis sets are constructed in a typical “solid state fashion” with a cutoff by energy, while the latter basis set is atomic-like and includes the primitives with very high angular momentum. In “standard” DFT calculations of molecules, there is no point in using more expensive basis sets; however, in the current study, we will use geometries far away from equilibrium conditions and with weak interactions between molecules. Therefore, the quality of the basis sets might be important.

OpenMX²¹ was used to explore atomic interactions between water and trehalose both in a vacuum (without periodic boundary conditions) and in an amorphous solid (with periodic boundary conditions). It is a DFT²² code with the basis functions in the form of LCPAO.²³ The generalized-gradient approximation was employed as parametrized by Perdew, Burke, and Ernzerhof (PBE).²⁴ In particular, the following basis functions within a confinement scheme^23,25 were employed: H7.0-s2p1, O7.0-s2p2d1, and C7.0-s2p2d1, where the principal symbol designates the chemical name, followed by the cutoff radius in the Bohr atomic units and the primitive orbitals. Henceforth, the employed basis set is referred to as LCPAO with the above-stated basis functions. When there is a difference in the basis functions, this is indicated. For instance, LCPAO (H5.0-s3p2, O5.0-s2p2d1) or LCPAO (H7.0-s3p2) that means that the basis functions for H and/or O are not identical to those specified earlier. An energy cutoff of 150 Ry was chosen to obtain the total energy precision of 10⁻⁶ Ry per cell within a typical real-space grid division of 72 × 105 × 96 for the solid state (orthorhombic symmetry). All configurations were relaxed at 0 K by minimizing interatomic forces. The binding energy was calculated by subtracting the total energy of water + trehalose (or water + water) from the individual counterparts. In vacuum, interactions within a water dimer as well as between water and a trehalose molecule were considered. Six configurations were explored for the latter, taking into account both hydrogen acceptor and hydrogen donor scenarios. All configurations were analyzed using the VESTA software.²⁶ 

ANO basis sets used in Molcas code can be classified based on the number of basis functions, from double-zeta (DZP) to triple (TZP), quadruple (QZP), and pentupal (5ZP). Molcas has implementations for several DFT functionals, but for consistency with OpenMX results, we used the same functional, such as PBE. The benchmark calculations with other functionals (36 combinations of basis sets and functionals) are presented in the supplementary material. Although D3 dispersion correction²⁷ is important for the computation of weak hydrogen bonds, it is less relevant for strong trehalose–water hydrogen bonds and hence was not used in this work.

To describe atomic interactions in a solid state, OpenMX²¹ was used with the same settings as in the case of molecular interactions in a vacuum (basis set LCPAO). An amorphous trehalose configuration was created using the rotation-annealing algorithm.²⁸ Four trehalose molecules (180 atoms) were initially stacked up in two layers and then randomly rotated with respect to each other. Four water molecules were then inserted at random positions. Such an initial configuration was held at 300 K for 2500 fs to mimic an amorphous configuration by allowing for a certain degree of molecular disorientation. DFT based molecular dynamics (canonical ensemble, velocity scaling, time step 1 fs) was run with the OpenMX code. At 0 K, this configuration was quenched and relaxed by minimizing the interatomic forces and adjusting the volume (constrained to the orthorhombic symmetry). Two more configurations were considered and derived from the original system, one with the thermal treatment at 200 K and another at 350 K for 2500 fs, followed by the same relaxation procedure at 0 K.

The molecular dynamics (MD) simulations were performed using the GROMACS²⁹ software package, version 2020. The CHARMM³⁰ force field with a modification proposed by Lay et al.³¹ and the TIP4P water model were used. The cutoffs for van der Waals and electrostatic interactions were 1.2 nm; the long-range electrostatic interactions were treated using the particle-mesh Ewald (PME) method. All simulations were performed using a Berendsen barostat with a reference pressure of 1 atm and a time constant of 5 ps. The temperature coupling for both isothermal and temperature scanning runs was performed using a Berendsen thermostat. For comparison, the Nose–Hoover thermostat in combination with the Parrinello–Rahman barostat was also tested. This combination worked well in isothermal experiments but caused strong oscillations in temperature scans, so it was not systematically used in this study. In temperature scanning experiments, five different scan rates from 3 to 10 K/ns were used. The absolute values of scan rates upon cooling and subsequent heating were equal. In temperature scanning runs, the simulation box consisted of 512 molecules of trehalose and 512 molecules of water.

In isothermal dehydration simulations, starting coordinates were taken from cooling scans at corresponding temperatures. In order to preserve the thermal history of the nonequilibrium glassy state, no further equilibration was performed. After a 1 ns run, dehydration was performed in 32 steps, where at each step, 16 molecules of water were removed, followed by a 1 ns isothermal MD run. Therefore, relatively gradual dehydration from 50 to 0 mol. % of water was simulated by 33 consecutive short runs.

The MD results were analyzed using GROMACS utilities and custom written MATLAB scripts. For visualization of MD data, the Savitzky–Golay filter was applied when appropriate.

For understanding trehalose–water interactions, hydrogen bonding in a water dimer is a natural reference. Therefore, the DFT method and the basis set used for calculations of trehalose–water interactions were first tested for water–water hydrogen bonding. Another reason to consider the water dimer is to test how variations in water molecules and hydrogen bonding geometry affect the binding energy in the simplest case.

Second, we considered the equilibrium geometry and interaction energy of the water dimer using various DFT methods and basis sets. For this, a single water molecule and water dimers were optimized with the ANO-L-VDZP basis set (double zeta + polarization for H-2s1p and O 3s2p1d) for each functional. Then, with fixed geometries, different basis sets were applied for each functional. Table S1 in the supplementary material shows calculated energies using six different methods and six different basis sets (36 energy values in total). It is not surprising that the largest and most expensive basis set, the ANO-L-5ZP, provides the lowest value of the total energy (data not shown), with the best agreement for interaction energy with the literature/experiment. For the calculation of relatively large molecular complexes such as trehalose–water, it is convenient to use a smaller basis set (e.g., ANO-L-VDZP or ANO-L-VTZP), which still provides a reasonable value of interaction energy (see Table S2).

Since in the amorphous glassy state, the hydrogen bonds are not in their optimal geometry, it is useful to quantify their energy dependence on the two main parameters: OO distance and HOO angle.

In order to characterize hydrogen bonding in multimolecular water complexes, we performed DFT calculations of the water tetramer using the PBE/ANO-L-VTZP combination (the use of the B3LYP functional gives numerically similar results). We should note that, in contrast to the water dimer, the considered tetramer is not stable. Therefore, the geometry optimization for the water tetramer, as shown in Fig. S3, was performed with certain constraints: a mirror plane has been used, so water molecules W2 and W3 are symmetry identical, the distance between molecule W1 and W4 has been fixed, as well as the angles between W2–W1–W3 and dihedral angle W2–W3–W1–W4.

The next step is to consider hydrogen bonds between trehalose and water molecules in a vacuum. We performed DFT calculations of 12 different complexes; two of them (with single hydrogen bonds) are presented in Fig. 1 (other complexes are illustrated in supplementary material, Fig. S5). When a complex has only one hydrogen bond, the binding energy is weaker when water donates the hydrogen bond, compared to the case when it accepts it (see Table I and Fig. 1, configurations 1 and 2, respectively). One should, however, keep in mind that different oxygen atoms of trehalose are involved in hydrogen bonding in these two cases.

When a complex had two trehalose–water hydrogen bonds, the strongest energies (around 45 kJ/mol) were observed when a water molecule donated one and accepted one hydrogen bond. We did not observe a significant difference in the binding energies when a water molecule bridges two glucose rings compared to the situation when both hydrogen bonds point to the same ring (even though there are data in the literature suggesting stronger binding in the case of bridging¹⁸). As one could expect, when a water molecule accepted one hydrogen bond and donated another one to a less polar ether oxygen, the binding energies were ∼10 kJ/mol lower (or about a half of the energy of a hydrogen bond between two hydroxyl groups, see Table II). When both hydrogen bonds were donated to trehalose OH groups by the same water molecule, the binding energies were also low, most probably due to the low values of the parameter n from Eq. (8). Overall, fitting of Eqs. (2)–(5) and (7) with the data from Table I resulted in reasonable values of fitting parameters (see Table II).

Moreover, a complex of two trehalose molecules with one water molecule forming three hydrogen bonds was also computed. Even though the geometry of the hydrogen bonds was apparently not optimal (see Table I), the binding energy was 75.5 kJ/mol, i.e., more than 25 kJ/mol per one hydrogen bond. This indicates that the local environment provided by trehalose molecules facilitates hydrogen bonding strength.

Three different structures that consist of four trehalose and four water molecules were generated as described in Sec. II B. Periodic boundary conditions were applied to all configurations. One of the structures is shown in Fig. 2.

To quantify the bonding energy of water molecules in the solid-state matrix, one water molecule at a time was removed from the structure; the remaining structure was not geometry-optimized in order to prevent the formation of new hydrogen bonds (which would strongly affect the system energy). The obtained energies shown in Table S4 range from about 52 to almost 100 kJ/mol, i.e., vary by a factor of two. The reason for a high variation can be found partly in the geometries of hydrogen bonds and partly in the properties of the local environment [quantified by the parameter n in Eq. (8)]. Still, these two factors do not fully explain the observed results. For example, the sum of n values is actually higher for the water molecule with the lowest binding energy than that for the water molecules with the two highest energies. The deviations from optimal geometries are not expected to provide the observed difference in energies (compare values in Table S4 and Figs. S2 and S3). This suggests that additional factors may affect the interaction energies in the solid state, which are not taken into account in the energy decomposition described by Eqs. (2)–(5) and (7)–(9). This can also be illustrated by the fact that configuration 12 from Table I described in vacuum is a part of solid-state structure N2, water molecule N6 (Table S4). The water molecule is bound by identical hydrogen bonds to trehalose molecules in those two structures; still, the energies are different (around 75 kJ/mol in vacuum vs around 100 kJ/mol in solid state). Although the basis sets used in these two calculations were different, this factor is unlikely to result in such a pronounced difference (the difference in energies of water dimers computed using these two basis sets was minor; see Table S2).

It is instructive to compare the obtained values of energies with the experimental data. The typically observed experimental hydration energy of −18 kJ/mol corresponds to a sorption energy of −62 kJ/mol. This is lower than the average water binding energy obtained in the solid-state computations (−79 kJ/mol) but falls within the range of energies presented in Table S4. The diversity of hydrogen bond energies reported here correlates with the results of infrared spectroscopy studies,³³ which suggest the existence of several populations of hydrogen bonds in the trehalose–water system. Additional factors, such as the temperature and annealing conditions, will affect the values of the energies of water–trehalose interactions. In the following sections, we analyze these issues using molecular dynamics simulations.

We performed MD simulations of the glass transition in a trehalose–water system with a composition of 1:1 on a mole basis, which corresponds to 5 wt. % of water. Moreover, simulations with a gradual dehydration were also performed; hence, a range of compositions from 0 to 1 mol/mol were studied.

Remarkably, the non-equilibrium heat capacities in the glassy state of the system are almost independent of the scan rate (and the direction of the temperature change). In contrast, the equilibrium heat capacities in the liquid state exhibit stronger deviations from the mean value, which further propagates into a variation of ΔC_p of the glass transition. The mean value of ΔC_p is 0.806 ± 0.029 J g⁻¹ K⁻¹, which is in excellent agreement with differential scanning calorimetry (DSC) data for pure trehalose of 0.81 J g⁻¹ K⁻¹ at 298 K.⁷ 

For a better comparison of the MD scans with different scan rates, the enthalpy as a function of temperature is presented in Fig. 3(b). In these coordinates, an enthalpy hysteresis is seen in the glass transition region. The origin of this hysteresis is an additional enthalpy relaxation during the heating scans, which further decreases the enthalpy values compared to the cooling curves. Moreover, in the low temperature plateau, a clear scan rate dependence of the enthalpy level is seen: the lower the scan rate, the lower the enthalpy. This implies that the glass transition calculated as a crossing point between the “glassy” and “liquid” baselines is lower for lower scan rates.

To understand the effect of temperature scans on intermolecular interactions, it is instructive to consider the change in hydrogen bonds between the components of the system. Table IV shows the number of hydrogen bonds between the components of the system at the two extreme temperatures of the simulation: 0 and 450 K. For comparison, the maximum possible number of hydrogen bonds is also shown. For trehalose–water hydrogen bonds, it was assumed that each water molecule could donate two and accept two hydrogen bonds from trehalose. In all other cases, it is assumed that the number of donors defines the maximum number of hydrogen bonds. The results show that at 0 K, hydrogen bonding is rather efficient: more than 90% of possible hydrogen bonds are formed in the system, while at 450 K, this is only 54%. At both temperatures, the highest percentage of the formed bonds is observed for trehalose–water interactions, while the lowest is observed for water–water interactions. We note, however, that this does not imply energetically stronger trehalose–water hydrogen bonds since trehalose molecules provide a much higher number of OH groups for potential hydrogen bonds.

The relative number of hydrogen bonds changes differently for different types of interactions (see Fig. 5). The strongest increase upon cooling is observed for trehalose–trehalose hydrogen bonds, while in the water–water case, the number of hydrogen bonds decreases. Interestingly, the main decrease of the water–water hydrogen bonds is observed in the glass transition region, while both in the liquid and glassy states, this number increases upon cooling. This, however, has a limited effect on the thermodynamic state of the system since the number of water–water hydrogen bonds is low (Table IV).

For comparison with sorption calorimetry data, we performed MD simulations of the gradual dehydration of the system. Starting coordinates for isothermal dehydration runs were taken from cooling scans at corresponding temperatures. The dehydration was performed in 32 steps, where at each step, 16 molecules of water were removed, followed by a short MD run. Therefore, the water content range from 1:1 to 0:1 was covered.

Immediately after each dehydration step, the energy of the system increased stepwise; see the example of data for 300 K in Fig. 6(a). After that, a quick energy relaxation on the time scale of 10 ps was observed. The new energy level was, however, higher than the previous one. This indicates that the removal of water molecules results in breaking a certain number of hydrogen bonds with a relatively fast formation of new hydrogen bonds. The overall change in hydrogen bonding is negative with every step.

For the gas phase, the energy of the water molecule $Ew0,gasNw0$ was taken to be 3RT, which is very close to the experimental values.³⁷ 

The MD results presented here reproduce the experimental value of the enthalpy of mixing of water $Hwm$ (−18 kJ/mol) observed in carbohydrates in the dry limit and 298 K. This parameter is dependent on temperature but was calculated only for 300 and 360 K since it requires data on pure liquid water and is not available for extreme temperatures. In contrast, the enthalpy and energy of sorption are not dependent on liquid water and can be calculated for the whole range of temperatures (Table V).

In order to interpret the observed exothermic effect of mixing from the point of view of hydrogen bonds, one needs to compare Eqs. (13) and (16). An exothermic effect of mixing arises when $ddrNhbtotNs$ is higher than $Nnb,w0Nw0$⁠, i.e., when the number of hydrogen bonds formed in the mixture as a result of the transfer of water is higher than the loss of hydrogen bonds in pure liquid water. Our MD results show that in pure liquid water at 300 K, a water molecule on average has around 1.8 hydrogen bonds (since each water molecule can both donate and accept hydrogen bonds, when considering every individual water molecule, this number should be multiplied by two, i.e., 3.6 hydrogen bonds per molecule). From Fig. 6(c), it is apparent that $ddrNhbtotNs$ is higher—around 2.5 hydrogen bonds. Hence, the transfer of one water molecule from a pure liquid to glassy trehalose adds 0.7 hydrogen bonds on average, producing the exothermic heat effect of −18 kJ/mol observed both in the experiments and simulations presented here. Interestingly, each water molecule adds three hydrogen bonds with trehalose and around 0.23 bonds with other water molecules (this is, however, concentration dependent) but removes 0.73 trehalose–trehalose hydrogen bonds, producing a total balance of 2.5 hydrogen bonds per transfer from the gas phase.

On the level of individual water molecules, an interesting and perhaps counterintuitive picture emerges. In pure water, each water molecule has 3.6 hydrogen bonds, while in glassy trehalose it has less—only 3.1. Upon the removal of a water molecule from a pure liquid, the remaining molecules readily form 1.8 new hydrogen bonds, producing a balance of −1.8 hydrogen bonds. In contrast, removing a water molecule with 3.1 hydrogen bonds results in the formation of only 0.73 new trehalose–trehalose hydrogen bonds (plus a small number of new trehalose–water bonds formed after breaking water–water bonds).

Hence, despite a lower number of hydrogen bonds for a water molecule in a trehalose–water system, the removal of such a molecule from the trehalose matrix is associated with higher energy compared to liquid water because the glassy matrix cannot effectively substitute for the missing hydrogen bonds. These findings can be relevant when considering water replacement and preferential hydration effects in protein–trehalose–water systems.

We performed an extensive study employing DFT and MD to investigate hydrogen bonding between trehalose and water, placing particular emphasis on interactions in the amorphous solid state. The observed hydrogen bonding energies, both in vacuum and in solid state, suggest that it is impossible to predict binding energy based solely on the geometrical parameters of hydrogen bonds. The local environment, such as the presence of other hydrogen bonds, plays an equally important role.

DFT quantum results reveal that the binding energy of an individual water molecule can be as high as 100 kJ/mol, but a broad distribution of energies reflecting different geometries and local environments exists. MD simulations show that the system energy, to a large extent arising from hydrogen bonding, is dependent on temperature and the presence of a glass transition. A strong increase in the total number of hydrogen bonds is observed upon cooling; however, this increase slows down with the glass transition. Simulation of the system dehydration reproduces the experimentally observed hydration enthalpy value of −18 kJ/mol in the glassy state. This effect arises from the limited ability of the glassy matrix to relax in response to the removal of water molecules.

See the supplementary material for additional plots, molecular structures, and tables illustrating the DFT calculations in vacuum and in solid state.

V.K. acknowledges support from the Swedish Governmental Agency for Innovation Systems (Vinnova), which was carried out within the competence center NextBioForm. Professor Daniel Harries and his group are acknowledged for their discussion of MD simulations. The computations using MOLCAS were enabled by resources provided by LUNARC. OpenMX modeling and MD simulations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) at the National Supercomputer Centre (NSC) in Linköping, Sweden, partially funded by the Swedish Research Council through Grant Agreement No. 2022-06725.

The authors have no conflicts to disclose.

Vitaly Kocherbitov: Conceptualization (lead); Data curation (equal); Investigation (equal); Project administration (lead); Visualization (equal); Writing – original draft (lead); Writing – review & editing (equal). Denis Music: Conceptualization (equal); Data curation (equal); Investigation (equal); Visualization (equal); Writing – review & editing (equal). Valera Veryazov: Conceptualization (equal); Data curation (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.