Internal water and microsecond dynamics in myoglobin.

Myoglobin (Mb) binds diatomic ligands, like O2, CO, and NO, in a cavity that is only transiently accessible. Crystallography and molecular simulations show that the ligands can migrate through an extensive network of transiently connected cavities but disagree on the locations and occupancy of internal hydration sites. Here, we use water (2)H and (17)O magnetic relaxation dispersion (MRD) to characterize the internal water molecules in Mb under physiological conditions. We find that equine carbonmonoxy Mb contains 4.5 ± 1.0 ordered internal water molecules with a mean survival time of 5.6 ± 0.5 μs at 25 °C. The likely locations of these water molecules are the four polar hydration sites, including one of the xenon-binding cavities, that are fully occupied in all high-resolution crystal structures of equine Mb. The finding that water escapes from these sites, located 17-31 Å apart in the protein, on the same μs time scale suggests a global exchange mechanism. We propose that this mechanism involves transient penetration of the protein by H-bonded water chains. Such a mechanism could play a functional role by eliminating trapped ligands. In addition, the MRD results indicate that 2 or 3 of the 11 histidine residues of equine Mb undergo intramolecular hydrogen exchange on a μs time scale.

Introduction

The 153-residue hemoprotein myoglobin (Mb) transports O2 from sarcolemma to mitochondria in cardiac and skeletal muscle tissue of most mammals.[1] In addition to O2, other diatomic ligands, such as CO and NO, bind to the heme iron and it has been suggested that Mb also acts as a NO scavenger.[2] Crystal structures[3−9] of Mb do not reveal the pathway(s) for ligand migration from the exterior to the buried binding site in the distal pocket (DP) at the heme-bound iron (Figure 1). Ligand access to the heme is therefore thought to be controlled by structural fluctuations.

Figure 1

Crystal structure[7]1DWR of equine MbCO in (a) ribbon and (b) cross-sectional surface (1.4 Å probe) representations, rendered with PyMol (Schrödinger, LLC). Color code: heme (orange), CO (cyan), internal water molecules (red), and xenon binding sites[52] (green).

Despite a large amount of experimental and computational work,[10−51] the mechanism of conformationally gated ligand migration through Mb has not been firmly established. In particular, while there is wide support for the direct pathway linking the DP to the exterior via the histidine (His-64) gate, the importance of ligand migration through the network of cavities in Mb, including the four Xe binding sites (Xe1–Xe4)[52] and the phantom (Ph) region (Figure 1), is still debated. Furthermore, the complex kinetics of ligand migration and escape under physiological conditions have not been fully elucidated, although the relevant time scale appears to be 0.1–10 μs at ambient temperature.

The permanent cavities in Mb, with a combined (probe-size-dependent[53]) volume approaching 300 Å3, could accommodate a dozen water molecules, but since most of them are lined by non-polar side chains, one expects a low water occupancy.[54,55] Indeed, high-resolution crystal structures of Mb have identified only four (equine Mb)[4,7,9] or five (sperm whale Mb)[5,6,8] internal water molecules outside the DP. The carbonmonoxy form[5−7] of Mb (MbCO) has no water in the DP, but a water molecule is found in the DP in the deoxy form[9] and is bound to the iron atom in the met form.[4−6,8] While ligand migration in Mb has been extensively studied, only a few studies have examined the involvement of internal water molecules in the mechanism and then the focus has been on the DP water.[17,30,34,41]

Information about internal water outside the DP has come from molecular dynamics (MD) simulations[44,51] and from magnetic relaxation dispersion (MRD) experiments.[56−58] The MD studies indicated several partly occupied internal hydration sites,[44,51] only some of which coincide with the crystallographic sites, but these simulations did not access the μs time scale of water and ligand migration.

The water MRD technique measures the magnetic field (or resonance frequency) dependence of the longitudinal spin relaxation rate R1 of one of the magnetic water nuclides 1H, 2H, or 17O.[59,60] Two of the MRD studies[56,57] employed the 1H nuclide to determine the paramagnetic relaxation enhancement for the DP water in met-Mb. The more recent MRD study,[58] which used the 2H and 17O nuclides, concluded that equine MbCO contains at least three internal water molecules with survival times exceeding the ∼10 ns tumbling time of Mb.

From MRD measurements on protein solutions, one can obtain the number of water molecules with a survival time longer than a lower bound set by the protein tumbling time (∼10 ns for Mb) and shorter than an upper bound set by the intrinsic spin relaxation rate for the internal water molecule (∼1 μs for 17O and ∼100 μs for 2H).[59,60] However, within this rather wide range, solution MRD experiments cannot differentiate water molecules with different survival times.

More information can be obtained if the protein is immobilized by a chemical cross-linking reagent such as glutaraldehyde (GA) so that molecular tumbling is inhibited. The MRD profile is then essentially the Fourier transform of the survival correlation function, and survival times up to ∼10 μs of individual internal water molecules can be determined directly from the 2H MRD profile.[61,62]

Here, we report an extensive set of 2H and 17O MRD data from equine MbCO. Water dynamics on the ns to μs time scale is characterized in detail by examining solutions with free-tumbling Mb as well as cross-linked gels with immobilized Mb, in both cases under fully hydrated conditions (∼1.6 mM Mb). To extract quantitative information about internal water molecules from 2H MRD data, it is necessary to isolate any contribution from rapidly exchanging labile deuterons in the protein. By analyzing the pD dependence of the 2H MRD and by comparing with the 17O MRD, we find that several histidine side chains in Mb exchange deuterons with water on the μs time scale, apparently by an intramolecular catalytic mechanism involving an adjacent deuteron exchange catalyst, such as a carboxylate.

The principal questions that we address here are the following: (1) How many internal water molecules occupy the internal cavities of Mb under physiological conditions? (2) Where are they located? In particular, are the Xe binding sites hydrated? (3) How long do the internal water molecules remain inside the cavity network? (4) By what mechanism do the internal water molecules exchange? (5) Does the internal water play a purely structural role, or are they involved more directly in the function of Mb?

Materials and Methods

Sample Preparation

Lyophilized equine skeletal muscle Mb was purchased from Sigma (≥95%) and further purified by cation-exchange chromatography (SP sepharose; GE Healthcare). The purified Mb was lyophilized following extensive dialysis against Milli-Q water. To prepare carbonmonoxy Mb (MbCO) samples, the purified protein was dissolved in 99.8% D2O (Cambridge Isotope Laboratories) with (gel samples) or without (solution samples) a buffer. For 17O MRD experiments, 20% 17O-enriched D2O was used. By using D2O for both 2H and 17O MRD experiments, correlation times can be compared directly without correcting for the (unknown) H/D kinetic isotope effect.

