Introduction

It is commonly accepted to state that native cellulose microfibrils have a dominant hydrophilic character. At the ultrastructural level, this statement is supported by recording of high-resolution details of large highly crystalline microfibrils such as those of Valonia, tunicin, Micrasterias etc. Indeed, the analysis of their cross-sections by electron diffraction and diffraction contrast electron microscopy, including lattice imaging, clearly indicates that their surface corresponds to either the (1 1 0) and (1 −1 0) planes of the monoclinic cellulose Iβ lattice or the (1 0 0) and (0 1 0) of its triclinic Iα counterpart (Sugiyama et al. 1991), which are rich in hydrophilic OH moieties (Sugiyama et al. 1985; Revol 1982; Helbert et al. 1998; Kim et al. 1996). Observations of Valonia surface at near atomic resolution by AFM is in agreement with the TEM results as they definitely show organized arrays of surface hydroxymethyl groups within the Iα phase of Valonia cellulose (Baker et al. 1997, 2000).

Despite the established dominance of the hydrophilic surfaces of cellulose, the amphiphilic character of cellulose, which has been underestimated in the past, is attracting more and more attention (Medronho et al. 2012). As it is believed that in some cases, it is the occurrence of minor hydrophobic surfaces that play a determinant role in the reactivity of cellulose. This phenomenon is best documented in the case of cellulose digestion, where the action of efficient modular cellulases is strongly dependent on the adsorption capability of the cellulose-binding module (CBM) of the cellulases, prior to the action of the enzymatic attack. It is now established that the adsorption of the CBM on crystalline cellulose relies on the affinity of a planar strip of amino acids, phenylalanine, tyrosine or tryptophan, that is systematically found within the various CBM three-dimensional structures (Reinikainen et al. 1995; Linder et al. 1995; Tormo et al. 1996). In a seminal paper, it was even proven by electron crystallography analysis, that it was indeed at the hydrophobic (1 1 0) surface of the Iα phase of Valonia that isolated CBM containing the appropriate tryptophan set-up could be selectively anchored (Lehtiö et al. 2003). In theory, this (1 1 0) surface should be negligible in cellulose, being reduced to an acute corner consisting of only one chain if the Valonia microfibrils had a perfect squarish section delineated by the hydrophilic (1 0 0) and (0 1 0) triclinic surfaces. In fact the ultrastructural observations on cross-sections have revealed that the corners of the Valonia microfibrils were frequently blunt (Sugiyama et al. 1985; Lehtiö et al. 2003; Xu et al. 2009; Mazeau 2011) and thus the occurrence of a significant secondary surface such as that of the hydrophobic triclinic (1 1 0) face appears to be far from negligible. This occurrence explains why it is a hydrophobic type of an interaction that is likely to account for the affinity of cellulase for crystalline cellulose. It may be useful to recall that this interaction can be taken as an example of the well documented stacking of aromatic amino acid residues at the hydrophobic surface or apolar “patches” of pyranose rings. The phenomenon has been extensively studied as it is recognized as a key factor in many sugar protein recognition systems that are central to a wealth of biological processes (Quiocho 1989; Weis and Drickamer 1996; Screen et al. 2007; Diaz et al. 2008).

