Hydrophobic Effect Drives Oxygen Uptake in Myoglobin via Histidine E7*

Background: The distal histidine in myoglobin is thought to act as a gate regulating oxygen uptake. Results: Neither the open nor closed conformation hinders oxygen uptake; a hydrophobic site is more apparent in the open conformation. Conclusion: The driving force for ligand uptake is the hydrophobic effect. Significance: This amplifies our understanding of the mechanisms through which proteins regulate ligand uptake. Since the elucidation of the myoglobin (Mb) structure, a histidine residue on the E helix (His-E7) has been proposed to act as a gate with an open or closed conformation controlling access to the active site. Although it is believed that at low pH, the His-E7 gate is in its open conformation, the full relationship between the His-E7 protonation state, its conformation, and ligand migration in Mb is hotly debated. We used molecular dynamics simulations to first address the effect of His-E7 protonation on its conformation. We observed the expected shift from the closed to the open conformation upon protonation, but more importantly, noted a significant difference between the conformations of the two neutral histidine tautomers. We further computed free energy profiles for oxygen migration in each of the possible His-E7 states as well as in two instructive Mb mutants: Ala-E7 and Trp-E7. Our results show that even in the closed conformation, the His-E7 gate does not create a large barrier to oxygen migration and permits oxygen entry with only a small rotation of the imidazole side chain and movement of the E helix. We identify, instead, a hydrophobic site in the E7 channel that can accommodate an apolar diatomic ligand and enhances ligand uptake particularly in the open His-E7 conformation. This rate enhancement is diminished in the closed conformation. Taken together, our results provide a new conceptual framework for the histidine gate hypothesis.

In the crystal structures of hemoglobin (Hb) 7 and myoglobin (Mb) (1,2), the histidine perched at the seventh position along the E helix resembles a gate guarding the entrance to the active site on the distal side of the heme group. Based on this observation, Perutz et al. (3) proposed that a swinging motion of this single amino acid side chain could control the migration rate for ligand entry into either Mb or Hb (Fig. 1), an idea that thus became known as the histidine gate hypothesis. More than 50 years of research have culminated in Mb becoming one of the most studied proteins in the field of structural biology (4). The histidine gate hypothesis and the role played by alternative ligand migration paths have been examined in detail by both experimental and theoretical methodologies but are still surrounded by controversy (5). It would be impractical to describe all the available data, and instead, we will attempt to summarize the key observations and discrepancies. Mb crystal structures have resolved the His-E7 side chain in the closed and open conformations (6). At neutral pH, the Mb crystal structures show the closed structure, which was originally observed by . Lowering the pH of Mb crystals resulted in His-E7 adopting an open conformation, which is attributed to His-E7 protonation and increased solvent exposure of imidazolium cation. Rabenstein et al. (7) analyzed the effect of pH on the His-E7 rotameric state theoretically and theirs results also indicated that when His-E7 is protonated the fraction of the open conformation increases.
The rates of ligand migration into and out of Mb are characterized experimentally by the overall ligand association (kЈ X ) and dissociation (k X ) rate constants, where X refers to the ligand being studied. The dissociation rate constant is determined by the speed of the iron-ligand bond breakage and the competition between the ligand migration to the solvent and internal rebinding to the iron atom. The limiting step is generally the iron-ligand bond breaking, with rates on the order of 0.01 to 100 s Ϫ1 , whereas rates of ligand migration out of the active site on the order of 10 6 to 10 7 s Ϫ1 , making the observed value of k X a relatively poor indicator of ligand trajectories. The rate of ligand association is determined by the speed of ligand entry into the distal site and the extent of iron-ligand bond formation versus ligand escape and is often a more direct indicator of ligand pathways. Thus, for protein variants or changes in pH, observed changes in kЈ X can be often be used to estimate changes in the initial barrier to the ligand entry (8 -10). If the fraction of geminate recombination, F geminate , is measured after photodissociation of the iron-ligand bond, the true bimolecular rate for ligand entry, kЈ entry , can be computed as kЈ entry ϭ kЈ X / F geminate because F geminate provides the fraction of ligand molecules that bind to the iron atom after they enter the active site (8).
Site-directed mutagenesis studies in Mb have shown that both kЈ entry and kЈ X increase up to 10-fold when His-E7 is replaced with smaller apolar residues, such as Ala or Gly, whereas both rate constants decrease ϳ3-fold when His-E7 is replaced with the larger Trp (11), highlighting the importance of the E7 pathway in ligand migration. These results also show that the Trp-E7 and Ala-E7 Mb mutants are two useful limiting cases for evaluating the E7 pathway (8,12). Association rate constants for ligand binding to Mb also suggest a pH dependence in the opening of the His-E7 gate because both kЈ O2 and kЈ CO increase 2-fold upon lowering the pH from 7.0 to 4.6 (13,14). This increase in the kЈ X is more dramatic upon further lowering the pH to 2.0, but under these conditions, the increase in rate is due to protonation of the proximal histidine and formation of a four-coordinate deoxyheme group (15). In addition to the mutant and pH kinetic work, numerous other studies have also lent support to the relevance of the His-E7 gate opening to ligand entry in Mb (16 -19).
Arguments against the His-E7 gate being the primary pathway for ligand entry come from studies showing the existence of multiple hydrophobic cavities (the Xe binding sites), which are thought to play a role in ligand migration, especially because internal pathways connecting the cavities have been identified (20 -28). In addition, all computational studies have even identified many possible ligand exit pathways involving the Xe cavities (29 -31), as well as ligand exit through the E7 gate (31). However, mutagenesis studies of Mb residues lining the Xe cavities only resulted in minor impacts on either kЈ O2 or kЈ entry (8,33), suggesting that the role of the Xe cavities in ligand uptake is likely secondary to that of the E7 pathway. Regardless of this controversy, there is a general consensus that the His-E7 gate is one of the principal pathways for ligand entry in myoglobin (11,29,30,34).
Although correlations between pH, His-E7 conformation, and ligand migration rates have been identified, there is still a need for a thermodynamic characterization of ligand migration through the His-E7 gate and its relationship to the conformation of the imidazole side chain. To fill this gap, we performed molecular dynamics simulations to understand the relationships between the following: 1) the His-E7 protonation state and the side chain conformation (open or closed), 2) the open and closed conformations and the free energy barriers to oxygen movement, and 3) Ala-E7 and Trp-E7 substitutions and rates of ligand entry. To our knowledge, this study is the first to explicitly attempt to characterize the relationship between distal histidine tautomeric states, side chain orientations, amino acid replacement, and the free energy barriers for oxygen migration through the E7 channel.

