INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

There exists ample experimental and theoretical evidence that a protein in an aqueous solution modifies both structure and dynamics of water around it (Teeter, 1991; Brunne et al., 1993; Lounnas and Pettitt, 1994a,b; Schoenborn et al., 1995; Phillips and Pettitt, 1995; Denisov and Halle, 1996; Burling et al., 1996). Previously, we explored the relationship between the solvent structure and dynamics at the protein-water interface in terms of diffusion rates and radial distribution functions (Makarov et al., 1998b).

Here we intend to study the protein-water interface from a different point of view, concentrating on those water molecules that are closely associated with the protein so that their motion can no longer be described simply as diffusive. The concept of hydration sites (Lounnas and Pettitt, 1994a,b; Hummer et al., 1995), as opposed to crystal waters of hydration, appears to be useful in this respect.

Hydration sites appear as local maxima in the time-averaged singlet solvent density (Lounnas and Pettitt, 1994a,b; Scanlon and Eisenberg, 1975, 1981). This definition immediately invites a comparison to the solvent structure observed in x-ray (Scanlon and Eisenberg, 1975, 1981; Lounnas and Pettitt, 1994a,b; Hummer et al., 1995) or neutron diffraction (Gu and Schoenborn, 1995) experiments. Although this comparison is useful for verification of the overall structure of the hydration shell observed in the simulation, a direct comparison can be problematic. There are two reasons for this. First, the simulated hydration sites, unlike experimental hydration sites, are not subject to normal refinement constraints, as they are simply points in space and not individual molecules. The conditions used to define the simulated sites may be chosen so as to best reproduce the complicated three-dimensional solvent distribution by a set of local maxima, rather than to comply with steric restrictions and hydrogen bonding requirements of individual molecules. Second, there are no crystal packing effects in the simulation, because we simulate a solution, albeit a very concentrated one. Therefore, we intend to use the hydration site concept not just to provide direct comparisons to crystal structures, but also as a tool for the analysis of solvent dynamics at the interface with the protein. Temporal ordering of solvent in the hydration shell of the protein may be characterized by the population analysis of hydration sites (Lounnas and Pettitt, 1994b; Gu and Schoenborn, 1995), from which the average water residence times may be derived (Brunne et al., 1993; Lounnas and Pettitt, 1994b; Schoenborn et al., 1995; Abseher et al., 1996; Rocchi et al., 1997). Hydration sites are ideal as volume elements for the molecular residence time calculation. Using the hydration sites, as opposed to shells (Abseher et al., 1996) or irregular volume elements proximal to specific solute residues (Brunne et al., 1993; Garcia and Stiller, 1993; Kovacs et al., 1997; Rocchi et al., 1997), allows us to relate the temporal order parameters (residence times) to the spatial structure of solvent in the most natural way.

In turn, this comparison enables us to compare our analysis to that of x-ray crystallography or nuclear magnetic resonance (NMR), even though the precise definitions of our densities and residence time parameters may differ in detail from those used in these experimental fields. Last, but not least, hydration sites offer a significant technical advantage over other definitions of volumes for the residence time calculation. All sites can be forced to be of the same shape and size, which makes all occupancy and temporal parameters consistent and enables a direct comparison of different sites.

Specifically, we seek to understand the correlation between the spatial and temporal structure of the protein hydration shells. We wish to quantify the extent to which it is possible to characterize the mobility of the interfacial water in terms of specific residence times and occupancies, as is commonly done both in simulation (Lounnas and Pettitt, 1994a,b; Gu and Schoenborn, 1995; Brunne et al., 1993; Schoenborn et al., 1995; Abseher et al., 1996; Rocchi et al., 1997) and experimental studies (Brunne et al., 1993; Schoenborn et al., 1995). In this respect, we wish to determine the factors that affect the lifetimes of water molecules on the protein surface.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Molecular dynamics (MD) simulation

