Water-assisted Proton Transfer in Ferredoxin I*

The role of water molecules in assisting proton transfer (PT) is investigated for the proton-pumping protein ferredoxin I (FdI) from Azotobacter vinelandii. It was shown previously that individual water molecules can stabilize between Asp15 and the buried [3Fe-4S]0 cluster and thus can potentially act as a proton relay in transferring H+ from the protein to the μ2 sulfur atom. Here, we generalize molecular mechanics with proton transfer to studying proton transfer reactions in the condensed phase. Both umbrella sampling simulations and electronic structure calculations suggest that the PT Asp15-COOH + H2O + [3Fe-4S]0 → Asp15-COO− + H2O + [3Fe-4S]0 H+ is concerted, and no stable intermediate hydronium ion (H3O+) is expected. The free energy difference of 11.7 kcal/mol for the forward reaction is in good agreement with the experimental value (13.3 kcal/mol). For the reverse reaction (Asp15-COO− + H2O + [3Fe-4S]0H+ → Asp15-COOH + H2O + [3Fe-4S]0), a larger barrier than for the forward reaction is correctly predicted, but it is quantitatively overestimated (23.1 kcal/mol from simulations versus 14.1 from experiment). Possible reasons for this discrepancy are discussed. Compared with the water-assisted process (ΔE ≈ 10 kcal/mol), water-unassisted proton transfer yields a considerably higher barrier of ΔE ≈ 35 kcal/mol.