MATERIALS AND METHODS
System Preparation-Classical molecular dynamics (MD) simulations were performed to sample the configurational space of the systems of study. Coordinates from the Protein Data Bank were used to create initial structural models (Table 1). Each conformation found in the crystal structures of the wild-type Mb (open and closed) was simulated in the three tautomeric states: the neutral ␦-tautomer (HID), the neutral ⑀-tautomer (HIE), and the positively charged protonated state (HIP), giving us six different initial structures for the wild-type myoglobin.
Aside from His-E7, all histidine protonation states were assigned based on the immediate environment in the crystal structures and cross-checked against the protonation prediction server Hϩϩ (35). Throughout the simulations, all of the histidines maintained the same immediate environments, validating the choice of protonation states.
All systems were built using the Amber ff99SB force field (36) and solvated in a box of TIP3P water molecules, where the minimum distance between the protein and the extreme of  the box was 12 Å. Crystallographic waters were left as found in the original crystal structures. The heme parameters were determined using restrained electrostatic potential charges (37) and HF/6 -31G(d) wave functions according to the Amber standard protocol, as described previously (38 -45). Molecular Dynamics Simulations-MD simulations were performed using the pmemd module within the AMBER11 simulation package. Periodic boundary conditions and Ewald sums were employed for the long range electrostatic interactions (46). With the SHAKE algorithm (47), we used a 2-fs time step for the propagation of coordinates in time in the isothermal-isobaric ensemble, using the Berendsen thermostat and barostat (48). The equilibration process was performed by slowly heating the system from 0 to 300 K over the course of 500 ps. We performed a total of 150 ns of MD simulations (isothermal-isobaric ensemble) for the HID and HIP system and 200 ns for the HIE system. To provide some evidence about the convergence the simulations, we used the last snapshot of each trajectory to perform (after an additional equilibration period) 50 ns of MD with switched protonation states: HIE to HID, HIE to HIP, HIP to HIE, and HID to HIE.
Multiple Steered Molecular Dynamics-In the multiple steered molecular dynamics method, the original potential is modified by adding an harmonic restraint called V add (t), where k is the spring constant, is an arbitrary reaction coordinate, and 0 (t) is the time-dependent center of the spring, which is modified based on an arbitrary velocity (v) as follows in Equation 2.
Using this approach and integrating the force applied in each step, it is possible to calculate the irreversible work (W) through the reaction coordinate chosen. From multiple irreversible works obtained starting from different initial microconfigurations it is possible to employ Jarzynski's equation (49) to estimate the free profile of the process. Jarzynski's equation relates the irreversible works with free energy change as follows, where k B is the Boltzmann constant, T is the temperature, and W is the calculated irreversible work for each of the independent non-equilibrium process. A total of 40 independent perturbations were performed. This approach was used in a variety of studies such as ligand entry, escape, dynamic control of tunnel opening, and various mutagenesis studies, yielding results that are in excellent agreement with the experimentally determined rates (42)(43)(44)(45)50). Because the reaction coordinate involves the formation of the O 2 -iron bond, a process not allowed in the standard MD force field, we introduced a Morse potential to describe this coordination (51), where D o corresponds to the well depth (10 kcal/mol), r 0 is the equilibrium bond length (1.8 Å), and the parameter ␤ was 4 Å Ϫ1 . Error analysis of the free energy profiles was done considering the bias inherent to the Jarzynski free energy estimator for a small number of independent simulations (n ϭ 40) as proposed by Gore et al. (52).