The analysis presented in this paper is based upon a MD simulation of sperm whale myoglobin that has been reported in detail in a separate publication (Andrews et al., 1998). Therefore, only the brief account of the simulation setup will be presented here. The all-atom CHARMM-23 parameter set (MacKerell et al., 1992) was used to model a single molecule of carbomonoxy myoglobin (Protein Data Bank entry 2mgk, Bernstein et al., 1977; crystal structure by Quillin et al., 1993) solvated by 3717 flexible TIP3P (Jorgensen et al., 1983) water molecules in a 60.4 Å × 54.7 Å × 40.7 Å box under periodic boundary conditions. The system was prepared through a series of energy minimization and heating steps. The heating stage was followed by equilibration at 294 K for 200 ps. Equations of motion were integrated using a 0.5-fs time step without constraints. Electrostatic interactions were treated with a 13 Å cutoff and a potential-based switching function beginning at 10 Å. The trajectory was continued for a total of 1.1 ns, of which the last 900 ps were chosen for analysis. The trajectory was recorded with 0.1-ps intervals.

Calculation of solvent density distribution

This part of our analysis followed procedures previously developed in this lab (Lounnas and Pettitt, 1994a,b; Lounnas et al., 1994; Rudnicki and Pettitt, 1997; Makarov et al., 1998a). For each step of the MD trajectory, the protein was fitted to a consistent frame of reference. Next, the same transformation was applied to the water molecule coordinates, taking the periodic boundaries into account. The coordinates of the water oxygen atoms were then mapped onto the three-dimensional rectangular grid with a 0.5 Å grid step, producing an average three-dimensional number density distribution. The particular choice of the grid step is a compromise between the uncertainty in location of the density features and the statistical error in the local density value that arises due to a lower number of counts in each grid cell. At the chosen grid step every cell in the regions corresponding to bulk solvent would have at least 50 counts over the entire trajectory. The density map was smoothed by averaging the value of each cell with six of its nearest neighbors before further manipulations.

Localization of hydration sites

For the purposes of this work, protein hydration sites are defined as local maxima in the water oxygen number density map that satisfy certain conditions. They should be no farther than 5 Å from any protein atom, no closer than at least 1 Å (two grid steps) from each other, and have a maximum density value no lower than 150% of the bulk water density in the simulation. For the normal water density of 1 g/mL, the last cutoff is equivalent to 0.05 particles/ų. The maxima were determined by comparison of the density in the current grid cell to each of its six nearest neighbors (the 7-point rule). If close peaks were found, their coordinates were weighted by their respective maximum density values and averaged together, resulting in a single peak.

Comparison with crystal water positions

To compare the hydration peaks from the above simulation with the positions of the crystal waters of hydration, we have selected a number of sperm whale myoglobin structures, all at the resolution of 2 Å or better, from the Protein Data Bank (Bernstein et al., 1977). There were 73 sperm whale myoglobin structures in the database as of March 1998. Every protein was oriented in space to fit the reference structure, and the positions of any hydration water molecules were then saved and compared with those from all other structures on the list. Sometimes this comparison would indicate a trivial dependency between a group of structures. In such cases, only one structure from the group was taken into consideration. The purpose of this elimination procedure was to avoid potential bias due to identical water positions that might have been copied from a different myoglobin structure in a molecular replacement during crystallographic refinement. This procedure resulted in the following selection of structure files: 1abs, 1bvc, 1bvd, 1fcs, 1mbc, 1mbd, 1mbi, 1mcy, 1mgn, 1mlf, 1moa, 1mti, 1mym, 1spe, 1swm, 1vxa, 1vxb, 1vxc, 1vxf, 1yog, 1yoi, 2cmm, 2mb5, 2mbw, 2mgb, 2mgc, 2mgd, 2mge, 2mgf, 2mgg, 2mgh, 2mgi, 2mgj, 2mgk, 2mgl, 2mgm, 2mya, 2spl, and 4mbn. We then compared the simulated hydration sites to a collection of crystallographic water molecules contained in the above list. Quantitatively, the agreement between the theoretical and experimental water positions is measured by the proximal radial distribution of the crystal waters around the simulated hydration sites (Hummer et al., 1995).

