When a proton attacks cellobiose in the gas phase: ab initio molecular dynamics simulations.

Investigations of reaction pathways between a proton and cellobiose (CB), a glucose disaccharide of importance, were carried out in cis and trans CB using Ab Initi o Molecular Dynamics (AIMD) simulations starting from optimized configurations where the proton is initially placed near groups with aﬃnity for it. Near and above 300 K, protonated CB (H + CB) undergoes several transient reactions including charge transfer to the sugar backbone, water formation and dehydration, ring breaking and glycosidic bond breaking events as well as mutarotation and ring puckering events, all on a 10 ps timescale. cis H + CB is energetically favoured over trans H + CB in vacuo , with an energy gap larger than for the neutral CB.


Introduction
Biological and industrial processes dealing with the conversion of polysaccharides to simpler sugars involve hydrolysis reactions, often mediated by enzymes and acids. [1][2][3] In many of these processes reactions of a proton with the sugar are assumed to take place. Nevertheless, while the importance of a proton interacting with a sugar is recognized, there is incomplete knowledge of the pathways of this interaction and its dynamics on a microscopic scale. The purpose of the present paper is to explore computationally possible elementary reactions that take place upon contact between H + and CB, an important saccharide that is the building block of cellulose.
As a context for our theoretical simulations we consider isolated CB in interaction with a proton, essentially under molecular beam conditions. Studies of sugars under such isolated conditions in the gas phase, pioneered by J. P. Simons and his school, have rendered a wealth of information on the structures of saccharide molecules, on their hydration and vibrational properties. 4,5 Such studies have not yet been extended to protonated sugar systems, with the exception of recent work on substituted b-galactose, combining experiment 6,7 and theoretical predictions from AIMD simulations. 8,9 However, it seems to the authors that such studies should be very desirable and may offer unique advantages to the understanding of proton-sugar interactions at the microscopic level.
In exploring reaction channels, the use of ab initio potentials seems to be essential, particularly in the context of dynamics and since standard force field molecular dynamics cannot account for reaction processes involving bond breaking. Ab initio molecular dynamics simulations (AIMD), in particular using density functional theory (DFT) with the BLYP functional, have been successfully utilized to describe peptide interactions and proton transfer processes, 10,11 for example in protonated peptides and in protonated complexes of amino acids. Qian et al. [12][13][14] were first to study protonated sugar monomers of glucose and xylose, by AIMD and recently by metadynamics. 15 The authors confirm a previously proposed hydrolysis mechanism for xylose 16,17 and propose several degradation pathways for glucose, involving water extraction. 15 Acidic dehydration and competing reaction pathways were studied computationally in b-glucose, and free energies for such reactions were characterized recently in great detail. 18 Glycosidic bond breaking was investigated by Liang et al. 19 in dilute HCl solution of trans CB (0.57 mol L À1 ), by means of metadynamics, and they report two possible pathways: one that first forms protonated CB (H + CB) and is followed by breaking of the glycosidic bond; the second in which protonation and dissociation of CB occur simultaneously.
In the present paper we study reactions between isolated cellobiose and a proton attacking several possible sites of that molecule. This choice of the saccharide is motivated by several considerations, First, cellobiose is a building block of cellulose and therefore may offer insights into corresponding processes in the extended polysaccharide. Secondly, computational studies of cellobiose in interaction with a proton were hitherto carried out only in solution and it is important to explore whether a range of processes have not been masked by the complexity of the environment. Indeed, we will be led to the conclusion that there is quite a number of different reaction channels between a proton and cellobiose, considerably more than that has so far been found in studies of protonated mono-saccharides, such as galactose and its derivatives. Last but not least, cellobiose offers the opportunity to explore the effect of different conformers on proton-sugar interaction. The influence of steric structure upon chemical reactivity has important implications for drug design and continues to be of fundamental interest for small systems as well as for macromolecules. 20,21 Sugars in general have very rich conformational structures, though close in energy, expressing different chemical properties. [22][23][24][25][26][27][28] The existence of the cis (or anti) and trans (or syn) cellobiose conformers and their relative stabilities has been a topic of considerable interest in the community of researchers focused on the study of saccharides. 22,23,25 Here we examine the effect of a proton on the relative stability of cis and trans CB and especially compare and contrast reaction channels for the two cases.
The structure of this article is as follows: the system description and computational approach are presented in Section 2; results of the calculations are brought and analyzed in Section 3 and concluding comments are found in Section 4.