Without the clear-cut demonstration obtained in the case of cellulose/cellulase interaction, further evidences of the partial hydrophobic character of cellulose can also be inferred from a series of other experiments. In a recent report, the amphiphilic character of cellulose was evidenced by observing the stabilization of oil/water interfaces by the addition of aqueous suspensions of cellulose nanocrystals (Kalashnikova et al. 2012). The action of direct dyes on cellulose is another area where hydrophobic interactions of the van der Waals type have been proposed to account for the affinity of the dye molecules on cellulose (Pielesz 2007). These azo dyes have in common an elongated flat geometry, based on a succession of aromatic rings and the presence of sulfonate moieties (Inglesby and Zeronian 2002; Bird et al. 2006). It is commonly accepted that the direct dyes will adsorb with their long axis parallel to the cellulose microfibrils (Morton 1946), but it not clear whether their adsorption is essentially based on hydrophobicity or whether the electrostatic interactions or hydrogen bonding also play an important role. So far, it has not been demonstrated whether such dyes were adsorbing on the dominant hydrophilic cellulose surfaces or on the minor hydrophobic ones. In order to address such question, the present paper describes an attempt to model the adsorption of the direct Congo red (CR) dye on model cellulose crystals. Even if the staining of cellulose with CR is not used industrially for toxicological and colour stability reasons, its adsorption on cellulose and other polysaccharides, which is used for analytical purpose, is well documented (Wood 1980, 1982; Yamaki et al. 2005). In a molecular mechanics modelling study, Woodcock et al. (1995), have investigated the docking of CR molecules on the hydrophilic faces of crystalline cellulose. These authors, however, did not consider the possibility of adsorption on the hydrophobic face of cellulose. In our continuing interest in modelling the surface interactions of cellulose with hydrophobic and hydrophilic chemicals (Mazeau and Vergelati 2002; Mazeau and Rivet 2008; Da Silva Perez et al. 2004), this work was therefore devised to see whether stable interactions of CR on hydrophobic planes could also be found, using molecular dynamics simulations.

Methods

Computational details

All calculations were performed with the modelling package Materials Studio 5.5 (Accelrys. Software Inc. San Diego, USA) using the Dreiding force field parameters (Mayo et al. 1990). The geometry minimizations were performed using an algorithm containing a cascade of the steepest descent, and quasi-Newton methods. The convergence criterions used were 0.5 (kcal/mol)/Å for force and 0.001 kcal/mol for energy. Molecular dynamics were used in the NVT ensemble. Temperature was controlled using Nosé algorithm (1984). The standard Verlet algorithm (1967) was used to integrate Newton’s law of motion with the time step of 1 fs. The length of the MD simulations was 2 ns if not stated differently. The charges were calculated using quantum equilibration method (Rappé and Goddard 1991).

System preparation

In this work, we have considered only the cellulose Iβ allomorph (Nishiyama et al. 2002), but the protocol could be transposed to the cellulose Iα phase after adequate transformation of the crystalline lattice parameters. We defined a model of cellulose microfibril (see Fig. 1) made of 18 glucan chains, each consisting of 8 glucosyl residues and organized in 4 layers.

Fig. 1
figure 1

A schematic representation of a native cellulose microfibril. Labelling of the four different surfaces is indicated according to Mazeau 2011; the system modelled in this study, represented by the black cube, corresponds to the hydrophobic corner of a segment of a cellulosic macrocrystal

This model was inserted in a periodic box having dimensions a, b and c of 50.0 × 52.0 × 41.52 Å (along the cellulose fiber axis, four times 10.38 Å, the unit cell dimension). Covalent bonds were created across the periodic box in order to model infinitely long chains. Following the models of Xu et al. (2009) and Fernandes et al. (2011) who described the surface of cellulose microfibrils as consisting of a mixture of dominant hydrophilic and minor hydrophobic components, our cellulose model exposes 3 chains of the hydrophobic (1 0 0) surface. Two types of cellulose chains were considered: (1) those unconstrained in the top two layers, which may interact with the CR molecules; (2) those in the two underlying layers, which are composed of chains that will not be in contact with the adsorbate; they represent the inner bulk chains of cellulose and were constrained. The CR-cellulose system was subjected to periodic boundary conditions. Regarding the CR in solution, it is known that the dye molecules may aggregate in solution (Hillson and McKay 1965; Wälchi 1947; Burdett 1983) and there is an equilibrium between isolated molecules and CR aggregates. In this work, we have not considered the adsorption of molecular aggregates but concentrated on the adsorption of single CR molecules.

A CR molecule (Fig. 2) having the two sulfonate groups on the same side was built and was first subjected to geometry minimization.

Fig. 2
figure 2

Top A molecular structure of Congo red, the Φ torsion angle of interest is shown. Bottom numbering of the atoms in a glucose residue of cellulose

The resulting structure was then subjected to molecular dynamics for 500 ps and the molecule in the last frame was selected for the adsorption studies. The conformation of the dye molecule was similar to the twisted one found in the CR crystal (Ojala et al. 1995).

In a typical simulation, a CR molecule was positioned about 12 Å away from the cellulose surface (Fig. 3) and the whole system was subjected to simulations of adsorption process with molecular dynamics after a preliminary energy minimization.