Water can affect chemical and biological systems at various levels, including the formation of specific interactions through hydrogen bonding, the screening of Coulomb interactions, the mediation of proton transfer, or as an intrinsic component in the secondary structure of proteins. Usually, directly probing this role is difficult because of the transient nature of the processes involved. The difficulty also extends to probing the role of single water molecules in chemical and biological catalysis. Particularly in proton and hydrogen transfer reactions, the role (or absence of it) of water is often postulated but cannot be unequivocally proven. Examples include water involvement in ribozyme catalysis (1) where a large body of data suggests water involvement in RNA backbone trans-esterification or proton transfer in ferredoxin where the experiments were explained as a direct protein-to-[3Fe-4S] cluster transfer without participation of water (2,3).
There are three profoundly different ways to characterize the role of water in a specific context. First, the water degrees of freedom can be essentially averaged out, and not much insight at the atomistic level can be gained. Second, a few individual water molecules are singled out, and their role can be analyzed in detail. Third, different types of water molecules can be distinguished, for example "surface-bound" versus bulk water molecules (4). Individual water molecules have been implicated in mediating proton transfer, and their role has been characterized by spectroscopic means (5). Structural studies, on the other hand, have established that individual water molecules play a central role in protein-ligand interactions such as in carbohydrate-protein binding (6) or in HIV-1 protease (7). In protein folding and protein-protein interactions, water has been found to act as a lubricant or "facilitator" in protein recognition (8,9), and structural waters can render proteins more flexible (10). Also, it has been argued that water may even be conserved evolutionarily as an integral part of a protein structure (10).
Water has also been suggested to play important roles in proton translocation (11,12). Two prominent examples are electron transfer processes (13) and ATP hydrolysis (14,15). In general, proteins performing such proton transfer (PT) 3 reactions are referred to as "proton pumps," among which cytochrome c oxidase found in the mitochondrial electron-transfer chain (16 -19) and bacteriorhodopsin found in photochemical reaction centers (20 -23) are prime examples. Experimental observation of an enzymatic proton transfer has been found by Fourier transform infrared difference spectroscopy (24 -26). However, direct experimental observation of PT at an atomistic level is generally difficult because it is a transient process.
Complementary to experimental work, theoretical and computational methods have been used to provide a more detailed understanding of PT processes (27)(28)(29)(30)(31)(32)(33)(34)(35)(36)(37)(38)(39). These methods are often based on empirical potentials (29 -32), Car-Parrinello, or related approaches (33)(34)(35)(36)(37)39). With standard force fields, it is not possible to examine in detail the dynamics of the proton transfer itself due to their inability to describe breaking and forming of chemical bonds. One way to circumvent this problem is to use mixed quantum mechanics/molecular mechanics calculations (QM/MM) (40,41). They decompose the system into a part that is directly involved in the reaction and treat it with quantum mechanics, whereas the rest of the system is treated with a molecular mechanics force field. Other methods * This work was supported by the Schweizerischer National Fonds through that use a QM/MM separation are empirical valence bond theory (42) or approximate valence bond theory (43). Recently, we introduced a new method for studying PT in protein-sized systems using molecular dynamics (MD) simulations (44,45). This approach, named molecular mechanics with proton transfer (MMPT), is inspired by QM/MM simulations but combines a potential energy surface (PES; the "QM" part), suitable for describing the proton transfer between an acceptor and a donor atom, with a force field (the MM part) for the remaining degrees of freedom. The MMPT potential includes a modified treatment of hydrogen bonding that allows for the formation and breaking of bonds involving the proton being transferred (45). The main advantage of this approach over QM/MM is the performance, which is comparable with a force field simulation and makes it applicable for long time simulations of large systems such as proteins. In the present work, we present a generalization of MMPT to condensed phase systems and apply it to proton transfer in ferredoxin (3,46,47).
Ferredoxins are a family of small iron-sulfur proteins that mediate electron transfer, which are involved in such fundamental biological roles as photosynthesis and nitrogen fixation (48). One ferredoxin for which considerable amounts of structural and experimental data are available is ferredoxin I (FdI) from Azotobacter vinelandii, a nitrogen-fixating soil bacterium. It is known that the one-electron reduction of the iron-sulfur cluster in FdI is immediately followed by the uptake of a proton from the solvent. The mechanism of this proton transfer, which can be used as a model of a redox-driven proton pump, has been the subject of particular attention (46,3,49,50). The kinetics of the electron-coupled proton transfer were probed experimentally using cyclic voltammetry, which, combined with site-directed mutagenesis, indicates that Asp 15 , a surface residue, plays an important role in catalyzing the proton transfer (51,3). The simplest possible mechanism would involve a direct (water-unassisted) transfer of the proton from Asp 15 to give the protonated [3Fe-4S] 0 H ϩ cluster. However, proton tunneling under such conditions is limited to distances of Ͻ0.25 Å (52, 53), and detailed MD and quantum chemical calculations show that the aspartic acid side chain in FdI is too far from the nearest 2 -sulfur for efficient proton transfer to occur (45,47). An alternative possibility is that proton transfer between Asp 15 and S is mediated by a water molecule. Although no crystallographic water molecules near the iron-sulfur cluster were reported, detailed atomistic simulations showed that the active site of FdI is water-accessible and that the active site water is stabilized over extended periods of time and could potentially act as a proton relay (47,54). It is known from experiment that when interior waters are mobile, they may not be detected in x-ray structures, and in many cases, the number of observed water molecules is smaller than that actually present (55)(56)(57). Thus, a water-mediated mechanism is plausible. The overall reaction studied with density functional theory (DFT) and free energy simulations in the present work is Asp 15 The forward barrier is found to be in good agreement with the experimentally determined rates, whereas the reverse barrier is somewhat overestimated. Contrary to this, the water-uncatalyzed reaction, which is also characterized by density functional theory methods, yields much higher barriers and is unlikely to occur.