The solution was centrifuged at 13 000 rpm for 3 min to remove any insoluble protein, flushed by CO gas for 3 min on ice, and equilibrated by stirring for 30 min at 6 °C. Mb was then reduced by adding ∼2 equiv of sodium dithionite (Sigma), and the process of flushing CO gas for 2–3 min and stirring for 30 min was repeated (usually five times) until no further change in the optical absorption spectrum[63] (350–650 nm) could be detected (Shimadzu UV-1800). The Mb concentration was determined from the absorbance at 542 nm with an extinction coefficient, 18.4 mM–1 cm–1, that had been calibrated using complete amino acid analysis (performed at Amino Acid Analysis Center, Dept. of Biochemistry and Organic Chemistry, Uppsala University, Sweden). MRD measurements were performed on seven solution samples with 2.0–2.6 mM Mb and on 26 gel samples, 23 of them with 1.6–1.7 mM Mb (Table S1, Supporting Information).

The protein was immobilized by adding 25% GA solution (Sigma) to the MbCO solution at 6 °C.[64] Except for two cases (Table S1, Supporting Information), the GA/Mb mole ratio, NGA, was in the range 29.7–31.1. After mixing by pipetting, an aliquot (∼50 μL) of the solution was removed for pH measurement and the remainder was transferred to a glass tube (8 mm o.d. × 12 mm height), where it was cross-linked overnight at 6 °C. For the solution samples, pH was adjusted by adding either HCl or NaOH. Reported pD values are corrected for the H/D isotope effect according to pD = pH* + 0.41, where pH* is the reading on a pH meter calibrated with H2O buffers.[65]

To pressurize Xe in a cross-linked MbCO sample, the protein was immobilized in a 10 mm heavy wall NMR tube with a pressure/vacuum valve (Wilmad). After cross-linking, the NMR tube was filled with 8 bar of Xe gas (AGA, Sweden) and equilibrated for 1 week before NMR experiments. To estimate the Xe concentration in the MbCO sample, two reference samples (cyclohexane and 50 mM PIPES pD 7.4 in D2O) were prepared. The reference solutions were placed in standard 10 mm NMR tubes, and Xe gas was bubbled through the liquid.

NMR Experiments

The water 2H or 17O longitudinal relaxation rate, R1, was recorded on six different NMR setups: (1) a Stelar Spinmaster 1 T fast field-cycling (FC) spectrometer (2H, 1.5 kHz to 5.4 MHz); (2) a Tecmag Apollo console equipped with a GMW field-variable (≤1.8 T) iron-core magnet (17O, 0.7–1.6 MHz); (3) a Tecmag Discovery console equipped with a Drusch field-variable (≤2.1 T) iron-core magnet (2H, 2.5–13.1 MHz; 17O, 2.2–11.5 MHz) or (4) with a Bruker 4.7 T superconducting magnet (2H, 30.7 MHz; 17O, 27.1 MHz); and (5) Varian DirectDrive 500 (17O, 67.8 MHz) and (6) 600 (2H, 92.1 MHz; 17O, 81.3 MHz) spectrometers. For the 2H FC measurements, pre-polarized single-pulse (≤2.6 MHz) or inversion recovery (>2.6 MHz) sequences were used with pre-polarization and detection at 5.4 and 4.8 MHz, respectively.[66] On other spectrometers, the standard inversion recovery pulse sequences were used. 129Xe spectra were recorded on a Varian DirectDrive 600 (129Xe, 166 MHz) using a single 90° pulse. All NMR experiments were performed at 25.0 ± 0.1 °C, maintained by a thermostatted air flow. The sample temperature was checked with a thermocouple referenced to an ice–water bath and by recording the bulk solvent relaxation rate on a reference water or buffer sample.

Analysis of MRD Data

Water 2H and 17O MRD profiles, R1(ω0), from free-tumbling Mb in solution samples were analyzed with the standard multicomponent exchange model[59,60]where ω0 is the Larmor frequency in rad s–1, R10 is the known (frequency-independent) bulk solvent relaxation rate, and NW is the known water/Mb mole ratio in the sample. The inverse proportionality of the excess relaxation rate, R1ex(ω0) ≡ R1(ω0) – R10, to NW was used to normalize all MRD data to the same Mb concentration, corresponding to NW = 3 × 104. In eq 1, a component k, as identified by MRD profile analysis, may include several sites with similar correlation times (see below).

The first, frequency-independent, term within brackets in eq 1 refers to the external hydration shell, comprising NH water molecules with average dynamic perturbation factor ξH ≡ τH/τ0.[67] We estimate NH = 727 by computing the solvent-accessible surface area for equine Mb (average of PDB structures 1DWR, 1WLA, and 2VIK) with a probe radius of 1.7 Å and dividing by the mean surface area, 10.75 Å2, occupied by one water molecule.[67] The second term within brackets in eq 1 arises from internal water molecules (and labile deuterons), with N water molecules (or labile-deuteron water equivalents) in component k with intrinsic relaxation rate R1(ω0). Formally, N is the sum of occupancies of the sites belonging to component k. N can be non-integer, since these sites may be partially occupied.

In MbCO, the iron is in the diamagnetic ferrous low-spin state so there is no paramagnetic contribution to the water 2H or 17O relaxation rate. The intrinsic relaxation rate, induced by the nuclear electric quadrupole coupling, is described in the model-free approach as[60,68]The contribution from internal motions in the hydration site, typically on a sub-ps time scale, is negligibly small and has been omitted in eq 2. Here, Siso, is the usual isotropic orientational order parameter[59,60,68] and ωQ is the rigid-lattice nuclear quadrupole frequency (8.70 × 105 rad s–1 for 2H and 7.61 × 106 rad s–1 for 17O).[59] The correlation time τC, is related to the tumbling time τR of the protein (assumed to undergo spherical-top rotational diffusion) and the mean survival time (MST) τS, of a water molecule in site k as[60,68] 1/τC, = 1/τR + 1/τS,. Note that the survival time, often referred to as the residence time in the literature, is the time interval from an arbitrary time point to the next water exchange event.[69]

The isotropic rank-2 rotational correlation time τR for Mb was obtained from molecular hydrodynamics calculations using the program[70] HYDROPRO v. 10 with the recommended[70,71] effective (non-hydrogen) atomic radius of 3.0 Å and the 1.45 Å resolution crystal structure 1DWR of equine MbCO.[7] Corrected to the viscosity of D2O at 25 °C, we thus obtained τR = 11.1 ns, close to the value, 11.9 ± 0.4 ns, deduced[72] (after correction to D2O viscosity) from the field-dependent paramagnetic transverse relaxation rate of the proximal histidine Hδ1 in equine deoxy-Mb at 25 °C. Because of its oblate-like shape, Mb actually undergoes anisotropic rotational diffusion, but the computed 17% span of the five asymmetric-top rotational correlation times cannot be resolved by our MRD data.