Fig. 3
figure 3

An initial model showing : the periodic limits, the (1 0 0) hydrophobic corner the Iβ cellulose elementary fibril and CR in one of its initial position. Left projection perpendicular to the cellulose chain axis, Right projection parallel to the cellulose chain axis

Several dynamics experiments were performed with the CR molecules initially oriented parallel and perpendicular with respect to the cellulose fiber axis and at different temperatures, ranging from 300 to 650 K. In order to improve the conformational sampling, the last structure of the trajectory at 300 K resulting from an initial parallel orientation of the CR was extracted. The CR was then manually translated by ±4 Å along the surface, in directions parallel and perpendicular to the fiber axis. The CR was also rotated by 45° clockwise and anti-clockwise around an axis normal to the cellulose surface. Each resulting model was minimized and then subjected to a short 50 ps molecular dynamics at 300 K.

We then considered the possibility of accounting for explicit water environment. We encountered difficulties in performing the MD in a box full of water molecules due to the computational demand of MD calculations with a very large number of atoms. To overcome such difficulty, we devised two strategies. At first, MDs were performed with a CR molecule placed 12 Å above a cellulose surface, both sulfonate groups were solvated with 10 water molecules (SPC model) each. Such semi-explicit approach was similar to that published by Fennell et al. (2011).

Water molecules were equally distributed around the sulfonate groups and the geometry was optimized according to the standard procedure. MD experiments at different temperatures and with different initial position of CR were run and the last frames were minimized. In the second strategy a single frame was taken from the end of the MD trajectory in vacuum where CR is fully adsorbed and then the system was solvated by inserting 2,824 water molecules into the computational box with the help of Monte Carlo algorithm. The solvated complex was then minimized and a 2 ns molecular dynamics simulation at 300 K was performed.

To quantify the energetic consequences of the adsorption process in the gas phase, the energy of interaction between CR and cellulose surface (ΔEi) was calculated. It was obtained by subtracting the energies of the isolated dye molecule (Edye) and that of cellulose (Ecell) from the total energy of the complete system (Etot):

$$ \Updelta {\text{E}}_{\text{i}} = {\text{E}}_{\text{tot}} -{\text{E}}_{\text{dye}} - {\text{E}}_{\text{cell}} $$

In an explicit water environment, the total energy of the CR—cellulose complex (Etot) is the sum of contributions of: the energy of the water molecules (Ewat), the energy of CR (ECR), the energy of cellulose (Ecell), the energy of interaction between CR and water (\( {\text{E}}^{\text{i}}_{\text{CR - wat}} \)), the energy of interaction between cellulose and water (\( {\text{E}}^{\text{i}}_{\text{wat - cell}} \)) and the energy of interaction between CR and cellulose (\( {\text{E}}^{\text{i}}_{{{\text{CR}} - {\text{cell}}}} \)). Therefore, the \( {\text{E}}^{\text{i}}_{{{\text{CR}} - {\text{cell}}}} \) was estimated according to:

$$ \Updelta {\text{E}}^{\text{i}}_{\text{CR - cell}} = {\text{E}}_{\text{tot}} -{\text{E}}_{\text{CR}} -{\text{E}}_{\text{cell}} - {\text{E}}^{\text{i}}_{\text{wat - cell}} - {\text{E}}^{\text{i}}_{\text{CR - wat}} - {\text{E}}_{\text{wat}} $$

The global orientation of the adsorbed CR molecule with respect to the cellulose microfibril surface can be quantified by three relative rotation parameters: roll, tilt, and twist and by three relative translation parameters: rise, shift, and slide (Dickerson 1989). Rotation parameters are defined as specific angles between principal axes (long, short and perpendicular) of two planes (Fig. 4): the first plane defining the CR molecule and the second plane including the (1 0 0) surface of cellulose.

Fig. 4
figure 4

Description of the rotation and translation parameters between two planes (adapted from Dickerson 1989)

Results and discussion

Adsorption process in the gas phase

First, to reveal the low energy intermolecular arrangements between CR and cellulose, we performed simulations in the gas phase, thus avoiding the solvation effects. Trajectories were run for different initial orientations of CR with respect to the fiber axis and at different temperatures.