Calculation of water residence times

We note that the definition of a "residence time" varies widely in the available simulation and NMR literature. Some of these definitions are discussed in a review by Schoenborn et al. (1995). The method used here is that of the time correlation function (Brunne et al., 1993), although the correlation function itself is modified and we offer a different interpretation. The correlation function in question is simply

(1)

where the function (x  y) takes the value of 1, when x = y (x  0 and y  0), or 0, when either x  y or both x and y = 0, meaning that the site is not occupied. N(t) is the index of the water molecule found in the hydration site at time t. The site was considered occupied if a water molecule was found in a spherical volume with radius of 1.5 Å (roughly equal to the van der Waals radius of water) centered on the coordinates of this site. The resulting time correlation function is usually fitted by a single exponential (Brunne et al., 1993), but we have noticed that the single exponential fit is often poor. Using a double exponential model resulted in a more accurate fit in all cases. The two new rate constants in case of solvent diffusion correspond to solvent molecules that stay in the hydration shell for prolonged periods of time or enter and then immediately leave. Both of these types of water behavior are observed in the simulation, and it is important to distinguish between them. In fact, it has been suggested (Garcia and Stiller, 1993; Rocchi et al., 1997) that in some cases, groups of water molecules may show different residence times with respect to the same atomic site. The process of solvent diffusion is clearly much more complex than that described by either single or double exponential model and probably should be characterized by a distribution of rate constants. Although a Kohlrausch-Williams-Watts stretched exponential is often used to describe phenomena governed by multiple relaxation rates (Abseher et al., 1996), its interpretation is not so straightforward as that of a simple biexponential model. The model used here is

(2)

where W₀ is the average site occupancy, defined as the fraction of time the site was occupied by any water molecule, and k₁, (1  w) and k₂, w are the rate constants and the weights for the fast and the slow components, respectively. The short and the long residence times are then defined simply as ₁ = 1/k₁ and ₂ = 1/k₂.

For comparison, we have also repeated this calculation for sites in the bulk region of solvent in the simulation. There are 8 such sites located in the corners of the rectangular simulation box. The occupancy and residence time parameters for these bulk sites were then averaged together.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Hydration site location and comparison with crystal data

Using the procedures described above, we have identified a total of 294 hydration sites around myoglobin. The number of hydration sites and their position in relation to the protein appear to be rather insensitive to the parameters of the localization procedure. For instance, in a lower resolution water oxygen number density map (at 1 Å grid step) the algorithm yields 288 sites at approximately same positions (there is, of course, a small error associated with the larger grid step). Thus, we conclude that our solvation sites are sufficiently reliable. For convenience, the sites are numbered in descending order depending on their maximum density value.

The spatial distribution of hydration sites is both highly anisotropic and asymmetric. In general, the hydration shell follows the shape of the protein, but there are regions where the hydration sites are tightly clustered and regions that are almost devoid of sites. In Fig. 1, the high density region is located on top of the protein, and the low density areas without hydration sites are to the right. The protein surface residues in the areas without high density sites (highlighted in yellow in Fig. 1) do not fall into any one particular residue class. Both charged and polar residues are found in that region along with the hydrophobic residues (Table 1). In contrast, the distribution of crystal hydration waters (Fig. 2) is much more uniform and contains few excessively populated or empty regions. To provide a quantitative measure for the comparison of crystal water positions to the MD hydration sites, we have constructed and displayed in Fig. 3 a cumulative distribution of nearest-neighbor distances between crystal waters and hydration sites (Hummer et al., 1995). The histogram achieves saturation after approximately 3 Å, meaning that even the worst positioned MD hydration sites would have at least one corresponding crystal water molecule at this distance. From this point of view, agreement of our data with the experiment is better than that found in the simulation used by Lounnas and Pettitt (1994a). On the other hand, even better reproduction of the experimental data is achieved in simulations of small hydrated molecules, such as nucleic acid fragments (Hummer et al., 1995; Feig and Pettitt, manuscript in preparation). The most probable reason for the discrepancy between the simulated and crystal hydration sites is the absence of the crystal packing effects in the simulation of the protein in solution.