As it stands, eq 1 is strictly valid only in the fast-exchange regime, where the MST in each site k that contributes significantly to R1(ω0) is sufficiently short that R1(0)τS, ≪ 1. In the dilute regime, where NW ≫ 1, however, eq 1 remains valid to an excellent approximation for arbitrarily long MST, provided that N and τC, in eqs 1 and 2 are reinterpreted as effective parameters according to[60,68]

Water 2H MRD profiles from immobilized Mb in gel samples were analyzed with the exchange-mediated orientational randomization (EMOR) model,[62] which has recently been quantitatively validated with the aid of an ultralong MD simulation of the protein BPTI.[69] Since the cross-linked protein cannot tumble, the water orientation is randomized by the exchange process itself, so that τC, = τS,. As a consequence, the conventional perturbation theory of spin relaxation is only valid in the fast-exchange regime, where ωQτS, ≪ 1. In this regime, the EMOR theory reduces to eqs 1 and 2 (with τC, = τS,). Outside this regime, meaning τS, > 1/ωQ ≈ 1 μs, the general form of the EMOR theory, based on the stochastic Liouville equation, must be used.[62] All fits to 2H MRD data from immobilized MbCO reported here are based on eq 4.7 of ref (62), which is exact in the dilute regime (NW ≫ 1). Whereas, in the solution case, each component k is modeled by the two parameters NSiso,2 and τC, (or their effective counterparts), there are four parameters in the gel case: N, S, η, and τS,.[62] Two of these parameters are related to the isotropic order parameter as[62]Siso, = S(1 + η2/3)1/2.

Since the 17O nuclide has spin quantum number I = 5/2, the relaxation behavior is in general more complex than that for 2H with I = 1.[59,61] In the dilute regime, however, eqs 1–3 remain valid to an excellent approximation in the solution case.[59,73] In the gel case, approximate expressions[61] are available for arbitrary τS,, while for τS, < 1/ωQ ≈ 100 ns the EMOR theory reduces to eqs 1 and 2, which are valid also for 17O.[59,73]

In fitting these models to the MRD data, we used the trust-region reflective non-linear optimization algorithm[74] with uniform 0.5% (17O solution MRD) or 1.0% (all other MRD profiles) R1 error, estimated from the scatter of the frequency independent bulk relaxation rate, R10. With this algorithm, the model parameters can be constrained to their physically admissible ranges.

Results and Discussion

Overview of MRD

Because 2H relaxation is relatively slow, the field-cycling technique[59,66] can be used to measure the MRD profile from immobilized Mb down to 1.5 kHz, allowing internal-water exchange to be monitored on time scales up to 10 μs. The drawback is that 2H relaxation does not probe water deuterons exclusively but also may contain a contribution from labile deuterons (LDs) in acidic protein side chains.[59,61,75,76] (This problem is even more severe for 1H MRD.[59,77,78]) An essential step in the interpretation of 2H MRD data is therefore to separate the contributions from D2O molecules and LDs. To do this, we use four strategies: (1) we vary pD to change the acid- and base-catalyzed LD exchange rate and to titrate the LD-bearing side chains; (2) we use different buffers to examine buffer-catalyzed LD exchange; (3) we compare 2H MRD profiles from immobilized Mb with profiles from solutions of free-tumbling Mb, where the LD effect is usually more pronounced; and (4) we record 17O MRD profiles (for both immobilized and free-tumbling Mb), which monitor water molecules exclusively but do not extend below 0.7 MHz.

In sections 3.2–3.4, we analyze the 2H and 17O MRD data. The interpretation of the data gives a consistent picture of the internal hydration of Mb but involves several technical details. As a guide to the reader, we summarize the key findings here.

2H MRD profiles from immobilized MbCO show three dispersion components with characteristic mean survival times (MSTs) of 5.6 μs, 121 ns, and 6 ns, respectively. The one with the slowest MST (5.6 μs) is the dominant component and is attributed to internal water molecules and labile 2H in histidine side chains. The other two components are assigned to internal motions of His side chains (121 ns) and to confined external water molecules (6 ns), as seen in the previous protein gel MRD study.[61] The pD dependence of the gel 2H MRD profiles demonstrates that 5.2 ± 0.6 internal water molecules and two or three His side chains are responsible for the 5.6 μs component.

Solution 2H MRD profiles from free-tumbling MbCO can be described by a single dispersion component. This component is attributed to the two slow components (5.6 μs and 121 ns), and the obtained parameters are consistent with the gel profiles.

The 17O MRD profile from immobilized MbCO reveals an additional internal water molecule with an MST of 32 ns. This water molecule is not evident in the 2H MRD data, presumably because 180° flips about the water dipole axis reduce the 2H order parameter. The solution 17O MRD profile is described by a single component just like the solution 2H MRD. The solution 17O MRD can be accounted for by the 5.6 μs and 32 ns components, whereas the solution 2H dispersion is produced by the 5.6 μs and 121 ns components. The gel and solution 17O MRD profiles indicate that 3.9 ± 0.6 internal water molecules are responsible for the 5.6 μs component, whereas the 2H analysis yields 5.2 ± 0.6. Our final estimate of 4.5 ± 1.0 internal water molecules with an MST of 5.6 μs is the average of the 2H and 17O derived results.

2H MRD from Immobilized MbCO

2H MRD profiles from MbCO immobilized by cross-linking with GA, which reacts primarily with lysine side chains,[64] were recorded at different Mb concentrations, GA/Mb mole ratios (NGA), pD values, buffer conditions, and temperatures, in all 25 profiles (Table S1, Supporting Information). Figure 2a shows a typical profile, measured at pD 7.0. The first thing to note is that the low-frequency dispersion has a large amplitude, ∼3 times (after NW normalization) that obtained from the protein BPTI with three contributing internal water molecules.[61] The dominant contribution to this large dispersion is a component with an MST of 5.6 ± 0.5 μs (Figure 2a, Table 1). To reveal any LD contribution to this large dispersion, we recorded MRD profiles at different pD values in the range 5.66–7.15. pD values outside this range were excluded by heme dissociation[63] and concomitant Mb aggregation (lower pD) or by the GA cross-linking chemistry[64] (higher pD).

Figure 2

(a) 2H MRD profile, R1(ω0), from immobilized MbCO at 25 °C, pD 7.0, and NGA = 29.9, scaled to NW = 3 × 104. The error bars do not exceed the symbol size. The solid curve is the fit that resulted in the parameter values in Table 1. The dashed curves are dispersion components 1 (magenta), 2 (blue), and 3 (red). The bulk water (black dotted) and external hydration shell (black dashed) contributions are also shown. (b) Zero-frequency excess relaxation rate, R1ex(0), obtained from fits to 22 MRD profiles, versus pD. The red curve is a fit based on the histidine pKa values in Mb, and the dashed line is the deduced pD-independent water contribution. The difference between the two curves hence represents the His LD contribution. Solid black symbols represent different buffer conditions as indicated. Colored symbols refer to different MbCO concentrations: 2.65 mM (green solid circle), 1.35 mM (green open circle), or 1.02 mM (blue solid circle); and NGA: 110 (blue solid circle) or 62.1 (blue open circle). In all other cases, CMb = 1.60–1.67 mM and NGA = 29.7–31.1.