MATERIALS AND METHODS
Electronic Structure Calculations-Following the previous study on FdI (54), electronic structure calculations were carried out at the spin-unrestricted level with the UB3LYP DFT functional and a 6 -31G(d,p) basis set, using the Gaussian 03 suite of programs (58). Three model systems differing in complexity (models A to C) were considered for scanning the potential energy surface. They are shown in Fig. 1, A-C. Model A (33 atoms) includes the [3Fe-4S] 0 cluster (total charge of the system q ϭ Ϫ3, total spin S ϭ 2) with thiomethoxy side chains replacing residues Cys 8 , Cys 16 , and Cys 49 from FdI; a water molecule coordinating to the S1 atom of the cluster; and acetic acid (AcOH) to model the Asp 15 side chain from the 7FDR x-ray structure recorded at 1.4 Å resolution (59). Internal bond and angle geometries of [3Fe-4S] 0 initially taken from 7FDR were optimized prior to proton transfer scans in this system. Model B (59 atoms, q ϭ Ϫ3, total spin S ϭ 2) is an extension of model A with a complete representation of residue Asp 15 and parts (terminated by hydrogen atoms) of residues Cys 16 , Thr 14 , and Tyr 13 to account for interactions of the nearby protein backbone with the proton transfer partners. The geometries of the cluster and amino acid backbone were taken from 7FDR, and the internal coordinates of the Asp 15 side chain and the water molecule were optimized before the potential energy scans. A version of this model with the water molecule deleted was also used to study water-unassisted proton transfer. For the MMPT parametrization (see below) model system C, including a protonated [3Fe-4S]H ϩ cluster with the same thiomethoxy side chains attached as in model A, and a water molecule was used (26 atoms, q ϭ Ϫ2, total spin S ϭ 2; see Fig. 1C). Internal bonds and angles of the cluster and all degrees of freedom of the water molecule were optimized prior to calculating the interaction potentials.
Density functional theory scans for the complete proton transfer including the initial transfer (PT1) from Asp 15 -OH to a water molecule and a second proton transfer (PT2) from the hydronium ion intermediate H 3 O ϩ to the [3Fe-4S] cluster were performed for models A and B of Fig. 1 along a reaction coordinate defined by the OD1 to HD1 distance r OD1-HD1 . For both models, internal and external coordinates of water were allowed to relax during the scans, whereas the coordinates of the cluster were kept fixed. This is a reasonable assumption as a comparison of the Fe-S distances between optimized protonated ([3Fe-4S]H ϩ ) and unprotonated ([3Fe-4S] 0 ) clusters show only small differences (0.05 Å) (59). In model A, all external coordinates of AcOH except the donor oxygen OD1 to proton acceptor atom S1 distance (R OD1-S1 ) were allowed to relax. R OD1-S1 was constrained to three different values: R OD1-S1 ϭ 4.4, 4.7 (the value found in the x-ray structure of FdI), and 5.0 Å, resulting in three individual scans for model A. For the energy scan of model system B, only the degrees of freedom of the water molecule and the angle and dihedral connecting HD1 to OD1 of Asp 15 , were optimized at each point of the scanning coordinate. Relaxing more degrees of freedom is computationally very demanding for this structure. Each scan in model sys-tems A and B in Fig. 1 was performed along r OD1-HD1 . Starting from r OD1-HD1 ϭ 0.8 Å, the distance was increased in steps of 0.1 Å until the transition state (TS) was reached, and atom H2 of the water molecule was transferred to S1 of the [3Fe-4S] cluster. Close to the transition state, the step size was reduced to improve the resolution of the potential scan in this region.
To parametrize the MMPT potential for PT2 a two-dimensional scan of the potential energy surface was calculated starting from the energy-optimized model C in Fig. 1C. The S1 to OH2 distance R 2 and the S1 to HS distance r 2 along the hydrogen bond were scanned, whereas all other degrees of freedom were kept rigid. A total of 347 points were calculated. The resulting energy surface (shown in Fig. 2A) has a single minimum that corresponds to having the proton closer to S1, forming a strong hydrogen bond to the water molecule (see Fig. 1C).
MD Simulations and Intermolecular Interactions-All atomistic simulations were carried out with CHARMM (60) and the CHARMM22 force field (61) with provisions for the MMPT potential (45). Further details of the simulations are given below and in the supplemental data. Conventional force fields cannot describe the breaking or formation of bonds because stretching potentials are typically parametrized as harmonic oscillators. In contrast, MMPT uses a functional form for V PT (R, , ), which is based on Morse potentials and allows for describing proton or hydrogen transfer between D and A. The total potential energy of the system with respect to all coordinates x ជ is decom-posed into a part for the proton transfer motif (V PT (R, , )) and the remaining degrees of freedom y ជ of the system. Thus, the total interaction is written as shown in Equation 1, where R is the distance between the heavy atoms, is the PT progression coordinate, defined as ϭ (r Ϫ 0.8)/(R Ϫ 1.6), and is the hydrogen bonding angle. In the following, PT1 takes place between Asp 15 -COOH and water (OD1-HD1⅐⅐⅐OH2) and PT2 between [3Fe-4S] and water (S1-H2⅐⅐⅐OH2) (see Fig. 1). According to this definition, PT1 and PT2 start at each end of the overall PT reaction and develop toward the hydronium ion. Thus, for PT1, is the hydrogen bonding angle, OH2-HD1-OD1, and for PT2, it is OH2-H2-S1. Subsequently, the variables R 1 , 1 , and 1 refer to PT1 and variables R 2 , 2 , and 2 to PT2.
A consequence of using a large model system such as the one in Fig. 1C as a reference for fitting is that nonbonded interactions beyond the immediate vicinity of the proton transfer play a potentially important role. This is unlike the case of the prototype systems used previously, such as H 5 O 2 ϩ , where these interactions are included in V PT2 (62). As a result, it is no longer appropriate to fit V PT directly to the DFT potential energy surface because that would result in double counting of the nonbonded interactions when V PT is added to the CHARMM force field. Instead, V PT has to be fitted to the difference between the DFT energy, and the energy computed using the force field

