Fluctuating hydrogen-bond networks govern anomalous electron transfer kinetics in a blue copper protein.

We combine experimental and computational methods to address the anomalous kinetics of long-range electron transfer (ET) in mutants of Pseudomonas aeruginosa azurin. ET rates and driving forces for wild type (WT) and three N47X mutants (X = L, S, and D) of Ru(2,2'-bipyridine)2 (imidazole)(His83) azurin are reported. An enhanced ET rate for the N47L mutant suggests either an increase of the donor-acceptor (DA) electronic coupling or a decrease in the reorganization energy for the reaction. The underlying atomistic features are investigated using a recently developed nonadiabatic molecular dynamics method to simulate ET in each of the azurin mutants, revealing unexpected aspects of DA electronic coupling. In particular, WT azurin and all studied mutants exhibit more DA compression during ET (>2 Å) than previously recognized. Moreover, it is found that DA compression involves an extended network of hydrogen bonds, the fluctuations of which gate the ET reaction, such that DA compression is facilitated by transiently rupturing hydrogen bonds. It is found that the N47L mutant intrinsically disrupts this hydrogen-bond network, enabling particularly facile DA compression. This work, which reveals the surprisingly fluctional nature of ET in azurin, suggests that hydrogen-bond networks can modulate the efficiency of long-range biological ET.

Electron transfer (ET) is of central importance in biological processes ranging from photosynthesis to respiration (1–6). Understanding the mechanisms by which biochemical systems facilitate efficient ET is essential for the elucidation of biochemical pathways and development of biomimetic catalysts.

Extensive investigations of Ru-modified derivatives of the blue copper protein azurin have shed light on the factors that control long-range biological ET (7, 8). Notably, the copper site in azurin is tuned for efficient ET with a low reaction activation energy (8, 9). Extensive mutagenesis studies have also been performed for this ET process, the results of which have been explained in terms of either the effect of individual residues on the Cu reduction potential (10–15) or the ET reorganization energy (16, 17). This previous work was consistently explained using the Marcus theory of ET subject to the assumption that the donor–acceptor (DA) electronic coupling for the ET reaction is unaffected by the mutations. Mutagenesis studies that account for changes in the DA coupling are less common and either assume static protein configurations or include only protein fluctuations in the reactant basin (18–20).

The current work explores mutations at the N47 site of Ru(2,2′-bipyridine)2 (imidazole)(His83) azurin (i.e., Ru-azurin). Both ET rates and Cu reduction potentials for wild type (WT) and three mutants (N47D, N47S, and N47L) are measured. The motivation for exploring residue 47 is that this asparagine participates in an extended hydrogen-bonding network that stabilizes the “rigid rack” for the Cu center in the protein (8–12, 21). A representative configuration of WT azurin in Fig. 1 shows the location of asparagine-47 between the two metal centers.

Fig. 1.

Representative configuration of WT azurin. ET occurs from the Cu center (large orange sphere) to the Ru center (large pink sphere). The Ru center is attached to the protein exterior at H83 (shown in cyan along with the Ru ligands); the Cu ligands are shown in red, and asparagine residue 47 is highlighted in green.

Unlike previous mutagenesis studies, the observed reactions cannot be explained using standard Marcus theory without accounting for the important role of protein fluctuations in the ET reaction mechanism. To understand these effects and to elucidate the protein motions that accompany ET in these azurin mutants, we use a recently developed nonadiabatic simulation technique, kinetically constrained ring polymer molecular dynamics (KC-RPMD). KC-RPMD is a Feynman path-integral method that enables the simulation of electronically nonadiabatic processes using classical trajectories; it accurately describes condensed-phase ET reaction dynamics and mechanisms across the normal, barrierless, and inverted regimes including the weak electronic-coupling regime associated with long-range ET (22, 23). The main benefit of the methodology is that it allows for the unbiased simulation of the ET dynamics, accounting for the full range of protein fluctuations that accompany the ET reaction.

Results and Discussion

Experimentally determined Cu reduction potentials and ET rates for WT azurin and the three mutants are shown in Table 1. Also reported is the driving force, , for each ET reaction, obtained from the measured Cu reduction potentials and the known Ru reduction potential, V vs. the normal hydrogen electrode (NHE) (8). The mutations at residue 47 significantly affect both the driving force and overall ET rate. Notably, while all three mutants show a decrease in the driving force for the reaction in comparison with WT, only N47D and N47S exhibit a decrease in the ET rate; the rate for N47L is enhanced relative to WT.

Table 1.

Experimentally determined copper reduction potentials (Cu II/I) vs. NHE, ET driving forces , and ET rates , for WT azurin and mutants

Mutant	E^0 (CuII/I)/V	ΔG^0/eV	k⋅s−1⋅10−5
WT	0.30(2)	−0.70(2)	11(2)
N47D	0.38(2)	−0.62(2)	5.2(7)
N47S	0.42(2)	−0.58(2)	2.5(2)
N47L	0.44(2)	−0.56(2)	21(4)

To disentangle the effect of the driving force on the ET reaction rate, the experimentally measured kinetics parameters are analyzed using the standard ET rate constant expression from Marcus theory (24–27),where is the DA electronic coupling between the donor and acceptor, is the reorganization energy, and is the driving force. Isolating the contribution from the driving force, we assume that all other parameters are unaffected by mutation, including and , such that the relative rate for a given mutant vs. WT iswhere eV (Table 1) and = 0.8 eV (28). Table 2 compares the ratio of ET rates measured from experiment, , to the ratio obtained using Eq. . For two of the mutations, N47D and N47S, the trend from this application of Marcus theory is consistent with the experimentally observed rates, with both mutations leading to lower rates. However, N47L shows clear enhancement of the observed ET rate compared with the prediction based solely on the change in driving force (Eq. ). This suggests that an additional aspect of protein motion plays an important role in the ET rates.

Table 2.

ET reaction rates from experiment, ; Marcus theory under the assumptions of Eq. 2, ; and KC-RPMD, , as well as the reorganization energy and average electronic coupling from KC-RPMD

Mutant	k/k^WT	k_MT/kMTWT	k_KC/kKCWT	λ_KC	⟨H_AB⟩_KC
WT	1.0	1.00	1.00	0.78	3.1
N47D	0.5(1)	0.76(8)	0.75(2)	0.78	3.1
N47S	0.23(5)	0.63(7)	0.88(4)	0.77	3.4
N47L	1.9(5)	0.56(7)	1.71(6)	0.78	5.3

Units for λKC and ⟨HAB⟩ are eV and 10−5 eV, respectively. Statistical error for these is smaller than the last reported digit.

To investigate the atomistic features governing the rate enhancement, KC-RPMD is used to simulate the nonadiabatic ET in both WT azurin and the investigated mutants. The KC-RPMD simulations use a fully atomistic representation of the metalloprotein, with over 15,000 atoms including explicit solvent. A two-state molecular mechanics force field, which is parameterized to fit WT experimental data, is used to describe the protein in the ET reactant and product diabatic states; the force field additionally describes the electronic coupling, , which explicitly depends on the distance between the Ru and Cu metal centers. All simulation details are provided in and .

Although KC-RPMD has been demonstrated to accurately describe ET across a broad range of electronic-coupling, driving-force, and solvent-coupling regimes (22, 23), the calculations converge with fewer trajectories when applied with a physically reasonable dividing surface (29). Throughout this work we use the “kink-pair” dividing surface in KC-RPMD (22, 23), which rapidly converges for ET reactions in the weak-coupling regime. This dividing surface corresponds to the ensemble of configurations for which the reactant and product electronic diabats are degenerate and weighted according to both the Boltzmann distribution and the magnitude of the DA electronic coupling, such that in the weak-coupling regimewhere is the probability of a nuclear configuration and and are the reactant and product potential energy surfaces, respectively. demonstrates that fewer than 8% of KC-RPMD trajectories undergo dynamical recrossing when initialized from the kink-pair dividing surface for WT azurin and for all mutants; this result confirms that the bottleneck for the ET reaction is well described by the kink-pair dividing surface, and it further indicates that the ET rate is accurately described by a Marcus-type transition-state rate theory of the formwhere is the average electronic coupling at the kink-pair dividing surface and is the outer-sphere reorganization energy calculated using KC-RPMD. A more detailed description of the rate-law analysis is provided in .

Table 2 gives the ET rates obtained using KC-RPMD, as well as the breakdown of the contributions to the rate in terms of the calculated outer-sphere reorganization energy and average electronic coupling. As anticipated, the outer-sphere reorganization energies computed from the KC-RPMD simulations confirm that this quantity is unaffected by the considered mutations. In contrast, while the WT, N47S, and N47D have almost identical values of the average electronic coupling, the value for N47L is markedly larger, fully accounting for the observed rate enhancement of N47L. Finally, we note that the relative KC-RPMD rates are in good agreement with experiment, suggesting that analysis of the contributions to the KC-RPMD rate will yield insight into the anomalous trend in the experimental kinetics upon mutation.

The effect of mutations on the electronic coupling can be seen from the distributions of DA distances calculated from KC-RPMD simulations in the reactant basin and at the dividing surface (Fig. 2). Surprisingly, WT and all three mutants show a large compression of the DA distance (> 2 Å) at the dividing surface compared with the reactant basin. In addition, N47L exhibits a stronger degree of compression in comparison with WT and the other mutants, displaying both a shifted peak and a fat tail at short DA distances for the dividing surface ensemble. The greater compression in N47L leads to an increased value of the average electronic coupling, , which in turn accounts for the anomalous N47L rate enhancement in Table 2.

Fig. 2.

Histogram of the Ru–Cu distance generated using KC-RPMD in the reactant basin and at the dividing surface for WT (red), N47D (purple), N47S (green), and N47L (blue). The dividing surface shows a clear compression of the Ru–Cu distance for WT and all mutants, with N47L showing the strongest degree of compression leading to an increased electronic coupling.

We now explore the detailed molecular rearrangements and interactions that lead to the DA compression trends in Fig. 2. For WT azurin, the color map in Fig. 3 shows the degree to which various residues are coupled to this compression during the ET reaction. Specifically, the crystal structure of the WT is depicted with the color and size of each atom scaled according to , where is the average distance of atom to the Cu center in the reactant ensemble and is the corresponding average from the dividing surface ensemble; blue indicates that the atom is farther from the Cu center in the dividing surface compared with the reactant basin, and red indicates that the atom is closer. Most of the protein is white, indicating that these regions do not compress toward the Cu center at the dividing surface. The main colored portion of the protein is the Ru moiety, which behaves as a rigid body bending toward the Cu center at the dividing surface. However, residue N47, highlighted in Inset, is observed to move away from the Cu center during the ET event, indicating that this residue is coupled in some way to the compression process.

Fig. 3.

The WT crystal structure with Cu (orange) and Ru (pink) indicated. The size and color of the remaining atoms are scaled according to the difference in distance from the atom to the Cu center in the reactant basin and dividing surface; blue/red indicates farther from/closer to Cu at the dividing surface than in the reactant basin. Inset expands the N47 region.

Fig. 4 compares representative configurations of the region around N47 obtained from the WT reactant basin and dividing surface ensembles. Fig. 4 shows that in the reactant basin, the N47 in WT forms a hydrogen bond to a nearby threonine (T113) and to a water molecule found in the pocket between the Ru moiety and N47; the water molecule additionally forms hydrogen bonds with backbone oxygens in the pocket. Fig. 4 illustrates that very different configurations are adopted in the dividing surface ensemble for WT. H83 compresses into the pocket, displacing the water and disrupting the hydrogen-bonding network around N47; four of the five reactant-basin hydrogen bonds are broken to allow for the motion of the Ru moiety into the pocket. As such, compression in WT during the ET process appears to be gated by a substantial fluctuation in which the local hydrogen-bonding network of the N47 residue is disrupted.

Fig. 4.

Representative WT configurations from KC-RPMD simulations illustrating the hydrogen-bonding network in the pocket surrounding residue 47. (A) The reactant basin exhibits a hydrogen-bond network between N47, the neighboring T113, a water molecule present in the pocket between N47 and H83, and nearby backbone oxygens. (B) The dividing surface is characterized by changes that include (1) breaking of the hydrogen bonds around N47, (2) displacement of the water molecule from the pocket, and (3) compression of H83 into the pocket. Hydrogen bonds are indicated by orange-dashed lines, the Cu center is in orange, residue 47 is highlighted in green, and the Ru center and H83 are highlighted in pink.

To quantify the observations from Fig. 4, and to compare the WT behavior during compression with that of the N47L mutant, Fig. 5 plots the 2D histogram of the Ru–Cu distance and the hydrogen-bonding distance between residue 47 and T113 in the reactant basin ensemble (red) and in the dividing surface ensemble (blue), for both WT (Fig. 5) and N47L (Fig. 5). The distance between residue 47 and T113 (indicated in Fig. 5 , Insets) reports on fluctuations in the hydrogen-bonding network and is defined as the distance between the backbone nitrogen on T113 and the carbonyl oxygen on N47 for WT and as the distance between the backbone nitrogen on T113 and one of the methyl carbons on L47 for N47L.

Fig. 5.

Two-dimensional histogram of the Ru–Cu distance and the hydrogen-bonding distance between residue 47 and the nearby T113 for (A) WT and (B) N47L. Reactant basin and dividing surface configurations are indicated by red and blue, respectively. The dashed lines indicate the average hydrogen-bond distance for the WT in the reactant basin. Insets indicate characteristic configurations for the two cases, with the hydrogen-bond distance indicated.

Fig. 5 shows strong correlation between the Ru–Cu distance and the N47–T113 distance in WT, confirming that the hydrogen bond in the reactant basin is broken during the compression that gates ET. In contrast, Fig. 5 shows that the corresponding distance in the N47 L mutant is uncorrelated with DA compression, reflecting that point mutation to a leucine residue intrinsically disrupts the hydrogen-bonding network around residue 47. These results emphasize that the N47L mutant is “primed” for an enhanced ET reaction rate relative to WT, since facile compression of the DA distance can be performed without the restraints of the hydrogen-bond network in the vicinity of residue 47. We further note that these observations are completely consistent with the reported DA distance distributions in Fig. 2; the lack of hydrogen-bonding network in the N47 L mutant leads to more facile compression in the dividing surface ensemble, which in turn leads to stronger electronic couplings and faster ET rates.

As may be anticipated from the relative similarity of the DA distance distributions in Fig. 2 for the WT and the N47S and N47D mutations, reveals that the hydrogen-bond fluctuations in the N47S and N47D mutations during ET are more similar to the WT results in Fig. 5 than to the N47L mutant results in Fig. 5. Specifically, like the WT, N47S exhibits an intact hydrogen bond in the reactant basin that becomes disrupted during the DA compression that gates the ET reaction. Interestingly, the D47–T113 hydrogen bond in the N47D mutant appears to be sufficiently strong that it remains preserved, even in the compressed configurations at the dividing surface.

Taken together, Figs. 4 and 5 reveal the strikingly fluctional nature of the hydrogen-bonding network that gates the ET reaction in WT azurin, and it provides a molecular basis for understanding the qualitatively different DA distance distributions in the N47L mutant that lead to anomalously fast ET kinetics. The key features of the ET conformational gating in WT azurin are illustrated in Fig. 6, which depicts the distribution of Ru–Cu distances [, black curve] and electronic coupling [, green curve]. The red portion of the distribution indicates the thermally accessible configurations in the reactant basin, which exhibit an extended Ru–Cu distance, a small electronic coupling, and a well-formed hydrogen-bonding network involving residue N47. Given the strong distance dependence of the electronic coupling on DA distance, the ET process is most favorable from configurations in which the system fluctuates to compressed DA distances, indicated by the blue region, which involves disrupting the hydrogen-bond network in the WT azurin. Because the N47L mutant intrinsically disrupts this hydrogen-bond network, DA compression is more facile, leading to the accessibility of configurations with substantially shorter DA distances (Fig. 2) and thus faster ET reaction rates. We note that although this mechanistic explanation was made possible with all-atom nonadiabatic simulations, the basic interpretation of the experimentally observed trends does not hinge on the details of the force field or the form of the electronic coupling, beyond the robust assumption that the electronic coupling is strongly dependent upon the DA distance.

Fig. 6.

Schematic illustration of the atomistic features governing ET in WT azurin. The black curve corresponds to the probability distribution, , of Ru–Cu distances, , while the green-dashed curve corresponds to the distance-dependent electronic coupling, . Reactant basin and dividing surface configurations are indicated in red and blue, respectively.

Finally, we show that the experimentally observed trends for the ET rates are poorly explained without explicit inclusion of the DA compression that gates the ET reaction, even when a more sophisticated description of the electronic coupling is used. Previous theoretical work illustrated the dynamical nature of azurin, but modeled ET on the assumption that the reaction proceeds from the ensemble of configurations in the reactant basin (i.e., without inclusion of the DA compression discussed here) (30–33). Table 3 presents the relative ET rates, , obtained using Eq. but with the average value of the electronic coupling calculated in the reactant basin instead of at the dividing surface. The results calculated in this way fail to capture the anomalous enhancement of the N47L rate. When the electronic coupling is calculated using the Pathways model (34–36), which is more sophisticated than the simple exponential form of the coupling (Eq. ) that is otherwise used in this study, the estimate of the N47L rate is in even greater disagreement with experiment. These results highlight that capturing the experimentally observed trends in azurin demands explicit inclusion of DA compression during the ET reaction, as in the KC-RPMD simulations reported here.

Table 3.

Comparison of the relative experimental rates, , to the those calculated while neglecting DA compression

Mutant	k/k^WT	k_A/kAWT	k_Path/kPathWT
WT	1.0	1.00	1.00
N47D	0.5(1)	0.64(1)	0.73(4)
N47S	0.23(5)	0.83(2)	0.75(5)
N47L	1.9(5)	0.93(2)	0.57(3)

uses Eq. for the electronic coupling, and uses the Pathways model.

Concluding Remarks

We present a combined experimental and computational study to elucidate the kinetics and conformational fluctuations associated with ET in Ru-modified Pseudomonas aeruginosa azurin. Experimental measurements of the ET rate and driving force of the WT Ru-azurin and three mutants (N47S, N47D, and N47L) show significant enhancement of the ET rate for N47L, which is inconsistent with behavior that is expected solely on the basis of the changes in the ET driving force among the mutants.

To understand this anomalous rate enhancement, the recently developed KC-RPMD method is used to directly simulate the nonadiabatic dynamics of the ET reaction in azurin and its mutants. Calculation of the ET rate constants for azurin and its mutants (Table 2) reveals good agreement with experiment, allowing for the analysis of the KC-RPMD trajectories to understand the anomalous kinetics of N47L. Analysis of the KC-RPMD trajectories reveals that, surprisingly, WT and all mutants exhibit a strong compression of the DA distance between the Cu (electron donor) and Ru (electron acceptor) species that gates the ET reaction (Fig. 2). The N47L mutant shows a particularly large degree of compression in comparison with WT and the other mutants that leads to stronger DA electronic coupling and the observed rate enhancement for N47L.

Analysis of the molecular fluctuations that accompany the ET reaction in the WT reveals an extended hydrogen-bond network that is disrupted during compression of the DA distance (Figs. 4 and 5). The nonpolar character of the leucine mutation in N47L intrinsically disrupts this hydrogen-bond network even in the reactant basin, giving rise to a more facile DA compression which leads to the higher electronic coupling and anomalously fast N47L ET rate.

This work reveals unexpected features of ET in protein systems, even for the extensively studied case of WT azurin, and it illustrates the importance of methods that naturally describe the fluctuations that accompany ET reactions. We emphasize that the results presented here do not indicate a breakdown in the Marcus theory for electron transfer—they simply demonstrate the importance of explicitly including the DA compression motions that gate the ET reaction, which may be accompanied by nontrivial fluctuations in extended hydrogen-bond networks. These effects are rigorously and conveniently captured via KC-RPMD, without the need for computationally costly multidimensional free-energy profile calculations. We expect KC-RPMD to prove useful in future studies of nonadiabatic chemical reactions in other protein systems, particularly those for which conformational fluctuations are thought to play an even more important role than in azurin, such as cytochromes (7, 18, 37).

Materials and Methods

Experimental Details.

Mutant azurin proteins were expressed, purified, and modified with [Ru(2,2′-bipyridine)2]2+, using literature protocols (38, 39). The CuII/I reduction potentials were determined using differential pulse voltammetry experiments with a standard three-electrode electrochemical cell. Reported potentials are referenced to the NHE. Transient absorption (TA) experiments were conducted in the Beckman Institute Laser Resource Center at California Institute of Technology. Reactions were initiated using 500 nm pump light and monitored using a white-light probe or 632.8 nm light from a HeNe continuous-wave laser. Established flash–quench conditions were used in all TA experiments (40). Full details of the experimental procedures can be found in .

Calculation Details.

All KC-RPMD simulations were performed at 300 K using a modified version of the Gromacs-5.0 molecular dynamics package (41). The potential energy surfaces used to define the ET reactant and product states are based on the GROMOS 53a6 force field (42) with additional terms to describe the interactions of the metal centers and the ET driving force. These additions are parameterized on the basis of WT experimental data; no additional potential energy fitting was performed for the mutants. Full details are provided in .

The KC-RPMD equations of motion used to simulate the nonadiabatic ET dynamics are (22, 23)which correspond to the classical limit for the nuclei. The position, velocity, and mass of nuclear coordinate are indicated by , , and , respectively, and the vector of nuclear positions is . The auxiliary electronic variable in KC-RPMD reports on nonadiabatic transitions between electronic states and is described by position , velocity , and mass . It is coupled to a Langevin bath with friction coefficient ; is a normal Gaussian random variable. Additional KC-RPMD details are provided in .

The electronic coupling, , is modeled using an exponential form that depends on the distance between the Ru and Cu metal centers, ,using the experimentally determined values of cm−1, d0=3 Å, and Å−1 for all mutants (8, 28).

Ensemble averages reported for the reactant basin and dividing surface are obtained from KC-RPMD sampling trajectories that are run for at least 20 ns and 4 ns with a time step of 1 fs and 0.5 fs, respectively. The reactant basin corresponds to configurations for which the auxiliary electronic variable is constrained to , while the kink-pair dividing surfaces configurations are constrained to (22, 23).