Table 1

Results of Fits to 2H and 17O MRD Profilesa

		immobilized Mb		free-tumbling Mb
kb parameter (unit)		2H	17O	2H	17O
1	τC,k (μs)	5.6 ± 0.5	[5.6]
1	Nk	6.7 ± 0.8	[5.2]
1	Sk	0.79 ± 0.02	[0.79]
1	ηk	0.3 ± 0.1	[0]
2	τC,k (ns)	121 ± 10	32 ± 3	9.5 ± 0.5	5.6 ± 0.2
2	NkSiso,k2	0.70 ± 0.03	1.0 ± 0.1	5.6 ± 0.3	3.0 ± 0.1
3	τC,k (ns)	6.0 ± 1.0	4.3 ± 1.0
3	NkSiso,k2	3.6 ± 0.6	4.0 ± 0.6
	ξH	7.3 ± 0.5	7.7 ± 0.4	5.7 ± 0.5	5.8 ± 0.1
	χred2	1.16	0.36	1.17	0.21

All four MRD profiles were measured at 25 °C and pD 7.0. Quoted errors correspond to one standard deviation, propagated from the uniform 0.5 or 1.0% uncertainty in R1. Parameter values within square brackets were fixed during the fit.

Component 2 has different physical origins in the four cases (see text).

Figure 2b shows the pD dependence of the zero-frequency excess relaxation rate R1ex(0) ≡ R1(0) – R10, extrapolated from the R1 data with the aid of the parameters derived from fits to 22 individual MRD profiles like the one in Figure 2a (Table S1, Supporting Information). Several conclusions can be drawn from these results. Since all data points fall on the same master curve despite having been recorded on samples with different Mb concentration (1.02–2.65 mM) and different NGA (29.7–110), we infer that the internal-water dynamics probed by the MRD profile do not depend significantly on parameters that might affect the structure of the cross-linked protein network. When the Mb concentration was raised to 5.14 mM, however, a 23% reduction of R1ex(0) was observed. This might be due to impaired exchange for some of the internal water molecules or labile deuterons caused by increased protein–protein contacts at high concentration. The subsequent analysis of MRD data from immobilized Mb is restricted to samples with CMb ≈ 1.6 mM and NGA ≈ 30.

The second major conclusion that follows from the data in Figure 2b is that R1ex(0) must have a significant LD contribution. Before attempting to separate the pD-dependent LD contribution from the pD-independent internal-water contribution, we shall identify the side chains and the exchange mechanism responsible for the LD contribution. In doing so, we are guided by two observational facts: (1) the LD contribution to R1ex(0) decreases with pD in the investigated pD range 5.66–7.15, and (2) the MST, τLD, of the LDs must be close to the 5.6 μs MST of the dominant dispersion component (Table 1). In the slow-exchange regime (τLD > 1 μs), we expect the LD contribution to R1ex(0) to be approximately proportional to NLD(pD)/τLD(pD).[62] The observed pD dependence can hence reflect changes in the number, NLD, of LDs and/or in τLD.

The LD exchange mechanism can be intermolecular (involving a single side chain and one or more solvent species including buffer ions) or intramolecular (involving two proximal side chains). We consider first the usual intermolecular mechanism. Due to the observed buffer independence of R1ex(0) (Figure 2b), buffer catalyzed LD exchange can be neglected. We therefore need to consider only the acid (D3O+) and base (OD–) catalyzed exchange mechanisms. LD exchange from basic groups (Arg, Lys, His, N-terminus) is base catalyzed in the examined pD range. Any LD contribution from these groups would vary with pD as NLD/τLD ∝ 1/(10–p + 10–pD); that is, it would increase monotonically with pD, in contrast to the observed decrease (Figure 2b). Furthermore, even if base catalyzed LD exchange proceeds at the maximum (diffusion controlled) rate, τLD at pD 7 is 3 orders of magnitude longer than the required μs time scale. LD exchange from hydroxyl-bearing groups (Ser, Thr, Tyr) is both acid and base catalyzed. In principle, the observed pD dependence (Figure 2b) could therefore reflect acid catalysis. The rate constant, ka, for acid catalysis is, however, typically 4 orders of magnitude smaller than the rate constant, kb, for base catalysis.[75,76,79] The minimum LD exchange rate should therefore occur close to pD* = (1/2) log[ka/(KWkb)] ≈ 5.5, where KW = 10–14.95. In the examined pD range, we would therefore expect R1ex(0) to increase with pD if hydroxyl LDs were responsible. Furthermore, even with the maximum conceivable rate constants[75,76,79] (ka ≈ 107 M–1 s–1 and kb ≈ 1011 M–1 s–1), τLD is too long, by at least 2 orders of magnitude, in the examined pD range. LD exchange from carboxyl groups (Asp, Glu, C-terminus) usually occurs directly to water without involvement of catalytic species, so that τLD is independent of pD.[75,76] In principle, the observed pD dependence (Figure 2b) could therefore reflect deprotonation of COOD groups. For immobilized BPTI, a small COOD contribution was seen at pD 4.4, but not at pD 6.5 where all carboxyl groups in BPTI are deprotonated.[61] For MbCO, carboxyl groups with strongly upshifted pKa are not indicated,[80] so carboxyl LDs cannot explain the pD dependence in Figure 2b.

On the basis of these considerations, we conclude that the LD contribution observed here is produced by intramolecular LD exchange with an essentially pD-independent rate. The observed pD dependence in R1ex(0) must then reflect the ionization equilibrium. With no pKa upshifted carboxyl group, only the imidazolium (His) groups are titrated in the examined pD range. Equine Mb contains 11 histidines (BPTI has none) with pKa values in the neutral range (Table S2, Supporting Information).[80−83] We then expect that R1ex(0) = a + NLD(pD)b, where a is the internal water contribution and NLD is the number of labile His deuterons (two in the acidic form, none in the basic form). Using published pKa values,[80−83] we find that this function can reproduce the observed pD dependence in R1ex(0). The fit shown in Figure 2b yields a pD-independent contribution (the parameter a) to R1ex(0) of 26.4 s–1, corresponding to 73.1% of R1ex(0) at pD 7.0. Although all the His residues were included in this analysis, the result only depends on the pKa distribution (in the examined pD range) and not on the number of His residues. Indeed, very similar results are obtained (26.2 s–1 and 73.3%) if we assume that the pD dependence is produced by a single His residue with a typical pKa value of 6.5.