RESULTS
Relationship between His-E7 Tautomeric and Conformational States-We performed extensive MD simulations to determine whether the His-E7 protonation state controls its conformation. To identify the protonation states of histidine, we use the following nomenclature: HIE, the neutral tautomer with the proton on the ⑀ nitrogen; HID, the neutral tautomer with the proton on the ␦ nitrogen; and HIP, the histidine in its positively charged protonated state. As discussed below, the fixed protonation simulations afford an interesting analysis of the nonequivalence of the two neutral histidine tautomers. Given the two possible conformations of the imidazole side chain (open or closed), we performed a total of 150-ns MD simulations for each tautomer (HIE, HIP, HID) starting from the two possible conformations to remove any bias due to the initial structure. Moreover, we perturbed the protonation state after each 50-ns block of simulation time (HIE to HID, HIE to HIP, HIP to HIE, and HID to HIE) to ensure causality between the change in protonation state and the observed conformation of the His-E7 side chain. We classified the conformation as open or closed according to the His-E7 side chain dihedral angle (C-C ␣ -C ␤ -C ␥ ) (ϳ160°open; ϳ60°closed) according to the observed values in the neutral and low pH Mb crystal structures as well as in previous MD studies (26). The results show that the protonated (HIP) and the neutral ⑀-tautomer (HIE) sample exclusively either the open or closed conformation, respectively (data not shown). The third form (HID), however, was found primarily in the open conformation, although not as exclusively as in the HIP form. In the HIP and HID forms, the open conformation is stabilized by a strong hydrogen bond interaction between the HN ␦ of His-E7 and the backbone carbonyl group of nearby Asp-60 (Fig. 2), as previously observed (7,53). Extending the gate analogy, this interaction serves as a "latch" for the His-E7 gate, maintaining it in the open conformation. Recall that the ␦ proton is present in both HIP and HID but not in the HIE state. These results suggest that His-E7 does not have to necessarily be in its charged (protonated) state to sample the open conformation. Thus, even at neutral pH, Mb should significantly populate the open state even if it cannot be "seen" in crystal structures.
Structure and Dynamics of the Ala-E7 and Trp-E7 Mb Mutants-To have a better understanding of the steric constraints in the E7 channel, we also simulated the two experimentally well characterized Ala-E7 and Trp-E7 Mb mutants, which provide the upper and lower bounds for the rates of ligand entry ( Table 2). Simulation of the Ala-E7 Mb mutant showed that the protein is stable with no observed changes in the distal site. The heme iron center is clearly accessible from the solvent, through the hole created by the absence of the E7 side chain (Fig. 3, top-right panel). The results for Trp-E7, on the other hand, were more surprising. The Trp-E7 side chain was drawn inward toward the distal site, as opposed to its initial outward conformation (Fig. 3, bottom panel, and Fig. 4). We observed a similar inward rotation of the hydrophobic side chain in the case of another mutant Phe-E7 Mb (data not shown). We should note that the new configuration adopted by Trp-E7 was previously observed in the crystal structures of the ␣ Trp-E7 mutant subunit of human recombinant deoxyHbA containing wild-type ␤ subunits (Fig. 4, right panel) (Protein Data Bank code 3NMM) (12). Both of the observed Trp-E7 conformations result in steric hindrance of the E7 pathway to the heme, as opposed to both Ala-E7 and the open His-E7 conformation (Fig. 3).
Oxygen Migration through the Conformation-dependent WT E7 Pathway-To characterize the oxygen entry process thermodynamically, we computed free energy profiles corresponding to the most populated conformation for each His-E7 state: the open conformation observed for HIP and HID and the closed conformation observed for HIE (Fig. 5). To compute each profile, we used multiple near-equilibrium perturbations and Jarzynski's equality (49,50) to take dioxygen in the solvent, placed just outside the protein, and bring it to the bound state at 1.8 Å from the iron. The free energy was set to a value of 0 at 10 Å, where the oxygen molecule is fully solvated. The iron-O 2 interaction was modeled using a Morse potential to avoid spurious steric repulsion when the oxygen is within bonding distance of the iron. This strategy has been successfully applied to the determination of ligand binding to heme proteins by our group and others (51, 54 -56).
The profiles for the open conformation (found in HIP and HID forms) present little or no barrier for oxygen entry into the E7 channel, which is characterized by a significant free energy