Model and the computational approach
Specific systems were chosen for investigation on the basis of qualitative considerations. In CB, OH groups and the glycosidic O are sites likely for a proton to attack. These were not statistically sampled because of the considerable computational effort involved. Rather, for each of the species (cis and trans) a proton was placed near (within 1.8-2 Å) a moiety with affinity for it. In total, 9 different configurations were constructed for each of the 2 species. All initial guess positions were optimized to the local minimum, to ensure that species selected for further investigation by AIMD have chances to be experimentally observed. From among these, structures lowest in energy and sufficiently different from already chosen structures, in each species, were selected for further investigation by AIMD. For the selected conformers, we ran 16 trajectories, most at B300 K and a few at B40 K and at B500 K. A subset of these trajectories starts from a local minimum representing a state in which the proton had overcome the potential energy barrier. For such trajectories, overcoming the barrier from an initial condition maybe a time-demanding process, but the reaction after the barrier is surpassed is of considerable interest, hence a suitable initial condition for exploration via AIMD.
AIMD were carried out using the DFT implementation within the CP2K simulation package, 24,29 on the gas-phase NVE micro-canonical ensemble (simulation carried at constant energy) that consisted of CB and a proton. Potentials were computed at the BLYP 30,31 level of theory, and using the triple-z valence with the three polarization (TZV2P) basis set for geometry optimization (convergence tolerance o 1.0 Â 10 À7 Hartree) and for AIMD (convergence tolerance 1.0 Â 10 À6 Hartree). To address recent findings 32-35 that cast doubt on the accuracy of BLYP functionals, with respect to their ability to predict accurate energetics of a system and kinetics of reactions involving proton transfer, we proceeded to validate (in a subset of our structures) the conformational energy rankings against those obtained with the B3LYP functional and using a comparable basis set, namely 6311++G** (3df, 3pd with diffuse sp). In all cases we used the Goedecker, Teter and Hutter (GTH) type pseudopotential (with plane waves as basis set), 36 Grimme's dispersion correction at the D3 level, 37,38 and a kinetic energy cutoff of 300 Ry, over a cubic cell with a 32 Å side. Trajectory calculations starting from geometry optimized structures were carried out for intervals of up to 12 ps at a time step interval of 0.5 fs. Mulliken 39 partial charges were recorded for all trajectories. Standard notation of CB atoms is illustrated in Fig. 1 (top), with OH groups numbered by the number of the C atom to which they are connected, and will be used throughout this paper, except that the glycosidic O will be referred to as Og, and the ring oxygen will be referred to as Or or Or 0 .

Structure and energetics
The set of optimized configurations selected for further investigation by AIMD are listed in Table 1. These were selected either because of the energy ranking within the groups (cis or trans) or because, as in trans H + CB(3), we suspected that an interesting reaction, namely glycosidic bond breaking, might be observed during AIMD simulations, see Fig. 1. The table shows the energy rankings relative to the configuration of lowest energy. The conformer column indicates, in parenthesis, the initial placement of the proton.
The optimization carried at the B3LYP level of theory serves as a semi-quantitative validation of the BLYP potential. We note that although there is a large difference (more than double) between relative energies obtained at the BLYP/TZV2P level and those obtained at the B3LYP/6311++G** level, the energetic ranking of the configurations is preserved and so is the conformational structure (see RMSD). Structures listed in Table 1 (obtained with BLYP and B3LYP) are provided in the ESI † for comparison purposes.
The configurations are shown in Fig. 1. It is worth noting that: in cis H + CB configuration 1 the excess proton is trapped between groups OH6 and OH6 0 for considerable parts of the trajectory, which is expected to affect the vibrational IR spectrum of these OH groups; cis in configuration 2 has already resulted in a mutarotation of OH1 groups from the equatorial orientation to an axial orientation above the plane of the ring; cis in configuration 3 has one H shared between OH3 and OH2 0 , i.e. an intra-molecular proton bound dimer, for good parts of the trajectory, and that too may affect the vibrational IR signature, but in addition, the bond C2 0 -O2 0 is elongated, hinting at possible dehydration events during AIMD. In the case of trans H + CB, configuration 1, the excess proton migrates to Or 0 , thus weakening the ring through elongation of the O5 0 -C1 0 bond; trans in configuration 2 traps the excess proton between OH4 0 and OH3 0 , in a motif similar to that present in cis H + CB(3) and in other configurations not shown here; finally, trans in configuration 3 was chosen because the excess proton migrates to within 1.1 Å from the glycosidic O, hinting that such bond breaking might be observed during dynamics at elevated temperature.
The most notable and surprising finding of this part of the study is that the protonated cis species of cellobiose is energetically more stable than protonated trans cellobiose, to an extent far greater (B18 kcal mol À1 ) than for the corresponding parent neutral species in vacuo (B3 kcal mol À1 , at similar levels of theory). 22,23,25,26 It is of course not possible to be confident, despite extensive searches, that we have found the global minimum. An observation worth noting is that when comparing the lowest energy conformations of the neutral species to lowest conformers found here, in the protonated cis species, the proton forms an intra-molecular proton bound dimer between OH6 0 and OH6, which renders the angle at Og more acute as compared to the neutral counterpart, while in trans H + CB(1), the proton binds to Or 0 opening the ring and flattening the angle at Og more than in the neutral species, perhaps enhancing the energetic difference. Nevertheless, the enhancement of cis over trans in the isolated cellobiose, by protonation, is of considerable interest in exploring the balance between the two species also in different environs. While it is known experimentally as well as theoretically 40-50 that cis is more stable than trans CB in the isolated species and also in small water clusters [25][26][27] in the gas phase, in neutral condensed phase environs, both crystal and in solution, the energetic balance favours the trans species. The reasons for these structural preferences are as yet not understood and continue to be a topic of considerable research interest. However, in light of our findings, cellobiose in the condensed state, when placed in a strongly acidic environment (perhaps consisting of a solvent other than water, which competes strongly for the proton), might present a different picture of the ratios of trans versus cis species than expected based on the neutral species abundance ratios, estimated by Larsson in water solution to be 93% trans CB and 7% cis CB conformations, at 300 K. 49