The LDs of His side chains can exchange via an intramolecular exchange mechanism, where a water molecule bridges a positive His residue and negative carboxylate group or a positive/neutral His pair (section S1, Supporting Information), similar to what has been proposed for certain low-molecular-weight His derivatives.[84] The crystal structure of equine MbCO reveals several His residues that could engage in such intramolecular exchange (section S1, Figure S1, Supporting Information). A more detailed discussion of different exchange mechanisms for His LDs can be found in section S1 (Supporting Information).

The 2H MRD profiles from immobilized MbCO are strongly dominated by the slowest kinetic component (Figure 2a), hereafter referred to as component 1. At pD 7.0, the fit yields N1 = 6.7 ± 0.8 for the number of water molecules and/or LD water equivalents (one acidic His residue equates to one D2O molecule) responsible for this component (Table 1). Hydration sites may have fractional occupancy and His residues may be partly titrated, so N1 need not be an integer. The analysis in Figure 2b indicates that 73.1% of R1ex(0) is produced by water deuterons at pD 7.0 and the MRD fit in Figure 2a shows that R1ex(0) at pD 7.0 can be broken down as follows: 91.3% component 1 (with contributions from water and LDs), 6.1% component 2 (entirely due to LDs; see below), and 2.6% from component 3 and the frequency-independent hydration-shell component (both entirely due to water; see below). We thus infer that (73.1 – 2.6)/0.913 = 77.2% of component 1 is produced by N1W = 0.772N1 = 5.2 ± 0.6 internal water molecules.

The remaining 22.8% of component 1 at pD 7.0 is attributed to LDs in His residues (see above). If the locally averaged nuclear quadrupole coupling is the same for the LDs as for the ∼5 internal water molecules, we can infer that N1LD = 3.0 ± 0.4 LDs contribute to component 1 at pD 7.0. Since the (rigid-lattice) 2H quadrupole coupling constant for the acidic imidazole deuterons[85] is expected to be smaller than for D2O, the actual value of N1LD may be as large as five. We therefore conclude that at least two, and possibly three, His side chains (with two LDs per imidazolium group) engage in water-mediated intramolecular LD exchange on the μs time scale.

According to the fit, the MST of the water molecules contributing to component 1 is τS,1 = 5.6 ± 0.5 μs (Table 1). That ∼5 internal water molecules, located at different sites within the MbCO molecule, have the same MST suggests that they exchange by a common mechanism (section 3.6). But it should be noted that even though the dominant low-frequency dispersion is reproduced by a single kinetic EMOR component, a modest MST distribution would not have been resolved by our data.

The value τS,1 = 5.6 μs is on the long side of the maximum of R1(0) versus τS, corresponding to τS ≈ 1/ωQ ≈ 1 μs.[61,62] In this slow-exchange regime, R1(0) decreases with increasing τS. For example, if τS,1 had been 1 μs, R1(0) would have been several-fold larger and the dispersion would have appeared at a higher frequency. If τS,1 had been much longer, say 50 μs, then the dispersion frequency would only have been slightly downshifted, but the amplitude would have been strongly reduced.[62] In other words, a much larger number of internal water molecules or LDs would have been required to account for the observed MRD profile if the MST were much longer. In fact, in the ultraslow motion regime, where τS ≫ 1/ωQ ≈ 1 μs, R1(ω0) is approximately proportional to N/τS.[62] Even though the observed 5.6 μs MST is not long enough to be in this regime, the parameters N1 and τS,1 have a sizable covariance. Reassuringly, constrained fits to the 2H MRD profile with different fixed values of N1 yield the best fit quality for N1 between six and seven (Figure S2, Supporting Information), consistent with the result, N1 = 6.7 ± 0.8, from the unconstrained fit.

To reproduce the 2H MRD profile also at higher frequencies, the model must include, in addition to the dominant 5 μs component, two components in the fast-exchange regime (Table 1). One of these components has τS,2 = 0.12 ± 0.01 μs and N2Siso,22 = 0.70 ± 0.03. As argued in section 3.4, this component can be attributed to internal motions of rapidly exchanging LDs in His side chains. The third dispersion component has τS,3 = 6 ± 1 ns and N3Siso,32 = 3.6 ± 0.6 (Table 1). Similar parameter values were obtained for the highest-frequency component of the 2H MRD profiles of immobilized BPTI and ubiquitin.[61] As before, we attribute this component to a small number of confined water molecules in the external hydration shell.[61] Finally, the dynamic perturbation factor, ξH = 7.3 ± 0.5, deduced from the frequency-independent excess relaxation, is somewhat larger than that for free-tumbling Mb (section 3.3). This was also the case for immobilized BPTI (ξH = 7.9 ± 0.2) and ubiquitin (6.0 ± 0.2),[61] both of which yield ξH ≈ 4 for the free-tumbling proteins.[67] This difference between gel and solution samples, as well as the fact that component 3 is not observed in solution (section 3.3 and 3.4), can be rationalized in terms of gel-induced confinement due to short intermolecular cross-links and/or intramolecular cross-links (equine Mb has 19 Lys residues).

In summary, the high-frequency tail of the 2H MRD profile from immobilized MbCO exhibits the expected features associated with the external hydration shell, whereas the low-frequency part (components 1 and 2) reflects the occupancy and dynamics of internal hydration sites and LDs in a few His residues. The observed pD dependence of the 2H MRD profile indicates that 5.2 ± 0.6 internal water molecules with MST 5.6 ± 0.5 μs contribute to the dominant dispersion component 1.

2H MRD from Free-Tumbling MbCO

As a check on the interpretation of the 2H MRD data from immobilized MbCO, we measured the 2H MRD profile from free-tumbling MbCO in solutions with pD in the range 5.0–10.0. (For the gel samples, the cross-linking chemistry[64] limited us to pD ≤ 7.2.) In the acidic pD range, a low-frequency MRD component was observed with a correlation time corresponding to tumbling of large protein aggregates, presumably induced by heme dissociation.[63] At pD 5.0, this component dominated the 2H MRD profile. We therefore restrict the analysis to R1 data measured at ω0/2π > 0.9 MHz on solutions with pD ≥ 6. Under these conditions, the fraction aggregated protein is so small (≪1%) that the associated R1 contribution can be neglected.