A visual inspection of the trajectories performed at room temperature shows that CR is quickly attracted by cellulose. It approaches the surface thanks to its polar sulfonate groups. The electrostatic interactions between negatively charged oxygen atoms of one of the sulfonate group and the positively charged hydrogen atoms of cellulose hydroxyl groups are obviously the driving force of the adsorption process. When the first sulfonate group grabs the cellulose surface, the process of adsorbing the second sulfonate group begins. Simultaneously, both naphthalene residues of CR come closer to the surface. Then, the entire molecule of the dye is in contact with the cellulose surface. Once CR is adsorbed on cellulose, a stable system is formed and no departure from the adsorbed conformation is observed.

A typical time evolution of the interaction energy, ΔEint, between CR and cellulose for the temperature 300 K is given in Fig. 5.

Fig. 5
figure 5

Evolution with time of two selected properties of the system (initial parallel CR position) during the 2 ns MD simulation at 300 K. Blue number of atomic close contacts (within a 5Å cutoff) between CR and the cellulose surface. Red interaction energy. (Color figure online)

As we expected, the decrease of the interaction energy at the beginning of the simulation revealed that the system reached a more stable state when CR left its free state to become adsorbed. The decrease occurred in steps: two of them are clearly visible in Fig. 5, where ΔEint takes the constant values of about −32 and −72 kcal/mol. Once CR is fully adsorbed, the interaction energy is stable at −105.4 ± 6.3 kcal/mol. For each trajectory, the average energies of interaction between CR and the cellulose microfibril were calculated and reported in Table 1.

Table 1 Energies of interaction (kcal/mol) at equilibrium between CR and cellulose calculated for different temperatures (K) and different initial positions in vacuum

The CR molecule was then translated along the y and z axes, rotated around the x axis, and the resulting complexes were simulated with molecular dynamics. Energies of interaction between CR and cellulose are also presented in Table 1. Changing the position and orientation of CR leads to the identification of different adsorption sites whose energies of interactions were of similar magnitude, the largest difference observed being of +12.0 kcal/mol. The minor difference between these energies suggests that there are no specific adsorption sites on the (1 0 0) cellulose surface. The selected examples of adsorption sites on the cellulose (1 0 0) surface are presented in Fig. 6.

Fig. 6
figure 6

Selected examples of adsorption sites of CR on cellulose. Adsorption on the (1 0 0) surface: Top left the last frame of the original trajectory, top right rotated around x axis (twist parameter = −45°), bottom left moved along z axis (shift parameter = + 4 Å), bottom right adsorption partly on the (1 0 0) and partly on the (1 1 0) surfaces

It should be noted that the absolute values of the interaction energies have no real meaning. The electrostatic interactions, estimated in vacuum, are about 100 times stronger than those in water solution. Consequently the energies reported in this section are strongly overestimated. However, the different adsorbed states of CR can be compared based on their energies, even though these are far too large. Such comparison is useful in the context of our study as it gives a rough estimate of the relative abundance of a given adsorption geometry in the equilibrium mixture.

Simulations were also performed at higher temperatures. Not surprisingly, when the temperature increased, the process of adsorption was more rapid. In addition, large modifications and disorder in the conformations of the glucan surface chains were observed, the largest changes being observed in the middle glucan chain, especially at its hydroxymethyl groups. It should be noticed that at a temperature as high as 600 K, CR did not adsorb as closely as in the lower temperatures: only the two sulfonate groups remained closely bound to the cellulose surface, whereas the rest of the CR molecule kept moving up and down.

In the majority of the recorded trajectories, CR adsorbed flat on the (1 0 0) surface of Iβ cellulose. However, in some trajectories, it laid partly on the expected (1 0 0) surface and partly on one of the two hydrophilic surfaces, (1 1 0) or (1 −1 0), with comparable energy of interaction of −105.1 (±10.4) kcal/mol.

To decipher the chemical groups involved in the interaction and thus the forces stabilizing the complex, we counted the number of close contacts between CR and cellulose occurring below 5Å. A time evolution of the close contacts during simulation of the adsorption process is presented at Fig. 5. As expected, the number of close contacts increased at the beginning of the simulation and remained steady when CR was adsorbed. The average number of close contacts at equilibrium was remarkably constant; they were between 400 and 500. The description of the groups involved in close contacts is summarized in Table 2. Names of the atoms in glucose residues are given in Fig. 1.