Reactions and other dynamical processes
In the following sections we discuss the most interesting events observed during the dynamics simulations: ring breaking events, glycosidic bond breaking and dehydration events. Puckering and mutarotation, which were also observed will be discussed in the context of the main events.
2.1 Ring breaking in cis H + CB(1). Dynamics at B300 K for the lowest energy protonated cis conformer reveals a fast reaction on a time-scale of 1 ps, the intermediate steps of which are captured in Fig. 2. In it, water initially forming at O6 transfers a hydrogen to the ring oxygen -Or 0 (see Fig. 3a), leading to the opening of the secondary ring at C1 0 -Or 0 (see Fig. 3b); this opening persists throughout the simulation, and is accompanied by formation of an intermittent double bond between C1 0 and  ) and the water at O6, which forms periodically throughout the dynamics, see Fig. 3c and d, indicating that the process has not settled.
2.2 Ring breaking in trans H + CB (1). In this, the lowest energy protonated trans conformer investigated (B18.12 kcal mol À1 higher in energy at the B3LYP level than the most stable cis conformer found, see Fig. 1), a barrierless and spontaneous reaction already took place at t = 0, leaving the secondary ring open at C1 0 -Or 0 (on average B1.7 Å) for the entire duration of the trajectory, see Fig. 4a. As in the previous case of cis H + CB(1), where C1 0 -Or 0 opened up within the first ps, the consequence here too is formation of a persistent double bond between C1 0 -Og (Fig. 4a, in blue, on average at 1.35 Å) and subsequent weakening (or elongation) of the other side of the glycosidic bond (Fig. 4a, red).
The partial charge trajectory at Or 0 (Fig. 4c, in red), on the other hand, exhibits fluctuations in a much narrower range than for the corresponding atom in the cis H + CB(1) case, which is consistent with this apparent barrierless reaction. It is also evident from the other partial charge trajectories in the panel that the charges on OH3, which remain loosely bound throughout the simulation (see Fig. 4b), contribute strongly to the stabilization of this conformer. Panel b of Fig. 4 also illustrates a conformational change taking place, at t B 4 ps, which brings OH6 0 and H3 closer, when the OH6 0 group rotates in a counter clockwise direction. Fig. 5 shows several snapshots along the trajectory: before and after the conformational rotation of OH6 0 and 2 snapshots taken near events of maximum separation between C1 0 -Or 0 , where a strong double bond forms at C1 0 -Og. The charge of both the excess proton and of H3 increases slightly after this transition (see Fig. 4c), with H3 retaining overall the most positive charge (among H + , H6 0 and H3), while correspondingly, Or 0 being the least negatively charged (among Or 0 , O6 0 and O3) throughout the trajectory.

