2. Results and Discussion
DFT Results. The properties of THAL + water mixtures stand on the hydrogen bonding abilities of both types of molecules, which should lead to an extensive network of interactions along the liquid structure. Therefore, the characteristics of these mixed fluids are produced by the nature of the possible THAL-water interactions, which can be accurately analysed by using DFT methods. For this purpose, clusters with 1 THAL + i water (i = 1 to 3) were optimized with DFT approach, thus providing information of water clustering around THAL molecules. Although interactions in liquid state will extend beyond these simple cluster models, the characterization provided by DFT calculations leads to information which can be extended through bulk liquid phased by MD simulations.
In a first stage, the structure of THAL monomers was optimized at the considered theoretical level (B3LYP-D3/6-311++g(d,p)),
Figure 2a. The α,α isomer reported in
Figure 2a is the naturally occurring form of THAL [
12], as well as the lowest energy conformer [
13], and thus, it was the only isomer considered for DFT calculations. Regarding the THAL structure, the presence of four hydroxyl groups in each hexose ring, which are fully exposed to external solvent molecules being not sterically hindered in contrast with the glycoside oxygen, which would be hindered for accepting hydrogen bonds, should lead to extensive interactions with surrounding water molecules, i.e., water–THAL hydrogen bonding can be developed through any of the eight hydroxyl groups with similar topological characteristics and strengths.
For the case of 1 THAL:1 water clusters (
Figure 2b,d), the most relevant interaction sites were considered: (i) water acting as a bridge between the two hexoses involving the hexose -CH2-OH group (
Figure 2b) and (ii) water bridging two hydroxyl groups in one of the hexoses at two different sites (
Figure 2c,d). The case of water bridging the two hexoses leads to larger binding energies than in the case of water molecules hydrogen bonding a single hexose ring (roughly 42% stronger interactions), with the water molecule developing hydrogen bonds with hydroxyl groups at roughly the same distances, thus pointing to strong hydrogen bonds with both rings;
Figure 2b. Likewise, the connection of the two hexose rings in
Figure 2b leads to large changes in THAL molecule, with changes in the relative disposition of the two hexose rings around the glycoside oxygen. The structure of THAL bonded to water in
Figure 2b is more compact than the one for isolated THAL in
Figure 2a, the calculated Connolly volumes being 295.24 and 301.04 Å
3, respectively. The development of water-THAL hydrogen bonding in a single hexose ring (
Figure 2c,d) leads to large binding energies but lower than those in
Figure 2b. In any case, water molecules develop two simultaneous hydrogen bonds with hydroxyl groups in hexose rings with similar H to O distances pointing in all the cases to very effective hydrogen bonding.
Results in
Figure 3 show AIM and RDG (Reduced Density Gradient) analysis for the 1:1 interacting pair in
Figure 2b. This interaction is characterized by the development of two binary critical points (BCPs, type (3, −1)) along the donor–acceptor line, closer to the corresponding H atoms, for which the electron density,
, and the corresponding Laplacian,
, are in the upper limit of the corresponding ranges considered for defining hydrogen bonding in the AIM framework (0.002 to 0.035 a.u., 0.024 to 0.139 a.u, for
and
, respectively) [
14]. Therefore, although some doubts have been raised in the literature about the quantification of hydrogen bonding by AIM parameters [
15], the results for BCPs point to very strong hydrogen bonds according to interaction energies in
Figure 2. Likewise, the interacting water-THAL region is characterized by the development of two ring critical critical points (RCPs, type (3, +1)), which confirm the reinforcement (cooperativity) of the binding by the presence of the two simultaneous interactions connecting both hexose rings [
16]. The analysis of the hydrogen bonding in terms of the Noncovalent interaction (NCI) analysis [
17] showing the Reduced Density Gradient (RDG) surfaces is also reported in
Figure 3a. The two large blue spots in the water -THAL hydrogen bonding bond path show strong attraction corresponding to the large
and negative λ
2 values (being λ
2 the largest eigenvalue of Hessian matrix of electron density). Additionally, the green spot around the RCPs confirms the reinforcement of the interaction by additional van der Waals interaction. The contour map of electron density in one of the water-THAL hydrogen bonds in
Figure 3b shows the densification along the hydrogen bonding bond path. For the case of 1 THAL:2 water and 1 THAL:3 water clusters (
Figure 2e–h), the reported results show that water molecules can be hydrogen bonded to neighbour sites in THAL molecule without weakening the THAL-water strength of interaction, the trend of energetics being maintained, thus leading to a cooperativity effect, which would lead to efficient THAL solvation in liquid water solutions.
MD simulations. The DFT results concluded strong hydrogen bonding between water and THAL molecules, which would compete in the liquid state with the strong term of THAL molecules to self-aggregate to form large clusters [
10]. These effects were analysed in this work by MD simulations as a function of mixture composition and temperature considering large ranges to study different dry preservation scenarios.
Water-THAL mixtures are characterized by the formation of glassy states characterized by the corresponding glass transition temperatures,
Tg [
18]. The formation of glassy structures was inferred from MD simulations by the mixture composition evolution of predicted density (
Figure 4a), which evolve through maxima with temperature for each mixture, thus allowing the prediction of
Tg (
Figure 4c). Predicted
Tg were compared with those from Couchman–Karasz (CK) model, which showed excellent agreement with experimental results [
19]. The MD predicted
Tg are in excellent agreement with those from CK model. The maxima of density vs.
T curves shifts towards lower temperatures as the water content increases. A non-linear evolution of
Tg vs. THAL mass fraction (ω) results, which is in agreement with the literature data [
10] showing the trend of THAL molecules to self-aggregate forming giant THAL clusters as the water content decreases (increasing ω in
Figure 4c), whereas the large
Tg point to small THAL clusters for water rich mixtures. The evolution of predicted MD density with mixture composition for isothermal conditions as reported in
Figure 4b also shows non-linear evolution, especially for low temperatures, which agrees with the clustering and thus glassy states for THAL rich mixtures.
The strength of intermolecular forces from MD simulations are reported in
Figure 5. The total potential energy,
E(total), is reported in
Figure 5a as a function of temperature and composition. The
E(total) of the system provides an estimate about the stability of the systems under investigation. In general, we observe the stability of the molecular system decreases upon increasing the concentration of THAL and upon increasing the temperature. Mixtures with low THAL content show negative
E(total) in the full temperature range, whereas as the THAL content increases, positive
E(total) are inferred, which can be related with phase segregation through THAL large clustering. Additionally, the composition evolution of
E(total) vs. THAL content under isothermal conditions is largely non-ideal for all the considered temperatures, which points to large changes in the aggregation schemes, i.e., type and size of homo and heteroassociated clusters, with composition evolution. The interaction energy was split in the different contributions: water-water, water-THAL and THAL-THAL;
Figure 5b–d. The water-THAL interactions (heteroassociations) lead to larger interaction energies than the developed homoassociations (water-water and THAL-THAL). The homoassociations are stronger for the corresponding mixture rich regions, i.e., water–water for water rich mixtures and THAL-THAL for THAL rich mixtures. It is especially remarkable for THAL-THAL interactions, which evolve non-linearly with mixture composition reaching very large values, equivalent to water-THAL ones, for THAL rich mixtures. This behaviour agrees with the increasing trend of THAL molecules to self-associate reaching large clusters for THAL rich mixtures, i.e., leading to very large THAL-THAL interaction energy with the decrease of water-THAL interactions.
The development of strong hydrogen bonding in the mixtures as well as its non-linear evolution with mixture composition should lead to different molecular mobility in the mixtures considering THAL content and temperature. This effect was analysed through the MD predicted self-diffusion coefficients,
D, which were calculated from mean square displacements and Einstein’s equation for fully diffusive regimes. The temperature evolution of THAL
D values reported in
Figure 6a shows a non-linear evolution corresponding to the formation of a glassy state for the corresponding
Tg for each mixture composition. The large deviations of Arrhenius behaviour for THAL self-diffusion as reported in
Figure 6a manifest cooperativity of molecular diffusion. The cooperativity is maintained in the whole concentration range, thus suggesting THAL trend to clustering. Likewise, the decrease of THAL
D values with increasing THAL content confirms the THAL aggregation.
The structuring of water–THAL mixtures is firstly analysed by site-site radials distribution functions, RDFs. The hydrogen bonding between water and THAL molecules considering the different donor/acceptor sites in THAL is analysed in
Figure 7a. The reported results discard the formation of hydrogen bonds through the oxygen linking the two hexose rings, Ob, and through the two ether oxygen atoms in the hexose rings, Oe (atom labelling in
Figure 1);
Figure 7a. Regarding the interaction with the four possible hydroxyl groups in each ring, O1 to O4 (
Figure 1), the corresponding oxygen (water) to oxygen (OH in THAL) RDFs are characterized by a strong and narrow peaks at 2.8 Å, showing strong hydrogen bonds in agreement with DFT results in
Figure 2. The intensity of these RDFs first peaks is almost the same for all the available hydroxyl THAL groups, thus confirming that water molecules are hydrogen bonded to all the available sites although with slightly different extension around each site. The differences arise in the water structuring beyond the first solvation shell (hydrogen bonds); the second water solvation shell is inferred from the second RDF peak, whose intensity increases from O1 to 4 OH groups, i.e., larger RDFs as the OH groups are further from the bridge oxygen.
The spatial distribution function of water molecules around a central THAL molecule is reported in
Figure 7b, which shows water molecules concentrating around the two hexose rings, although the distribution is not homogeneous with larger water concentration around one of the hexose ring in comparison with the other one, which can be justified considering that most of the water molecules will be placed in the outer cores of THAL-THAL clusters, thus with a part of the THAL molecule involved in THAL-THAL interactions and the other part (second hexose ring) exposed to water and thus leading to water-THAL hydrogen bonding in larger extension. The composition effect on RDFs is analysed (at 310 K) in
Figure 8a, where Ow (water)–O1 (THAL) as a function of THAL content; the reported results show that RDF first peaks is maintained in the whole composition range with very minor changes in RDFs beyond the peak corresponding to hydrogen bonding. The integration of RDFs leads to the corresponding running integrals,
N, i.e., the number of water molecules around each THAL hydroxyl site;
Figure 8b. The
N values in
Figure 8b follows the ordering O1 > O4 > O2 > O4 > O3, i.e., larger concentration of water molecules directly hydrogen bonded to THAL for OH groups closer to the O bridge site. Nevertheless, the differences are minor, and for low THAL content,
N values are close to two for all the considered sites, thus confirming that water molecules are hydrogen bonded to O1–O4 sites. The increasing THAL concentration leads to a decrease of
N following a non-linear evolution with THAL mass fraction. This behaviour can be justified with the increase of THAL-THAL clusters, which decrease the free available OH sites to be hydrogen bonded to water molecules. The non-linear evolution of
N shows that THAL clusters do not increase their size just because the increase of THAL content, with the size increasing abruptly from a certain THAL content [
10]. The temperature effect on water–THAL aggregation is reported in
Figure 9 in the 100 to 400 K range (for ω = 0.112). The reported results show that although the position of the first RDF peak is maintained in the whole temperature range, i.e., water and THAL are hydrogen bonded in the whole range; the intensity of this first peak decreases upon heating with a large change through the
Tg and large differences in the ordination beyond the first solvation shell, which correspond to different structural arrangements in the glassy and liquid states. Regarding the THAL-THAL interactions, RDFs for Ob-Ob (as centre-of-mass to centre-of-mass interactions) are reported in
Figure 10a for isothermal conditions (
T = 310 K) as a function of THAL content. The reported RDFs show peaks in the 5 to 10 Å range, with narrower peaks for water rich mixtures evolving to wider bands as THAL concentration increases. Therefore, it shows THAL clustering, which is confirmed by the corresponding
Ns in
Figure 10b, the non-linear increase of THAL molecules around a central THAL shows increase of THAL-THAL clusters sizes but with an abrupt increase beyond certain water THAL content as showed in
Figure 10c.
The THAL clustering is analysed as visualized in
Figure 11. RDFs for Ob-Ob sites (THAL-THAL self-association) are plotted for two THAL contents as a function of temperature as well as the corresponding snapshots showing THAL clusters. For low THAL content (
Figure 11a) RDFs show poor ordering for low temperatures and thus only once the
Tg is surpassed, the clustering is inferred. Analogous results are inferred as the THAL content is increased, although the RDFs show wider bands corresponding to larger THAL-THAL aggregates;
Figure 11b. The snapshots reported show poorly defined THAL clusters for low temperatures below the
Tg whereas clusters are properly defined once
Tg is surpassed. The increase in cluster size with increasing THAL content is clearly inferred from results in
Figure 11e,f; the massive clusters (even propagating through neighbour simulated cells) at ω = 0.240 contrast with the smaller clusters for ω = 0.112, which confirm that THAL-THAL homoassociations are very favoured, and massive [
10] clusters are formed thus leading to water molecules solvating, by hydrogen bonding, the external shells of these large clusters.
The orientation of THAL molecules in THAL-THAL large clusters was quantified by combined distribution functions, CDFs (plotted using TRAVIS software [
20]), for the vectors joining two opposite carbon atoms in the hexose rings and the Ob-Ob intermolecular distance. It should be remarked that, as results in
Figure 7b show, the two hexose rings are not coplanar, which agrees with DFT results in
Figure 2. Results in
Figure 12 show a strong composition dependence on CDFs, as it may be expected considering the different geometries of THAL clusters reported in
Figure 11e,f. As a rule, for most of the mixture compositions no large spots for the studied angle are inferred, thus discarding parallel orientation of neighbour THAL molecules. Therefore, stacking is not inferred and THAL molecules seem to be mostly randomly oriented for leading to the clusters reported in
Figure 11e. The temperature effect on CDFs is reported in
Figure 13 showing the poor THAL clustering below the
Tg, in agreement with the results in
Figure 11 for cluster properties. Additionally, large temperatures (e.g., 400 K in
Figure 13d) lead to less ordered THAL clusters.
The kinetic properties of the different hydrogen bonds were analysed by using the reactive flux method by Gehrke et al. [
21] using the forward (hydrogen bond dissociation) and backward (hydrogen bond reformation) processes characterized by the corresponding lifetimes,
and
. Therefore, the
can be considered to quantify the duration of the corresponding hydrogen bonds as well
to measure the difficulties to form a new hydrogen bond after breaking. Hydrogen bonds were defined with 3.5 Å and 60° as criteria for donor-acceptor separation and angle. Results in
Figure 14a for
show that the hydrogen bonds developed in water–THAL mixtures have very different lifespans. THAL-THAL hydrogen bonds, independently of the considered hydroxyl group involved in the interaction, have
an order of magnitude larger than any other hydrogen bond in the full composition range. The lifespan for all the considered hydrogen bonds increases with increasing THAL content in a non-linear way, which can be related with the development of larger THAL clusters, which decreases molecular mobility;
Figure 6. It should be remarked that water-water hydrogen bonds have the lower lifespans, but they are also quickly reformed,
Figure 14b. THAL-THAL hydrogen bonds have large lifespans, but the reformation time is even larger than the
; thus, once THAL-THAL hydrogen bonds are broken, it is difficult to find a new suitable position to develop a new hydrogen bond, which is again related with the giant clusters formed for THAL rich mixtures, which decrease molecular mobility and also molecular reorientation forming hydrogen bonding reformation. THAL-water hydrogen bonds show similar properties considering THAL as donor or as acceptor, with lifespans in between the THAL-THAL and water-water interactions. The non-linear increase of kinetic properties for THAL-water interactions is again related with the THAL clustering considering that for THAL rich mixtures, once the large THAL clusters are formed, most of the water–THAL interactions will be developed with THAL molecules in the external shells of the THAL clusters. Regarding the temperature effect (
Figure 15), the effect of
Tg is inferred with extremely large lifespans for all the considered interactions below
Tg, which increase upon heating above
Tg. Likewise, reformation times are very large for low temperatures (it was not possible to measure it for 100 K).
The trend of molecules to self-aggregate forming large clusters should lead to changes in the properties of THAL molecule structure. Results in
Figure 16a compare the average structure of THAL in water-THAL mixtures (ω = 0.112,
T = 110 K) with the structure obtained in gas phase from DFT results. This comparison shows that THAL molecule in water-THAL mixtures is very different to that in gas phase, whereas in gas phase, both hexose rings are close to coplanar (151° for the dihedral angle around the glycoside oxygen,
Figure 16a,b); in the case of water-THAL, one of the hexose rings suffers a rotation remaining almost perpendicular to the other ring;
Figure 16a. Additionally, MD simulations of isolated THAL molecules (i.e., gas phase MD simulations) were carried out leading to 140 ± 9.5° for the dihedral angle around the glycoside oxygen, which are similar to those from DFT results and points to rotation of the THAL rings as produced by water presence. This rotation upon mixing with water is quantified in
Figure 16b; the distribution of the dihedral angle is almost independent of the THAL content; thus, hexose rotation is inferred even for low THAL content. The reason of this THAL configuration in water mixtures could be the formation of a conformer, which allows more efficient THAL packing for the formation of THAL large clusters and allowing THAL-THAL hydrogen bonding in a more efficient way, thus fulfilling the THAL trend to self-associate.
The extension of hydrogen bonding in water-THAL mixtures is quantified in
Figure 17, using the same hydrogen bonding criteria as for the kinetic properties in
Figure 14 and
Figure 15. Results in
Figure 17a show the evolution of total number of hydrogen bonds in the simulation cell for the different interacting pairs. The non-linear evolution of hydrogen bonds with mixture composition stands on the large non-ideal behaviour of the water-THAL mixtures. Likewise, the strong affinity of water molecules for THAL molecules is also quantified as well as the extension of THAL-THAL hydrogen bonding. The average number of hydrogen bonds per molecule are reported in
Figure 17b,c. Fort THAL-THAL interactions, a low number of hydrogen bonds per THAL molecule is inferred, starting from one hydrogen bond per molecule at low THAL content and evolving to two with increasing THAL content. Therefore, the clustering of THAL molecules is not characterized by a large extension of THAL-THAL hydrogen bonds; each THAL molecule is connected to another THAL molecule by a small number of hydrogen bonds, which are forming a network of hydrogen bonds leading to the large THAL clusters. The THAL molecules are also largely hydrogen bonded to surrounding water molecules both acting as donor and acceptor;
Figure 17b. The increasing THAL concentration decreases water-THAL hydrogen molecules; as the THAL cluster size increases, the THAL hydroxyl sites are partially hidden into the cluster structure, thus hindering hydrogen bonding with water molecules. It is confirmed that for rich THAL mixtures, THAL-THAL interactions are just two per molecule whereas the number of water-THAL interactions has decreased by half when compared with water rich mixtures. In the case of water-water hydrogen bonding, it is decreased by half with increasing THAL content in comparison with water rich mixtures; water molecules tend to be hydrogen bonded with THAL molecules, which hinders water-water interactions with increasing THAL content. Therefore, the considered mixtures are characterized by hydrogen bonding networks involving THAL molecules, leading to THAL cluster, with water molecules hydrogen bonded to the remaining free THAL hydroxyl sites, developing water-THAL hydrogen bonds mainly with THAL molecules in the external regions of cluster, thus hindering water-water interactions.
3. Methods
Density functional theory (DFT) calculations were carried out to analyse THAL-water hydrogen bonding using around central THAL molecules considering small 1 THAL +
i water molecular clusters (
i = 1 to 3). These simulations allow accurate characterization of hydrogen bonding between the considered molecules in terms of interaction energy and the topology of the interacting sites as a function of the possible donor and acceptor sites. DFT modelling was carried out using ORCA program [
22] at B3LYP [
23,
24,
25] functional and 6-311++G(d, p) basis set, with dispersion correction from semiempirical Grimme’s method (DFT-D3) [
26]. Interaction energies, Δ
E, for the studied 1 THAL +
i water clusters were calculated as the difference of the energy for the cluster and the sum of the energies of the corresponding monomers, at the same theoretical level, with Basis Set Superposition Error (BSSE) corrected through counterpoise method [
27].
Classical MD simulations were carried out with MDynaMix v.5.2 [
28] program. The forcefield parameterization for THAL is reported in
Table S2 (Supplementary Material); charges and internal terms (bonds, angles and dihedrals) were obtained from Merck Molecular ForceField (MMFF) [
29] whereas Lennard–Jones terms were obtained from CHARMM22 [
30] forcefield as included in SwissParam database [
31]. Water was described according to SPW-flexible water model [
32]. The systems used for MD simulations, as described in
Table S1 (Supplementary Material) were prepared considering a fixed number of water molecules (1500) and variable number of THAL molecules covering the 0 to 0.760 THAL mass fraction (ω). Initial simulation boxes were built with Packmol [
33] program. Periodic boundary conditions were used for all MD simulations. Simulations were done in the NPT ensemble at 1 bar and temperatures in the 100 to 400 K range (100, 150, 180, 200, 230, 270, 310, 340, 370 and 400 K). Equilibration steps (2 ns) were followed by 20 ns long production runs. The equations of motion were treated using the Tuckerman–Berne double time step algorithm [
34] (1 and 0.1 fs for long- and short-time steps). Coulombic terms were handled with the Ewald method [
35], the number of terms in the reciprocal part was selected considering that the reciprocal series were cut when expression in the
exp of the reciprocal part exceeds 9. For Lennard–Jones interactions, 15 Å cut-off was considered, and cross-terms were calculated with Lorentz–Berthelot mixing rules for all interacting pairs. No scaling factor was used for 1-4 Lennard–Jones interactions. This approach allowed covering both different Trehalose weight percentage (0% to 100%) and the temperature range. The objective of this approach was to mimic the drying process and the structural changes upon water content reduction. Although Lammert et al. [
36] reported THAL solubility in water being 46.6 g THAL per 100 g solution, studies in the full composition range are required to study dehydration process. In the literature, previous MD studies have carried out MD simulations for THAL concentration well above 46.6 wt % [
19,
37]. Olgenblum et al. [
29] used an analogous procedure to the one reported in this work to study Trehalose–water mixtures in the 0 to 100 wt%.
A relevant property for the study of THAL-water mixtures is the glass transition temperature,
Tg, as a function of solution concentration. The MD simulation of
Tg, is a quite challenging task [
38], and it may be largely dependent on the considered cooling rate. Therefore, to avoid this artefact, the considered computational procedure, analogous to Olgenblum et al. [
29], does not consider a particular cooling rate, generally employed in previous MD studies, but simulations for each of ten different temperatures (ranging from 100 K–470 K) and for each of the trehalose-water mixture composition. This procedure allowed enough equilibration for each temperature then followed by simulations at lower temperatures. Subsequently, the evolution of the predicted density for the different trehalose-water mixture at temperature range between 100–400 K, allows to infer
Tg values for each solution concentration. The
Tg values were determined from the temperature evolution of density for each mixture composition, although the most common approaches for the determination of
Tg stand on the evolution of properties such as heat capacity, thermal expansion coefficient, diffusion coefficient or specific volume [
29]; the use of density has also been considered in the literature [
39].