Table 2 Description of the close contacts (below 5Å) between structural groups of CR molecule and those of cellulose

The vast majority of close contacts involves aromatic groups of CR and hydrophobic aliphatic groups of cellulose, in particular both faces of the glucose rings interact with the aromatic moieties of CR. This indicates that van der Waals forces, in addition to the electrostatic ones, also account for the adsorption strength. In agreement with the experiments, adsorption of dyes on cellulose is driven by van der Waals forces (Diaz et al. 2008). The order of magnitude of the strength of these van der Waals interactions is given by density functional theory (DFT) of l-tyrosine interacting with a methyl β-D-glucopyranoside, a system that can be considered as a model for the cellulose–CBM complex (Mohamed et al. 2010). Such calculations have been performed in vacuum, as we did, where the interaction energies values were estimated to be between −5 and −8 kcal/mol. Despite this main feature, the presence of short contacts between hydrogen atoms of cellulose and electronegative nitrogen and oxygen atoms of CR suggests the existence of hydrogen bonds. It is important to emphasize that most of the hydrogen bonds occur via the hydroxyl groups associated with the C2 and C6 carbon atoms of each glucosyl rings. Sulfonate, amino, and azo groups of CR form hydrogen bonds with the hydroxymethyl groups of cellulose. Even if the (1 0 0) surface has a dominant hydrophobic character, still a significant number of hydrogen bonds do occur between CR and cellulose. In addition to these “classical” hydrogen bonds involving polar groups, a significant number of close contacts involving aliphatic C–H groups interacting with an electronegative atom are also found. We found three kinds of them: (1) C–H···O, with lengths of 2.80 and 2.65 Å, also found experimentally (Bird et al. 2006), (2) C–H···N, with lengths of 2.47 and 2.95 Å, (3) N–H···O, length of 2.59 Å.

The forces involved in the interaction revealed by the present simulations are supported by the analogy with those measured in saturation transfer difference and line broadening NMR studies of the binding of cellohexaose to the CBM of Clostridium thermocellum (Viegas et al. 2008). In such a system, it was observed that in addition to the specific hydrogen bonds, the aliphatic protons of the glucosyl units were close to the tyrosine aromatic rings of the CBM (Viegas et al. 2008).

The description of the position and the orientation of adsorbed CR at equilibrium with respect to the cellulose surface, quantified by the three rotation and three translation parameters, are presented in Table 3.

Table 3 Equilibrium values of the orientation parameters between CR molecule and the (1 0 0) surface of cellulose at different temperatures and for different initial positions of a dye molecule in vacuum

At low temperatures (300 and 350 K), the roll and tilt angle parameters are small, suggesting that both planes of the CR molecule and of the top layer of the cellulose chains become more coplanar. In fact, the values of these two parameters were not rigorously zero, since as aforementioned the CR molecule was not fully planar and since the cellulose surface possessed a residual roughness. As for the twist angle, it became equilibrated to either the low values of 0–5° or the large ones of 46–48°. In addition, the shift and slide translational parameters were equilibrated at different values, suggesting that there were several adsorption sites on cellulose, in agreement with a recent study describing the adsorption of the CBM from Cel6A on crystalline cellulose (Shiiba et al. 2012). When adsorbed on the hydrophobic surface, the rise parameter, which describes the altitude of CR with respect to the cellulose surface is about 4 Å, a value similar to the one found for the adsorption of the CBM of Clostridium thermocellulm on cellohexaose (Viegas et al. 2008). The large values of the roll and twist angles observed for some of the high temperature simulations together with the negative values of the rise parameter illustrate this particular adsorption state of CR on both hydrophobic and hydrophilic surfaces. It occurs in this particular situation where the gravity centre of CR is lower than that of the surface chains of cellulose.