Glycosidic bond breaking in trans H + CB(3).
Glycosidic bond breaking is one of the most important reactions and is of interest to a large community of scientists bound on unlocking the vast amounts of energy stored in cellulose -the main building block of plant cell walls. 51,52 In cellobiose, the fundamental building block of cellulose, a dimer of two glucose rings is linked by the glycosidic oxygen. Several recent computational studies have explored the glycosidic bond breaking for cellobiose in weak acidic solutions. 19,53 Not much is known about this process for sugars in the gas phase, yet recent experimental advances 6-8 point to the possibility that this process can be explored and confirmed in the gas phase, by mass spectroscopic means. In this conformer, trans H + CB(3), resulting from positioning a proton B1.8 Å from O6, the excess proton exchanges places with H6, while the later migrates to bind to the glycosidic oxygen in the final optimized structure. This structure is B24 kcal mol À1 higher (at the b3lyp level of theory) than the lowest cis structure.
In the first 3 ps of the trajectory, but also at later times, H6 binds to Og, a seemingly barrierless process. Snapshots presented in Fig. 6 describe significant events occurring in the first ps of the trajectory: (a) the glycosidic bond C1 0 -Og breaks; (b) as the C1 0 -Og bond reforms, the other side, Og-C4, opens; (c) while C1 0 -Og is open, the secondary ring rotates 901 relative to its original configuration, with the planes of the rings reaching perpendicular orientation to each other, on a transition path to cis, and reaching a maximum separation of 1.76 Å at t B 0.55 ps.
During this period, in which the rotation takes place, the partial charge at Og is least negative, on average À0.29 a.u. (see Fig. 6b, black), while the total charge over the primary ring fluctuates around 1 a.u. throughout the trajectory (Fig. 6b,  blue). The total charge on the secondary ring plus the charge on Og, on the other hand fluctuates around neutral (Fig. 6b,  red). Could it be that the neutral charge over the secondary ring plus Og permits the ring rotation (shown in Fig. 7)? Is the charge at Og responsible for keeping the pyranose rings together?
The charge and bond trajectories presented in Fig. 6 reveal other dramatic transient events taking place. At t E 3 ps, the charge at Og drops abruptly to À0.55 a.u, where it remains, except for 2 events centered at t B 4.4 ps and t B 6.9 ps, when the Og partial charge again reaches maximum; echoes of these events are also reflected in the bond trajectories of Fig. 6a. Two snapshots taken at the transition where the Og charge drops abruptly (before and after), indicated by the arrow labelled 1, are shown top right in Fig. 6, indicating this to be the result of Og losing the bond to H6. Conversely, the 2 snapshots taken during the short events when the Og partial charge raises again to values it had at the start of the trajectory (at t B 4.5 ps and t B 7 ps), indicated by arrows labelled 2 and 3, and shown under the ''transition state'' snapshots in Fig. 6, confirm that it is the bond to H6, which reforms during these two brief periods, that is responsible for lowering the Og charge.

Dihedral angle dynamics.
More detailed examination of the trajectory shows that, for example, as the secondary ring rotates back (1 o t o 3 ps) to some semblance of the initial structure, significant puckering 54 events take place. As an example, shown in Fig. 8 is the dihedral formed by C2 0 C1 0 OgC4. It shows two significant sharp reconfigurations, at t B 2.3 ps (lasting B80 fs) and at t B 4 ps. Snapshots for these events are presented in the figure above the plot, but although the snapshots were selected to represent equivalent dihedral angles, they are very different: note the Og bond to H6 at t E 2.3 ps and its absence in the other snapshot, and also note the different puckering configurations of the secondary ring.
In summary, it appears that in the gas phase (1) glycosidic bond breaking can be very fast in trans species, which after all is less favoured energetically; (2) the reaction is associated with a trans to cis conformational transformation that (3) may be enabled by the very different charge distribution over the two pyranose rings of CB.
2.4 Dehydration processes. It turns out from this work that water formation processes are some of the most important features of chemical attack on sugars. For the most part, the water molecules do not leave the sugar within the simulation time-scale. Rather, the water formation is transient: formed water  This journal is c the Owner Societies 2013 Phys. Chem. Chem. Phys., 2013, 15, 15382--15391 15387 reacts again with the sugar backbone, which is positively charged, then starts to break away and the process oscillates between these states. We shall see later, however, that events exist in which the formed water molecule departs from the sugar, leaving an isolated ion, presumably analogous to the carboxonium, observed in other sugar monomers in an acidic environment. 8,[12][13][14][15]19,53 Following are examples of such processes. (2). The initial state of this conformer can be characterized by a rigid bridge of hydrogen bonds (HB) connecting OH6 0 to OH6 to OH1, through the sharing of the excess proton -a proton wire. Fig. 9 describes a dehydration reaction taking place at B4.9 ps, as the excess proton transfers to O1 (Fig. 9a, black) and the bond C1-O1 is severed (Fig. 9a,  blue), reaching a maximum separation of 2.9 Å, at t B 5 ps (and averaging B2.2 Å thereafter). The total charge on the water forming at O1 drops on average to B0.1 a.u. but continues to fluctuate between B0 a.u and 0.3 a.u (Fig. 9c). The fact that the total charge at the water periodically drops to neutral indicates that the charge transfer reaction had completed, but water is  still held in place by dipole-charged dipole interactions. The process is accompanied by the charging of the sugar backbone, leaving it in a much more reactive state: for example, the link between C1-Or decreases to a distance characteristic of the double bond, less than 1.3 Å (Fig. 9b, red), while the bond on the other side, namely Or-C5, gets elongated (average of B1.6 Å Fig. 9b, black), a motif observed in the ring breaking events noted earlier and observed by others in acid catalyzed water abstraction in b-glucose 12,18,19 and in protonated galactose. [6][7][8][9] The snapshot labelled 2 (corresponding to the arrow in panel Fig. 9a), taken at t B 5 ps, illustrates water reaching maximum separation from the ring (C1).