View larger version (91K): [in this window] [in a new window]   FIGURE 1   Three-dimensional distribution of high-density solvent peaks around myoglobin. Exposed protein residues that do not have solvent peaks associated with them are colored yellow. These residues form the bald spots on the surface of the protein. The heme is in red.

View larger version (85K): [in this window] [in a new window]   FIGURE 2   A view of the myoglobin hydration sites generated in the MD simulation (yellow dots) compared to the collection of the hydration sites contained in the Protein Data Bank (blue dots).

View larger version (26K): [in this window] [in a new window]   FIGURE 3   Cumulative histogram of nearest neighbor distances between the MD and the Protein Data Bank water sites.

The simulated hydration sites form an elaborate hydrogen bonding network that encircles the protein and involves participation of many of its surface-accessible residues. A small fragment of this network is shown in Fig. 4. This polygonal structure is located in the high density region and is actually a part of an even more complex network of sites. Similar arrangements have been observed in the earlier study of myoglobin hydration (Lounnas and Pettitt, 1994a,b) and in high-resolution crystal structures (Teeter, 1991). Fig. 5 shows the site-site radial distribution function (RDF). It has a broad peak centered at 2.6 Å (the histogram step is 0.2 Å), which corresponds to water-to-water hydrogen bonding distance. The peaks at 1 Å and 1.4 Å are induced by the grid. In the previous study of myoglobin hydration sites (Lounnas and Pettitt, 1994a), the site-site RDF had a peak at 2.1 Å. Although the hydration sites do not represent individual water molecules and thus are not subject to usual steric constraints, our present site-site RDF is more consistent with chemical intuition and data from protein crystallography (Schoenborn et al., 1995).

View larger version (39K): [in this window] [in a new window]   FIGURE 4   A fragment of the hydrogen-bonding network of hydration sites in the vicinity of site number 2 and Gln-129. Only the participating protein residues and the highest density solvent sites are labeled. The numbers of lower-density sites are omitted for clarity. Hydration sites were treated as atoms that may participate in a hydrogen bond as both a donor and an acceptor at the same time. The definition of a hydrogen bond followed a general convention that requires a distance between donor and acceptor no more than 3.3 Å.

View larger version (30K): [in this window] [in a new window]   FIGURE 5   Site-site radial distribution function.

Temporal ordering of water around the protein

Population analysis of the hydration sites shows that all of them are never occupied simultaneously. In fact, the number of sites occupied at the same time never exceeds 164 out of 294, and there are no cases of double occupancy. Fig. 6 shows the population time correlation functions for two distinct hydration sites (numbers 20 and 36) along with that for the bulk solvent. We begin with the description of these bulk sites to establish a baseline for future comparisons. Average occupancy (W₀) in the bulk sites, which are positioned in the corners of the simulation box as described above, is 49%, and the residence time parameters are t₁ = 0.34 ps, t₂ = 4.1 ps, and w = 0.39. We note that both bulk residence times are short and the weight of the slow component (w) is low.

View larger version (13K): [in this window] [in a new window]   FIGURE 6   Examples of the population time correlation functions for different hydration sites: solid line, site number 20; dotted line, site 36; dashed line, bulk water. Thin solid lines represent the double exponential fits (Eq. 2) to the correlation functions. The corresponding fit parameters are t1(20) = 0.35 ps, t2(20) = 15.0 ps, w(20) = 0.55 and t1(36) = 0.6 ps, t2(36) = 54.6 ps, w(36) = 0.64.