Water-assisted Proton Transfer in Ferredoxin I
. This relation is shown graphically in Fig. 2B for PT2. The difference between the potential including doublecounted interactions (V PT2 , black curve in Fig. 2B) and that without it (gray curve in Fig. 2B) remains nearly constant in the region where PT takes place. V PT2 is shifted by Ϸ 6 kcal/mol after subtraction of V MM (y ជ). For the problem studied here, this means that total energies are moderately influenced by double counting, whereas the PT barrier remains nearly unaffected. The procedure of refitting the surface to account for the double counting of nonbonded interactions is described in detail in the supplemental data.
To distinguish between the different protonation states, the following nomenclature is used: "state 1" has Asp 15 protonated (-COOH), "state 2" has the water molecule protonated (H 3 O ϩ ), and "state 3" has the iron-sulfur cluster protonated ([3Fe-4S]H ϩ ). As a first approximation, the transition between state 1 and state 2, and between state 2 and state 3 may be assumed to lie at 1 ϭ 2 Ϸ 0.5. However, this assumption can be refined based on the shape of the free energy surface obtained from the simulations.
Potentials of Mean Force-One-dimensional potentials of mean force were calculated from MD simulations using umbrella sampling (63). The simulations were started from an equilibrated structure of the reduced protein (Protein Data Bank code 7FDR) at 300 K in a 62 Å cubic TIP3P (64) water box with periodic boundary conditions. The time step in the simulations was 0.5 fs, and all bonds to hydrogen atoms except for those in the water molecule involved in the proton transfer were constrained using the SHAKE algorithm (65). The potentials of mean force are built from individual 50-ps simulations (after 2 ps of further equilibration) with harmonic biasing potentials ranging from k umb ϭ 30 to 600 kcal mol Ϫ1 along the reaction coordinates. To follow PT1, the driving coordinate was defined as 1 ϭ r 1 /R 1 , where r 1 is the OD1-HD1 distance, and R 1 is the OD1-OH2 distance (see Fig. 1A). The correspond-ing coordinate for PT2 is 2 ϭ r 2 /R 2 , where r 2 is the S1-H2 distance, and R 2 is the S1-OH2 distance already used above. 1 and 2 were modified in intervals of 0.05 along the umbrella sampling simulation. The -coordinate was chosen as a proxy for for practical reasons. A graphical representation for the relationship between and is given in the supplemental data. To account for some of the charge transfer that accompanies proton transfer, atomic charges appropriate for each value of the reaction coordinates were modified (see supplemental data). The data from individual simulations was combined with a weighted histogram analysis (66).
Minimum Energy Pathway-The potential energies along 1 and 2 were also calculated starting from the MMPT-equilibrated protein structures in states 1 and 3, respectively, with both proton transfers occurring toward state 2. The absolute minimum along the PT scan was determined prior to constraining the system at different values of by running 1000 steps of steepest-descent minimization followed by multiple steps of adopted basis Newton-Raphson minimization with a gradient-based cut-off of 10 Ϫ6 . Subsequent adopted basis Newton-Raphson minimization with the same average gradient tolerance were carried out starting from the minimum energy structures by constraining the system to the average -geometries obtained from each individual umbrella sampling window.