Figure 3a shows the 2H MRD profile from free-tumbling MbCO at pD 7.0. At all pD values in the range 6.0–10.0, the MRD data (for ω0/2π > 0.9 MHz) can be adequately represented by a single dispersion component. The correlation time, τC, deduced from the single-component fit is shorter than the Mb tumbling time, τR = 11.1 ns (section 2.3) and it decreases with increasing pD: from 9.6 ± 0.4 ns at pD 6.0 to 7.4 ± 0.5 ns at pD 10.0 (Figure 3b). These findings can be explained if there is a pD-dependent contribution from LDs exchanging on the ∼100 μs time scale, as expected for base-catalyzed LD exchange in hydroxyl, ammonium, and guanidinium groups.[61,75,76,79] In contrast to the rapidly exchanging His LDs (section 3.2), these more slowly exchanging LDs do not contribute to the 2H MRD profile from immobilized MbCO because of the more restrictive fast-exchange criterion ωQτS < 1. For a free-tumbling protein, the fast-exchange criterion is ωQ(τRτS)1/2 < 1 (section 2.3), so LDs with τS on the order of 100 μs make a significant contribution, albeit with reduced effective weight Neff and correlation time τCeff according to eq 3. As explained in more detail in section S2, τC decreases with pD (Figure 3b) because the relative contribution to the MRD profile from LDs with reduced effective correlation time increases with pD.

Figure 3

(a) 2H MRD profile, R1(ω0), from free-tumbling MbCO at 25 °C and pD 7.0, scaled to NW = 3 × 104. The error bars are up to twice the symbol size. The solid curve is a single-component fit to the solid data points (ω0/2π > 0.9 MHz), with the parameters in Table 1. The bulk water (dotted line) and external hydration shell (dashed line) contributions are shown. (b) Amplitude NSiso2 (red squares) and correlation time τC (blue circles), obtained from single-component fits to five MRD profiles, versus pD. The dotted lines are visual guides.

The pD dependence of the amplitude parameter NSiso2 has a broad minimum of 5.0 ± 0.3 at pD 8–9 (Figure 3b). It decreases first because the His residues involved in intramolecular LD exchange are titrated (section 3.2), consistent with the decrease of R1ex(0) with pD in Figure 2b. Above pD ∼ 9, where all His residues in equine Mb are titrated,[80−83]NSiso2 increases because base-catalyzed exchange brings an increasing number of hydroxyl-bearing or basic (mainly Arg) side chains in Mb into the fast-exchange regime.[75]

The dynamic perturbation factor ξH for the external hydration shell does not vary significantly in the pD range 6.0–8.0, where ξH = 5.5 ± 0.2 (mean ± standard deviation). This value is within the range found for other proteins.[67,86] At pD 10, ξH = 6.5 ± 0.6 is slightly larger, most likely due to a frequency-independent R1 contribution from sub-ns internal motions of rapidly exchanging LDs.[75]

The single component in the solution MRD profile (Figure 3a) must include both components 1 and 2 in the gel MRD profile (Figure 2a). In solution at pD 7.0, these two components should appear as a single (unresolved) component with amplitude (see section 2.3) NSiso2 = N1S12(1 + η12/3) + N2Siso,22 = 5.0 ± 0.9 (Table 1), consistent with the value NSiso2 = 5.6 ± 0.3 deduced from the solution profile in Figure 3a. A slightly different way of showing this consistency is to compute NSiso2 as a function of N1 using parameters from constrained gel MRD fits with different fixed N1 values. The intersection of the resulting NSiso2 = f(N1) curve with the solution value 5.6 ± 0.3 yields N1 = 7.2 ± 0.6 (Figure S3, Supporting Information). In this analysis, we ignored the minor contribution from slowly exchanging LDs to the solution NSiso2 at pD 7.0. Using the minimum (at pD 8.0) solution value NSiso2 = 5.0 ± 0.2, we obtain instead N1 = 6.6 ± 0.6. In either case, the result of this analysis provides further support for the value N1 = 6.7 ± 0.8 obtained from the unconstrained fit (Table 1).

In summary, the mutually consistent 2H MRD profiles from immobilized and free-tumbling protein indicate that MbCO contains 5.2 ± 0.6 highly ordered (S = 0.79 ± 0.02) internal water molecules with a common MST of 5.6 ± 0.5 μs. In addition, the MRD profiles contain pD-dependent contributions from LD exchange and internal motions in His (neutral pD) and other (basic pD) side chains.

17O MRD from Immobilized and Free-Tumbling MbCO

While 17O MRD is rigorously water selective, the fast 17O relaxation precludes the use of field-cycling to access the low-frequency range of the 5 μs dispersion.[59] However, even if this were possible, a water molecule with MST of 5 μs would be far into the slow-exchange regime (since ωQ2 is 2 orders of magnitude larger for 17O than for 2H) and it would therefore contribute very little to the 17O MRD profile from immobilized MbCO.

Nevertheless, the 17O MRD profile from immobilized MbCO at pD 7.0 exhibits a large dispersion in the examined frequency range (Figure 4a). In analyzing this MRD profile, we include the small contribution from component 1, taking τS,1 = 5.6 μs and S1 = 0.79 from the 2H MRD fit (Table 1) and using the estimate N1 = 5.2 for the number of water molecules associated with this component (section 3.2). Even at the lowest frequency, this contribution hardly exceeds the measurement error in R1 (Figure 4a). Beyond this minor contribution, two fast-exchange components are required to account for the 17O MRD profile (Table 1).

Figure 4

(a) 17O MRD profile, R1(ω0), from MbCO gel at 25 °C, pD 7.0, and NGA = 29.9, scaled to NW = 3 × 104. The error bars are up to twice the symbol size. The solid curve is the fit that resulted in the parameter values in Table 1. The dashed curves are dispersion components 1 (magenta), 2 (blue), and 3 (red). The bulk water (dotted line) and external hydration shell (dashed line) contributions are shown in both panels. (b) 17O MRD profile, R1(ω0), from MbCO solution at 25 °C and pD 7.0, scaled to NW = 3 × 104. The error bars are comparable to the symbol size. The solid curve is a single-component fit, with the parameters in Table 1.

The dynamic perturbation factor ξH and the parameters τS,3 and N3Siso,32 of the fastest dispersion component do not differ significantly from the corresponding 2H parameters (Table 1), as expected if these contributions reflect weakly anisotropic water motions in the external hydration shell in the gel. The amplitude parameter, Niso,22 = 1.0 ± 0.1, of the slower 17O component is similar to that of 2H component 2 (Table 1), but the factor 4 difference in the associated correlation times shows that these components have different physical origins. A water molecule with 120 ns MST would give rise to a very large 17O dispersion, which is not observed, so we must conclude that this 2H component is produced by LDs. Since 120 ns is too short to be an MST for a LD (all carboxyl groups are titrated at pD 7.0), we assign 2H component 2 to internal motions of the His residues that contribute to 2H component 1. Internal motions of protein side chains on this time scale have previously been inferred from 2H MRD studies of several proteins.[61,75,76]