Distributions of the site occupancy, t₁, t₂, and w for all 294 sites are shown in Fig. 7. The distributions of the residence times and occupancy peak at their respective bulk values, which indicates that many of the hydration sites are rather bulk-like. Typical residence times (t₂) vary between 5 and 35 ps, and the average site occupancy varies from 15% (site 113) to 59% (site 2). Long tails on the right hand side of the residence time distributions indicate that the protein may prolong the life of the bound water molecules by up to one order of magnitude, but such sites are few. A residence time on the order of hundreds of picoseconds in a 1-ns simulation is likely to be caused by a single visit of a water molecule to the particular site, occurring only once over the course of the simulation, and this may make the calculated residence time statistically unreliable. Nonetheless, examples of these long-lived sites are listed in Table 2. The influence of the protein on the surrounding solvent is manifested to an even greater degree in the distribution of the slow component weight parameter, w (Fig. 7 d). It is shifted to the right of its respective bulk value in agreement with the fast component dominating the motion in the bulk. In terms of water mobility this means that there are two alternative scenarios possible for every water molecule that enters a site: one may stay there for a time period of t₂, on the average, or enter and immediately exit with a characteristic short residence time t₁. The w parameter may be thought of as simply a fraction of water molecules that follow one of these two alternative scenarios. Thus, the distribution of residence times at the interface as opposed to bulk appears to be skewed toward the longer components.

View larger version (20K): [in this window] [in a new window]   FIGURE 7   Distributions of the various parameters of the population time correlation functions fitted with the double exponential. (a) Long residence time t2. (b) Short residence time t1. (c) Site occupancy. (d) w parameter. In each plot, the thin vertical line denotes the corresponding bulk value.

Our present observations with regard to the water residence times agree very well with available experimental and simulation data. Short residence times for water molecules in protein hydration shells were observed in many other simulations, without respect to the nature of the protein studied or to the exact method for the calculation of the residence times (Abseher et al., 1996; Brunne et al., 1993; Garcia and Stiller, 1993; Kovacs et al., 1997; Lounnas and Pettitt, 1994b; Muegge and Knapp, 1995; Rocchi et al., 1997). Water residence times on a 10- to 50-ps time scale were also observed in protein NMR experiments (Brunne et al., 1993; Denisov and Halle, 1996). Long-tail residence time distributions with a peak below 10 ps were found in the analysis of a MD simulation of hydrated ubiquitin (Abseher et al., 1996). These distributions are of essentially the same shape as those computed here. It has been suggested (Abseher et al., 1996) that the existence of multiple relaxation rates and residence time distributions in the protein hydration shell is due to the structural and chemical heterogeneity of the protein surface.

In this regard, it is important to note that there is no simple correlation between the maximum site density, occupancy, and residence times. Table 3 illustrates the absence of a relationship between the site maximum density and the temporal order parameters. Further, lack of correlation between density and residence time is evident from Table 2. The sites are numbered depending on their maximum density, starting with the highest density site (number 1), but the sites with longest residence times are 42, 1, 214, and 6. 

This lack of correlation between the spatial (density) and temporal (residence times) order in the surrounding solvent is hardly surprising, given the observations in NMR experiments. A vast majority of water molecules that appear ordered in protein crystal structures cannot be resolved in NMR, partly because of the lack of crystal packing in solution, but also due to their rapid exchange with bulk solvent (Denisov and Halle, 1996). On the other hand, many protein crystal structures appear to have rather large empty internal cavities, where the presence of positionally disordered water is indicated only by NMR dispersion experiments (Denisov and Halle, 1996).

Yet another comparison is possible to the xenon binding sites observed experimentally in the structure of metmyoglobin (Tilton et al., 1984). Of the four xenon binding sites, we observe water in three, namely the two xenon sites in the proximal heme pocket (overlapping hydration sites 113, 115, and 19) and the xenon site in the EF/H helix corner. The latter xenon binding site overlaps with many hydration sites, of which the site number 34 has the longest residence time. The last xenon binding site, on the distal side of the myoglobin molecule, apparently does not contain spatially or temporally ordered water. This may be correlated to its low occupancy by xenon in experiment, which is the lowest of all four xenon sites at 0.46 (Tilton et al., 1984).