RESULTS
Proton Transfer from DFT Calculations on Model Systems-Potential energy scans as described in "Materials and Methods" were carried out for model systems A and B (Fig. 1, A and B) along the OD1-HD1 distance r 1 . This yields the minimum energy path (MEP) for double proton transfer and allows to identify the geometries of formation and decomposition of the intermediate H 3 O ϩ species. Furthermore, it is possible to locate the transition state. The corresponding MEPs are shown in Fig. 3. The initial transfer (PT1) from acetate to the water molecule, which leads to AcO Ϫ ⅐⅐⅐H 3 O ϩ , completes at r 1 ϭ 1.40 For model B, only one scan was carried out for R OD1-S1 ϭ 4.8 Å. The TS for the forward double-PT reaction is at r 1 ϭ 1.58 Å with a barrier of 30.2 kcal/mol. The difference between models A and B is related to the fewer degrees of freedom that are allowed to relax in model B. Thus, the system has more strain, which is also indicated by the much steeper slope of the surface. On the other hand, the back transfer barrier of this model system is only 8.1 kcal/mol. Overall, the electronic structure calculations consistently find no stabilized intermediate (H 3 O ϩ ), and the barrier for proton transfer is sensitive to rearrangements in the Asp 15 side chain, which are possible only for the relaxed model A scans.
It is also possible that such a concerted PT reaction involves electronic coupling. To test this, a reverse scan was performed for R OD1-S1 ϭ 4.4 Å starting in state 3. The energy of the transition state differs by only 0.2 kcal/mol compared with the forward scan but because the potential energy surface is flat around the transition state, the S1-H2 distance differs by Ϸ 0.25 Å between the two structures. Despite this considerable geometrical difference the largest deviation between the Mulliken charges of the atoms involved in PT does not exceed 0.07 e. Thus, the different states are not coupled electronically, and the process is ground-state proton transfer.
Reactions PT1 and PT2 from Umbrella Sampling-To validate the MMPT force field in the protein environment, constant energy MD simulations were carried out. For this, 100 ps were run starting from an equilibrated structure of solvated 7FDR. All hydrogen atoms, except for the hydrogen atoms in the water molecule as well as the proton being transferred, were constrained using SHAKE (65). A time step of 0.5 fs was used to follow the rapid H-motion explicitly. Fig. 4 reports the histogram of E tot (t) and temperature T(t), which show no drift over the time interval studied here and establish that the computational procedure is reliable and meaningful.
The experimental rate constant describing proton transfer from the protein (Asp 15 ) to the [3Fe-4S] 0 cluster was determined from fast scan protein film voltammetry at pH 8.34 as k on hop ϭ 1294 Ϯ 100 s Ϫ1 and the back transfer rate from [3Fe-4S] 0 H ϩ to the protein (Asp 15Ϫ ) as k off hop ϭ 332 Ϯ 25 s Ϫ1 (2). According to transition state theory, such a rate k TST is related to the free energy of activation ⌬G ‡ by Equation 2, where k B is the Boltzmann constant, T is the temperature in K, h is the Planck's constant, and R is the gas constant (67). From the experimental rate constants (k on hop and k off hop ) free energies of activation ⌬G on ‡ ϭ 13.3 Ϯ 0.05 and ⌬G off ‡ ϭ 14.1 Ϯ 0.05 kcal/ mol for the forward and backward reaction are obtained, respectively. Such barriers are too large to be accessible to unconstrained molecular dynamics simulations. Therefore, the reactive steps were in the following investigated by using umbrella sampling (63).
The first proton transfer (PT1) Asp 15 -COOH ϩ H 2 O 3 Asp 15 -COO Ϫ ϩ H 3 O ϩ was reinvestigated with the asymmetric potential energy surface (see Ref. 45) but now included fluctuating charges on the atoms involved in the proton transfer motif (see "Materials and Methods" and supplemental Fig. S3). The potential of mean force from umbrella sampling simulations is shown in Fig. 5A (solid line). The free energy curve does not show a barrier but a characteristic flattening around 1 ϭ 0.56,  which is related to formation of H 3 O ϩ . This configuration is ⌬G PT1 ‡ ϭ 11.7 kcal/mol above the Asp 15 -COOH ϩ H 2 O state, which is somewhat higher than in the earlier work, which was 8.2 kcal/mol (45). This difference is related to the use of fluctuating charges in the present work, which takes into account partial charge transfer from OD1 of Asp 15 to OH2 of water (for PT1) and subsequently to S1 of the cluster (for PT2), see supplemental data. The dotted and dashed curves in Fig. 5A correspond to the gas-phase minimum energy path from PES 1 and the path calculated in the condensed phase (protein) environment. The gas-phase minimum is shifted from 1 ϭ 0.34 to 1 ϭ 0.36 compared with the minimum energy path (from condensed phase minimizations) and the potential of mean force (from MD simulations), respectively. Fig. 5B reports a projection of the average PT coordinate from umbrella sampling onto the corresponding potential energy surface from DFT calculations from Ref. 45. In the minimum of state 1, R 1 lies between 2.8 to 2.9 Å. Decreasing 1 increases R 1 and thus separates the water from Asp 15 . As the reaction progresses, the water approaches R 1 ϭ 2.5 Å, and as soon as the proton is completely transferred (in state 2; 1 Ϸ 0.7), R 1 rises again with increasing 1 , which corresponds to the positively charged hydronium ion moving away from the deprotonated Asp 15 residue.
For PT2, the potential of mean force is shown in Fig. 5C as a solid line. The minimum at 2 ϭ 0.36 corresponds to the proton closer to the S1 sulfur atom (although forming a strong hydrogen bond with the water molecule). For comparison, the gasphase minimum energy path from the potential energy surface in Fig. 2A is also shown (dotted line in Fig. 5C). The two curves are qualitatively similar, although the MEP rises more steeply at higher 2 . The dashed line represents the condensed phase minimum energy path. Between 2 ϭ 0.4 and 0.5, the curve is Ϸ 3 kcal/mol higher than the gas-phase MEP but for 2 Ͼ 0.5, it lies between the gas-phase MEP and the condensed phase PMF. This suggests that the protein environment provides additional stabilization of the hydronium ion. State 2 appears around 2 ϭ 0.65 in Fig. 5C, where the gradient of the potential of mean force flattens out. The corresponding free energy difference relative to Fig. 5D shows the average trajectory of the PMF projected onto the R 2 , 2 surface. At small 2 , the water molecule is far from the S1 sulfur (R 2 ϭ 4.0 Å). The initial increase of 2 is primarily related to a reduction of R 2 , i.e. the water molecule shifts toward the acceptor S1. The actual PT2 takes place at R 2 Ϸ 2.5 Å. Beyond 2 ϭ 0.8, R 2 increases again, which corresponds to the completed formation of the hydronium ion. Subsequently, it is detached from the [3Fe-4S] cluster without the formation of a stable minimum. Thus, the overall proton transfer reaction is most likely concerted.
Having identified proton transfers PT1 and PT2 separately with both exhibiting a flattening of the PMF around state 2, it is also possible to connect the two processes. This amounts to a diabatic treatment around the transition state and provides an overarching description of H ϩ transport from the protein via the water molecule toward the buried iron-sulfur cluster. Fig. 6, A and B, shows typical structures from the umbrella sampling simulation for PT1 and PT2 in state 2. PT1 occurs with the proton-accepting water molecule preferentially arranged near the bulk. During PT2, the same water molecule coordinates to the S1 atom of [3Fe-4S] and the OD1 atom of Asp 15 . The structural transition from state 2 of PT1 (Fig. 6A) to the corresponding state for PT2 (Fig. 6B) therefore consists mainly of a rotation of the hydronium ion. The two separate potentials of mean force from Fig. 5, A and C, describing PT1 and PT2 are shown together in Fig. 6C. It should be noted that there is no a priori relationship between the two progression coordinates 1 and 2 , which is illustrated by the dashed line that symbolizes the uncharacterized parts of the free energy surface. However, as the two structures of Fig. 6, A and B, demonstrate, the PES is probably very flat in this region. Because both potentials of mean force do not have a barrier for reaching state 2 (H 3 O ϩ ) from their corresponding minima, the merged PMF also consists of only one transition state region with the previously determined forward (⌬G PT1 ‡ ϭ 11.7 kcal/mol) and backward (⌬G PT2 ‡ ϭ 23.1 kcal/mol) barriers. In each of the density functional theory scans (models A and B), the initial transfer from OD1 to the water (PT1) is completed at an OD1-HD1 distance of r 1 ϭ 1.4 Å and a OD1-OH2 distance of R 1 ϭ 2.50 Å. The OH2-HD1 distance at this point therefore is 1.1 Å ( ϭ 180°), and PT1 is regarded to be completed. Nearly the same average geometry was sampled in the potential of mean force of PT1 using MMPT at 1 ϭ 0.56. This is the region where the surface in Fig. 5A is flattening out and suggests completed formation of H 3 O ϩ . Geometries determining state 2 of PT2 were extracted from scans of model A with density functional theory at a S1-OD1 distance of 5.0 Å (dotted line in Fig. 3). In this scan, the H 3 O ϩ species remains stable up to r 1 ϭ 1.8 Å. The PT2-related coordinates at this point are r 2 ϭ 1.96 Å and R 2 ϭ 3.00 Å corresponding to 2 ϭ 0.65, which is the region on the PMF of PT2 (Fig. 5C) that exhibits flattening.