17O component 2, on the other hand, must reflect water motions. In principle, this could be local motions of the 5 μs internal waters. The only local water motion that might be as slow as 32 ns, however, is a 180° flip about the dipole axis, which, for symmetry reasons, cannot induce 17O relaxation.[59] We therefore conclude that the dominant 17O component 2 represents a single highly ordered water molecule with an MST of 32 ± 3 ns. This water molecule is unlikely to be deeply buried, and it must exchange by a different mechanism than that used by the 5 μs waters. But why is this water molecule not evident as a 32 ns component in the 2H profile? If it undergoes 180° flips on the time scale 0.1–10 μs, the 2H parameter Siso2 is reduced by a factor of 3.[59] The contribution of this water molecule would then not be resolved because it would only add a small dispersion amplitude (0.3 s–1 in R1) between the 6 and 121 ns dispersion steps (Figure 2a).

As in the 2H case, the 17O MRD profile from free-tumbling MbCO is well described by a single dispersion component (Figure 4b). The dynamic perturbation factor, ξH = 5.8 ± 0.1, is the same as that obtained from the 2H solution profile (Table 1), as expected since water motions in the external hydration shell are nearly isotropic.[59] The correlation time, τC = 5.6 ± 0.2 ns, is shorter than the Mb tumbling time, τR = 11.1 ns, because the single 17O solution MRD component is an unresolved superposition of two components, both with correlation times shorter than τR. The first component, due to water molecules with τS = 5.6 μs as deduced from the gel 2H MRD profile, is in the intermediate-exchange regime for 17O, with the effective correlation time shortened by a factor of ∼2 according to eq 3. The contribution of this component to the total amplitude parameter, NSiso2 = 3.0 ± 0.1 (Table 1), is reduced by the same factor. The second component, due to the single water molecule with τS = 32 ns as deduced from the gel 17O MRD profile, is in the fast-exchange regime with effective correlation time (1/τR + 1/τS)−1 = 8.3 ± 0.2 ns (section 2.3).

A quantitative analysis of the two contributions to NSiso2τC, which is proportional to the magnitude of the dispersion, shows that the solution and gel 17O MRD profiles are consistent, within the experimental error, with the gel 2H MRD profile (section S3, Figure S4, Supporting Information). Of particular interest is the number N1W of internal water molecules with an MST of 5.6 μs. Our analysis of the pD dependence of the gel MRD profile (section 3.2) indicates that N1W = 5.2, but the 17O MRD data suggest a somewhat smaller value (section S3, Figure S4, Supporting Information). Averaging these results, we arrive at the conservative estimate N1W = 4.5 ± 1.0.

In summary, all the 2H and 17O MRD profiles measured on MbCO gels and solutions can be interpreted in a self-consistent way. The gel 2H MRD profile and the solution 17O MRD profile are both consistent with N1W = 4.5 ± 1.0 ordered water molecules with an MST of 5.6 ± 0.5 μs. The solution and gel 17O MRD profiles indicate an additional component, not resolved in the 2H profiles, comprising a single highly ordered water molecule with an MST of 32 ± 3 ns.

Internal Hydration Sites

MRD experiments can provide the number of internal water molecules and their mean survival times, but the location of a hydration site can only be determined by MRD if the water molecule can be displaced by a ligand or a substituted side chain.[59] Mb contains four cavities, denoted Xe1–Xe4 (Figure 1), that have been shown to bind xenon at elevated pressure.[52] To determine if the internal water molecules identified here occupy any of the Xe sites, we recorded the 2H MRD profile from an MbCO gel that had been equilibrated with 8 bar of Xe gas. We determined the Xe occupancies of these sites from 129Xe NMR spectra of the MbCO gel and of two reference samples (Figure S5, Supporting Information), using also the reported[87] Xe binding constants for equine Mb at pH 7.0 and 25 °C. In this way, we obtained a Xe occupancy of 0.80 ± 0.09 for site Xe1 and a combined occupancy of 0.10 ± 0.05 for sites Xe2–Xe4 (section S4, Tables S3 and S4, Supporting Information).

As seen from Figure 5 and Table S5 (Supporting Information), Xe binding at this level has no significant effect on the 2H MRD profile. With 80% Xe occupancy in Xe1, we would have expected 100 × 0.8/6.3 = 13% reduction of R1 at low frequency if Xe1 had been fully occupied by a long-lived water molecule. We can therefore rule out internal water in this site. This is hardly surprising, since Xe1 has by far the highest Xe affinity[87] among the four Xe sites and thus is expected to be the least polar site. For the other sites, the combined 10 ± 5% Xe occupancy is probably not sufficient to produce a detectable reduction of the MRD amplitude even if all these sites were fully occupied by water in the absence of Xe. If one Xe atom displaces one water molecule, MRD component 1 would be reduced by merely 1.5% (0.1/6.7), comparable to the measurement error.

Figure 5

2H MRD profile, R1(ω0), from immobilized MbCO (25 °C, pD 7.16, scaled to NW = 3 × 104) in the presence (blue symbols) or absence (red symbols) of 8 bar of xenon. The error bars do not exceed the symbol size. The parameter values obtained from the fits (solid red and dashed blue curves) are collected in Table S5 (Supporting Information).

All high-resolution crystal structures of equine and sperm whale Mb, in the CO, deoxy, or met forms,[4−9] identify two internal hydration sites in the Xe3 cavity, with full (equine Mb) or nearly full (sperm whale Mb) occupancy (Figures 1 and S6, Tables S6 and S7, Supporting Information). The water molecule in subsite Xe3a has a small B factor and is engaged in four H-bonds, donating to Ile-75 O and His-82 Nδ1 and accepting from Gly-80 N and the adjacent water molecule. The water molecule in subsite Xe3b has a much larger B factor and makes only one H-bond with the protein (Ala-134 O). In addition, all but one of the high-resolution crystal structures listed in Tables S6 and S7 (Supporting Information) locate an internal water molecule in the so-called apical site (Figures 1 and S6, Tables S6 and S7, Supporting Information). The apical water makes three H-bonds with the protein. The crystal structures of equine Mb (but not sperm whale Mb) in the CO and deoxy (but not met) forms also locate a water molecule in the “bottom” site (Figures 1 and S6, Tables S6 and S7, Supporting Information). The bottom water also makes three H-bonds with the protein, but it is less deeply buried than the apical water.