TABLE 2 Experimentally reported ligand entry rates of oxygen migration into Sperm whale myoglobin at pH 7.0
The bimolecular rate constant for ligand entry, kЈ entry , was calculated as kЈ O2 /F geminate (8). The errors for wild-type (WT) myoglobin were estimated from averages of multiple geminate analyses (Ͼ10) (8).
well at ϳ5 Å from the heme iron atom (Fig. 5). This well (ϳ3 kcal/mol) corresponds to an apolar cavity just above the edge of the porphyrin ring, but still away from the iron-O 2 bond formation region. The ligand in this free energy well is in contact with apolar portions of the Thr-E10, His-E7 and Val-E11 side chains (Fig. 6). From this well, the oxygen does have to overcome a small barrier to be oriented just above the iron center for bond formation. The origin of this small barrier is the presence of Val-E11 (Fig. 6, right panel). Surprisingly, not even the closed conformation exhibits a significant barrier for oxygen entry into the E7 channel free energy well from the solvent (Fig. 5, black curve). As in the open conformation, there is only a small barrier to reach the iron, which in this case is due to the presence of Val-E11 and His-E7 (Fig. 6, left panel). The only notable difference between the closed and the open conformations is the depth of the free energy well, which is shallower in the closed conformation (ϳ1.5 versus 3 kcal/mol). In the closed conformation of the His-E7 side chain, the imidazole ring occu-pies part of apolar channel that corresponds to the free energy well, sterically restricting O 2 in this pocket (Fig. 6, left panel). Oxygen Entry through E7 Pathway in Trp-E7 and Ala-E7 Mb Mutants-The free energy profiles for Ala-E7 and Trp-E7 Mb provide an internal control of the methodology as well as a deeper understanding of the microscopic basis of the different ligand entry rate constants listed in Table 2. As expected, the profile for the Ala-E7 mutant exhibits absolutely no barrier for oxygen entry through the E7 channel. Surprisingly, the same free energy well, ϳ5 Å from the iron atom, is found in the channel of Ala-E7 Mb and is considerably deeper than in the wild type case, ϳ4 kcal/mol more favorable than bulk solvent (Fig. 7, magenta curve). The much smaller Ala side chain creates a larger and completely apolar cavity above the edge of the heme ring in which the oxygen molecule is more readily accommodated.
As described in the last section, the Trp-E7 Mb mutant was observed to sample two conformations in our simulations, including the one found in the Mb crystal structure (outward), and the one similar to that observed in the Trp-E7 ␣ subunit of Hb mutant hybrid (inward). The resulting free energy profile for the outward but still blocked conformation of Trp-E7 shows a large ϳ3 kcal/mol barrier at 8 Å from the iron center (Fig. 7, blue curve). This barrier is directly related to the presence of the bulky Trp-E7 hindering access to the heme group (Figs. 3 and  4). In this case, the free energy well, which is 5 Å from the iron atom, is completely blocked. The inward conformation of Trp-E7 presents an even higher (ϳ6 kcal/mol) barrier, which in this case is only 3 Å from the heme because the indole side chain has rotated into the distal pocket completely hindering access to the iron atom (Fig. 7, orange curve). The location of both barriers (outward and inward) corresponds to the shifted position of the Trp-E7 side chain in each case: at 8 Å from the iron center in the outward conformation and at 3 Å from the iron center in the inward conformation. The differences in heights are due to the mobility of the Trp-E7, which has more flexibility in the outward conformation, allowing the oxygen to circumvent the bulky side chain more easily than when it is in the inward conformation.  (right panel). Trp-E7, Phe-CD1, Thr-E10, and Val-E11 are depicted in cyan. In both panels, the E7 pathway is blocked by the bulky Trp side chain (see Fig. 3). Position of Trp-E7 in the crystal structure of the Trp-E7 mutant ␣ subunit of recombinant human hemoglobin with wild-type ␤ subunit was superimposed in violet (Protein Data Bank code 3NMM) (12).