DISCUSSION
In the present work, atomistic simulations with provisions to describe hydrogen-or proton-transfer reactions and electronic structure calculations were combined to characterize waterassisted proton transfer in a protein. A concrete procedure for using the MMPT potential in extended and fully solvated condensed-phase systems is presented. The strategy is based on model ab initio calculations around critical points of the interaction potential and avoids double counting of nonbonded interactions.
For proton transfer between the Asp 15 side chain and the buried [3Fe-4S] cluster in FdI, both approaches (MMPT and DFT) support a concerted double PT in which PT1 is the ratedetermining step. Previous computational studies also concluded that PT reactions on or near a protein surface proceed semi-to fully concerted rather than stepwise (38,39,68). However, these studies were carried out either at the semiempirical level with variational transition state theory (which neglects explicit dynamics) (68) or with calculation of minimum energy pathways on an ensemble of structures (39), or by using density functional theory without dynamics (38).
Computed values for ⌬G PT1 ‡ and ⌬G PT2 ‡ can be directly related to the free energies of activation ⌬G on ‡ and ⌬G off ‡ evaluated from experimentally measured rate constants. ⌬G PT1 ‡ (Fig.  5A) underestimates ⌬G on ‡ slightly by 1.6 kcal/mol, whereas ⌬G PT2 ‡ overestimates the experimental barrier ⌬G off ‡ by 9.0 kcal/mol and is only in qualitative agreement with experiment. Taking zero-point energy corrections into account, both barriers will be further lowered by Ϸ 2 kcal/mol. Furthermore, the experimental rate k on hop ϭ 1294 s Ϫ1 was measured at a pH of 8.34. At a physiologically more relevant acidic pH (69,70), the experimental rate for PT2 is smaller (k off hop ϭ 36 Ϯ 5 s Ϫ1 at pH 5.0 (2)), which, according to transition state theory, corresponds to a barrier increase by ⌬G PT2 ‡ ϭ 1.4 kcal/mol compared with pH 8.34. In the simulations for PT2, the carboxyl group of Asp 15 is deprotonated, which covers the situation down to a pH of Ϸ 4.1 (the pK a of Asp lies at 4.1 (71)). Because the measured barriers for proton transfer are pH-dependent, this has to be taken into account in comparing with the computed value. Applying both zero-point corrections (2.0 kcal/mol) and corrections due to the pH dependence (1.4 kcal/mol), the difference between the computed and the experimentally determined barrier for PT2 reduces to 5.6 kcal/mol.
Similarly to this, DFT scans (model compound A) at different R OD1-S1 show larger variations of barrier heights for PT2 compared with PT1. A possible effect in umbrella sampling is local relaxation of the environment which can lead to an overstabilization of one state over the other (72). To assess such relaxation effects, 100 individual and unconstrained MMPT MD simulations starting in state 2 ( 1 ϭ 0.56 and 2 ϭ 0.65; i.e. proton-transferred state) were run and the potential energy along 1 and 2 was collected until they reached the potential minima. The fluctuations ␦V of the potential energy (indicated as bars in Fig. 6C) remain moderate for PT2 down to 2 Ϸ 0.55 but increase sharply up to Ϸ 10 kcal/mol when approaching geometries of the [3Fe-4S]H ϩ state. Corresponding fluctuations for PT1 are of similar size for large 1 but remain moderate (␦V Ϸ 5 kcal/mol) for structures around the minimum of PT1. It is important to note that these fluctuations are not "errors;" rather, they reflect the variation in potential energy of downhill dynamics from the TS toward the respective local minimum, and thus, they also exclude entropic contributions. Similar amplitudes for instantaneous relaxation in large condensed phase systems have previously been found to lead to significant changes in barrier heights due to long range forces and global perturbations (73). In conclusion, both DFT and MMPT MD simulations suggest stronger fluctuations for PT2 compared with PT1. An additional reason for the different degrees of agreement of the forward and reverse barriers for PT1 and PT2 is the increased complexity of PT2. PT1 involves essentially PT between COOH and water, which is a relatively simple process to capture with electronic structure calculations. On the other hand, proton transfer between [3Fe-4S]H ϩ and water is much more challenging and quantitative information is likely to be more difficult to obtain.
The transition state energies calculated with models A (smaller, more degrees of freedom relaxed) and B (larger, fewer degrees of freedom relaxed) differ considerably (10.1 kcal/mol versus 30.2 kcal/mol). This can be primarily related to the increased flexibility allowed for in model A, which is also supported by the umbrella sampling simulations that find a barrier energy of ⌬G on ‡ ϭ 13.3 kcal/mol. Comparing the structures of model A obtained in the minimum of state 1 with the one of state 3 shows that AcO Ϫ reorients considerably during proton transfer. The RMSD of AcO Ϫ between the two structures aligned to the [3Fe-4S] atoms is 5.8 Å for the scan with R OD1-S1 ϭ 5.0 Å. The important structural differences leading to such a large rearrangement are a distance increase of 1 Å between the CG atom of AcOH and the S1 atom of [3Fe-4S] in going from state 1 to state 3 and more noticeably a near inversion of AcOH along the CB-CG vector. This finding corroborates the swinging-arm conduction mechanism suggested from the rapid-scan voltammetry on FdI (2), which is also referred to as a "piggy back" mechanism. Other long range proton-shuttling systems for which the swinging-arm conduction was proposed are NADH or NADPH (74) and bacteriorhodopsin (22). Structures from umbrella sampling simulations in state 1 and state 3 also show an elongation of the CG-S1 distance by Ϸ 0.5 Å. On the other hand, a complete inversion of the Asp 15 side chain along the CB-CG axis in the simulations could not be observed as it is expected to induce major reorientation in the protein backbone. Electronic effects involved in reducing the barrier from 30 kcal/mol to 10 kcal/mol between models A and B can be excluded. For example, Mulliken charges on the atoms involved in PT from structures around the transition state show no significant differences (Ͻ0.05 e) between models A and B. The DFT profile from model A with R OD1-S1 ϭ 5.0 Å yields a forward barrier of 10.1 kcal/mol and a reverse barrier of 13.7 kcal/ mol, both of which are close to the barriers derived from the experimental rates. All other DFT scans performed in this work (especially those with short R OD1-S1 ) have their global minimum in state 1. These observations suggest that state 3 in FdI must have larger R OD1-S1 compared with the 7FDR x-ray structure. These general conclusions are supported by the umbrella sampling simulations with MMPT, which find ͳR OD1-S1 ʹ ϭ 4.8 Å in state 1, ͳR OD1-S1 ʹ ϭ 4.5 Å around the TS, and ͳR OD1-S1 ʹ ϭ 5.3 Å after formation of [3Fe-4S]H ϩ , all consistent with the DFT calculations.
It is worthwhile to mention that relaxed DFT scans for model B at the UB3LYP/6 -31G(d,p) level for water-unassisted proton transfer from the ϪCOOH group to the accepting sulfur atom of [3Fe-4S] 0 yield a barrier of ⌬E Ϸ 35 kcal/mol, which makes such a process considerably less likely. Compared with this, water-assisted proton transfer reduces this barrier to ⌬E Յ 10 kcal/mol and brings it within the range of experimentally observed rates (2).
In summary, generalized MMPT has been applied to PT in FdI, which involves transport of one proton from the solvent-exposed Asp 15 residue via a water molecule to a buried [3Fe-4S] cluster. The PMFs from MD simulations show that the hydronium ion (H 3 O ϩ ) is not a stable intermediate in the protein cavity. Therefore, the postulated water-mediated PT from Asp 15 to [3Fe-4S] 0 is most likely a concerted process. Barriers from the present PMF curves qualitatively agree with experimentally determined rate constants of the PT process in FdI. Quantitatively, the forward PMF barrier is in good agreement with experimental data. Our DFT scans using different model systems either over-or underestimated this barrier mainly due to the difficulty of finding a model structure that balances flexibility and protein backbone constraints that affect the conformational dynamics of Asp 15 (cf. swinging-arm mechanism). In addition, the DFT models include only the immediate protein environment without explicit solvent. From such models, no quantitative information can be expected. The present study establishes that it is possible to investigate PT reactions in condensed phase environments, including explicit solvation and at atomistic detail by using generalized MMPT. By studying PT in FdI, MMPT has proven to be a practical tool to simulate multistep PT processes in proteins.