QM/MM Molecular Dynamics Study of the Galactopyranose → Galactofuranose Reaction Catalysed by Trypanosoma cruzi UDP-Galactopyranose Mutase

The enzyme UDP-Galactopyranose Mutase (UGM) catalyses the conversion of galactopyranose into galactofuranose. It is known to be critical for the survival and proliferation of several pathogenic agents, both prokaryotic and eukaryotic. Among them is Trypanosoma cruzi, the parasite responsible for Chagas' disease. Since the enzyme is not present in mammals, it appears as a promising target for the design of drugs to treat this illness. A precise knowledge of the mechanism of the catalysed reaction would be crucial to assist in such design. In this article we present a detailed study of all the putative steps of the mechanism. The study is based on QM/MM free energy calculations along properly selected reaction coordinates, and on the analysis of the main structural changes and interactions taking place at every step. The results are discussed in connection with the experimental evidence and previous theoretical studies.


Introduction
Chagas' disease, also known as American trypanosomiasis, affects approximately 8 million people worldwide. It is endemic in Latin America but in the last decades it has also spread towards North America and Europe [1]. Its pathogenic agent is the flagellate protozoan Trypanosoma cruzi (T. cruzi), which is transmitted to humans by the faeces of triatomine insects. The disease was first described by Dr. Carlos Chagas in Brazil in 1909. Despite this early discovery there are still no drugs capable of curing it. Nifurtimox and Benznidazole are used in the acute phase of the disease. However none of them are efficient and both have strong side effects [2][3][4]. Most patients discontinue the treatment when the side effects become too severe. For these reasons new and more efficient drugs are needed.
Galactose is a common monosaccharide. In mammals it is exclusively found as galactopyranose (Galp), the six-membered ring hemiacetal form. On the other hand, in T. cruzi and many other human pathogens such as Mycobacterium tuberculosis, Escherichia coli, Leishmania major, Aspergillus fumigatus, Salmonella typhimurium and Klebsiella pneumoniae [5][6][7][8], it is found as galactofuranose (Galf), the five-membered ring hemiacetal form [5,6,[8][9][10][11][12][13][14][15][16][17]. The sole source of Galf in these species is the enzyme UDP-Galactopryranose Mutase (UGM), which catalyses the isomerization between UDP-Galp and UDP-Galf, the precursor of Galf [18,19]. It is known that Galf is an essential component of the cell wall and the extracellular matrix of these pathogens [13,19]. Suppression of the UGM gene in many of them caused attenuated virulence and increased sensitivity to drugs [20][21][22]. In T. cruzi, particularly, Galf is attached to glycoinositolphospholipids and glycosylphosphatidylinositol anchor proteins [23,24], which are highly expressed throughout the life cycle of the parasite and are essential for its survival and proliferation [25][26][27]. When T. cruzi is incubated with specific antibodies against Galf, the binding of the parasite to the mammalian cells is blocked, leading to an 80% decrease in infectivity [13]. Since neither Galf nor UGM have ever been found in mammals, UGM has gathered significant interest as a target for drugs design [28]. Due to this interest, it has been subjected to several structural and mechanistic studies [28][29][30].
In 2001 was presented the first known crystallographic structure of a UGM. It corresponded to E. coli, [31]. After that, other bacterial structures were also determined [31][32][33][34][35][36]. Eukaryotic UGMs received less attention. The first structure of that kind, corresponding to Aspargillus fumigatus, was published in 2012 [37]. Shortly after, the one of T. cruzi (TcUGM) became also available [38]. The comparison between eukaryotic and prokaryotic UGMs revealed that they share a common folding and a GxGxxG motif, necessary to bind the cofactor, flavin adenine dinucleotide (FAD) [39]. Moreover, the cofactor conformation and its interaction with the enzyme environment is highly conserved in both groups. However, the interactions with the substrate differ significantly and the sequence identity is pretty low (15%) [38]. In the active site, only 5 out of 13 residues are shared. Besides eukaryotic UGMs are approximately 100 residues longer than prokaryotic ones. This additional part of the chain forms extra secondary structures, modifying the active site flexibility and the oligomerization state of the enzyme [39]. Fig. 1 shows the main species of the catalysed reaction. The transformations between these species we will be denoted as ''stages'' of the mechanism. The first and last stages consist of just one reaction step while the second and third stages involve two. All the steps of the mechanism under analysis are presented in Fig. 2. According to different experimental studies the reaction initiates with the formation of a flavin-galactose adduct (conversion from I to II in Fig. 1) [28,34,40,41]. This requires the rupture of the Galp-UDP bond and the creation of a bond between Galp and the nitrogen at position 5 of the reduced flavin adenine dinucleotide (FADH), N5 FADH [40,[42][43][44].
It was experimentally found that no conversion between Galp and Galf occurred when the native cofactor was replaced by 5deaza-FAD [45]. Since this modified cofactor can only participate in two-electron transfers, it was argued that the mechanism in UGM should involved a one electron transfer. In particular, it was suggested that an oxocarbenium ion was first formed, followed by a single electron transfer, and that the recombination of the radicals so formed would produce the flavin-galactose adduct. However, it was then argued that the evidence presented does not exclude the possibility of a nucleophilic attack of N5 FADH onto the anomeric C of Galp, C1X GAL , with a S N 2 type mechanism [46]. Positional isotope effects experiments, together with studies that employed FAD analogues with different electron density on N5 FADH , uphold this hypothesis [46]. Besides, the analysis of the crystallographic structures, as well as recent investigations on TcUGM, give further support to this mechanism [28,34,40,43]. The next stage (conversion from II to III in Fig. 1), involves the opening of the ring to form an iminium ion [40,43]. This intermediate species has been trapped using NaCNBH 3 in two independent studies [28,40]. Naively, one would suggest that the iminium is formed by a direct proton transfer from N5 FADH to the cyclic oxygen of galactose, O5X GAL . However, as noted by Huang et. al., such transference involves the passage through a fourmembered ring structure which is rather high in energy. As an alternative, the same authors proposed that the proton is first passed from N5 FADH to O4 FADH , and then transferred to Galp to initiate the opening of the ring [41]. Once the iminium intermediate is formed, two stages are needed to complete the reaction. They can be considered as the reverse of the two previous stages, except for the fact that galactose is now in the furanose form. Thus, stage three involves the sugar ring closure to form Galf (conversion from III to IV in Fig 1). Sobrado et. al. indicated that this is the stage that determines the rate of the whole process [28]. Stage four consists of the breaking of the flavin-substrate bond along with the binding of UDP to the sugar (conversion from IV to V in Fig. 1).
Huang et. al. performed a theoretical study on the mechanism of the reaction catalysed by UGM [41]. They carried out electronic structure computations on active site models built from the PDB structure of Klebsiella pneumoniae UGM (KpUGM). The largest of their models contained 26 active site residues plus the substrate, the cofactor and several crystallographic water molecules. A quantum mechanics-molecular mechanics level theory (QM/MM) was employed to characterize the structures of reactants, products, intermediate species and transitions states appearing in the mechanism. More recently, the involvement of several active site residues on the catalytic activity of TcUGM was evaluated through site directed mutageneis experiments [29].
In this article we present a QM/MM molecular dynamics study of the reaction catalysed by TcUGM. We applied the umbrella sampling technique to obtain the free energy profiles along different reaction coordinates, conveniently defined to describe every step of the mechanism. QM/MM free energy computations have become a widely employed tool to gain information on the atomistic details of enzymatic reactions. One of their main assets is the ability to reveal both, energetic and dynamical contributions to catalysis. We also analysed the most significant conformational changes and interactions taking place at each step. This includes the monitoring of bond distances, dihedral angles, H-bonds, partial charges, bond orders as well as the Cremer-Pople angles that describe the conformations of the pyranose and furanose rings [47]. Finally, we implemented an energy decomposition method to evaluate the contribution of the active site residues to the lowering of the barriers at every step. The results of the simulations are discussed in connection with previous experimental findings, as well as with the theoretical analysis of Huang et. al.

Results and Discussion
In Fig. 3 we present a sketch of the free energy changes (DG 0 r ) and free energy barriers (DG { r ) for the successive steps of the mechanism presented in Fig. 2. The profile shows that the barrier for the ring opening (step 3) is sensibly smaller than that of the ring closure (step 4). In fact, the barrier for step 4 is the highest. This is in agreement with the experimental findings of Sobrado et. al. [28]. The profile also indicates that products are more stable than reactants. The same result was found in the computations of Huang et. al. [41]. For the reverse reaction the largest barrier corresponds to the tautomerization of FADH. We also note that for both, forward and backward reactions, the appearance of the iminium ion species presents a small barrier.
In the following sections we describe in detail the outcome of the QM/MM computations for all the stages of the catalysed reaction. When pertinent, the results are compared with those recently reported for KpUGM [41]. We note, however, that a meaningful  . Detailed mechanism for the reaction catalysed by Tc UGM. The mechanism includes the intermediates detected by experiments as well as those whose existence was inferred from theoretical considerations. Red color is used to denote the bonds being broken (solid line) or formed (dashed line), as well as the atoms involved. The distances between these atoms are labelled because they are used to define the reaction coordinates. doi:10.1371/journal.pone.0109559.g002 comparison between these computations requires keeping in mind the aspects in which they differ. Among the differences we have: (1) that KpUGM and TcUGM bear low sequence homology: 18.3% for the whole protein and only 5 out of 13 residues for the active site [39]; (2) that free energy computations include dynamical effects that are not considered in electronic structure computations; (3) that residue His62 was protonated in the present work but was set as neutral in the work of Huang et. al.; (4) that we modelled the whole TcUGM crystal structure in explicit solvent while Huang et. al. considered an active site model consisting of the cofactor, the substrate, 26 active site residues and 8 water molecules; (5) that the quantum subsystems and the DFT levels of theory employed in each computation were different.
In what follows the numbers of the residues and the names of the atoms correspond to the crystal structure of TcUGM taken from the Protein Data Bank (PDB code 4DSH). Figure 2 should serve as a guide for the reading of the following sections. It describes every step of the mechanism, highlighting the distances involved in the definition of the corresponding reaction coordinates. On the other hand, at the supplementary information section, we provide pdb files with the average structures of reactants, products and all the intermediates of the reaction (Text S1 to Text S7). These structures can be used to obtain information that was not included in the main text in order to keep the article at a reasonable length. For the same reason, several pictures depicting the evolution of important distances and angles along the different steps of the reaction are given at the supplementary information section.

Stage 1: Formation of the flavin-Galp adduct
This stage consists of just one concerted step in which N5 FADH bonds covalently to the anomeric C of Galp while the UDP moiety detaches from it. In panel a) of Figure S1 the evolution of the distances involved in the definition of the reaction coordinate (z 1~d2 {d 1 ) along this step, are shown. At the transition state, TS1, we obtained d 1 = 2.52 + 0.05 Å and d 2 = 2.15 + 0.05 Å . This corresponds to bond orders of 0.30 and 0.19 for the C1X GAL -O3B UDP and C1X GAL -N5 FADH bonds, respectively. These orders support the hypothesis that the reaction proceeds via a dissociative S N 2 mechanism, as has been suggested by several independent experimental studies [28,34,46]. The calculated DG  Table 1. It is observed that the cyclic oxygen, O5X GAL , losses considerable electron density in going from reactants to TS1, but it partially recovers it when the adduct is finally reached. C1X GAL , on the other hand, gains substantial electron density along the whole process. Finally, the partial charge of N5 FADH increases from 20.18 to 0.14 while its configuration changes from planar to tetrahedral. We note that the substantial loss of electron density of the nucleophile nitrogen in this step was predicted by the experiments in which FAD analogues with different electron-withdrawing/donating groups where used to determine the S N 2 character of this step [46]. These changes weaken the N5 FADH -H bond facilitating the transference of the proton during the next step. The evolution of the Cremer-Pople angles is shown in Fig. S1 panel b). At the Michaelis complex, the pyranose ring shows a 4 C 1 conformation and its h angle oscillates around 0 0 . However, after surmounting the transition state, this angle increases by ,50 0 while w diminishes to -60 0 . This correspond to a 2 H 3 six membered ring conforma- Table 1. Partial charges of key atoms for the first two steps of the mechanism. tion. This conformation persists until the formation of the flavingalactose adduct, when the pyranose ring returns to 4 C 1 .
The phosphate group of UDP bears strong H-bond interactions with Tyr395 and Tyr429 during the whole step. It also forms a Hbond with Arg327 but, at the Michaelis complex, this interaction is rather weak. However, once the covalent bond between UDP and the sugar is broken, the interaction gains strength because of the negative charge acquired by the reactive oxygen of the phosphate (see Table 1). Thus, while only 32.3% of the structures sampled before TS1 present a H-bond between Arg327 and the phosphate, the percentage raises to 69.7% for those sampled between TS1 and products. This indicates that Arg327 plays an important role in stabilizing TS1, as well as the products of the current step. Further support for this conclusion comes from Table 2 which shows that Arg327 has the most negative DE i R?TS value. The location of this arginine within the active site can be seen in Figure 4. Arg327 is conserved in UGMs of both, prokaryotic and eukaryotic organisms [39]. For KpUGM, it was found that its substitution by Ala completely abolishes the enzyme activity [48]. For TcUGM the same substitution was found to reduce the catalytic activity, measured by k cat =K M , in 50% [29]. This effect is due to a reduction in k cat since the mutated enzyme has a larger affinity by the substrate, as indicated by the decrease in K M [29]. Our results agree with these experimental findings and suggest that the effect could be due to the role played by Arg327 in the catalysis of the flavin-galactose adduct formation. Finally, we note that the H-bonds between the phosphate and Tyr395, Tyr429 and Arg327 were found to be present in all of the remaining steps. Thus, they restrain the mobility of the phosphate group so that it is ready to re-bind the sugar moiety once the furanose ring is formed. This fact could explain the detrimental effect on k cat observed when Tyr395 and Tyr429 are substituted by Ala [29], and could also contribute to the negative effect observed upon the substitution of Arg327 with Ala.
Another interaction that is worth mentioning is the H-bond between the H atom bonded to N5 FADH and the carbonyl oxygen of Gly61. This interaction already exists at the reactants configuration but becomes stronger once TS1 is reached. The distance between the carbonyl oxygen of Gly61 and the H atom is 2.30 + 0.39 Å before TS1, but decreases to 1.89 + 0.11 Å after it (see Fig. S1 panel b). Because of this, the H-bond is present in 26.4% of the configurations sampled before TS1 and 59.2% of those sampled afterwards. Since the interaction is stronger for TS1 than for reactants, it certainly helps to reduce the barrier to reaction. Unfortunately, the stabilizing effect of glycine residues cannot be evaluated with the energy decomposition method employed in this work. Because of this, Gly61 is not mentioned in Table 2. It has been found that the replacement of Gly61 with Ala or Pro has a profound detrimental effect on the activity of TcUGM (,90%) [38]. A putative explanation for this fact would be that these alternative residues reduce the flexibility of the backbone chain, hindering its ability to locate the carbonyl group in an appropriate position for the H-bond interaction. However, the evaluation of this hypothesis requires further MD computations that are outside the reach of this work.

Stage 2: Formation of the iminium ion
This stage consists of two steps. The first one involves an intramolecular proton transfer from N5 FADH to O4 FADH (see Fig. 2 Both calculations agree to indicate that the products of this step are more stable than the reactants. However, the energy difference is smaller in our computations and our estimated barrier is sensibly higher. The large barrier is not surprising since the initial distance between the H atom being transferred and the acceptor is rather long, 2.46 + 0.8 Å , while the donor-hydrogen-acceptor angle is far from collinear, 98.7 + 7.1 0 . Because of this, a large configurational change needs to take place to enable the transference. The probability of H-tunneling for this step was found to be negligibly, with a transmission coefficient of 0.002 calculated at the zero point energy of the N-H bond. This results is not surprising since the barrier width is rather large. In Fig. S2a we show the variations of the distances involved in the definition of the reaction coordinate z 2~d3 {d 4 during step 2. At the transition state TS2 the distance between O4 FADH and the H atom has decreased by ,1.3 Å but the N5 FADH -H distance has not increased substantially. On the other hand the donor-acceptor distance, presented in panel (b), decreases by 0.36 Å in going from reactants to TS2. The analysis of the structures reveals that this shortening is mainly caused by the twisting of the C5X-N5-C4X-C4 torsion of FADH, also shown panel (b). This observation is confirmed by a correlation coefficient of 0.87 between the donoracceptor distance and the corresponding dihedral angle. Thus, the ability of the isoalloxazine ring to bear this distortion contributes to make possible the transfer. In addition, Fig. S2b shows the distance between the H atom being transferred and O5X GAL . This distance reaches a minimum of 1.5 Å when the transference is complete generating a new H-bond between O4 FADH and O5X GAL . As can be seen in Table 1, the partial charge of O5X GAL decreases along this step while that of the H atom increases. Both, the new H-bond and the variations in the partial charges make the scenario more prone to the attack of the proton onto O5X GAL that initiates the ring opening at the following step. The partial charges of the N5 FADH , C4X FADH , C5X FADH , C4 FADH and O4 FADH also suffer noticeable variations along this step. Basically, the donor and the acceptor of the proton gain electronic density in going from reactants to TS2, while C4X FADH , C5X FADH and C4 FADH loss it.
Residue His62 lies pretty close to the isoalloxazine ring of FADH (See Fig. 4). They attract to each other since His62 has a positive charge while FADH has a negative charge. Along step 2 His62 moves closer to FADH and its side chain rotates so that the  Table 2. It is observed that His62 exerts the largest stabilizing effect on TS2, 216.35 kcal/mol, being the main responsible for the acceleration of this step. The involvement of His62 in the catalysis of step 2 could explain why the substitution of this residue by Ala diminishes the catalytic activity of TcUGM in 98% [38].
The following step (step 3 of Fig. 2) involves the transfer of the H atom from O4 FADH to O5X GAL to open the sugar ring and form the iminium ion. The calculated DG { r and DG 0 r for this process are 13.7 + 0.2 and 21.1 + 0.3 kcal/mol, respectively. Experimental results suggested that the opening of the ring, detected as the appearance of the iminium ion, is fast. [28] In agreement with the experimental evidence we found that this step presents a relatively small barrier. Huang et. al. estimated an energy barrier rather similar to our DG { r , but their DE 0 r (8.12 kcal/ mol) is quite distinct to our DG 0 r estimation. Fig. S3 shows the evolution of the distances involved in the definition of the reaction coordinate z 3~d5 zd 6 {d 7 . It is observed that at TS3 the proton is halfway between the donor and the acceptor, while the bond between O5X GAL and C1X GAL is partially broken. When this bond gets completely broken and the ring opens, the iminium ion is formed. During this process the distance between C1X GAL and N5 FADH decreases from 1.50 + 0.04 Å to 1.34 + 0.03 Å because the bond between them acquires a partial double bond character. In agreement with experimental findings [46], it was found that N5 FADH losses considerably electron density during this step (See Table 1). Besides two out the three atoms bound to this nitrogen (C1X GAL and C4X FADH ) gain electron density. The analysis of the TS3 stabilization pattern, presented in Table 2, shows that none of the active site residues has a significant influence on the energetic of this step. However, it should be kept in mind that the values of DE i R?TS only provide information on static contributions to the interaction energy between a given residue and the quantum subsystem. Any effect of the residue on the conformational freedom of the active site will no be spotted by this analysis.
In that regard, it is interesting to analyse the movement of the hydroxyl groups of the sugar moiety once the chain is open. Fig. S4 shows the probability distributions for the dihedral angles C2X GAL -C3X GAL -C4X GAL -O4X GAL , C3X GAL -C4X GAL -C5X GAL -O5X GAL and C4X GAL -C5X GAL -C6X GAL -O6X GAL . They describe the ability to rotate of the hydroxyl groups formed by O4X GAL , O5X GAL and O6X GAL , respectively. The widest distribution corresponds to the hydroxyl group formed by O6X GAL , which rotates almost freely. On the contrary, the groups involving O4X GAL and O5X GAL liberate around to their average values. These two atoms participate in the bonds that close the ring in Galf and Galp, respectively. In both cases the impediment to rotate is mainly caused by the strong H-bonds that these hydroxyl groups maintain with O4 FADH . The O4X GAL -O4 FADH H-bond is present in 95% of the configurations. The O5X GAL -O4 FADH H-bond is present in 96%. While the sugar chain remains open, these groups do not participate in any other hydrogen bonding interaction. The remaining hydroxyl groups of the sugar also form H-bonds. O3X GAL interacts with Asn201 and O2X GAL with the phosphate group.

Stage 3: Formation of the flavin-Galf adduct
This stage also has two steps. The first one is the cyclization of the sugar into the furanose form (step 4 of Fig. 2). It occurs accompanied of the proton transfer from O4X GAL to O4 FADH . We found that this step supports the highest barrier of the whole mechanism: 23.4 + 0.4 kcal/mol. This agrees with the measurements of Sobrado et. al. who determined that this formation of the furanose ring is sensible slower than the opening process TcUGM [28]. The free energy change of the step is 2.9 + 0.2 kcal/mol. In Fig. S5 panel a) we show the evolution of the distances involved in the definition of the reaction coordinate z 4~d8 {d 9 {d 10 , while panel (b) displays the dihedral angles that determine the orientation of the hydroxyl groups at positions 4 and 5. It is observed that, at the beginning of the process, the two hydroxyl groups change their orientation in a concerted way, while O4X GAL and C1X GAL approximate to each other. These movements destroy the H-bond interaction between O5X GAL and O4 FADH and, initially, also drive the H atom to be transferred away from O4 FADH . However, once the O4X GAL -C1X GAL distance gets short enough, a fine tuning in the orientation of the hydroxyl group at position 4 is observed. This reorientation takes its H atom closer to O4 FADH . Initially, while all these rearrangements take place, the H-O4X GAL bond hardly stretches. Only when the O4X GAL -C1X GAL distance gets smaller than 2.3 Å , the H-O4X GAL bond starts to weaken. At the transition state, which appears rather late, the H-O4X GAL and H-O4 FADH distances are almost the same. Finally, when the products configuration is reached and the sugar is in the furanose form, O5X GAL and O6X GAL present no H-bond interactions. In contrast, O3X GAL and O2X GAL keep their interactions with Asn201 and the phosphate group, respectively. The stabilization pattern presented in Table 2 shows that none of the active site residues lowers the energy of TS4 with respect to the reactants in a significant amount. However we note that entropy could play an important role for this step. In general, the conformational freedom of a sugar moiety is larger when the chain is open than when it is closed. Accordingly, it is expected that the cyclization process occurs accompanied by a continuous loss of entropy in the sugar molecule. This would make the reaction slower and the equilibrium position more favourable to the open form. In water solution this effect is compensated by an increase in the entropy of the solvent, but the situation is different within an enzyme. The results of the MD simulation for the open form, discussed at the end of the previous section, suggest that for the reaction under analysis this deleterious entropic effect is ameliorated by the interactions between the carbonyl oxygen of the cofactor and the -OH groups at positions 4 and 5 of the sugar moiety. These interactions attenuate the mobility of the open form and therefore reduce its entropy. Consequently, the entropy change in going from reactants to TS4 or to products is not so adverse. We note, however, that our MD runs are not long enough to allow for an accurate estimation of entropic effects. Thus, the hypothesis put forth in this paragraph needs to be evaluated by additional simulations, especially tailored to that end.
A direct comparison with the results of Huang et. al. cannot be done for this step because the authors split the process in three parts. First, a rotation around the C4X GAL -C5X GAL bond to take O5X GAL away from O4 FADH ; second a rotation around the C4X GAL -C3X GAL bond to place O4X GAL close to C1X GAL ; third, the attack of O4X GAL onto C1X GAL . We found that the computed free energy curve along reaction coordinate z 4 presents no intermediate minimum, indicating that the ring closure takes place in a single concerted step. In other words, the energy minima corresponding to the intermediates detected by Huang et. al. do not appear in our free energy computations. This difference could be attributed to dynamical effects which, as explained above, are expected to be large for this cyclization process but are not considered in electronic structure computations.
In the second step of this stage (step 5 of Fig. 2) the H atom is transferred back from O4 FADH to N5 FADH . The calculated DG { r and DG 0 r for this process are 19.8 + 0.4 and 1.0 + 0.3 kcal/mol, respectively. These results are similar to those of Huang et. al. It should also be noted that the barrier and free energy change of this step are pretty similar to those of the reverse of step 2. The result is not surprising since the only difference between them is whether the sugar moiety is in the furanose or pyranose form. The transmission coefficient for H-tunneling for this step is more than fifteen times smaller than that of step 2. This is due to the fact that this step is slightly endothermic while step 2 is exothermic. Nevertheless, in both cases, the probability for H-tunneling is negligibly.
In Fig. S6a we present the evolution of the distances involved in the definition of reaction coordinate z 5~d4 {d 3 , while in panel (b) we show the evolution of the C1X GAL -N5 FADH distance and the Cremer-Pople angle h of the furanose ring. The curves for the H-O4 FADH and H-N5 FADH distances evolve according to the expectations for a direct proton transfer. However, the curve of the C1X GAL -N5 FADH distance shows an unexpected increase after TS5. The final value, ,1.85 Å , is significantly larger than a typical C-N single bond. Our results thus indicate that the adduct between Galf and FADH becomes rather weak when N5 FADH adopts the sp3 hybridization. The enlargement of the C1X GAL -N5 FADH distance was also described by Huang et. al. However, it that case, the final value was somewhat smaller than in our calculations (1.70 Å ). In order to check the final distance between C1X GAL and N5 FADH we re-simulated the transference using longer simulation times for each window, as well as employing larger QM subsystems. However, we consistently obtained the same result. Moreover, the H transference was simulated applying a restriction on the C1X GAL -N5 FADH distance, so that it was forced to get values smaller than 1.65 Å . These calculations provided higher DG { r and DG 0 r than those obtained without the restriction. Besides, when an unrestricted MD was performed on the products of the restricted transfer, the system spontaneously relaxed to a stable conformation with a C1X GAL -N5 FADH distance of ,1.85 Å . Fig. S6b indicates that the increase in the C1X GAL -N5 FADH distance is accompanied with an increase in the Cremer-Pople angle h. This takes the configuration of the sugar ring from 2 T 3 , for reactants, to E 3 for products. It has to be noted that both, the enlargement of the C1X GAL -N5 FADH distance and the change in the conformation of the furanose ring, are required to avoid the steric clash between the cofactor and the substrate. The values of DE i R?TS presented in Table 2 show that His62 and Arg423 destabilize TS5, but the effect is more than compensated by the stabilization produced by Arg176 and Asn201. Site directed mutagenesis experiments determined that the mutation of Arg176 by Ala produces an impressive reduction in k cat [29]. The involvement of this residue in the catalysis of step 5 could be one of the reasons of this finding.

Stage 4: Formation of UDP-Galf
At first sight, this stage could be considered as the reverse of stage 1, except for the fact that the sugar is now in the funarose form. However, as stated in the previous section, the flavin-Galf bond is already very weak when the transference of the proton from O4 FADH to N5 FADH is completed. Because of this, the barrier for this step is pretty low, 5.8 + 0.2 kcal/mol, and DG 0 r is quite negative, 214.8 + 0.1 kcal/mol. In Fig. S7 we show the evolution of the distances involved in the definition of reaction coordinate z 6~d1 {d 2 . At the transition state, the N5 FADH -C1X GAL and C1X GAL -O3B UDP distances are both * 2.18 Å . This corresponds to a bond order of ,0.32 for the two bonds. As in the case of the flavin-Galp adduct formation, this is consistent with a dissociative S N 2 mechanism. The evolution of the Cremer-Pople angle h of the furanose ring is also shown in Fig. S7. It changes from 255 0 to 240 0 , indicating that the ring conformation returns from E 3 to 2 T 3 . Once the products are formed, the O5X GAL hydroxyl group establishes a new H-bond with the phosphate group. Table 2 shows that the stabilization pattern of this step is similar to that of step 1, with Arg327 being the most stabilizing residue of the TS.

Conclusions
We have presented a detailed description of the energetic and structural changes that take place during the entire mechanism of the reaction catalysed by TcUGM. The results confirm and explain several previous experimental findings, and they also provide new insights on the dynamics of the active site along the reaction.
In agreement with experiments, our results confirm that the first stage of the reaction (formation of the Galp-flavin adduct) proceeds on a single step via a S N 2 dissociative mechanism [28,34,46]. Moreover, the computations indicate that Arg327 is the main responsible for the selective stabilization of the TS of this step. This could explain why the substitution of this Arg by Ala reduces the k cat of TcUGM by 69% [29]. The carbonyl oxygen of Gly61 also plays a role in that regard. However, the stabilization energy of this residue could not be quantified.
The second stage of the reaction (formation of the iminium ion) occurs in two steps, as predicted by the electronic structure calculations of Huang et. al. First, a proton passes from N5 FADH to O4 FADH . This is followed by the transfer of the proton from O4 FADH to O5X GAL that triggers the opening of the ring. The energy decomposition analysis indicated that residue His62 is the main responsible for the catalysis of the first transference. This result could explain why mutations introduced on that position produced a high detrimental effect on the activity of TcUGM. The transfer of the proton from O4 FADH to the cyclic oxygen, on the other hand, presents a relatively low barrier and none of the active site residues is particularly relevant to stabilize or destabilize its TS with respect to the reactants.
At the end of stage 2, galactose is in the open-chain form. We analysed, the interactions and dynamics of the hydroxyl groups of the sugar in such situation. We found that the group at position 6 moves freely, without interacting with any active site residue, the one at position 3 forms a H-bond with Asn201 while the one at position 2 forms a H-bond with the phosphate group of UDP. More importantly, the hydroxyl groups at positions 4 and 5 strongly interact with O4 FADH via H-bonds. These interactions significantly reduce the conformational freedom of the sugar moiety.
During the closing of the sugar ring to form Galf, the dihedral angles around the C2X GAL -C3X GAL and C3X GAL -C4X GAL bonds rotate in a concerted way. These rotations take O5X GAL away from O4 FADH and locate O4X GAL close to C1X GAL . It was found that this step presents the highest free energy barrier of the whole mechanism, in agreement with the proposal of Sobrado et. al. [28]. The energy decomposition analysis showed that none of the active site residues is particularly important to reduce the energy of the TS with respect to the reactants, but the MD results discussed in the previous paragraph suggested that entropy could play an important role. In general, the cyclization of a sugar chain occurs with a reduction of its entropy, a fact that hampers the reaction. In this case, by reducing the conformational freedom of the open-chain form, the active site of TcUGM could make the entropy change and the activation entropy of this step less adverse. Unfortunately, the characteristics of our simulations do not allow to quantify this effect. We note, however, that since this step has the largest free energy barrier, any small reduction on that barrier can be significant.
Once Galf is formed, the next step involves the transference of the proton bound to O4 FADH towards N5 FADH . We observed that something unexpected occurs during this process. Once the system has passed over the TS, the furanose ring changes its conformation from 2 T 3 to E 3 while the distance between C1X GAL and N5 FADH increases to get a final value of ,1.85 Å . The visual inspection of the structures reveals that these modifications are required to avoid the steric clash between the substrate and the cofactor. Huang et. al., who used a different level of theory, different quantum subsystem and different model for the active site, also found a rather long C1X GAL -N5 FADH distance at the end of this transference. Residues Arg176 and Asn201 make the main contributions to the lowering of the barrier. This role of Arg176 is in line with recent experiments which found that the mutation of this residue by Ala reduce the k cat of TcUGM [29]. During the last step of the reaction, the sugar in the furanose form re-binds to UDP as it detaches from the cofactor. Since the C1X GAL -N5 FADH bond is already rather weak at the end of the previous step, this last transformation presents a small barrier and a very negative energy change.
Tyr395 and Tyr429 also play an important role in the reaction. Both residues bear strong H-bond interactions with the phosphate group of the cofactor. These bonds are stable throughout the whole catalysed mechanism. Since these interactions are always present, they do not modify the energy of the barriers found along the reaction. Instead, they facilitate the process by keeping the phosphate group at a relatively fixed position, close to the sugar moiety. Thus, UDP is ready to re-bind to the sugar once it adopts the furanose form. Not surprisingly, experiments determined that the substitution of any of these tyrosines by phenylalanine reduced the k cat of TcUGM [29].
Summarizing, the QM/MM molecular dynamics computations presented in this article determined that residues His62, Arg176, Asn201 and Arg327 contribute to the catalytic activity of TcUGM by reducing the barriers of different steps of the mechanism. Tyr385 and Tyr429, on the other hand, play a role by keeping UDP always close to the sugar moiety. Also, the results highlight the participation of the carbonylic oxygen at position 4 of the cofactor. As predicted by Huang et. al. this atom provides an alternative route for the transference of the proton between N5 FADH and the cyclic oxygen of the substrate. Without this route the barrier for the transference would be prohibitively high. Besides this oxygen restricts the mobility of the open-chain form of the sugar facilitating the ciclyzation process. We hope that the insights obtained from this computational study can contribute to the design of efficient inhibitors of TcUGM.

Initial settings
The crystallographic structure of reduced TcUGM with UDP was taken from the Protein Data Bank, entry 4DSH. To determine the coordinates of Galp within UGM we superimposed the UDP-Galp molecule, taken from the crystal structure of Asparragilus fumigatus UGM (PDB code 3UKF), with the crystallographic UDP of TcUGM. The resultant coordinates of UDP-Galp, together with those of TcUGM, were used as the starting geometry of TcUGM in its holo form. In the initial configuration the nucleophilic group and the leaving group laid on opposite sides of the sugar ring. The distance between C1X GAL and N5 FADH was 3.78 Å . The angle between N5 FADH , C1X GAL and the oxygen atom of UDP, O3B UDP , was 144.2u. The flavin cofactor was set in the reduced deprotonated state since it was recently shown that this form augments the nucleophilic character of N5 FADH [40]. Besides, since experiments indicate that the pKa of N1 FADH is , 6.7 while that of N5 FADH is w 20, the proton of the reduced flavin was located on N5 FADH [49]. The protonation state of the enzyme residues was assigned according with the standard rules except for His62, since recent experiments showed that this residue is protonated when the cofactor is in the reduced state [38]. The resulting file was fed into the Leap module of AMBER and the system was solvated in a 10.0 Å truncated octahedral cell of TIP3P explicit water molecules [50], including the crystallographic water molecules.
The QM/MM molecular dynamics and free energy simulations were performed with the AMBER12 package [51], using periodic boundary conditions with a cutoff distance of 10.0 Å and a time step of 1.0 fs. The potential energy of the classical region was computed with the Amber99SB force field [52] while the selfconsistent charge Density Functional Tight Binding method (scc-DFTB) [51] was employed for the QM subsystem. The DFTB method has proved to be appropriate to describe the energetics of many chemical [53] and biochemical reactions [54][55][56]. More recently, it was shown to provide the best semiempirical description for six-membered carbohydrate rings deformation [57,58]. The QM subsystem was formed with the flavin cofactor, the substrate, Gly61, His62, Val63, as well as the lateral chains of Arg176, Arg327 and Arg423. This adds up to 232 atoms with a net charge of -1.
The initial structure was first minimized at constant volume and then heated at NVT conditions from 0 K to 310 K by a simulated annealing technique. A weak harmonic restraint on the C a atoms was implemented during this period. This was followed by 200 ps of equilibration at NPT conditions at 310 K and 1 bar. No restrains were applied in this case. The Pauling Bond Orders, n x , were determined when galactose either attaches or detaches from the flavin cofactor (process a?b and f?g of Fig. 2). In both cases, the bonds involved are C-O and C-N. The equation used to calculate the orders was, n x~n0 e rx{r 0 ð Þ=0:6 : ð1Þ Here n 0 denotes the bond order of the fully formed bond while r 0 is the equilibrium distance, which was considered equal to 1.5 Å for the two bonds involved in these reactions. The value of r x was computed as the average distance among the structures sampled in the umbrella simulations at the transition state. The presence of Hbonds was monitored considering that a H-bond exists if the distance between the donor and the acceptor is v 3.15 Å and the donor-H-acceptor angle is w 145u.
When relevant (steps 2 and 5), the probability of H-tunneling was estimated employing the expression for the microcanonical transmission coefficient given at equation 14a of reference [59]. This expression corresponds to tunneling through a one-dimensional barrier whose shape, height and exothermicity are determined by three adjustable parameters. In our estimations these parameters were obtained by fitting the free energy curves of the corresponding proton transfer steps. Since the coordinate in these curves is not the proton coordinate but a difference between the distances of the bond being broken and formed, the effective mass of the particle being transferred was set as m H =4, where m H is the mass of the proton [60]. The energies employed in these estimations were 1700 cm {1 for step 2 and 1800 cm {1 for step 5. These are approximately the zero point energies of the N-H and O-H bonds.

Umbrella sampling calculations
The umbrella sampling technique was employed to analyse all the steps involved in the conversion between Galp and Galf within TcUGM (see Fig. 2). Free energy profiles were computed along different reaction coordinates, conveniently defined for each transformation. Harmonic restraints were applied in order to force the system to wander around the selected values of the reaction coordinate. A restraining force of 350.0 kcal/molÅ 2 was employed in all cases and the reaction coordinate was sampled considering windows of 0.08 Å wide. Within each window, an equilibration phase of 65 ps was followed by a production phase of 0.2 ns. The actual values of the reaction coordinate were recorded every 2 fs. Snapshots of the structures were downloaded every 3 ps. The last 30000 data of each window were used to compute the unbiased probability by means of the weighted histogram analysis method (WHAM) [61]. When following each reaction coordinate, the last structure of a given window was used as the starting point for the next one. Simulations of 0.5 ns without any restraint were performed for reactants, products and for each intermediate species in order to check their stability. To check the convergence of the free energy computations, several tests were performed. First, for each step, we compared the free energy profiles obtained using the first half of the data (i.e. the first 15000 values selected at each specific value of the reaction coordinate) with the second half. Besides, every reaction coordinate was sampled both, forward and backwards. Finally, computations for steps 3 and 4 (the opening and closure of the sugar ring) were repeated three times using different initial configurations. Summing all the steps of the reaction, with their corresponding convergence tests, the total length of the QMMM-MD simulations was 187.5 ns. Below we provide the numerical details of the umbrella sampling calculations for each stage.
Stage 1: Formation of the flavin-Galp adduct. This stage consists of just one concerted step in which the bond between Galp and UDP breaks while Galp joins FADH 2 (step 1 in Fig. 2). Accordingly, the reaction coordinate was defined as z 1~d2 {d 1 , where d 1 is the distance between C1X GAL and N5 FADH , while d 2 is the distance between C1X GAL and O3B UDP . This coordinate was sampled from 22.03 to 1.89 Å .
Stage 2: Formation of the iminium ion. This stage was proposed to occur in two consecutive steps. The first one is the tautomerization of FADH -via the transfer of the H atom bonded to N5 FADH towards O4 FADH (step 2 in Fig. 2). The reaction coordinate for this tautomerization was defined as z 2~d3 {d 4 , where d 3 denotes the N5 FADH -H distance while d 4 is the O4 FADH -H distance. Coordinate z 2 was sampled from 21.63 to 1.73 Å .
During the following step, the opening of the ring is initiated by the transfer of the proton linked to O4 FADH towards the oxygen atom of the Galp ring, O5X GAL , and proceeds with the breakage of the bond between C1X GAL and O5X GAL (step 3 in Fig. 2). We found that a correct description of this process required a reaction coordinate defined as z 3~d5 zd 6 {d 7 , where d 5 is the distance between O4 FADH and the proton being transferred, d 6 is the distance between C1X GAL and O5X GAL and d 7 is the distance between the proton and O5X GAL . Coordinate z 3 was sampled from 0.82 to 4.18 Å .
We also analysed the possibility of a direct proton transfer from N5 FADH to O5X GAL . In order to do that we defined a reaction coordinate z' 3~d ' 5 zd' 6 {d' 7 , where d' 5 is the N5 FADH -H distance, d' 6 the one between O5X GAL and C1X GAL and d' 7 the O5X GAL -H distance. This coordinate was sampled from 21.85 to 0.71 Å . In agreement with previous results of Huang et. al. [41] we found that this direct proton transfer is very unlikely since it has a barrier significantly higher than the alternative path.
Stage 3: Formation of the flavin-Galf adduct. This stage also occurs in two steps. First, the hydrogen attached to O4X GAL is transferred to O4 FADH while a bond between O4X GAL and C1X GAL is formed (step 4 in Fig. 2). The reaction coordinate for this step was defined as z 4~d8 {d 9 {d 10 , where d 8 is the distance between the H atom being transferred and O4X GAL , d 9 the distance between the H atom and O4 FADH and d 10 is the distance between O4X GAL and the C1X GAL . Coordinate z 4 was sampled from 24.85 to 21.01 Å .
The following step consists of a proton transfer from O4 FADH to N5 FADH (step 5 in Fig. 2). This can be seen as the reverse of step 2, except for the fact that galactose is now in the furanose form. Therefore, the reaction coordinate was defined as the reverse of step 2 (z 5~d4 {d 3 ) and it was scanned from 21.63 to 1.65 Å .
Stage 4: Formation of UDP-Galf. This last step corresponds to the breakage of the bond between FADH 2 and Galf along with the formation of a bond between Galf and UDP (step 6 of Fig. 2).
Since this process is analogous to step 1 but occurs in reverse sense we defined z 6~d1 {d 2~{ z 1 and scanned it from 21.98 to 1.38 Å .

Energy decomposition
An energy decomposition analysis was performed to evaluate how the active site residues stabilize or destabilize the transition states of the successive steps with respect to their correspondent reactants. Different variations of this idea have been implemented to study enzymatic reactions [62][63][64][65][66][67][68][69][70][71][72]. In this case we followed the approach recently employed to compare the catalytic mechanisms of T. cruzi transialidase and T. rangeli sialidase [54]. Since the approach has been discussed in detail elsewhere we only present here the most relevant equations.
In the QM/MM study of an enzymatic reaction the influence of the classical environment on the activation energy of a given step, DDE env R?TS , can be evaluated as, Here DE QM=MM R?TS is the activation energy within the enzyme (computed with the QM/MM approach) while DE QM R?TS is the activation energy of the isolated quantum subsystem (computed at the same QM level). The terms appearing in the summation, DE i R?TS , measure the influence of each individual residue on the reaction barrier. They are strictly given by, where DY TS T and DY R T are the wave functions of the quantum subsystem at the transition state and reactants configurations, respectively, while V i is the non-bonded interaction energy of classical residue i with the quantum subsystem. The evaluation of SY X DV i DY X T, with X = TS or R, is not trivial since the AMBER code does not compute these values. Instead it provides the energy of the whole system, which accounts for the quantum hamiltonian, H QM , plus the sum of all the non-bonded interactions between the QM subsystem and the classical environment. Thus we estimated each SY X DV i DY X T as, Here the first term on the right side gives the actual energy of the system at the given configuration. The second one is a fictitious energy calculated with the same wave function by setting the classical environment at exactly the same configuration except for the i-th residue which is transformed into Gly. Average values of SY X DV i DY X T, with X = TS or R, were computed employing 100 snapshots taken from the umbrella sampling calculations with the reaction coordinate set at the TS or reactants configurations, respectively. For these calculations we defined the QM subsystem as the substrate plus the cofactor, while the active site residues under analysis were His62, Val63, Arg176, Asn201, Tyr317, Arg327, Tyr395, Arg423, Tyr429 and Asn433.
The DE i R?TS computed in this way measures the difference between the actual barrier to reaction and the barrier that would be observed if the interaction between the side chain of residue i and the QM subsystem were turned off. Because of this, neither can it provide information about the effect of the backbone atoms or the effect of Gly residues. Moreover, since no dynamics is run when the i-th residue is replaced by Gly, DE i R?TS does not take into account dynamic effects arising from changes in the conformational freedom of the enzyme upon replacement. Finally we note that positive/negative values of DE i R?TS provide a strong indication about a deleterious/beneficial effect of residue i for the reaction step under consideration. However, they cannot be used to quantitatively estimate changes in k cat produced by the mutation of the i-th residue by Gly because such changes depend on variations in the activation free energy, DG R?TS . Text S1 PDB file for UDP-Galp bound to UGM (Michaelis complex). This species is labelled as a in Fig. 2

. (PDB)
Text S2 PDB file for the flavin-Galp adduct in UGM. This species if labelled as b in Fig. 2

. (PDB)
Text S3 PDB file for the third species of the mechanism proposed for the reaction catalysed by UGM. This species is labelled as c in Fig. 2

. (PDB)
Text S4 PDB file for the iminium ion in UGM. This species is labelled as d in Fig. 2

. (PDB)
Text S5 PDB file for the fifth species of the mechanism proposed for the reaction catalysed by UGM. This species is labelled as e in Fig. 2

. (PDB)
Text S6 PDB file for the flavin-Galf adduct in UGM. This species is labelled as f in Fig. 2

. (PDB)
Text S7 PDB file for UDP-Galf bound to UGM. This species is labelled as g in Fig. 2. (PDB)