DISCUSSION
The goal of the present study was to understand the relationship between the rate of ligand entry into Mb and the His-E7 side chain, specifically the E7 protonation, tautomerization, and conformation, as well as its replacement with larger and smaller functional groups. Simulation of the limiting case mutants Trp-E7 and Ala-E7 Mbs provided internal tests for the validity of our methodology and a clearer understanding of the ligand entry process. The pH-dependent conformational change of the His-E7 side chain and the open conformation of the HID tautomer highlights the contribution of the hydrogen bond between the ␦ proton of the imidazole group and the nearby Asp-60 carbonyl, which stabilizes the open conformation. Our theoretical results are in agreement with experimen-tal crystallographic studies showing an increase in the open conformation upon lowering the pH. These simulations also underscore the importance of the amphipathic nature of the histidine residue, which allows it to adopt the open conformation, even in a neutral tautomeric state (HID) or, alternatively, serve as an interface between the hydrophobic site and solvent molecules in the HIE tautomeric state.
The Impact of the Apolar Free Energy Well in the E7 Channel-The most significant, and novel characteristic seen in the free energy profiles is the presence of a well located near the heme at ϳ5 Å from the iron. The well depth, relative to the free energy of O 2 in water, increases in the order: Trp-E7 (no distinguishable well) Ͻ His-E7-closed (ϳϪ1.5 kcal/mol) Ͻ His-E7-open (ϳϪ3 kcal/mol) Ͻ Ala-E7 (ϳϪ4 kcal/mol). The measured bimolecular rate constants for ligand entry, kЈ entry , shown in Table 2 follow the opposite trend: Ala-E7 Mb Ͼ His-E7-open (low pH) Mb Ͼ His-E7-closed (neutral pH) Mb Ͼ Trp-E7 Mb. Thus, the depth of the free energy well is directly related to the ligand entry rate. The deeper free energy well can enhance the rate of ligand entry given that solvated oxygen crosses no large barrier to reach the apolar cavity or from the cavity to the distal site itself. One possible mechanistic interpretation involves pre-concentration of molecular oxygen (or other small ligands) in the apolar cavity, which would be more favorable when His-E7 is in the open conformation or replaced with smaller amino acids. In this interpretation, the overall bimolecular rate of entry into the distal pocket is given by the first order rate of ligand movement from the apolar cavity into the distal site times the fraction of O 2 molecules in the well. The latter fraction is very small and would be determined by [ In a previous work studying oxygen uptake in the ␣ chain of an Antarctic fish Trematomus bernacchi hemoglobin, using a FIGURE 6. Schematic representation of the apolar channel site corresponding to the free energy well for oxygen migration, which is ϳ5 Å from the iron in wt Mb. In the left panel, His-E7 is in the closed conformation, perturbing the well, whereas in the right panel, it is open, leaving the well unrestricted for oxygen entry. Heme group (yellow) with iron atom (green) and relevant E7, E10 and E11 residues (cyan) are shown explicitly. The arrow indicates the small lateral movement that the His-E7 side chain and the E helix undergo to allow oxygen to migrate toward the iron in the closed conformation. FIGURE 7. Free energy profile of ligand migration through the E7 pathway for mutant proteins: Ala-E7 Mb (magenta curve), Trp-E7 Mb in the outward conformation as found in the crystal structure (blue curve) and Trp-E7 Mb in the inward conformation adopted during the simulation (orange curve). The free energy was set to a value of 0 at 10 Å, where the oxygen molecule is fully solvated. Error bars represent mean square error considering the bias inherent to the Jarzynski free energy estimator for a small number of independent simulations (n ϭ 40) (52). The free energy profile for WT HIE (dashed black curve) and HID (dashed green curve) from Fig. 5 are included as a reference. similar methodological approach, we also found that the closed conformation of His-E7 presents a very small free energy barrier to oxygen entry and that a free energy well at 5 Å is more prominent in the open conformation (38). Thus, although more exhaustive studies are needed, the conclusions from our simulations for ligand movement through the E7 channel seem to be valid for other globins with histidine in the E7 position as well.
The E7 Channel and Alternative Pathways-As mentioned in the Introduction, multiple pathways that do not involve the His-E7 gate have been proposed for ligand entry and exit in Mb. Such alternative migration pathways have been observed in other globins, such as in the Mycobacterium tuberculosis truncated Hb N (58 -60) and in the Cerebratulus lacteus neuronal mini-globin (18), in which very clear apolar tunnels appear to be the main routes for ligand migration and not the E7 channel. In mammalian Mbs and Hbs, however, extensive experimental evidence has pointed to the E7 channel as the principal avenue for ligand migration (8,11,12,18). Our calculated free energy profiles for both the open and closed His-E7 conformational states also lend considerable support to the E7 channel being a principal ligand migration route in Mb.
Comparison with Other Computational Studies-Previous simulation results of Mb have been interpreted to either support or discard the histidine gate hypothesis, based on the observed dynamics of His-E7. Bossa et al. (61) and Maragliano et al. (29) report the suspected opening of the His-E7 supporting the histidine gate hypothesis. On the other hand, some studies in both Hb and Mb report no opening motion of the His-E7 and raise doubt about the histidine gate hypothesis (32,62). We believe that the difference between the results that lead to opposing conclusions may depend only on the tautomeric form of His-E7 used. The first set of authors used HID (mostly open in our simulations) and the second set used either HIE tautomer (exclusively closed) or did not comment on the protonation state. By considering all three forms of His-E7, we could differentiate between the conformational preference of each neutral tautomer (HID and HIE), which may reconcile some of the seemingly contradictory results above. The same exhaustive approach was used in earlier works by Jewsbury et al. (53) and more recently by Rabenstein et al. 7) leading to similar conclusions.
Regarding the ligand migration through the E7 channel, Cohen et al. report a free energy well at ϳ5-6 Å from the iron center (30), which may be the same that we identify in the HIP/ HID profiles. A high barrier for ligand exit through the E7 channel was later reported by Maragliano et al. (29), where the HID state was used, supporting our computed free energy profile. The counterintuitive result where ligand escaped through the E7 pathway when the histidine is in the closed conformation has recently been observed for Hb, although it was interpreted as evidence against the histidine gate hypothesis (32). By considering free energy profiles for all the protonation/conformational states of His-E7, we show that although ligand entry may occur in the closed His-E7 conformation, the presence of a deeper free energy well in the open conformation enhances the uptake of oxygen, lending new support for the histidine gate hypothesis.
Conclusions-The natural image of a gate controlling the entry to a given place is that of a steric barrier when the gate is closed and then the absence of one when the gate is open. Our analysis showed little or no barrier for oxygen entry even in the closed His-E7 conformation. The opening of His-E7 does not significantly lower a barrier; rather, it vacates a hydrophobic cavity and results in a free energy well for small apolar ligands such as molecular oxygen. This analysis suggests that movement to the open conformation creates a favorable space for the entering ligand. Thus, the His-E7 conformation influences the thermodynamics of ligand migration through this pathway by creating a lower free energy pathway rather than altering a steric barrier.
The driving force for ligand uptake into this cavity is the hydrophobic effect, in which solvent water expels the apolar diatomic ligand toward the apolar cavity, which is more pronounced by outward movement of the imidazole side chain.