Studies reveal the mechanisms behind tumor metastasis and how to stymie it. But primary tumors still get the lion’s share of researchers’ attention.
Image courtesy of Science Source/Phanie.
New Research In
Physical Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
Edited by J. Andrew McCammon, University of California at San Diego, La Jolla, CA, and approved February 24, 2012 (received for review December 22, 2011)
Abstract
[NiFe]-hydrogenases are fascinating biological catalysts with potential application in biofuel cells. However, a severe problem in practical application is the strong sensitivity of hydrogenase to gaseous inhibitor molecules such as CO and O2. Recently, a number of successful protein engineering studies have been reported that aimed at lowering the access of diatomic inhibitors to the active site pocket, but the molecular mechanism conferring increased resistance remained unclear. Here we use a multiscale simulation approach combining molecular dynamics with a master equation formalism to explain the steady drop in CO diffusion rate observed for the mutants V74M L122A, V74M L122M, and V74M of Desulfovibrio fructosovorans [NiFe]-hydrogenase. We find that diffusion in these variants is controlled by two gates, one between residues 74 and 476 and the other between residues 74 and 122. The existence of two control points in different locations explains why the reduction in the experimental diffusion rate does not simply correlate with the width of the main gas channel. We also find that in the more effective mutation (V74M) CO molecules are still able to reach the active site through transitions that are gated by the microsecond dihedral motions of the side chain of R476 and the thermal fluctuations of the width of the gas channel defined by M74 and L122. Reflecting on the molecular information gained from simulation, we discuss future mutation experiments that could further lower the diffusion rates of small ligands inhibiting [NiFe]-hydrogenase.
There is currently a strong research effort in developing biology-inspired approaches for the production of renewable energy sources. Examples include the use of natural or artificial photosystems for light harvesting and exciton generation (1, 2), hydrogenases for hydrogen production (3, 4) and oxidation in biofuel cells (5, 6), and carbon-monoxide dehydrogenase for CO2 reduction and synthesis of precursors of biofuels (7, 8). Yet, a general problem associated with the use of biocatalysts is their sensitivity to changing environmental conditions. A good example is the inhibition of hydrogenases by diatomic molecules. Naturally evolved under anaerobic conditions, the enzyme becomes inhibited by molecular oxygen that is present in the atmosphere and only few species can resist oxidative damage. Thus, an important aspect of biotechnological approaches is the ability to design and engineer enzyme variants that are more robust and durable than their natural counterparts.
In the field of hydrogenase research, three strategies have been proposed to reduce the effect of inhibition or damage (5): (i) restricting the access of gas molecules to the active site, (ii) lowering the binding affinity of the gas molecules to the metal centers, and (iii) facilitating the removal of undesired oxidation products. Recently, evidence was found that oxygen-tolerance in the membrane-bound [NiFe]-hydrogenase from the Knallgas bacterium Ralstonia eutropha is conferred by the presence of two additional cysteine residues facilitating four-electron reduction of bound oxygen to water (strategy iii) (9, 10). However, these variants are typically less active than the oxygen sensitive enzymes.
Important progress toward aim (i) has been recently reported by Leger and coworkers (11⇓⇓–14). In a series of elegant mutation experiments, the authors showed using protein-film voltammetry that the diffusion rate of inhibitors of [NiFe]-hydrogenase, CO and O2, could be reduced by up to four orders of magnitude relative to the WT enzyme, whereas the turnover rate of substrate (H2) did not decrease by more than a factor of four (14). All effective mutations involved a change of V74 (V), which forms together with L122 (L) a structural motif (denoted 74–122 motif in the following) that is speculated to act as a control point for access of gas molecules to the active site (see Fig. 1 A and B). Particularly effective were mutations where V74 was changed into residues with longer side chains, such as glutamine (Q), methionine (M), and glutamate (E), or bulkier residues such as tryptophan (W). Increase of the side chain by one CH2 unit from V74D to V74E or V74D to V74N led to an approximately 10-fold decrease in diffusion rate.
Fig. 1.
(A) WT [NiFe]-hydrogenase from D. fructosovorans [Protein Data Bank (PDB) 1YQW] and coarse-grained clusters used to describe the kinetics of CO diffusion. The clusters displayed as colored spheres were identified in previous simulations of the enzyme (33). Red spheres indicate the pathway relevant for CO molecules (pathway 1) (34), blue and yellow spheres indicate additional pathways for H2 and O2 diffusion (33). The clusters are labeled by numerals and by lowercase characters as reported in the crystal structure analysis of ref. 25. The upper part of the gas channel and the active site pocket are encircled and shown enlarged in two different perspectives in B and C. (B) V74 and L122 form the 74–122 motif. (C) V74 and arginine 476 (R476) form the 74–476 motif. (D–F) Active site pockets of mutant enzymes V74M (PDB ID 3H3X) (D), V74M L122M (MM, PDB ID 3CUR) (E), and V74M L122A (MA, snapshot from MD simulation) (F), respectively. The color code for atoms is C, green; N, blue; O,red; S, yellow in all panels.
A first explanation of these observations is that the mutations narrow the 74–122 motif and thereby hinder the diffusion of the large diatomics (CO, O2) toward the active site. However, there are a few observations that cannot be explained in this way. For instance, there is an absence of correlation between diffusion rate and width of the 74–122 motif for WT, V74M L122A (MA), V74M L122M (MM), and V74M mutant enzymes. Gas diffusion in the MA mutant is slowed by a factor of 42 relative to the WT enzyme, even though the gas channel diameter is not significantly reduced in this mutant (15). Another puzzling observation is that the MM double mutant is less effective than the V74M single mutant even though the gas channel diameter is similar to the one for the single mutant according to crystal structure (12, 13). A fact that is often overlooked is that the side chains can exhibit large thermal fluctuations at room temperature, which render interpretations based on static crystal structures questionable.
Adequately validated computer simulations based on molecular dynamics are a powerful tool to probe the effect of thermal fluctuations on the transport of molecules in proteins. Previous studies have characterized gas diffusion in myoglobin (16⇓⇓–19), (reviewed in ref. 20), lipoxygenase (21), flavin-containing monooxygenase (22) and oxidases (22, 23) (reviewed in ref. 24), and hydrogenase (25⇓⇓–28). Moreover, the diffusion mechanism of larger ligands in acetylcholinesterase has been investigated using molecular dynamics (MD) simulation (29) in combination with a gating model (29⇓–31) (see also ref. 32). These studies have provided valuable theoretical insights, but quantitative kinetic data that would allow for direct experimental validation were only very rarely reported. To fill this gap, we have devised in recent work a multiscale methodology for the calculation of gas diffusion rates in proteins (33, 34). In this approach, gas transport is described by a Markov process, in which gas molecules “jump” between sites (or cavities) in a protein. Using molecular dynamics simulation, a master equation is constructed that describes the time-dependent gas population of the sites. A fit of the populations to a phenomenological rate equation gives rate constants that can be compared to experimental measurements.
A recent application of this methodology to WT [NiFe]-hydrogenase gave diffusion rates for CO and O2 that were in good agreement with experimental results (12, 14). We also found that CO molecules have a strong preference for diffusion through the main tunnel of the enzyme as identified by Xe-binding X-ray crystallography (25), whereas O2 and particularly H2 can access the active site by multiple pathways. Diffusion of the latter two gases was found to exhibit similar characteristics to O2 diffusion in flavin-containing monooxygenase and oxidase as reported previously by Baron et al. (22). Their study indicated that flavoenzymes employ multiple funnel-shaped diffusion pathways that guide O2 molecules from the solvent to the enzyme active site. Another striking similarity was that diffusion from the solvent to the active site occurs in a stepwise fashion, a behavior that we used for the definition of the states in our kinetic model.
In the present investigation, we employ the combined MD/master equation approach to obtain a rationale for the impact that some mutants have on gas diffusion in [NiFe]-hydrogenase (14). For this purpose, we have chosen the three aforementioned mutant enzymes, MA, MM, and V74M (Fig. 1 D–F), which exhibit a small, medium, and large effect on CO diffusion, with a diffusion rate 1.5, 2.5, and 3.5 orders of magnitude lower than in the WT enzyme. In contrast to previous interpretations, we find that gas diffusion in enzymes that contain the V74M motif is controlled by two dynamical bottlenecks, one between V74M and R476 (Fig. 1C), hereafter referred to as 74–476 motif and the other between V74M and L122X (X = A,M,L, 74–122 motif, Fig. 1B). In the MA mutant, only the 74–476 passage provides a significant barrier to CO diffusion, whereas in the MM mutant and increasingly in the V74M mutant, the 74–122 passage gives rise to a second barrier of similar height. Yet we show that, even in the most effective mutation, the two motifs are subject to CO leakage due to the microsecond protein motions that temporarily widen the two passages. The molecular picture obtained is supported by the calculation of CO diffusion rates for the three mutant enzymes, which compare very favorably with experimental measurements.
Results
Thermal Fluctuations of the Gas Channel.
A key structural feature for gas diffusion in [NiFe]-hydrogenase is the width of the main gas channel at its narrowest point. We have taken the distance between the two terminal C atoms of the side chains of the residues at positions 74 and 122 to define the width, although the qualitative trends are not very sensitive to this particular choice. The probability distributions of this distance obtained from 10 ns MD simulation at 300 K are shown in Fig. 2 for WT and the three mutant enzymes. The distributions are rather broad (rmsd = 0.7–0.8 Å), approximately symmetric for WT and MA, and asymmetric for V74M and MM. As expected, the distribution for the WT enzyme peaks at the largest distance (7.5 Å). Yet the distribution for the MA mutant is almost superimposable on the one for WT. Going from the MA to the MM mutant the peak position shifts significantly to smaller values (4.7 Å), and a shoulder appears for MM at 6.5 Å. The peak position for V74M is only sightly smaller than for the MM mutant (4.3 Å).
Fig. 2.
Distribution of the width of the gas channel in WT and mutant [NiFe]-hydrogenases obtained from molecular dynamics simulation at 300 K. For WT the width was defined by the distance between the first γ-carbon atom of V74 and the first δ-carbon atom of L122, for V74M by the distance between the ε-carbon atom of M74 and the first δ-carbon atom of L122, for MM by the ε-carbon atom of M74 and the other ε-carbon atom of M122, and for MA by the distance between the ε-carbon atom of M74 and the β-carbon atom of A122. The corresponding distances in the crystal structures (PDB ID 1YQW, 3H3X, and 3CUR for WT, V74M, and MM, respectively) are shown as spikes. Notice that for MM two distances are given corresponding to the two conformations of M122 in the crystal structure.
In Fig. 2, we have also indicated the corresponding distances in the crystal structures by vertical dashed lines. Each computed distribution includes the crystal structure distance, confirming that the equilibrium simulations preserve the structural features found in the crystal structures while enabling the sampling of thermal fluctuations. The same is true for the fluctuations of the secondary structure rmsd, which are 1.48 and 1.67 Å relative to the crystal structure for V74M and MM mutants (Fig. S1). The peak positions for MA and MM are only slightly shifted to larger distances relative to the crystal structure, indicating that the narrowing of the 74–122 motif in these mutants is stable against thermal fluctuations, but only in a time-averaged sense. The shoulder observed for MM is due to rapid (nanosecond) exchange between two side-chain conformations, which is consistent with the report of two conformers in the crystal structure (dashed lines in red).
We notice that there is no correlation between the center of the distributions and the steady drop in experimental diffusion rates that are summarized in Table 1. Evidently, the 42-fold decrease in diffusion rate observed for the MA mutant is not due to a narrowing of the 74–122 passage. The significant drop in diffusion rate observed for the V74M and MM mutants may be interpreted in this way, though it is difficult to rationalize the large difference in diffusion rate between V74M and MM mutants by the distributions displayed in Fig. 2. It seems that the observed mutation effects cannot be simply explained by the width of the 74–122 motif, implying that other structural effects or the different physiochemical nature of the mutations are likely to play a role in the lowering of the diffusion rate.
View this table:
Table 1.
Computed (comp) rates for CO diffusion from the solvent to the [NiFe]-hydrogenase active site cluster G, k+1, and for the reverse direction, k-1, at 300 K, and the experimentally (exp) determined on-rate for CO binding, kin
Transition and Diffusion Rates.
The methodology for the calculation of diffusion rates is described in Materials and Methods. For the mutant enzymes we use the coarse-grained states previously defined for the WT enzyme and depicted as spheres in Fig. 1. The rate matrix is the same as for WT (34) except for transitions that are critically affected by the mutations. These are the transitions that occur in the vicinity of the mutations, one crossing the 74–122 motif, denoted E←68, and the other crossing the 74–476 motif and ending in the active site pocket denoted G←E. The mean first passage times (MFPTs) and the corresponding transition rates for these transitions are computed using nonequilibrium pulling simulations. To make the computations feasible, we have used a reduced system where only residues within 15 Å from the active site cluster G are allowed to move freely and the rest of the system is frozen (see SI Text for further details). Thus, only the motions of the residues in the wider vicinity of the mutation sites are taken into account, which should be a sufficiently good approximation because the mutation effects are expected to be local. Indeed, comparing the MFPT obtained for the reduced system with the one obtained for the full system, where all parts of the protein are allowed to move freely, we find only minor differences.
The results of the pulling simulations for the reduced model are shown in Fig. 3. In Fig. 3A, we show the MFPT as a function of pulling force for the E←68 transition, τE68(F). The data points fit the Dudko–Hummer–Szabo (DHS) model for force-dependent kinetics (35) very well (inverse of Eq. 5), with correlation coefficients of 0.98, 0.97, and 0.96 for V74M, MM, and MA mutant enzymes, respectively. Extrapolation of the fits to zero force gives τE68 = τE68(0). For the MA mutant, we obtain only a small increase to 110 ns compared to 30 ns for the WT enzyme. By contrast, τE68 for MM and V74M are significantly longer than for WT, 2.4 and 15 μs, respectively, in line with the smaller channel diameter in these mutants.
Fig. 3.
MFPTs for transition of CO molecules between clusters as obtained from constant-force pulling simulation. The MFPT is shown as a function of pulling force (A) for the transition E←68 and (B) for the transition G←E. See Fig. 1 for the location of the clusters within the protein. Statistical error bars for each pulling simulation are indicated, but may be smaller than the symbol size. Fits to Eq. 5 are shown in solid lines, data for WT are taken from ref. 34.
In Fig. 3B, we show the computed MFPT for the second transition affected by the mutations, G←E. The data for the V74M mutant fit again the DHS model very well with a correlation coefficient of 0.98. Compared to the WT enzyme the MFPT increases more sharply as the pulling force is lowered, giving a zero force MFPT τGE = 3.5 μs compared to 16 ns for WT. Evidently, the mutation from valine to methionine at position 74 affects the kinetics of this transition almost as much as the one for the transition E←68. We have not carried out pulling simulations for G←E for the double mutants MA and MM because MD simulations showed that the structure and dynamics of the 74–476 motif is hardly affected by the additional mutation at position 122. Thus, we assume here that the G←E transition rate in MA and MM mutants are the same as in V74M.
Converting the MFPTs to transition rates, kij = 1/τij, inserting the rates into the rate matrix K and solving Eq. 4 for initial conditions where all gas molecules are in the solvent, we obtain the time-dependent population of the active site cluster G, pG(t). The results are shown in Fig. 4. As for WT enzyme, the kinetics for all three mutants is monoexponential fitting well the phenomenological rate equation, Eq. 7. The mutations clearly impact the onset of the populations but do not change the equilibrium probability, pG(t → ∞). The rate constants for diffusion from the solvent to cluster G, k+1, and for the reverse direction, k-1, were extracted from the fit and are summarized in Table 1. For CO, k+1 can be directly compared with the experimentally determined on-rate kin (see refs. 14 and 34), which is also given in Table 1. We find that the present computations predict the correct order of magnitude for k+1 for all three mutant enzymes, and hence the correct order of the mutation effect, MA < MM < V74M.
Fig. 4.
Diffusion kinetics for CO in WT and mutant [NiFe]-hydrogenases. Data points are obtained by solving the master equation (Eq. 4) for the population of the active site cluster G, pG. At t = 0, the probability for the cluster S (solvent) was set equal to one, and the probabilities for all other clusters were set equal to zero. Fits to the phenomenological first-order equation (Eq. 7) are shown in lines. Data for WT are taken from ref. 34.
Free Energy Landscape.
The data obtained from nonequilibrium pulling can be used to construct a free energy profile along a typical diffusion path of the gas molecule (Fig. 5). The gas molecule starts in the solvent (S) and diffuses via clusters of pathway 1 (18,15,14,68) to the 74–122 motif and further on to the 74–476 motif to reach the active site (G). The MFPT for each transition from cluster i to j was converted into effective activation free energies (equal to barrier heights) using transition state theory, 
. The corresponding free energy differences were obtained from the ratio of MFPTs for forward and back transitions, ΔAji = -kBT ln(kji/kij) = -kBT ln(τij/τji). Computed free energies and time constants are indicated in Fig. 5 and summarized in Tables S1 and S2.
Fig. 5.
Free energy landscape along a typical diffusion path in WT [NiFe]-hydrogenase and mutant enzymes. The states or clusters corresponding to the minima of the landscape are denoted at the bottom of the figure. Their location within the protein is shown in Fig. 1. Computed MFPTs between states i and j (τji) are indicated on top of the barriers for the forward direction. MFPTs for the reverse direction are on the same order of magnitude except for the first step, for which the MFPT for the reverse direction is three orders of magnitude higher than for the forward direction. Time constants for CO diffusion from the solvent to the active site (τ+1) and for the reverse direction (τ-1) are indicated at the top of the figure.
We find that only the first transition from the solvent to a protein cluster is uphill in free energy (ΔA = 4.0 kcal/mol at 1 mM CO), which can be explained by the smaller volume of a given protein cavity relative to the volume per molecule in solution (entropic penalty). All other transitions are approximately thermoneutral. The time it takes for a gas molecule to enter a protein cluster (0.1–1 μs) is about three orders of magnitude longer than the time for subsequent transitions up to cluster 68 at the entrance to the 74–122 motif. Here mutation of V74 to M and L122 to A (MA) leads to a minor increase for the E←68 transition but to a strong increase for the G←E transition from 16 ns for WT to 3.5 μs. Consequently, the total diffusion time (τ+1) increases significantly from 94 μs for WT to 1.7 ms for MA. Thus, it is the transition across the 74–476 motif that is the main bottleneck to diffusion in the MA mutant, which explains why in experiment the diffusion rate is 42-fold lower than for WT (14), even though the width of the proposed “gate” (74–122 motif) remains virtually unchanged (see Fig. 2).
Changing the alanine residue at position 122 to methionine (MM mutant), the transition time across the 74–122 motif increases sharply from 0.11 to 2.4 μs. The latter is close to τGE, suggesting that the sevenfold reduction in the diffusion rate relative to MA, is due to the buildup of a second barrier across the 74–122 motif, which is of similar height as the one across the 74–476 motif. This finding is in line with the drop in the channel diameter distribution from MA to MM illustrated in Fig. 2. Mutating the methionine residue at position 122 back to the original leucine residue (V74M mutant), the transition time further increases from 2.4 to 15 μs. Thus, the 12-fold reduction in diffusion rate relative to MM is due to a “tightening” of the 74–122 motif, as indicated in Fig. 2 by the slight shift of the distributions to shorter distances and the disappearance of the shoulder that is observed for MM.
Discussion
Molecular Descriptors for Gas Diffusion.
Experiments and the present calculations show that the mutants exhibit a significantly smaller diffusion rate for the inhibitor CO than for the WT enzyme. However, even in the most effective mutation studied here (V74M), one CO molecule will reach the active site approximately every 50 ms, or every 25 turnover steps (assuming a turnover rate of 460 s-1 at 1 atm H2, ref. 14, and a solution 1 mM in CO). A particularly interesting question and the starting point for any “rational design” strategy concerns the nature of the protein motions that allow gas molecules to cross or circumvent barrier(s) to gas diffusion. More specifically, what are the protein motions that facilitate the crossing of the 74–122 and the 74–476 motifs, and can these motions be described by a suitable parameter, a diffusion coordinate?
The slowest transitions between protein cavities can be considered as rare events. The gas molecule diffuses within the initial cluster and will have to wait until the protein residues lining the transition path adopt a favorable position so that the gas molecule can cross to the final state cluster (“diffusive jumps”). Clearly, a good coordinate describing this process would have distinct values for the initial state (when the gas molecule cannot cross) and for the transition state (when the molecule is in the process of crossing). Thus, one would seek a coordinate σ for which the thermal distributions in the initial (IS) and transition state (TS), pIS(σ) and pTS(σ), respectively, have as little overlap as possible. 
[1]Here we choose σ as a linear combination of collective variables χn (e.g., distances, dihedrals) 
[2]and optimize the coefficients cn so as to minimize the integral in Eq. 1. The transition state distribution is obtained from pulling simulation at the lowest force bias, using only the part of the reactive trajectory where the molecule transits between initial and final state cluster (see SI Text for further details). The initial state distribution is obtained from equilibrium MD simulations.
We discuss at first the slowest transition in the mutant most effective in lowering CO diffusion, E←68 in V74M. Considering different linear combinations of distances and side chain dihedrals, we find that the best descriptor for this transition (the one with lowest overlap) is simply the channel width defined in Fig. 2. The distribution for the transition state (Fig. S2A, dashed line) peaks at about 7 Å, which is very similar to the peak position for WT in Fig. 2, showing that the crossing of the 74–122 motif is most likely in configurations where the channel width in the V74M mutant is similar to the average width in the WT enzyme. Such events, induced by thermal fluctuations of the protein backbone, are very rare and hence the low diffusion rate in this mutant. We have also further analyzed the nature of the interactions that cause the mutation effect in the V74M mutant. We found that electrostatic effects are of minor importance due to the small dipole moment of CO (see SI Text for details), implying that the mutation effect is a consequence of increased steric interactions.
The situation is more complex for the MM mutant. In contrast to V74M, the IS and TS distributions for MM exhibit a strong overlap implying that CO molecules can transit from cluster 68 to E even if the channel at the 74–122 motif is so narrow that no molecule would fit through (peak at 4 Å in Fig. S2B). Analyzing the pulling trajectories at lowest force bias, we observe that CO molecules can indeed circumvent the 74–122 motif. In one reactive path, the hydrogen network formed around a nearby water molecule rearranges and the CO molecule diffuses sideways via Asp541 and Pro542 to cluster E. In another path the backbone atoms of Glu285 depart from its original position and move away from the two methionines, thereby facilitating the transit of a CO molecule. We also observe that the conformational flexibility of the two methionine residues facilitates a direct transit through the 74–122 motif. If the second dihedral angle of either M74 or M122 is -60°, the channel width can increase to about 7 Å and CO molecules can cross the motif, which explains the second peak of the transition state distribution in Fig. S2B. In summary, the methionine residue at position 122 allows for a larger degree of fluctuations of the 74–122 motif than the Leu122 residue in the V74M mutant. Consequently, the diffusion rate for MM is higher than for V74M although the time-averaged channel width is similar. For a discussion of the results obtained for the G←E transition, we refer to the SI Text.
Comparison with Conformational Gating Mechanism in Actelycholinesterase (ACE).
In previous work Zhou et al. proposed the existence of a gating mechanism for ligand access to the active site of this enzyme (29). Rapid reorientation of five aromatic side chains led to an opening and closing of the passage to the active site, merely lowering the diffusion rate of the natural substrate (acetylcholine) by a factor of two, but significantly reducing the diffusion of bulkier ligands. This finding led to the conclusion that the gate formed by the side chains allowed for enzyme specificity without sacrificing efficiency. There is some similarity between ACE and hydrogenase in that active site access is controlled by the side-chain motions of a few specific residues. However, there is an important difference. In ACE, the gate was reported to open and close on the picosecond timescale (29), whereas in hydrogenase mutants we did not observe a single “open” conformation for CO within 10 ns of equilibrium simulation, which is in accord with the computed mean first passage times which are on the microsecond timescale except for WT (nanoseconds). Thus, although gating in ACE is fast compared to the characteristic gas diffusion time of the substrate (≪ τD = R2/D = 10 ns, R the enzyme radius and D the diffusion constant of the ligand) (29, 30), the relevant side-chain motions in hydrogenase mutants are likely to be at the other, slow extreme (≫τD = 3 ns). The WT hydrogenase can be considered to be in the intermediate regime. Another feature that ACE and hydrogenase share is that the diffusion rate does not correlate with the width of the main access tunnel. In ACE, this observation was explained by fast conformational gating (29). For hydrogenase, we found that the lack of correlation is due to mutations not only changing the diameter of the main channel but also creating an additional dynamical bottleneck and/or opening up alternative pathways bypassing the main diffusion path.
Outlook.
Can active site access of the inhibitor gases CO and O2 be further lowered? It would be interesting to probe the potential of the 74–476 motif in this regard by carrying out, e.g., mutations of R476. Unfortunately, this residue is strictly conserved within [NiFe]-hydrogenases, which may lead to complications in actual experimental mutations. Alternatively, one may speculate about the possibility of restricting the dihedral motion of R476 that gates the final transition to the active site cluster, for instance by engineering a stable salt bridge between a negatively charged residue at position 74 and R476. Similarly, the engineering of strong electrostatic interactions or the introduction of a disulfide bridge between residues 74 and 122 could lead to a further reduction in the thermal fluctuations of the channel width that are responsible for the residual CO leakage in existing mutants.
Materials and Methods
Computation of Diffusion Rates.
We describe the gas diffusive dynamics in the protein and solvent using a Markovian model (33, 34). Each gas molecule is assumed to stay either in a discrete state within the protein (also called “cluster” or “site” and denoted by i or j in the following), or in the solvent. The dynamics between these states can be approximated by a master equation, 
[3]
[4]where the kij are the rate constants for the transitions i←j, representing the constant elements of the rate matrix K, and pi is the population of cluster i. The rate matrix K is obtained from equilibrium simulations, where one assumes that the simulations are sufficiently long enough so that all relevant transitions are observed on the simulation timescale. Rates for transitions that are not adequately sampled on the timescale of the equilibrium simulations are obtained from nonequilibrium simulations in the presence of a pulling force. The molecule is pulled from the initial (j) to the final cluster (i) for a series of external forces of constant magnitude 
. Averaging over initial conditions, the simulations yield the force-dependent MFPT, τij(F), and corresponding transition rates, kij(F) = 1/τij(F). We use the one-dimensional DHS model for force-dependent kinetics (35), derived from Kramers theory (36), in order to extrapolate the transition rate constant from the pulling simulations to zero force, 
[5]where 
is the rate at zero force, which is equal to the desired rate matrix element kij. In Eq. 5, ΔG‡ is the height of the barrier, and x‡ is the distance to the transition state. The parameter ν is set equal to two-thirds, corresponding to a linear cubic form of the energy surface. In the next step, the time-dependent populations are solved according to Eq. 4 for given initial conditions. Because we are interested in the diffusion kinetics of gas molecules from the solvent (state S) to the protein active site (state G, for “geminate”), we set the initial population of statea, pS(t = 0), equal to one and the ones of all other states equal to zero as would be the case in experiment. The time-dependent population of site G, pG(t), is then obtained by solving Eq. 4. In the final step, pG(t) is fit to a phenomenological rate equation which yields phenomenological rate constants that can be compared to experiment. For hydrogenase, the association reaction between the free enzyme (Enz) and gas molecules (Gas) has been described in terms of a two step kinetic scheme (14), 
[6]The first step describes the diffusion of the gas into the active site of Enz, represented in the current calculations by the state where the gas molecule occupies cluster G. The corresponding diffusion rates are k+1 and k-1, respectively. The second step, occurring with forward and reverse rates k+2 and k-2, results in the chemical binding of the gas molecule to the active site (Enz–Gas). Using a classical force field, this second step is not represented in present calculations, meaning that the computed probability pG(t) is only due to the kinetics of the first reaction step. Hence, k+2 and k-2 can be set equal to zero in the kinetic model for pG(t). Assuming that the rate matrix K in Eq. 4 has been constructed for constant gas concentration, the population of the geminate cluster is given by the following pseudo-first-order kinetic equation, 
[7]Thus, a fit of pG(t) to Eq. 7 gives the desired phenomenological rate constants for diffusion from the solvent to cluster G, k+1, and for the reverse process, k-1.
Simulation Details.
Starting structures for the mutant enzymes were obtained from an equilibrated configuration of the WT enzyme (33) by mutating residues so that their side-chain conformation was the same as reported in the crystal structure where available (Protein Data Bank ID 3H3X, ref. 13, for V74M and 3CUR, ref. 12, for MM). The GROMOS96 43a1 protein force field (37), and the extended simple point charge model for liquid water (38) were used. For CO, we used a previously published three site model (16) (see Table S3 for parameters) that reproduces the experimental diffusion constant in water and n-hexane very well. The MFPT for transitions E←68 and G←E were obtained from pulling simulation at constant external force. For each force, it was necessary to average the MFPT over 50–200 trajectories of varying length (0.1–1 ns). The total simulation time for the slowest transition (E←68 in V74M) was 1.1 μs. Full simulation details are given in the SI Text.
Acknowledgments
We would like to thank Dr. Robert Best for fruitful discussions at the initial stage of this work and Dr. Christophe Leger for helpful comments. P.-h.W. acknowledges the Ministry of Education, Republic of China (Taiwan) for a PhD scholarship, and J. B. would like to thank the Royal Society for a University Research Fellowship. The University of Cambridge, Department of Chemistry, is acknowledged for computer support, and the United Kingdom’s High Performance Computing Materials Chemistry Consortium (funded by Engineering and Physical Sciences Research Council, EP/F067496) for access to the high-performance computing facility HECToR.
Footnotes
- ↵1To whom correspondence should be addressed. E-mail: j.blumberger{at}ucl.ac.uk.
Author contributions: J.B. designed research; P.-h.W. performed research; P.-h.W. contributed new reagents/analytic tools; P.-h.W. and J.B. analyzed data; and P.-h.W. and J.B. wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission.
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1121176109/-/DCSupplemental.
References
Simulation of gas diffusion in hydrogenase mutants
More Articles of This Classification
Physical Sciences
Chemistry
Biological Sciences
Related Content
Cited by...
Similar Articles
You May Also be Interested in
Both social and engineering factors must be further explored and scrutinized across the globe—as should notions of justice related to flooding impacts and responses.
Image courtesy of Shutterstock.com/Ekaterina Pokrovsky.
Blair Hedges discusses the state of Haiti's remaining forests and the consequences to biodiversity.
A study uncovers genetic evidence of early colonization of Jamaica by primates and suggests that adaptation can shape species morphology in novel environments, even when species come from morphologically conservative families.
Image courtesy of Lorraine Meeker (American Museum of Natural History, New York).
Researchers report that some moth species use sound absorbing-scales to avoid detection by bats and suggest that the findings could aid the design of biologically inspired sound absorbers for an array of applications, including noise mitigation and building acoustics.
Image courtesy of Thomas R. Neil.