The molecule of CR is expected to be rigid due to the low number of rotatable bonds. However, the bond linking the central biphenyl moieties (the torsion angle φ in Fig. 2) is a source of flexibility. Systematic rotation of this torsion angle on a biphenyl model molecule (not shown) in the 0°–90° range revealed the existence of an energy minimum located at 48°. This result is in good agreement with that found using quantum mechanical calculations (40° using the semi empirical AM1 method, 51° estimated using HF with 6-31G* basis set (Pan et al. 2000) and 45° with MP2/6-31G basis set of DFT (Alves-Santos et al. 2006). This good concordance suggests that the DREIDING force field adequately models the conformations of such molecule.

The distributions of the φ angle for the isolated and the adsorbed states of CR are given in Fig. 7, the most frequently observed structures of CR are given in Fig. 8.

Fig. 7
figure 7

A torsion angle distribution of the central biphenyl angle (Φ) of CR from the MD simulations. The isolated state is displayed in red, the adsorbed state in vacuum in green and the adsorbed state in water environment in black. (Color figure online)

Fig. 8
figure 8

Selected examples of the conformations of CR adsorbed on cellulose. The structure in red has the φ angle at +138°, that in green at −40° and that in black at +40°. (Color figure online)

When isolated, CR explored two equally populated conformational states, at about +60°, and +140°. According to the rotational profile of biphenyl, the distribution of this angle should in principle display four equally populated states, symmetrically distributed around φ = 0°. However, the considered initial conformation of CR had its sulfonate groups on the same side of the molecule (see above) and the thermally activated transition was not observed at room temperature: only the positive values of the φ angle were consequently observed. Importantly, CR is far from planar as flat conformations are obtained when the φ angle takes the exact value of 0°.

When CR was adsorbed on cellulose, only one conformation, at about φ = + 38°, is massively populated. Its symmetrical conformation, at −38°, is only weakly populated and the conformation at +140° almost disappears. The most important feature is that CR is flattened when adsorbed on cellulose, but it is still not rigorously flat. Steric interactions between the hydrogen atoms of adjacent phenyl rings prevail over the effect of the overlap of their π-orbitals and a fully planar conformation is not accessible (Lye et al. 2000).

The sulfonate groups are on the same side of CR when adsorbed. The only exception is adsorption at 600 K, when in the final adsorbed state, both sulfonate groups are on opposite sides of the CR molecule. Interestingly, it has been experimentally measured that dyes in which the sulfonate groups are on one side of the molecule have a larger binding enthalpy for adsorption on cotton cellulose than those where the sulfonate groups are on the opposite side of the dye molecule (Bird et al. 2006). To check if modelling successfully reproduces this experimental result, we calculated the interaction energies between cellulose surface and CR for two different conformations of the dye molecule: one where both sulfonate groups are on one side of the molecule and one where they are on opposite sides. In this modelling, the dye molecule was positioned 4 Å from the cellulose surface and a short 50 ps dynamics was run: the interaction between the CR and the cellulose surface was stronger when the sulfonate groups were on one side, as opposed to the situation where they were on both sides, the difference being of 20.7 kcal/mol.

Adsorption process in aqueous environment

Since in reality the adsorption occurs mostly in aqueous environment, the effect of adding water molecules was investigated. To account for both the screening of electrostatic interactions and their potential structural role, we undertook a two-step strategy. In the first step, we studied the adsorption process of CR on cellulose with both sulfonate groups being solvated with 10 water molecules (SPC model). This procedure, called the semi-explicit water approach, was drawn from a recent treatment of the hydration, in which only the important water molecules were treated explicitly and the rest of the solvent phase was treated implicitly by a continuum approach (Fennell et al. 2011). Basically the same MD experiments were run as in the vacuum case, at a 300-650 K temperature range and in two different initial positions of CR, i.e., parallel and perpendicular to the cellulose fiber axis.

At room temperature, the process of adsorption differed from that in vacuum. At 300 K, the hydrated sulfonate groups were not the first chemical groups of CR to touch the surface of cellulose, but it was the aromatic rings of the CR molecule that created the first contacts. In contrast, at higher temperatures, water left CR quickly and thus the sulfonate groups were the first to approach the cellulose surface. The highest temperature, 650 K, was clearly not adapted in this context since the water molecules escaped from the sulfonate groups of CR. The water desorption temperature was then assumed to be at about 650 K.

The geometrical features of the adsorption are however preserved when water molecules are explicitly present in the simulated medium. At low temperature, CR adsorbed entirely on the (1 0 0) surface, and adsorption on the hydrophilic surfaces occurred at high temperatures. The equilibrium values of the orientation parameters between CR molecule and the (1 0 0) surface of cellulose at different temperatures and for different initial positions of the dye molecule, shown in Table 4, indicates also a strong similarity with that recorded for the vacuum trajectories.

Table 4 Equilibrium values of the orientation parameters between CR molecule and (1 0 0) surface of cellulose at different temperatures and for different initial positions of a dye molecule when both sulfonate groups were hydrated

The small values explored by the roll and tilt angles suggest that CR adsorbed flat on the cellulose surface. The twist angle adopts different values, revealing that CR can accommodate different stable orientations with respect to the cellulose fiber axis. A noticeable difference between trajectories performed in vacuum and when only sulfonate groups were hydrated is that the adsorbed dye molecule is more rigid in water compared to its relative flexibility in the vacuum phase (see the black curve in Fig. 7). In the explicit water environment, only the φ angle at 60° is populated.

The energies of interaction between CR and cellulose, when both sulfonate groups of CR are hydrated with 20 water molecules, were globally 20–30 kcal/mol higher than those observed under vacuum conditions (Table 5).

Table 5 Energies of interaction (kcal/mol) between CR molecule and Iβ cellulose in CR-cellulose complexes, when both sulfonate groups are hydrated with 20 water molecules, calculated for different temperatures (K) and different initial positions

Such lowering of the interaction strength between the substrate and the surface induced by the explicit presence of water molecules have been recently documented and discussed in the case of the interaction between xylan and cellulose (Mazeau and Charlier 2012).

In a following experiment, we studied the behaviour of a stable pre-formed complex established in vacuum under a full explicit water environment. At a temperature of 300 K for a parallel initial position of CR, the interaction energy between a dye molecule and cellulose moiety slightly increased to −92.5 (±10.2) kcal/mol, thus lowering the difference between the energies from vacuum and hydrated systems by a small number of water molecules. Additionally, no noticeable modification of the geometry of the complex could be observed.

These results suggest that although the interaction energy between CR and cellulose is estimated to be less efficient when explicit hydration is considered (either in a partial or a full manner), the basic associative features remain the same.

Concluding remarks

Our work has revealed the basic associative features of the interaction of CR with the hydrophobic surface (1 0 0) of cellulose Iβ. The results clearly show that CR indeed presents a real affinity for the hydrophobic surface of cellulose but unlike in the case of the well documented adsorption of aromatic amino acids, the adsorbed aromatic parts of the CR molecules are not laid completely flat on the cellulose substrate, but are slightly distorted. In fact, the CR molecules consist of two symmetric planar aromatic parts, connected by a central biphenyl bond, which acts as a rotatable joint. With isolated CR, the two parts become substantially rotated with respect to one another under the repulsion of the hydrogen atoms located in the ortho positions of the two adjacent central phenyl rings (Takei et al. 1988). In the present modelling study, we observe that the torsion angle along the biphenyl bond does not go to zero, but is drastically reduced when the CR molecules become adsorbed on the hydrophobic surface of cellulose. This conformational adaptation clearly indicates that in addition to the electrostatic and hydrogen-bonding interactions, the CR adsorption also results from a substantial van-der-Waals type attraction of its two aromatic moieties on the corresponding hydrophobic cellulose surface.

In the previous study, Woodcock et al. (1995) used a molecular mechanics approach to model the adsorption of CR on the hydrophilic surfaces of crystalline cellulose. In similarity with the present study, they also obtained a series of adsorption sites of minimum energy, indicating that possibilities existed for CR to also adsorb on such surfaces. To explain their results, they indicated that “the adsorption phenomenon was dominated the electrostatic interaction between the highly polar groups of the CR and the hydroxyl groups and exposed acetal oxygen on the cellulose surfaces” (Woodcock et al. 1995). Our modelling, which, in addition to the electrostatic phenomena, also takes advantage of the van der Waals interactions, is nevertheless different from their approach. At present it is not possible to compare their results with ours to decide whether CR will prefer our adsorption sites or theirs or will cover both. Work is in progress to answer such question.