Because there is no obvious relation between the number density in the hydration shell and the water residence times, one may rightfully question what factors determine the residence times. At least three potential answers have been offered in the literature: (i) the chemical nature of the protein residues close to the hydration shells for which the residence times are measured (Schoenborn et al., 1995); (ii) number of hydrogen bonding opportunities available at a given site (Denisov and Halle, 1996); and (iii) the local geometry of the protein surface, i.e., whether the site is buried or exposed to bulk (Denisov and Halle, 1996).

The first hypothesis has been investigated in many protein simulation studies (Brunne et al., 1993; Garcia and Stiller, 1993; Muegge and Knapp, 1995; Kovacs et al., 1997; Rocchi et al., 1997) that attempted to produce ranking relations for the residence times of water molecules in hydration shells of charged, polar and non-polar protein residues. The results of these studies have been inconclusive or conflicting (Schoenborn et al., 1995; Denisov and Halle, 1996). Our own attempt in this direction was also unsuccessful.

The second hypothesis has been dismissed by Denisov and Halle (1996) upon examination of the hydration waters in various protein crystal structures. While ordered water molecules make, on average, 2 to 4 hydrogen bonds to the protein, this number is comparable to the average number of hydrogen bonds per water molecule in the bulk (3 to 4), and nearly all of this ordered water is short-lived (Denisov and Halle, 1996).

The third hypothesis may be tested immediately by mapping the sites with different residence times onto the protein structure. Without exception, all sites with very long residence times (defined here as t₂  80 ps) are trapped either in the cavities inside the protein or in the grooves and concave regions, whereas the low residence time sites (t₂  10 ps) are in convex areas and appear to be only loosely associated with the protein (Fig. 8). The opposite is true for the occupancy: high occupancy sites (W₀  0.5) are exposed, and low occupancy sites (W₀  0.3) are more buried. Interestingly, all sites with a high slow exponent weight (w  0.6) are also buried, and all sites with a low w (w  0.4) are exposed. Although the three classes of buried sites (high t₂, low W₀, and high w) do overlap, they are not identical, with only 3 sites in common.

View larger version (90K): [in this window] [in a new window]   FIGURE 8   Spatial distributions of myoglobin hydration sites. (a) Sites with high residence time t2. (b) Sites with low residence time t2. (c) Sites with high occupancy. (d) Sites with low occupancy. (e) Sites with high w parameter (f) Sites with low w parameter. The gray dotted area represents the molecular surface of the protein.

Therefore, following numerous experimental observations (reviewed by Schoenborn et al., 1995; Denisov and Halle, 1996), we find that the water residence times are determined almost exclusively by the position of the water molecules on the protein surface: cleft and buried waters have significantly longer lifetimes than those that are exposed. Although other factors mentioned above may also affect the residence times, they appear to be secondary to the geometry with respect to local topography. We hypothesize that these factors do influence the residence time of the interfacial water, but on different time scale. For instance, an ability to form hydrogen bonds to the protein or to the neighboring sites may determine the length of the water residence in a surface site with t₂ 50 ps, but it will not be as important in a cleft site with t₂ on the order of 500 ps. This would explain the failure of the ranking relation analysis mentioned above. Additional analysis, preferably using even longer and more accurate computer simulations, should shed more light on this issue.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

In conclusion, we confirm that there is no direct correlation between the spatial (number density) and temporal (residence time) order of solvent in the protein hydration shell. We find that there is no single specific residence time for any particular location in the protein-solvent interface. Instead, there exists a distribution of residence times, that may have both extremely short and extremely long components. In our study, the long and the short residence times t₂ and t₁ play this role. Close to the protein surface the residence time distribution is biased toward the longer components. The residence times of the interfacial water molecules appear to depend most strongly upon the degree of exposure of the hydration site to bulk solvent. Buried and concave sites have much longer residence times than those that are convex and exposed to bulk. In addition, hydrogen bonding opportunities and the hydrophobic character of the underlying protein residues may influence the water residence times, but these factors only appear to be important for the exposed sites with inherently short residence times.

Supplementary material

Hydration site coordinates and their occupancy and residence time parameters are available in electronic form on request from the authors.