It is tempting to identify these four crystallographic water molecules with the 4.5 ± 1.0 most long-lived (5 μs) water molecules indicated by the MRD data. The Xe difference MRD experiment does not rule out the presence of two water molecules in the Xe3 cavity because Xe appears to displace only one of these water molecules (at least for sperm whale Mb; see PDB entries 1J52 and 2W6Y), and then, as noted above, the effect would be too small to detect. Furthermore, the 10 ± 5% combined Xe occupancy is presumably associated mainly with the less polar Xe2 and Xe4 cavities. The 32 ns water inferred from the 17O MRD data might be located in a deep surface pocket. Alternatively, it could refer to one of the two least deeply buried among the four internal hydration sites, the bottom site and subsite Xe3b. In that case, another internal hydration site, not identified in the crystal structures, must contribute to the 5 μs MRD component. For example, the phantom (Ph) cavity (Figure 1) is large enough to accommodate four (mutually H-bonded) water molecules, but this cavity is empty in all the crystal structures of equine Mb.

The two MD simulation studies that have focused on internal water in Mb both find multiple partially occupied hydration sites[44,51] most of which have not been identified in crystal structures. The first of these studies,[44] a conventional MD simulation of equine met-Mb at 300 K, found, in addition to the Fe3+-bound water in the DP, eight internal hydration sites with a combined occupancy of 2.75. The apical site (occupancy 0.8) was identified but not the bottom site. Both the Xe3 and Xe4 cavities were assigned an occupancy of 0.3. Except for the iron-bound water, all internal water molecules identified from this 48 ns MD trajectory exchanged on a nanosecond time scale. Therefore, this MD study supports neither the crystallographically inferred locations and occupancies of internal hydration sites in Mb nor the MRD-inferred microsecond time scale of internal water exchange.

The second simulation study[51] examined sperm whale deoxy-Mb at 310 K and used temperature-accelerated MD for enhanced sampling. Numerous internal hydration sites (but no occupancies) were reported, including the Xe2, Xe3, and Xe4 (but not Xe1) cavities. On the basis of computed minimum free energy pathways for water migration, it was suggested that all internal water molecules enter the protein via the His-64 gate and the DP, from where water molecules can migrate further to the Xe cavities. With the aid of transition state theory, the time scale for water escape from the DP was estimated at 300 ns (at 310 K), an order of magnitude faster than the 5 μs MST deduced from our MRD data.

In summary, crystallography[4−9] supports a scenario where most or all of the ∼4 long-lived internal water molecules inferred by MRD reside in small polar cavities (Xe3, apical, and bottom sites). A larger number of partially occupied sites, as suggested by MD simulations, is not likely to account for the MRD data. X-ray crystallography may miss hydration sites with low occupancy and/or positional disorder, but water molecules in less polar sites (say, with less than three H-bonds) are expected to have a small order parameter and would then not contribute much to the MRD profile.

Mechanism of Internal-Water Exchange

The common MST, τS = 5.6 ± 0.5 μs, of the ∼4 water internal molecules inferred from our MRD data suggests that they exchange with external solvent via a common mechanism. Moreover, since the likely internal hydration sites are located far apart in the Mb molecule (the apical, bottom, and Xe3 sites are separated by 17, 27, and 31 Å), the common exchange mechanism must be global in spatial extent.

In an effort to gain more insight into the exchange mechanism, we studied the temperature dependence of the 2H MRD profile from immobilized MbCO. Irreversible changes in the MRD profile, possibly caused by gel restructuring, were noted above 35 °C. Therefore, the analysis was restricted to the rather narrow temperature range 5–25 °C, where fits like the one in Figure 2a yielded a surprisingly small Arrhenius activation energy of 17.4 ± 0.4 kJ mol–1 for the MST τS,1 at pD 7.0 (Figure S7, Supporting Information). Such a small activation energy might suggest that water exchange is rate-limited by an entropic bottleneck. The narrow temperature range and the strong covariance between N1 and τS,1 (section 3.2), however, preclude an unambiguous interpretation of the temperature dependence. The MSTs used for the Arrhenius fit in Figure S7 (Supporting Information) were obtained from constrained fits where τS,1 was the only adjustable parameter for component 1 even though the relative contributions of internal water and LDs to component 1 are expected to vary with temperature.

The ∼5 μs MST deduced from our MRD analysis implies that internal-water exchange in Mb is a rare event that has not yet been captured by MD simulations, which so far do not extend beyond 100 ns for Mb. Neither the average protein structure seen in crystal structures nor the fluctuations accessed by sub-μs simulations tell us much about the transient structures involved in the water exchange event. An ultralong MD simulation[88] was recently used to demonstrate that the internal water molecules in the protein BPTI exchange by a transient aqueduct mechanism,[69] where single-file water chains penetrate the protein through transiently formed tunnels or pores. It is conceivable that a similar mechanism operates in Mb, with the difference that the transient water-filled tunnels are partly built from preexisting cavities. Whereas a single water molecule is unlikely to migrate from a polar site into a largely non-polar cavity network,[55] an H-bonded water chain can move rapidly through a transiently formed tunnel.[54] The peripherally located internal hydration sites might then act both as portals for solvent penetration and as seeds for transient water chains.

MD simulations indicate that the ligands can migrate from the DP into the extensive network of permanent and transient cavities that permeates the Mb molecule.[29,33,35,42,46] The ligand migration pathways seen in these simulations include the apical and Xe3 hydration sites.[32,35] A transient aqueduct mechanism might act to intermittently “flush” the cavity network, thereby removing internally trapped ligands. Experiments[15,18,23,30,34,39] and simulations[35] indicate that photodissociated CO can escape from the DP via the His-64 gate on a 100 ns time scale. The half-life of CO that has migrated from the DP to the Xe1 cavity was, however, reported to be ∼10 μs (at 20 °C).[19] Since this process occurs on the same time scale as internal-water exchange, it may, in fact, be concomitant with the proposed “flushing” of the cavity system by transient water chains.

Conclusions

On the basis of 2H and 17O MRD measurements on gel and solution samples, we have arrived at the following answers to the first three of the five questions posed in section 1:In addition, the MRD results indicate that 2 or 3 of the 11 His residues of equine Mb undergo intramolecular hydrogen exchange on a μs time scale.

Under physiological conditions, equine MbCO contains 4.5 ± 1.0 internal water molecules that exchange on a μs time scale. The root-mean-square orientational order parameter of these water molecules is 0.8. Another ordered water molecule, possibly located in a deep surface pocket, exchanges on a time scale of 30 ns.

The likely locations of the ∼4 long-lived water molecules are the crystallographically identified Xe3, apical, and bottom hydration sites. In agreement with crystallography and MD simulations, the MRD results indicate that the Xe1 cavity is devoid of water.

Despite being located as much as 30 Å apart, the ∼4 long-lived internal water molecules have the same (or very similar) MST, τS = 5.6 ± 0.5 μs.

As regards the last two questions in the Introduction, we speculate the following:

The internal water molecules exchange by a global mechanism, that may involve penetration of the protein by transient H-bonded water chains, entering the protein at the peripherally located permanent hydration sites.

The proposed water exchange mechanism suggests a functional role, where intermittent “flushing” of the cavity system removes trapped ligands.