article: The dissociation mechanism of the thioredoxin (Trx) mixed disulfide complexes is unknown and has been debated for more than twenty years. Specifically, opposing arguments for the activation of the nucleophilic cysteine as a thiolate during the dissociation of the complex have been put forward. As a key model, the complex between Trx and its endogenous substrate, arsenate reductase (ArsC), was used. In this structure, a Cys29Trx-Cys89ArsC intermediate disulfide is formed by the nucleophilic attack of Cys29Trx on the exposed Cys82ArsC-Cys89ArsC in oxidized ArsC. With theoretical reactivity analysis, molecular dynamics simulations, and biochemical complex formation experiments with Cys-mutants, Trx mixed disulfide dissociation was studied. We observed that the conformational changes around the intermediate disulfide bring Cys32Trx in contact with Cys29Trx. Cys32Trx is activated for its nucleophilic attack by hydrogen bonds, and Cys32Trx is found to be more reactive than Cys82ArsC. Additionally, Cys32Trx directs its nucleophilic attack on the more susceptible Cys29Trx and not on Cys89ArsC. This multidisciplinary approach provides fresh insights into a universal thiol/disulfide exchange reaction mechanism that results in reduced substrate and oxidized Trx.
  Thioredoxin (Trx) is a powerful and universal protein disulfide bond oxido-reductase with a very low redox potential [1]–[3]. All thioredoxins have a similar three-dimensional fold comprising a central core of four β-strands surrounded by three α-helices [4]. All feature a conserved active-site loop containing two redox-active cysteine residues in the sequence Trp-Cys-Gly-Pro-Cys [5] numbered as Trp28 to Cys32 in both Bacillus subtilis (Bs_Trx) and S. aureus (Sa_Trx) Trx. This numbering will be used throughout. The active site is at the end of an alpha helix, the α1-helix, extending from Lys33 to Glu45. The pKa of the N-terminal active-site cysteine [6]–[9] is significantly lower than the pKa of a cysteine in the absence of a structured protein environment. Under physiological conditions, this low pKa enables thioredoxin to act as a nucleophile, attacking disulfides in proteins [1], [8], [10]. Well-documented endogenous substrates of Trx are the arsenate reductases (ArsC) from the thioredoxin-coupled ArsC family, including B. subtilis (Bs_ArsC) and S. aureus (Sa_ArsC) arsenate reductase [11]–[14]. ArsC catalyzes the reduction of arsenate [As (V) ] to arsenite [As (III) ] and is a key enzyme involved in arsenic detoxification [15]. For the reduction of arsenate, Bs_ArsC and Sa_ArsC combine a phosphatase-like nucleophilic displacement reaction in the active site with a distinct intra-molecular disulfide-bond cascade [13], [16]–[18]. Three redox active cysteines are involved (Cys10, Cys82 and Cys89) [19]. After a single catalytic arsenate reduction event, oxidized ArsC exposes a disulfide between Cys82 and Cys89 on a looped-out redox helix [17], [18]. Thioredoxin converts oxidized ArsC back to its initial reduced state [1]. Cys29Trx nucleophilically attacks Cys89ArsC of the Cys82ArsC-Cys89ArsC disulfide, leading to the reduction of Cys82ArsC and the formation of the Trx-ArsC mixed disulfide intermediate complex between Cys29Trx and Cys89ArsC (Figure 1 and Figure 2A) [1], [20]. In this complex, Cys32Trx performs a nucleophilic attack on Cys29Trx of the Cys29Trx-Cys89ArsC disulfide (Figure 2B). Accordingly, the Trx-ArsC complex dissociates, releasing reduced ArsC and oxidized Trx (Figure 2C). The mechanism behind the experimentally observed regioselectivity (attack on Cys89ArsC instead of on Cys82ArsC during complex formation and on Cys29Trx instead of on Cys89ArsC during complex dissociation) is not known. Experimentally, it has been observed that Cys32Trx had performed a nucleophilic attack on Cys29Trx, but why and how this can occur was not known. A priori, there is another possibility, i. e. a nucleophilic attack by Cys82Trx leading to an unproductive cycle. The regioselectivity analysis gives explanatory insights in the observed mechanism, which could not be understood on the basis of prior information including the experimental structure of the Trx-ArsC complex. This regioselectivity is critical for the function of the involved proteins, and is a reoccurring question within thiol-oxidoreductases. Also, it is not known how Cys32Trx is activated as a thiolate in the Trx-ArsC mixed disulfide complex, resulting in a productive cycle. Especially the deprotonation mechanism of the Cys32Trx thiol is extensively debated [21]–[26]. Although structures of Trx-peptide mixed disulfide complexes have been described, i. e. Trx-Ref-1 [27] and Trx-NFκβ [28], today, the only structure of a Trx-protein mixed disulfide complex is that recently obtained between Trx and Bs_ ArsC [20]. This presents an opportunity to study Trx mixed disulfide complex formation and dissociation. Here, molecular dynamics (MD), pKa calculations, reactivity analysis with conceptual density functional theory (DFT) and biochemical complex formation experiments with Cys-mutants are applied to probe different aspects of the reaction mechanisms of interest: differential reactivity of the cysteines, regioselectivity of the nucleophilic attacks, structure of the Trx-ArsC activated complex, and the role of specific residues such as Asp23Trx. These various elements are combined to propose a new and consistent view of the mechanisms, firmly grounded in structural information. A computational approach is well suited here in particular because it allows the direct investigation of the wild-type sequence with Cys32Trx. Also, computation offers the possibility to probe directly the structure, dynamics and energetics of the reacting species in atomic details. That is why computational approaches have made profound contributions to the study of enzymatic reaction mechanisms in recent years, typically combining MD simulations and quantum-mechanical calculations [29]. Nowadays, molecular dynamics (MD) simulations have become an important tool to explore the structural dynamics of proteins and their implications for function [30]. MD simulations have been applied to a number of thiol disulfide oxidoreductases in the Trx superfamilly to determine the structural factors that control the pKa of the thiols [31]–[33]. These studies and others [34]–[37] showed that hydrogen-bonding to cysteine sulphur atoms is crucial to stabilize the thiolate form and influence the reactivity. We have exploited MD simulations in explicit solvent to uncover conformational changes at the Trx-ArsC interface, which provide new insights in the structural factors underpinning the chemistry of the Trx-ArsC complex dissociation. The NMR structure of the mixed disulfide was used as a starting point for MD simulations, pKa and reactivity calculations. The conformational changes at the Trx-ArsC interface during the MD simulations revealed how Cys32Trx can move close enough to Cys29Trx to bring the two sulphur atoms in contact, primed for reaction. This local conformational change is a functionally new conformation, not present in the NMR structure and clearly different from the conformers in the NMR structure of the Trx-ArsC mixed disulfide. In this conformation, the sulphur of Cys32Trx can form two hydrogen-bonds stabilizing its thiolate form. In the following we refer to this conformation as the “activated complex”. Knowledge of the protonation state of the residues involved in the reaction mechanism is highly important, but the respective pKa values are not always easy to determine experimentally [38]. Therefore, much effort is devoted to develop methods for theoretical pKa computation [38], [39]. Several approaches exist, including classical electrostatics in the context of the Poisson-Boltzmann formalism [36], [40] and more empirically trained models [34]. In the present work, the linear correlation between the natural population analysis (NPA) charge [41] calculated on the sulphur atom of thiols and their experimentally determined pKa [37], [42] is used to quantitatively calculate the cysteine pKa' s in thioredoxin, arsenate reductase, and their mixed disulfide complex, including the activated complex. The more negative the NPA charge on the sulphur atom, the higher the tendency to bind a proton, and the higher the thiol pKa is. Calculation of the NPA charge with DFT means that the electronic environment of the sulphur atom is treated explicitly. To account for solvent effects, a continuum solvent model is applied. As such, the NPA-pKa correlation method is better rooted in first principles than the more empirical approaches [34]. The NPA-pKa method was successfully used in the study of the activation of the Cys82 and the Cys89 thiolates in ArsC [42] and in the study of the origin of the pKa perturbation of N-terminal cysteines in α- and 310-helices [37]. The NPA-pKa method remains tractable for proteins, because models including the relevant protein environment for the considered cysteine were designed. Another aspect of reaction mechanisms related to structure is their regioselectivity. The regioselectivity of the disulfide exchange reactions during Trx-ArsC mixed disulfide complex formation and dissociation can be addressed by the soft acids an bases principle (HSAB) [43], [44]. This principle is defined in the conceptual DFT [45] context and states that hard acids prefer to react with hard bases whereas soft acids prefer to interact with soft bases. Disulfide exchange reactions are soft-soft interactions; as such the softness is used as a reactivity descriptor. The smaller the difference in the local softness (s) between the sulphur atoms the more preferred the reactivity between the attacking nucleophilic cysteine (Cys29Trx, Cys32Trx and Cys82ArsC) and the accepting electrophilic disulfide (Cys82ArsC-Cys89ArsC and Cys29Trx-Cys89ArsC). This might result in several possible reaction paths. Similar reaction paths do not cross according to Klopman' s rule [46], and as such the relative energies of the reagents at the ground state (at the beginning of the reaction) correlate with the relative energies at the transition state. An advantage of this conceptual framework is that no activation energies need to be calculated to predict the reaction path. Only reactivities on structures at the beginning of the reaction were calculated. When investigating complex formation, the reactivity analysis was performed on free Trx and ArsC. When addressing complex dissociation, conformers of the Trx-ArsC complex were used. In keeping with the conceptual framework, snapshots of MD simulations representing activated complexes along the reaction path were not used in reactivity analysis. The local softness (s) results from the multiplication of the global softness (S) with the Fukui function (f) [47]. S is a global property and correlates with the system polarizability [48]. S is the inverse of the chemical hardness, which is the second derivative of the energy with respect to the numbers of electrons, evaluated at a fixed molecular geometry. This derivative can be approximated as the difference between the vertical ionization energy (IE) and the electron affinity (EA). Differentiation of the energy with respect to the external potential υ (r) (i. e. the potential felt by the electrons due to the nuclei) introduces a local character into the global reactivity descriptors resulting in the Fukui function f. f is a local descriptor and indicates regions where the molecule preferentially reacts (regioselectivity): i. e. which nucleophilic sulphur atom (f−) will attack and which sulphur (f+) of the disulfide will receive the electrons. The − and + sign indicate the reactivity towards nucleophilic and electrophlic attacks respectively. f− is related to the electron density of the highest occupied molecular orbital and f+ to the density of the lowest unoccupied molecular orbital when electrons are received [49]. The global and local softness and the Fukui function are well established and correlate with reactivity data [45]. In sum, the local softness descriptors are combined with pKa calculations, MD simulations on the Trx-ArsC complex and with biochemical experiments to give fresh insights in the mechanism behind mixed disulfide complex dissociation.
  The NPA charge has been shown to be an effective descriptor for the pKa. In a series of small substituted thiolate molecules, a linear relationship was obtained between the NPA-charge of the sulphur atom and the experimental pKa value (Figure 3A) [37]. The more negative the NPA charge on the sulphur atom, the higher the tendency to bind a proton and, the more basic (i. e. higher pKa) the compound is. This linear relationship can be used as a calibration curve to quantify the pKa perturbing effect. Remarkably, this initial NPA-pKa correlation [37] obtained for a series of small molecules was found to be directly transferable to the pKa of the cysteine residues in Trx and ArsC (Figure 3A, Table S1 and Figure S4 in Supplemental Data). The protein environments are represented by models which include all residues interacting with the considered cysteines. The models of all Trxs include the WCPGC active site and the adjacent α1-helix. ArsC is represented by the redox helix region (Cys82-Cys89), Arg16 and Thr11. The calculated cysteine pKa values are in agreement with the experimentally obtained pKa' s with a maximum deviation of 0. 5 pKa units (Figure 3A, red diamonds). This confirms that the model systems used for the calculations are appropriate to represent the cysteines in their protein environment. So far, the only outliner is the non-nucleophilic cysteine of E. coli Trx1 (Ec_Trx1; data point ‘g’ in Figure 3A), which has a severely underestimated pKa. This is likely due to the presence of the nucleophilic cysteine thiolate at 4. 7 Å (sulphur-sulphur distance). In the Ec_Trx1 model system used for pKa analysis of this non-nucleophilic cysteine, the negative charge of the neighbouring thiolates, is equally distributed among the Sγ atoms of both residues. As such, no negative charge can build up on the non-nucleophilic cysteine, leading to the underestimation of its pKa. The theoretical pKa estimation via the NPA charge is not suited for two neighbouring negatively charged thiolates, but that is not a situation which dominates or influences particularly the object of our study. The linear correlation between the NPA charge on the sulphur atom and the experimental pKa was adjusted specifically for the data points a–f of the Trx and ArsC systems (Figure 3B, Table S1 in Supplemental Data). This correlation is particularly suited for the active site thiols of the enzymes Trx and ArsC and is used throughout this work to quantitatively calculate the pKa' s of the cysteine residues involved in Trx-ArsC complex formation and dissociation. Details on the models used for pKa and reactivity analysis (vide infra) of the cysteine residues involved in Trx-ArsC complex formation and dissociation are given in Figure 4. With oxidized Bs_ArsC as a substrate, the Cys29Trx nucleophile of reduced Bs_Trx attacks Cys89ArsC of the Cys82ArsC-Cys89ArsC disulfide (Figure 2A). This section addresses the reasons behind the experimentally observed regioselectivity of this nucleophilic attack of Cys29Trx on the Cys82ArsC-Cys89ArsC disulfide. In this initial reaction, Cys29Trx is stabilized as a nucleophilic thiolate, in contrast with Cys32Trx. Based on the NPA-pKa correlation (Figure 3B), the pKa' s of Cys29Trx and Cys32Trx (in respectively the Trx_red_Cys29 and Trx_red_Cys32 models; Figure 4A) in reduced Bs_Trx (2GZY) are calculated to be 5. 5 and 8. 2 respectively (Table 1). The minimal difference in the local softness between the sulphur atoms of electrophile and nucleophile favors the nucleophilic attack of Cys29Trx on Cys89ArsC (Table 2). In the Cys82ArsC-Cys89ArsC disulfide of oxidized ArsC (1Z2D, Bs_ox in Figure 4B), the Fukui function value of Cys89ArsC is clearly higher than that of Cys82ArsC (Table 3). As such, Cys89ArsC is intrinsically more susceptible to a nucleophilic attack than Cys82ArsC. This regioselectivity is largely determined by the side chains of the looped-out redox helix in which the Cys82ArsC-Cys89ArsC disulfide is embedded. Indeed, removing the side chains of the redox helix (‘Bs_ox without side chains’ model system) blurs this regioselectivity (Table 3). In reduced Bs_Trx (2GZY, Trx_red_Cys29 and Trx_red_Cys32 in Figure 4A), Cys29Trx and Cys32Trx have comparable Fukui function values, corresponding to the same intrinsic reactivity. However, here, Cys29Trx has a lower softness than Cys32Trx, explaining (in addition to the pKa values) its higher reactivity towards the less soft Cys82ArsC-Cys89ArsC disulfide (Table 3). So, pKa calculations combined with reactivity analysis explain the regioselectivity of the first reaction in the catalytic cycle leading to reduction of ArsC by Trx. The NMR structure of the Bs_Trx-ArsC complex (PDB code 2IPA) [1] and the derived MD simulations (75 ns) have been used to investigate the dynamics of the hydrogen bond interactions formed with the sulphur atom of Cys32Trx and of Cys82ArsC. Indeed, hydrogen bonds to the sulphur can influence its reactivity by stabilizing its thiolate form (see Introduction). The NMR structure of the Trx-ArsC mixed disulfide was obtained after mutation of Cys32 and Cys82 to serines. These serines were converted to the wild-type cysteines by straightforward modeling in the 21 conformers of the NMR structure. In seven out of the 21 conformers, Cys82ArsC forms a hydrogen bond with Arg16ArsC, and in eight of them Cys82ArsC, forms a hydrogen bond with Thr11ArsC. Cys32Trx is hydrogen-bonded to the N-atom of the Cys29Trx backbone (Figure 4D) in eleven conformers. In eight conformers, no hydrogen bonds are formed with Cys82ArsC or with Cys32Trx (Figure 4C). This simple analysis suggested that both protonation states of Cys82ArsC, neutral and thiolate, should be considered, since in the 21 conformers, Cys82ArsC forms either no or two hydrogen bonds. Starting from a conformer of the NMR Trx-ArsC structure devoid of hydrogen bonds to Cys32Trx and Cys82ArsC, two MD simulations were set up, with Cys82ArsC modeled as neutral or as a thiolate. When Cys82ArsC is deprotonated, the following hydrogen bonds are formed during the simulation (Figure 5 and Figure S1 in Supplemental Data): Cys82ArsCSγ—Arg16ArsCNη1 (15% of the time), Cys82ArsCSγ—Arg16ArsCNε (1% of the time), Cys82ArsCSγ—Thr11ArsCN (94% of the time), Cys82ArsCSγ—Thr11ArsCOγ (74% of the time), Cys82ArsCSγ—Arg108ArsCNη2 (48% of the time), Cys82ArsCSγ—Arg108ArsCNη1 (40% of the time), Cys32TrxSγ—Cys29TrxN (60% of the time) and Cys32TrxSγ—Trp28TrxN (23% of the time). Cys32Trx is simultaneously hydrogen bonded to Cys29TrxN and Trp28TrxN 20% of the time. Not surprisingly, hydrogen bonds to Cys82ArsC were formed significantly less frequently in the MD simulation where it is neutral (Cys82ArsCSγ—Arg16ArsCNη1 (19% of the time), Cys82ArsCSγ—Arg16ArsCNε (1% of the time), Cys82ArsCSγ—Thr11ArsCN (14% of the time), Cys82ArsCSγ—Thr11ArsCOγ (1% of the time), Cys82ArsCSγ—Arg108ArsCNη2 (0% of the time), Cys82ArsCSγ—Arg108ArsCNη1 (0% of the time). In addition, Cys32Trx was hydrogen bonded to Cys29TrxN and to Trp28TrxN for respectively 1. 6% and 0. 5% of the time, much less frequently than when Cys82ArsC is treated as thiolate. This is a strong indication that the protonation state of Cys82ArsC influences the behavior of Cys32Trx. Of course, the exact numbers presenting the statistics from the simulations, e. g. regarding how frequently a hydrogen-bond is formed, are subject to the length of a simulation and its particular starting point, and are thus expected to vary without changing the main trends. The formed hydrogen bonds influence the pKa of both Cys32Trx and Cys82ArsC. In the Bs_Trx-ArsC complex, both the Cys32TrxSγ—Cys29TrxN hydrogen bond and the α1-helix of Trx stabilize Cys32Trx as a thiolate. When Cys32Trx is not hydrogen bonded, its calculated pKa is 8. 9 (Trx_ArsC_1_Cys32, Figure 4C) (Table 1). In the presence of the Cys32TrxSγ—Cys29TrxN hydrogen bond (Trx_ArsC_2_Cys32, Figure 4D), the pKa drops to 8. 3. To account for the effect of the α1-helix, we designed the Trx_ArsC_1_trunc and Trx_ArsC_2_trunc (Figure 4E–F), in which only Cys89ArsC of the ArsC part is included. Here, the α1-helix additionally decreases the pKa of Cys32Trx with 0. 4 and 0. 8 pKa units respectively. Cys82ArsC has a pKa of 8. 1 (Table 1) in the absence of hydrogen bond interactions (Trx_ArsC_1_Cys82, Figure 4C). Its pKa drops to 6. 3 when Arg16ArsC comes within hydrogen bonding distance (Trx_ArsC_2_Cys82, Figure 4D). In sum, in the Bs_Trx-ArsC complex, hydrogen bonds formed with the sulphurs of both Cys32Trx and Cys82ArsC are functionally important and stabilize the thiolate form of these cysteines. In the Bs_Trx-ArsC complex, dissociation takes place via the nucleophilic attack of Cys32Trx on Cys29Trx of the mixed Cys29Trx-Cys89ArsC disulfide (Figure 2B). Given that both Cys82ArsC and Cys32Trx can be stabilized as a thiolate, both are considered as possible nucleophiles. The nucleophilic attack of Cys32Trx on Cys29Trx is the favored reaction. Comparing the differences in local softness between Cys32Trx/Cys82ArsC (potentially attacking nucleophiles) and Cys29Trx/Cys89ArsC (potentially attacked electrophiles) shows a minimal local softness difference between Cys32Trx and Cys29Trx (Table 4). The Fukui function values of Cys29Trx-Cys89ArsC obtained when Cys32Trx is deprotonated (Trx_ArsC_2_Cys32, Figure 4D) are higher than those obtained when Cys32Trx is protonated (Trx_ArsC_2_Cys82, Table 4). Thus, the deprotonation of Cys32Trx in the Trx_ArsC_2_Cys32 model increases the susceptibility of the Cys29Trx-Cys89ArsC disulfide to a nucleophilic attack. This is consistent with the need for an activated nucleophile in order to perform the reaction. When Cys32Trx is deprotonated, Cys29Trx is the preferred target for nucleophilic attack (Figure 2, Table 5), since the Fukui function of Cys29Trx in the mixed disulfide is clearly higher than that of Cys89ArsC. The ArsC redox helix is essential for the regioselectivity, since in the absence of this redox helix (Trx_ArsC_2_trunc), the difference in the Fukui function values between Cys29Trx and Cys89ArsC is fading. The Trx α1-helix (Trx_ArsC_2_trunc+helix) is not essential for the regioselectivity since the Fukui function values are not affected by adding this helix to the model system (Table 5). In the Trx_ArsC_2_Cys32 model (Figure 4D), Cys32Trx has only a slightly higher Fukui function value than Cys82ArsC. The higher reactivity of Cys32Trx compared to Cys82ArsC towards the less soft Cys29Trx-Cys89ArsC disulfide is consistent with the lower softness of Cys32Trx compared to Cys82Trx (Table 5). This regioselectivity analysis was repeated with the Trx_ArsC_1 model (Figure 4C), confirming the obtained results (results not shown), indicating that our conclusions are not particularly sensitive on structural models, if they include enough relevant protein environment. In addition, geometric factors are at work. During both MD simulations of the Trx-ArsC complex, Cys82ArsC (and its Ser82ArsC equivalent in the NMR structure with PDB code 2IPA) is more than 15 Å away from Cys29Trx, while a local conformational change brings Cys32Trx up to 3. 7 Å of Cys29Trx in the simulation with ionized Cys82ArsC (Figure 6A). In this simulation, the sulphur atoms of Cys32Trx and Cys29Trx are within 4. 0 Å of each other 50. 7% of the time. We refer to this new conformation where Cys32Trx comes in contact with Cys29Trx (not present in the PDB structure 2IPA) as the activated complex. Indeed, a geometric proximity between sulphurs is a pre-requisite to a reaction between them. Such proximity was not observed between Cys29Trx and Cys82ArsC in the MD simulations. The conformational change bringing Cys32Trx in contact with Cys29Trx is associated with the formation of the key hydrogen bonds (Figure 6 B–C), stabilizing Cys32Trx as a thiolate, primed for a nucleophilic attack on Cys29Trx. The pKa of Cys32Trx was calculated in selected snapshots of the activated complex where both the Cys32TrxSγ—Cys29TrxN and the Cys32TrxSγ—Trp28TrxN hydrogen bonds are present and particularly strong (based on geometric criteria, see Figure 7). These calculations were performed with models similar to the Trx_ArsC_2_Cys32 model (Table 6). The snapshots of interest were observed at various stages during the simulation, including towards its end, well past what could be considered an equilibration phase. The results indicate that the introduction of the Cys32TrxSγ—Trp28TrxN hydrogen bond allows for an extra pKa decrease to 7. 4 compared to the Trx_ArsC_2_Cys32 model, in which only the Cys32TrxSγ—Cys29TrxN hydrogen bond is present. These results indicate that the pKa of Cys32Trx can be instantaneously being low. From a mechanistic point of view, it is enough to deprotonate Cys32Tx for a split second for the subsequent reaction of interest to take place. Indeed, reactivity may not be determined by the average value of the pKa, but may be gated by the occasional low values. Therefore, it is legitimated to select some snapshots of special interest. Although not within hydrogen bonding distance of Cys32Trx (>6 Å), Asp23Trx has been suggested to be the key residue for the activation of Cys32Trx as a nucleophile [26]. Hydrogen bonding via a structurally conserved water molecule located between Asp23Trx and Cys32Trx has been put forward for such activation (Figure 8A) [23], [26]. We revisited the suggested role of Asp23Trx by biochemical experiments and MD simulations. Complex formation of wild type Sa_Trx with Sa_ArsC C10S/C15A/C82A triple mutant (Sa_ArsCTrip) via a 2-nitro-5-mercaptobenzoic acid (TNB) -mixed disulfide with Cys89ArsCtrip (see Experimental procedure section) showed the release of TNB−. This means that Cys29Trx attacks Cys89ArsCtrip, indicating complex formation (Figure S2 in Supplemental data). Evaluation on non-reducing SDS-PAGE showed two bands, corresponding to ArsC and Trx. No Trx-ArsC complex could be detected (Figure 8B), indicating the dissociation of the mixed disulfide. Reaction with a Trx D23A mutant gave a similar result. With Trx D23A, we observed TNB− release, indicating complex formation, and only the two bands corresponding to Trx and ArsC were detected on a non-reducing SDS-PAGE (Figure 8B). So, in the absence of Asp23Trx as putative nucleophile activator, the Trx-ArsC complex was formed but still able to dissociate. We also evaluated the oxidation state of Sa_Trx D23A on reversed phase chromatography. Before the reaction, Sa_Trx D23A was in its reduced form, while after the reaction Sa_Trx D23A was oxidized (Figure 8C). This proves that Sa_Trx D23A reduces its substrate, consistent with Cys32Trx being activated as a nucleophilic thiolate and this independently of Asp23Trx. The assumed role of Asp23Trx was also investigated by MD experiments on Bs_Trx-ArsC. When Cys82ArsC is ionized, the distance of Asp23Oδ1 (or Asp23Oδ2) to Cys32Sγ increased to 10 Å and this increase was maintained during the remainder of the simulation up to 75 ns. Associated with this movement, the Cys32TrxSγ—Cys29TrxN and Cys32TrxSγ—Trp28TrxN hydrogen bonds are introduced, forming the activated complex (Figure 8 D, E). This adds to the view that Cys32Sγ is activated by hydrogen bonds to the protein backbone rather than interactions with Asp23, consistent with the results from the complex formation experiments. Recently it was suggested that when the N-terminal thiolate of Trx attacks its substrate disulfide to form a mixed disulfide complex, the leaving thiol group deprotonates the C-terminal active site thiol of Trx [21]. With ArsC as substrate, Cys82ArsC would deprotonate Cys32Trx when Cys29Trx attacks the Cys82ArsC-Cys89ArsC disulfide. In this hypothesis, mutating Cys82ArsC to Ala would trap the mixed disulfide complex as a stable entity, which we did not observe. In complexation experiments between wild type Sa_Trx, Sa_Trx D23A and Sa_ArsCTrip in which Cys82ArsC is mutated to alanine, Trx-ArsC complex formation was observed via TNB− release. Neither wild type Sa_Trx nor Sa_Trx D23A formed a stable complex with Sa_ArsCTrip as visualized on non-reducing SDS-PAGE (Figure 7B–C). Our experiments do not support the mechanism proposed in Carvalho et al [21], [32], [33], [42], [50].
  The intention of this work is to explain the result of the Trx-ArsC mixed disulfide dissociation seen experimentally, since this was not understood on the basis of the experimental structure of the Trx-ArsC complex. The used methods to calculate pKa, regioselectivity and MD are robust and of state-of-the-art quality in the field [31]–[33], [42], [50]. Before forming a mixed disulfide between Trx and ArsC, Cys29Trx has to be stabilized as a thiolate to nucleophilically attack the exposed Cys82ArsC-Cys89ArsC disulfide in oxidized ArsC. In the active site of reduced Bs_Trx, both the α1-helix and the Cys29TrxSγ—Cys32TrxNH hydrogen bond stabilize Cys29Trx as a thiolate. In the presence of the Cys29TrxSγ—Cys32TrxNH hydrogen bond the calculated pKa of Cys29Trx (5. 5) is remarkably lower than the experimentally determined pKa of the corresponding residue in Ec_Trx1 (7. 1) [24] or Sa_Trx (7. 1) [1], in which the Cys29TrxSγ—Cys32TrxNH hydrogen bond is not present. The calculated pKa of the C-terminal cysteine of the WCGPC catalytic site (Cys32Trx) is 8. 1 (Table 1), and only slightly higher in Sa_Trx (pKa∼9) (24). The pKa difference between the Trx active site N-terminal and C-terminal cysteines implies that Cys29Trx and not Cys32Trx will attack the Cys82ArsC-Cys89ArsC disulfide in oxidized ArsC. These pKa arguments are concordant with the reactivity analysis: being less soft than Cys32Trx, Cys29Trx is more reactive towards the less soft Cys82ArsC-Cys89ArsC disulfide. The regioselectivity of Cys29Trx for Cys89ArsC is difficult to rationalize based on the pKa of the attacking nucleophile, or on the negative correlation between the rate of thiol-disulfide exchange reactions and the pKa of the leaving thiol group [51]–[54]. Upon formation of the mixed disulfide complex, either Cys82ArsC or Cys89ArsC can be the leaving thiol (Figure 2). Initially in ArsC, these thiols are flanking a flexible short helix – the so-called redox helix -, but during a single catalytic cycle, this redox helix unrolls, exposing the Cys82ArsC-Cys89ArsC disulfide [17]. As the chemical environment of Cys82ArsC and Cys89ArsC changes the moment the helix starts to unfold, their respective pKa values will change. We were not able to determine the pKa values of the leaving Cys82ArsC and Cys89ArsC thiols, neither experimentally nor by calculation since structural information of ArsC with an unfolded redox helix and with both cysteines reduced is lacking. However, we obtained information regarding the selectivity of the nucleophilic attack with a DFT reactivity analysis. The unfolded redox helix in oxidized ArsC sets a higher Fukui function value for Cys89ArsC compared to Cys82ArsC (Table 3). Thus, Cys89ArsC is intrinsically more reactive towards a nucleophilic attack than Cys82ArsC, leading to the formation of the Cys29Trx-Cys89ArsC mixed disulfide and the release of Cys82ArsC. In a productive cycle, the Cys29Trx-Cys89ArsC mixed disulfide is dissociated by the nucleophilic attack of Cys32Trx, leading to oxidized Trx and reduced ArsC [13], [16], [17]. The deprotonation mechanism of the Cys32Trx thiol has been debated for more than twenty years [21]–[26]. Structural analyses, kinetic assays, site directed mutagenesis and pKa titrations have led to the suggestion that Asp23Trx may deprotonate Cys32Trx [23]–[26]. Recently it was suggested that the leaving thiol group (here Cys82Trx) of the disulfide attacked during the first step deprotonates the C-terminal active site thiol of Trx [21]. Yet, our experiments strongly indicate that the Trx-ArsC mixed disulfide still dissociates when Cys82ArsC is mutated to Ala. In addition, our MD simulations strongly suggest that the local conformational changes around the mixed disulfide (Figure 6C) are responsible for the introduction of extra hydrogen bonds with the Cys32Trx sulphur (Figure 5C). Experimentally, it is not possible to measure the pKa of Cys32Trx in the Trx-ArsC complex, since the complexes formed between wild type Sa_Trx and Sa_ArsCTrip, and between the Sa_Trx D23A mutant and Sa_ArsCTrip were only transient. Due to the transient nature of the Trx-ArsC complex experimental measurements using the wild type sequence is a true challenge. This is in particular the case regarding pKa measurements and probing conformational changes leading to the dissociation. Further, the pKa of Cys32Trx is lowered by hydrogen bonds to the protein backbone. Thus, the role of these hydrogen bonds cannot be probed by additional mutagenesis. Therefore, in the present context, calculations fill a gap on questions that are currently experimentally out of reach, providing a more detailed picture of complex dissociation. The Cys32TrxSγ—Cys29TrxN hydrogen bond (Trx_ArsC_2 model, Figure 4D), present during 60. 0% of the time of the MD simulation, decreases the Cys32Trx pKa from 8. 9 to 8. 3 (Table 1). This was calculated from the Trx_ArsC_1_Cys32 and the Trx_ArsC_2_Cys32 model respectively. In selected snapshots of the activated complex where both the Cys32TrxSγ—Cys29TrxN and the Cys32TrxSγ—Trp28TrxN hydrogen bonds are present (and particularly strong based on geometric criteria), the pKa of Cys32TrxSγ drops to 7. 4, sufficient for the dissociation to proceed (Table 6, Figure 7). For the mixed-disulfide complex to dissociate, it is sufficient that Cys32Trx is transiently stabilized as a thiolate in the activated complex. Thus, instantaneous low values of the pKa of Cys32Trx are expected to be more relevant than the corresponding average value. Concomitant with the introduction of these hydrogen bonds, Cys32TrxSγ approaches the sulphur of Cys29Trx. We refer to this conformation where Cys32Trx is primed for its nucleophilic attack onto Cys29TrxSγ as the “activated complex”, as it is functionally relevant, and different from any conformation present in the NMR structure of the complex (PDB code: 2IPA). In sum, the pKa of Cys32Trx clearly depends on the presence of hydrogen bonds. These hydrogen bonds stabilize the thiolate form of Cys32Trx. Additionally, Cys32Trx is part of the α1-helix in Trx accounting for another decrease of the pKa up to 0. 8 pKa units. This α1-helix in Trx is solvent exposed and as such the effective dipole moment influencing the pKa of Cys32Trx is relatively low. The effect of the Cys32TrxSγ—Cys29TrxN hydrogen bond and the α1-helix on the calculated pKa is rather small, but can be considered as significant seen the very good correlation between NPA and pKa (Figure 3B). MD snapshots clearly strongly indicate the activation of Cys32Trx as a nucleophile via a conformational change, which brings Cys32Trx up to 3. 7 Å to Cys29Trx (Figure 6A) concurrently with the formation of the Cys32TrxSγ—Cys29TrxN and Cys32TrxSγ—Trp28TrxN hydrogen bonds. These hydrogen-bonds form even though Cys32TrxSγ was simulated in its neutral form. Simulating Cys32TrxSγ as a thiolate would presumably provide an even stronger driving force for the formation of such hydrogen bonds. Concomitantly with the formation of the activated complex, Asp23Trx moves 10 Å away from Cys32Trx, arguably too far away to deprotonate Cys32Trx. Indeed, mutating Asp23Trx to Ala cannot prevent complex dissociation (Figure 8B). This experiment demonstrates that Asp23Trx is not critical for activation of Cys32Trx. The same was found for the putative deprotonation role of Cys82ArsC. In the absence of Cys82ArsC, the mixed disulfide complex still dissociates. So, the transient hydrogen bonding of Cys32TrxSγ on the timescale explored during the MD simulations is enough to stabilize Cys32Trx as a nucleophilic thiolate, leading to the reduction of the disulfide bond in its substrate. Hydrogen bonding by Arg16, Arg108, Thr11 and Gly12 will drop the pKa of Cys82ArsC to 6. 3 (Figure S1 in Supplemental data). Under physiological conditions, not only Cys82ArsC but also Cys32Trx (pKa = 7. 4) predominates as a thiolate, due to two hydrogen bonds to its sulphur. Thiols with pKa values close to but lower than the pH of the solution react most rapidly [55], indicating that Cys82ArsC would be more reactive than Cys32Trx. However, Cys82ArsC never approaches the mixed disulfide during the MD simulations. Further, reactivity analysis predicts that Cys32Trx is less soft than Cys82ArsC and thus more prone to nucleophilically attack the less soft Cys29Trx-Cys89ArsC disulfide. During the Cys29Trx-Cys89ArsC complex dissociation, Cys29Trx or Cys89ArsC might be the leaving group of the thiol/disulfide exchange reaction (Figure 2). The calculated pKa of Cys89ArsC at the C-terminal of the partially unfolded redox helix in Sa_ArsC is 6. 7 [42], while the calculated pKa of Cys29Trx in Bs_Trx is 5. 5 (Table 1). Based solely on the negative correlation between the rate of a thiol-disulfide exchange and the pKa of the leaving thiol [51]–[54], the rate with Cys29Trx is expected to be higher than with Cys89ArsC as leaving group. However, Cys89ArsC is the leaving group. More important is the structure argument in which a thiolate attacks a disulfide bond preferentially along the sulphur-sulphur axis [56], [57] favoring the attack of Cys32Trx on Cys29Trx. Further, we calculated that Cys29Trx in the Cys29Trx-Cys89ArsC disulfide has a higher intrinsic reactivity than Cys89ArsC (Table 5). As such, the regioselectivity analysis emerged as consistent with other, independent, structural arguments, further strengthening the reactivity analysis. Thus, the correlation between the rate of a thiol-disulfide exchange and the pKa' s of the reacting and leaving thiols is only one aspect to consider and is insufficient to interpret the regioselectivity of a nucleophilic attack. The used models for the pKa and reactivity analysis include all hydrogen bond interacting residues with the considered cysteine. The MD, pKa and reactivity calculations are performed with state-of-the-art protocols which have proved their worth in a variety of contexts [29]–[33]. The theory of chemical reactivity is well established and widely used to study generalized acid/base reactions - including most of the organic reactions (additions, substitutions, eliminations) and the inorganic complexation reactions - of which thiol-exchange reactions are an example. The application of MD simulations, pKa and reactivity calculations and experimental studies yields a consistent picture of the studied mechanism, which at the very least indicates that the choice of methods is sound. Therefore, we argue that the results presented are stronger than just a plausible hypothesis. Proving chemical mechanisms is always very difficult and proceeds via successively improved working hypotheses. Thus, the present work cannot claim to have achieved a final proof. Yet, our proposed mechanism is grounded in a credible new structural model. It allows the exclusion of some currently proposed hypotheses, and crucially, it has explanatory power regarding regioselectivity, the activation of Cys32Trx and the lack of direct role for Asp23Trx. Therefore, the evidence gathered here arguably offers currently the most convincing and operational working model of the thiol/disulfide exchange reactions catalyzed by thioredoxin. To give fresh insights in the experimentally observed regioselectivity of Trx-ArsC complex dissociation, information from the MD structural analysis combined with the local softness used in the context of the HSAB principle, has provided a more complete and accurate picture. Equally important, we have demonstrated a methodology of general interest. MD simulations combined with pKa calculations, reactivity analysis and biochemical experiments offer a general and powerful strategy to study thiol/disulfide exchange reactions.
  The linear correlation between the pKa of cysteines of different Trx and ArsC systems and the Natural Population Analysis (NPA) charge on its Sγ-atom is obtained as described [37] (Table S1). For the used Trx and ArsC models, see Supplementary Data (legend Table S1 and Figure S3). The DFT calculations are performed with Bs_Trx and Bs_ArsC, while all experimental work is executed with Sa_Trx and Sa_ArsC, since no structures of reduced Sa_Trx or Sa_Trx-ArsC complex are available. Bs_ArsC and Sa_ArsC are structurally very similar [16], [18], [58] and they use a similar thioredoxin-coupled reaction mechanism [13], [16]–[18]. Model systems of reduced Bs_Trx (2GZY) [58], oxidized Bs_ArsC (1Z2E) [18] and of the Bs_Trx C32S – Bs_ArsC C10S/C15A/C82S complex (2IPA) [34] are made (Figure 4) for the pKa and reactivity analysis. They contain all residues that hydrogen bond interact with the considered cysteine residues. The protein environment was modeled by a continuum solvent model with a di-electric constant of 20 [59]. See Supplemental data (Text S1) for details. To calculate the pKa of the activated Trx-ArsC complex, different MD snapshots in which Cys32-SG makes “good” hydrogen bonds to Trp28-NH and Cys29-NH were used. Good hydrogen bonds were defined as follows: SG…N distance<4 Å and SG…H-N angle>150 deg, with both Trp28-NH and Cys29-NH. A model similar to the Trx_ArsC_2_Cys32 model was made. In all models, hydrogen atoms are placed and optimized together with the Sγ atoms of the reduced cysteine residues at the B3LYP/6-31G* level. Starting from these optimized geometries which are shown in Figure 4, NPA [41] charges for pKa estimations are calculated at the B3LYP/6-31+G** level in a continuum solvent model (PCM) with a di-electric constant of 20 [59]. Bondi radii are used to construct the cavity. An electrostatic model in which the optimized structures are represented as ChelpG point charges [42], [60], [61] was used as an approximation for the enzymatic environment, when the Fukui function and the (local) softness of Cys82, Cys32 and the Cys29-Cys89 mixed disulfide in Trx_ArsC_1 and Trx_ArsC_2, Cys29 in Trx_red and the Cys82-Cys89 disulfide in ArsC_ox are calculated. This procedure is necessary since the effect of taking one electron out or adding one electron to the system (needed to obtain the Fukui function) on the NPA charge is not sizeable in a hundred atoms large system. The ChelpG charges are calculated at B3LYP/6-31G**. Fukui functions and softnesses are calculated at the B3LYP/6-31+G** level as described in the Supplemental data (Text S2). All calculations are performed using the GAUSSIAN 03 package [62]. Energy minimizations and MD simulations of the Bs_Trx-ArsC mixed-disulfide complex were performed with the program CHARMM [63], the version 22 of its protein force-field [64], [65], a dielectric constant of 1. 0, and atom-based non-bonded interactions truncated beyond 12 Å with force-shift [66]. There is ample evidence that the 12 Å force-shift spherical cut-off used in the present work performs well, and as well as alternative PME schemes also in use currently. This is evidenced by detailed studies, which compared spherical cutoffs to PME and concluded that both schemes work equally well [67]–[69]. In addition, a recent study of a DNA oligomer in solution obtained very similar results with either CHARMM and the force-shift cutoff or AMBER and PME [70]. Also, the 12 Å force-shift spherical cut-off was used across a series of studies on proteins in the Trx superfamily. These studies characterized the structure, dynamics and pKa' s of the ionized active-site cysteines, alongside detailed comparisons to experimental observables [31]–[33]. In every case, excellent agreement was found between simulations and experiment, even at predictive level. This is a very strong validation of the 12 Å force-shift spherical cut-off. Therefore, we continue to use this cutoff in the present work for reasons of consistency and to allow for comparisons across different proteins of the Trx superfamily. Non-bonded lists were maintained to 14 Å and updated heuristically. The initial coordinates were from conformer 7 of the NMR structure (2IPA). To stabilize this complex for NMR work, Cys32Trx, Cys10ArsC and Cys82ArsC had been mutated to Ser residues, and Cys15ArsC to Ala. For the simulations, these residues were restored to their wild type sequence by a straightforward modeling operation. Then, two complexes were setup, with Cys82ArsC neutral or as a thiolate. Force-field parameters for the thiolate were as described previously [31]. Conventional protonation states were assigned to other titratable residues. The complexes were overlaid with a rhombic dodecahedron of pre-equilibrated CHARMM TIP3P water [71], [72], with a normal distance of 88 Å between opposing faces of the dodecahedron. The protein did not see its image during the simulation since the shortest distance between the protein in the primary cell and its image was initially 21 Å, which is much larger than the longest cut-off used in the calculations (12 Å). The water molecules overlapping with the protein were removed. The systems were neutralized by the addition of sodium ions. Periodic boundary conditions were applied, and all covalent bonds involving a hydrogen were constrained with SHAKE [73]. The system was energy minimized in several stages, with the solute initially fixed and finally allowed to fully relax. The system was then submitted to MD simulations, using the leap-frog integrator, a 0. 002 ps timestep, and the NPT ensemble. One MD simulation per complex (Cys82ArsC neutral or thiolate) was performed. Heating was performed from 50 K to 300 K in 5ps by 5 K increments, with the protein atoms harmonically constrained to their initial position with a force constant of 2. 0 kcal/mol/Å2. The constraints on the solute were kept for a further 20 ps of equilibration. Then, the simulations were pursued at 300 K for at least 75 ns with each system. Details on the Site-directed mutagenesis, expression and purification of Sa_Trx and Sa_ArsC can be found in Supplemental Data (Text S3). An aliquot of purified Sa_ArsCTrip was incubated with 20 mM DTT for 30 min at room temperature to assure that the thiol was fully reduced. The excess of DTT was removed on a Superdex75 HR (10/30) column equilibrated in poly buffer A, containing 10 mM Na-borate, 10 mM Na-phosphate, 10 mM Na-citrate. The fractions containing ArsC were pooled, concentrated and 10 mM DTNB was added. The reaction was monitored at 412 nm for the completion of the reaction. The reaction mixture was applied a second time to the gel filtration column to remove excess DTNB and released TNB−. The fraction of mixed disulfide Sa_ArsCTrip–TNB was concentrated and used directly in the reaction with wild type or Sa_Trx D23A. Wild type and D23A Sa_Trx were incubated with 20 mM fresh DTT at room temperature for 30 minutes to assure that the remaining active site thiol was in the reduced state. The excess of DTT was removed on a superdex75 HR (10/30) column equilibrated in poly buffer A. Wild type and Sa_Trx C32A were reacted with an equimolar concentration of Sa_ArsCTrip at room temperature. Upon completion of the reaction, the mixture was loaded on a Superdex75 HR (10/30) column equilibrated in poly buffer A. Complex formation was evaluated on a non-reducing SDS-PAGE and on reversed phase. The reaction mixture was injected on a Grace Vydac C18 column (4. 6 mm×250 mm) equilibrated in 15% (v/v) acetonitrile, 0. 1% (v/v) trifluoroacetic acid (TFA) at 0. 8 ml/min. The column was eluted with a 30 min linear gradient from 15% to 60% acetonitrile at room temperature. Absorption data collection at 280 nm was performed under Empower (Waters, Milford Massachusetts).
summary: Thioredoxins are found in all types of cells and control several essential functions of life, including promotion of cell growth, inhibition of apoptosis, and modulation of inflammation. Thioredoxin has two ‘free’ cysteines in its active site, which are used to break disulfide bonds in oxidized substrate proteins. In the first step, an intermediate thioredoxin-protein complex is formed, which is broken in the second step, releasing the substrate protein in its reduced state. In other words, the disulfide is being transferred from the substrate protein to thioredoxin, or the electrons coming from thioredoxin are shuttled to the protein substrate. The exact reaction mechanism, i. e. , the detailed succession of steps in which the reaction takes place, of how this mixed disulfide is broken is not known and has been debated over the last twenty years. With a multidisciplinary approach, combining computational and experimental work, we provide fresh insights into how conformational changes activate the catalytic cysteines with which this universal reduction mechanism is completed with the correct regioselectivity. This work illustrates the strengths of computational approaches in probing phenomena which are otherwise very difficult to investigate experimentally.
keywords: chemical biology
  biochemistry/protein chemistry
  biophysics/theory and simulation
title: How Thioredoxin Dissociates Its Mixed Disulfide