cis H + CB
2.4.1.1 Mutarotational event. The charge trajectories, in particular for Or, O1 and O6, and the bond length trajectories for O1-H + and O6-H + (Fig. 9a), show an additional transient event occurring at t B 0.88 ps (indicated by the arrow labelled 1 in panel a of Fig. 9), which upon inspection, corresponds to a mutarotational conversion, shown in the snapshot labelled 1, in Fig. 9. The event is characterized by a maximum distance between the excess proton and the OH1 group, of 2.63 Å, which is localized 0.98 Å from O6; this is related to a rotation in which the dihedral angle for O1-C1-C2-O2 undergoes a mutarotation, with OH1 above the glucose ring and perpendicular to it, and with OH2 aligned in the opposite direction, perpendicular to the glucose ring but underneath it. This event is transient, lasting B80 fs, during which the OH1 and C1 are closest. In fact, a similar transformation was observed in b-glucose, exposed to hydronium in the gas phase (at B300 K), 12 which converted to a-glucose along an alternate pathway to the primary path -formation of carboxonium ions (and with a free energy barrier of about 10 kcal mol À1 higher than that for the dominant step). During the course of the dynamics ring puckering transitions are frequently observed. In view of the work of B. J. Smith, 59 such events should be expected, as they were shown to reduce the energy of the glycosidic bondstretch and as a result, lower the transition state energy for bond cleavage in protonated glycosides.
Dehydration events can also be observed in cis H + CB(3) but these are transient events and no definitive transitions occur during the 9 ps course of the trajectory. It may be of some significance to note that acid catalyzed dehydration studies carried out in b-glucose, in water 18 and in the gas phase 12,15 identify the anomeric glucose site as the most reactive or most  Phys. Chem. Chem. Phys., 2013, 15, 15382--15391 15389 likely to undergo dehydration, with a free energy of transition of B3.5 kcal mol À1 at C1 versus 18.4 kcal mol À1 for a transition at C3. 15 2.4.2 trans H + CB (2). A similar situation is obtained in the case of trans H + CB(2), where the excess proton was initially placed near OH4 0 , but in this case, there is less indication that water formed may be departing in the course of the trajectory, though the bond O4 0 -C4 0 is considerably weaker throughout the simulation, see Fig. 10c, red. The charge distribution at the water differs from that in the previous case: here, water retains a considerable fraction of positive charge (see Fig. 10a). In addition, the excess proton is shared by OH4 0 and OH6 0 groups, forming a proton bound complex. This motif has appeared in other situations [55][56][57][58] and was identified spectroscopically. [6][7][8] It is analogous to proton wires observed in protein and peptides and we expect it to produce a spectroscopic manifestation similar to that observed in the case of protonated galactose. 7 Other events, similar to those observed in cis H + CB(2), such as the intermittent breaking of the Or 0 -C1 0 bond and the shortening of the bond on the other side of Or 0 (Or 0 -C5 0 ) to lengths characteristic of the double bond (1.34 Å at minimum) are also observed. Puckering changes of the pyranose rings occur during the length of the trajectory as well as angle changes between the planes of the pyranose rings, which render the charged sugar more planar and the glycosidic bond more exposed.
2.5 Water evaporation and formation of an isolated carboxonium-like sugar ion. This is one case where we could observe water evaporation on the time scale of the simulation. As a consequence of placement of a proton near OH3 0 , water forms and separates from the ring at C3 0 to a distance of 2.24 Å. The departing water in this case has 0 charge and can therefore