J. Phys. Chem. 1993,97, 11339-1 1343
11339
Heterogeneous Diffusion of Water at Protein Surfaces: Application to BPTI E. W.Knapp' and I. Muegge Fachbereich Chemie, Institut f k r Kristallographie, Freie Universitiit Berlin, Takustrasse 6, D 141 95 Berlin, Germany Received: June 8, 1993; In Final Form: August 13, 19938
The influence of the protein surface on the dynamics of water is investigated by a computer simulation of the dynamics of bovine pancreatic trypsin inhibitor (BPTI) in water for 500 ps. The diffusion constant for water moving away from the protein surface has been determined by comparing computer simulation data with those of an analytical model. For 75% of the surface waters, the diffusion off the protein surface is a factor of 4 slower, but for the other 25% the diffusion is a factor of 15 faster than in bulk water.
1. Introduction
diffusion model to the data of computer simulation is performed in section 5. A summary is given in section 6.
It is now generally agreed that water in protein-water systems plays more than a passive role.' Often dry proteins do not function properly. Dry proteins differ also in their dynamical behavior as compared to their wet counterparts.2 Water is necessary for folding into the native conformation.' To investigate structure and dynamics of water in proteinwater systems in detail, coordinates of the water molecules are necessary. Coordinates of protein-water systems are usually obtained by X-ray structure analysis. For several years, protein structure determination has also been performed by NMR and neutron scattering data.&' With a resolution better than 2.0 A, the oxygen positions of ordered crystal water can be determined by X-ray structure analysis of protein crystals. Only recently it has become possible to distinguish surface water and internal water of protein-water systems by NMR method^.^.^ Ordered water molecules in protein crystals are immobile. However, in solution, practically all surface waters of proteinwater systems are mobile on a time scale of several hundreds of picoseconds.8 The purpose of the present contribution is to investigatethe mobility ofwater at the protein surface by computer simulation of the dynamics of a protein-water system. Bovine pancreatic trypsineinhibitor (BPTI) is a suitablemodel system to investigate the dynamics of water on protein surfaces. BPTI is a typical globular protein with one a-helix and an
[email protected] is small enough to do computer simulations of long-time dynamics. Several high-resolution X-ray structures with up to 60 crystals waters are available.lOJ1 Only very few crystal waters from different structures are equivalent. Residence times of specific surface waters are known from NMR data.899 Diffusion of surface water at BPTI has been studied several times by computer simulation of the dynamics.12-lS Since water molecules are moving away from the protein surface soon, the diffusion constant was evaluated by analyzing residence times or a small number of selected water trajectories only. So far, there has been no possibility to discriminate between diffusion along the protein surface and diffusion away from the protein surface. In the present treatment, diffusion away from the protein surface is considered by comparing simulation data of the BPTI-water system with those of an analytically solvable model. This comparison provides reliable estimates of the diffusion constant and is quite generally applicable for protein-water systems. Section 2 provides details of the computer simulation; section 3 gives a survey of simulation data. The diffusion model is introduced in section 4 and the Appendix. The application of the Abatract published in Adounce ACS Abstructs. October 1, 1993.
0022-3654/93/2097-11339$04.00/0
2. Computational Methods The dynamics of BPTI in water is simulated by solving the classical equations of motion for all atoms except nonpolar hydrogens by using the program CHARMm.16 Nonpolar hydrogen atoms are accounted for by suitable defined extended heavy atoms, resulting in 568 atoms for BPTI. With charged residues Arg+, Lys+, Glu-, and Asp-, the total charge of BPTI amounts to +6e. The TIP3P water model of Jorgensen17 is used. A cutoff radius of 8 A with shift boundary conditions16 for the electrostatic interactions is used. The initial coordinates for the computer simulation are taken from the X-ray structure of BPTI.l0J1 To saturate the protein with water, the solvent of the protein-water system is exchanged several times by using an overlay technique,'* which allows the placement of water molecules inside the protein matrix. Nonetheless, besides the four water molecules which are known to,be inside the protein, no other water could be placed. After this procedure, the protein-water system is placed in a water sphere whose radius is adjusted to 22.5 A to provide the proper bulk water density. The boundary of the sphere is given by a mean force field potential.19 There are 1319 water molecules in the system. Water molecules in the outer shell of 2.5-A thickness are coupled to a heat bath20 at 300 K with stochastic dynamics setting the friction constant at 50 ps-l. This coupling is used for a gentle heating of the whole system, which is completed within 30 ps. The hydrogen bond length is fixed by using the SHAKE algorithm.21 3. Survey of Simulation Data
Thediffusion constant of theTIP3P water model17is evaluated by applying the Einstein relation D = (?)/6At for a water box of length 31.10 A with periodic boundary conditions containing 1000 water molecules. A 30-ps simulation of the dynamics of water from which 10 ps are used for analysis is sufficient to obtain reliable results. In order to calculate the diffusion constant by interpolation at an arbitrary temperature, the simulation of thedynamics of the water box has been performed at four different temperatures, 286 f 4, 298 f 4, 308 f 4, and 320 f 4 K. To check the spherical model of diffusion also for the box geometry, a longer simulation of 100-ps dynamics was made at 308 K. At the reference temperature T = 294 K of the computer simulation of the BPTI-water system, the interpolatedvalue of thediffusion constant is D = 0.20 A2/ps as compared to the experimental value of Dcxp= 0.21 A2/ps.21 0 1993 American Chemical Society
Knapp and Muegge
11340 The Journal of Physical Chemistry, Vol. 97, No. 43, 1993
0.05
'
0.00 0.0
100.0
200.0
300.0
400.0
I
500.0
time (ps)
Figure 1. Time dependenceof rmsdeviationsbetweenthe X-ray structure of BPTI and the instantaneous structure from computer simulations of the dynamics of BPTI. The upper curve refers to the average of all non-hydrogen atoms of BPTI, and the lower curve refers to the average of protein backbone atoms only.
To demonstrate the usefulness of the spherical model for diffusion, a computer simulation of 200-ps dynamics of TIP3P17 water at 294 K in a sphere of radius 16.0 A is performed. Also here the walls of the sphere are modeled by a mean force field potential19 and stochastic dynamics with a friction constant of 50 ps-1 are applied to water molecules within 2.5 A of the wall of the sphere.20 To avoid a slowing down of the diffusion process by Langevin dynamics, only 13%of the water oxygens are coupled to the heat bath in the wall region. Computer simulations of the dynamics of water have demonstrated that the water dipoles are not isotropically oriented if spherical boundary conditions are used.23J4 In the present work it is shown that the translational diffusion constant is independent of the boundary conditions (spherical or periodic). For the BPTI-water system a total 500 ps have been simulated.25 Sixty picoseconds are used for heating and equilibrating of the system. The remaining 440 ps are suitable for analysis. Only 7% of water oxygens in the 2.5-A wall region are subject to the stochastic dynamics to avoid slowing down of the diffusion process. Every 0.2 ps, atomic coordinates are stored. The time average rms deviation of the BPTI structure from the X-ray structure is 1.9 A for all 455 non-hydrogen atoms of BPTI and 1.4 A if only main chain atoms are considered. These values are similar to rms deviations obtained with other force fields and for other protein-water systems. The time evolution of the rms deviation exhibitsavery steep increasein the first few picoseconds, followed by a gradual increase which merges into a plateau value after 200 ps (Figure 1). Smaller rms deviations are obtained for shorter simulation times.13J4 Water molecules may stay in the bulk, at the protein surface or in the inner part of the protein. In agreement with NMR data: the four internal waters 11 1, 112, 1 13, and 122 common to all three X-ray structuresIOJ1 stick to their places during the whole simulation (Figure 2 (top)). The majority of water molecules stays in the bulk during the whole simulation, as for instance water E393 (Figure 2 (bottom)), a water which has been added by theoverlay technique. Waters 121 and 143 remain at the same spot on the protein surface. Water 121 is most tightly bound, and water 143 is one of the two protein surface waters which has been found in all X-ray structures and the NMR Structure? Waters 105 and 140 move along the protein surface. Water 105 visits different parts of the a-helixof the BPTI protein, whereas water 140 jumps between two different surface sites. Analyzing each water molecule with respect to its protein and water environment using a finger and cluster algorithm, respectively, water molecules in protein pores can be identified.25 Even a protein as simple as BPTI possesses several water channels with
Figure 2. Projections of selected water trajectories shown on the background of the BPTI structure. The trajectories are drift corrected by least-square fitting of the protein structures at different times using the Kabsch algorithm.26 The X-ray protein structure is represented schematically with M~lscript.~'Water molecules marked with "E" are from the overlays. The other water molecules belong to the 60 X-ray waters.IoJ1 In the top projection, the trajectories of the four internal waters 1 1 1 , 112, 113, and 122 are displayed. The bottom projection shows selected surface and bulk water molecules. More explanations are given in the text.
up to 15 water molecules. The number of water channels fluctuates with a characteristic time of 5 ps typical for the lifetime of hydrogen bridges. At the end of the dynamics of the BPTIwater system, the cluster of the three internal waters (1 1 1, 1 12, and 113) forms a water channel with the help of two additional bulk water clusters.25 4. The Diffusion Model
For the analytically solvable diffusion model, the globular protein is considered to be inpenetrable and to have the shape of a sphere with radius a. The total protein-water system is contained in a sphere of radius c. The volume of surface water is represented by a spherical shell with outer radius b, whose value fulfills a < b < c (see Figure 3). The water molecules
The Journal of Physical Chemistry, Vol. 97, No. 43, 1993 11341
Heterogeneous Diffusion of Water at Protein Surfaces
k
1[
Q)
u
B
55 Q
ICI
0
c
0 .rl
[I)
5
ICI W
.d
a h
u c
I
0
30
60
90
120
150
180
time(ps)
geometry of diffusion model Figure 3. Geometry of the diffusion model. The inner sphere of radius represents the protein molecule. The solvent molecules are situated between the inner sphere and the outer sphere, whose radius is c. The walls of the two spheres are inpenetrable. The volume of the surface waters is represented by the layer between the inner sphere and the sphere of radius b. D
which are initially in the surface volume are represented by the radial density
which is normalized to unity
socdr3po(r) = 1 The fraction of solvent particles after time t remaining in the surface volume between the spheres of radii a and b is given by
n(t) = sabdrh ( r , t )
dr,O) = po(r)
(3)
By using the general solution of p(r,t),eq A6 from the Appendix, one finally obtains
for the fraction of solvent particles remaining in the surface volume. The coefficients sn2 in eq 4 are given by
=
3 -
2
1
+ C%,z
b 3 - a 3 (c-a)a!,6 a2+ c 2 + ac+ a2c2a!, 2 [ ( l + a,,%b) sin a!,,(b- a ) - a,(b - a) cos a,(b - a)I2 ( 5 )
The values a, are obtained by solving the transcendental eq A4 from the Appendix with a Newton-Raphson algorithm. 5. Application The usefulness of the spherical model for diffusion can be demonstrated by evaluating the diffusion constant for pure bulk water. The application is based on a 200-ps computer simulation of the dynamics of TIP3PL7water at 294 K in a sphere of radius c = 16.0 A (180 ps are used for analysis). The walls of the sphere are modeled by a mean force field potential. Stochastic dynamics with a friction constant of 50 ps-1 are applied to water molecules within 2.5 A of the sphere boundary.20 The absence of a protein molecule can be accounted for by a vanishing inner sphere radius a = 0.0 A. Diffusion across the surface of a sphere with radius b is considered. The value of b = 8.7 A is taken such that the
Figure 4. Fraction of water molecules moving out of a sphere of radius b = 8.7 A in the center part of a sphere of radius c = 16.0 A surrounding the whole system. There is no protein molecule in the system, i.e. the radius D = 0.0 A. The dashed line represents data from computer simulations. The solid line refers to the analytical model, yielding a diffusion constant of D = 0.23 A2/p for bulk water as compared to the experimentalvalue, whichis D , = 0.21 Az/ps”at thesame temperature.
asymptotic value of the fraction of water molecules remaining in the inner sphere agrees with the asymptotic value of the considered BPTI-water system. The fraction of water molecules remaining in the inner sphere after an elapsed time t is plotted in Figure 4. The quality of the statistics of the simulation data is considerably improved at short delay times by using multiple time origins. The simulation data exhibit a 5% component of instantaneous time decay of water moleculesremainingintheinner sphere. Thisisdueto fastrattling motions of water molecules in solvent cages close to the boundary of the inner sphere. No correct is made for this effect. Since the geometric parameters, the radii a, 6, and c, are determined by the system size, the diffusion constant D is the only remaining fit parameter of the analytical model. The diffusion constant can be varied by changing the time scale. The fit with the simulation data provides a diffusion constant of D = 0.23 f 0.02 A2/ps. This value is about 15% larger than the diffusion constant obtained by applying the Einstein relation D = ( r 2 ) / 6 A t to a water box and agrees quite well with the experimental value Dexp= 0.21 Az/psZ2 at the same temperature T = 294 K. This close agreement demonstrates that solvent polarizations induced by spherical boundary conditions do not spoil the value of the translational diffusion c o n ~ t a n t . ~ 3 ~ ~ ~ The dependence of the diffusion constant on the mean force field potential at the sphere boundary is tested by considering the diffusion of water from the spherical boundary of 1.%Athickness into the inner sphere of radius 14.2A. Using the same simulation data, the diffusion constant derived with the analytical model adopts the same value. From this agreement, one can conclude that a water molecule escaping from the Lennart Jones well of the mean force field at the boundary of the sphere is immediately replaced by a bulk water molecule such that no activation barrier occurs. The sensitivity of the evaluation of the diffusion constant on deviationsfrom spherical geometryis also tested. For this purpose, simulation data of 100-ps dynamics of water at T = 308 K in a box of length 3 1.10A and periodic boundary conditions are taken. These data have also been used to determine the diffusion constant with the Einstein relation yielding 0.28 f 0.02 A2/ps. The simulation data of the time-dependent fraction of water molecules remaining in the inner sphere of radius b = 10.5 A are fitted with the analytical model. Even though the system does not have a spherical geometry, the diffusion constant adopts the same value. On the other hand, small deviations in the volume ratio between theinitially tagged water moleculesand the total amount of water in the system result in significant discrepancies.
11342 The Journal of Physical Chemistry, Vol. 97, No. 43, 1993
To test deviations due to the nonspherical shape of BPTI, diffusion of water from an ellipsoid centered in the water box is studied. The half-axes of the water ellipsoid considered, x = 13.0A and y = z = 7.4A, have the same ratios as the half-axes, x = 18.5 A and y = z = 10.5 A, obtained by fitting the BPTI backbone to an ellipsoid. Also in this case a fit of the simulation data with the analytical model provides a diffusion constant of 0.28 f 0.02 A2/ps. From these results, we deduce that the spherical diffusion model can be applied also to protein-water systems as long as the shape of the protein does not deviate too strongly from a sphere and thevolume ratios of protein and surface to bulk water are considered properly. To apply the analytical model of diffusion to the BPTI-water system, spherical geometry is assumed. The true shape of BPTI is closer to an ellipsoid with half-axes x = 18.5 A and y = z = 10.5 A. The protein radius a is determined from the protein volume V, by assuming that BPTI has the shape of a sphere. The protein volume is calculated by a grid algorithm yielding the radius a = 14.0A. The surface water is defined by a distance criterium. A water molecule belongs to the protein surface if the distance of the water oxygen atom to the nearest non-hydrogen atom of the protein is less than do = 3.4A. The radius b of the sphere containing the surface water can be estimated from the volume of surface water V, as the volume difference of spheres with radii b and o according to
47r v,= -(b3 3
- a3)
The radius of the surface water sphere is calculated as b = 15.93A. The net difference between the two radii band a is 1.93 A. This is 1.47A less than the distance doused for the definition of the surface water. The major portion of this difference is due to the excluded volume of protein atoms at the surface. The corresponding average van der Waals radius is about 1.8 A. The discrepancy of 1.47 A indicates that the density of water in the surface layer is higher than that in the bulk. The radius c of the total protein-water system is determined from the volume of all waters in the protein-water system V,, which can be represented as the volume of the sphere containing the whole system minus the volume of the sphere representing the protein
4a
V, = -(c3 - a3) (7) 3 The radius of the sphere containing the total system amounts to c = 23.3 A, a value slightly larger than that for the starting structure. With this procedure, the geometry parameters of the analytical diffusion model are fixed. Figure 5 shows the time decay of the fraction of water molecules which are in the surfacevolume of the protein-water system. The simulation data are represented by the dashed line. Also here an average by multiple time origins is performed. Here the fraction of water molecules exhibits a 20% component of instantaneous decay due to the rattling motions of water molecules in their solvent cages. For the BPTI-water system, the cage effect is not negligible because the ratio of surface to volume is much larger. By considering a water molecule to change from surface to bulk and vice versa only if it remains in the new part at least for two subsequent time points, this component is reduced to 8%. The nonmonotoneous decay of the simulation data is due to water molecules which are initially in the surface volume and enter it again at later times. This reentry effect is also responsible for a nonvanishing water concentration at longer times. With a single diffusion process according to the analytical model, the simulation data cannot be fitted. Even with combination of the diffusion process with a single exponential decay term, no satisfactory approximation can be obtained. However, the simulation data can be fitted by using a superposition of two diffusionprocesses with diffusion constants D1 = 3.5 f 0.1 AZ/ps and 0 2 = 0.06 0.005 AZ/ps and relative weights of a1 = 0.25
*
Knapp and Muegge 1
0.8 0.6
0.4 0.2 0
100
200
300
400
time ( p s )
Figure 5. Fraction of water molecules remaining in the surface volume of the protein depicted as a function of time. The dashed line represents the computer simulation data, and the solid line refer to the analytical model using the geometry parameters a = 14.0 A,b = 15.93 A, and c = 23.3 A as explained in the text. The fit with the analytical model corresponds to 25% of the water molecules with a diffusion constant of DI = 3.5 A 2 / p and 75% of the water molecules with a diffusion constant of D2 0.06 A1f p.
and a2 = 0.75,respectively (solid line). This can be interpreted as follows: 25% of the water molecules encounter a hydrophobic repulsive protein surface potential from which they bounce back into the bulkwater immediately. Theother 75% areat hydrophilic surface areas where the diffusion is slowed down by a factor of 4 as compared to bulk water conditions. In future studies it has to be investigated whether details of the hydrogen bonding of water molecules at the protein surface are influenced by the polarization of the water molecules due to the spherical boundary conditions.24 6. Conclusions
Analyzing the computer simulation data of a protein-water system with a spherical model for diffusion,the diffusionconstant for motion of water molecules at the protein surface can be derived. So far, only a few hand-selected trajectories of surface waters could be used to derive a diffusion constant at the protein surface. Hence the statistical uncertainty for estimating the diffusion constant was large. The present method considers all water molecules which are at different timeorigins initiallyin the surface volume, yielding estimates of diffusionconstants which are highly reliable. In agreement with former results of computer simulation of BPTI-water sy~tems,~~-IS it has been found that the diffusion process of the majority of surface waters is slowed down by a factor of 4. This is also in agreement with recent NMR experiments on BPTI, where residence times of surface water molecules in the time regime of 100-300 ps have been measured.s.9 Assuming a distance of 3 A for the sensitivity of NMR measurements, the corresponding diffusion constant varies between 0. l and 0.03 R2/ps. A new aspect of the present work is the fraction of 25% surface waters which encounter hydrophobic areas of the protein surface such that they move away from the surface much faster than they move in the bulk. Without using the analytical model for diffusion it is not possible to obtain this result. It is likely that the diffusion along the rough protein surface is much slower than the diffusion away for the surface. Then it would be difficult to measure the surface diffusion seperately since water molecules would reach other places at the protein surface by diffusion in the bulk rather than at the surface. Acknowledgment. The authors are thankful for discussions with Daniel Hoffmann and for the programming skills of Katrin Blumhagen. The CHARMm source code has been provided by
The Journal of Physical Chemistry, Vol. 97, No. 43, 1993 11343
Heterogeneous Diffusion of Water at Protein Surfaces Prof. Martin Karplus and Molecular Simulations Inc. This work is supported by the Deutsche Forschungsgemeinschaft SFB 3 12, Projekt D7.
Appendix In the Appendix, the general time-dependent solution for diffusionwith spherical boundary conditions according to Figure 3 is given by using ref 28. The radial part of the diffusionequation can be written as
where u(r,t) = r*p(r,t), D is the diffusion constant, and p(r,t) is the particle density of the solvent. With reflecting boundary conditions
($)
r=a,e
=o
The time-dependent radial part of the diffusion equation (Al) can be solved with the following expression2*
u(r,t) = exp(a2Dt)[sin a(r - a ) + aa cos a(r - a ) ] (A3) where a fulfills the transcendental equation z
tan z = 1 +pz2 with and z = a ( c - a ) ( c - a)2 The complete solution is given by
-‘P
-s
3
c3-03
(A5)
c dr’rl2p0(r’) (A6) a
where po(r) is the initial particle distribution. The last term in expression A6 accounts for the stationary value at infinite time. According to the expression A3, the normalized basis functions R,(r) are given by Nn
R,(r) = -[sin r
a,(r - a )
+ anacos a,(r - a ) ]
(A7)
with normalization constant
References and Notes (1) Saenger, W. Ann. Rev. Biophys. Chem. 1987,16, 93. (2) Frauenfelder, H.; Parak, F.; Young, R. D. Ann. Rev. Biophys. Chem. 1988, 17,451. (3) Dill, K. A.Biochemistry 1990,29,7133. (4) Wagner, G.; Wfithrich, K. J . Magn. Reson. 1975,20, 435. (5) Wagner, G. Q. Rev. Biophys. 1983, 16, 1. (6) Norvell, J. C.; Nunes, A. C.; Schoenborn, 9. P. Science 1975,190, 568. (7) Bee, M. Applications of Quasielastic Neutron Scattering to SolidState Chemistry, Biology, and Material Science; Adam Hilger: London, 1988. (8) Otting, G.; Wiithrich, K. J . Am. Chem. Soc. 1989, Ill, 1871. (9) Otting, G.; Liepinsh, E.; WUthrich, K. Science 1991, 251, 974;J. Am. Chem. Soc. 1991,113,4363. (10) Deisenhofer, J.; Steigemann, W. Acta Crystallogr. 1975,831,238. (11) Wlodawer, A,; Walter, Y.; Huber, R.; Sj6lin. L. J . Mol. Biol. 1984, 180, 301. (12) Van Gunsteren, W. F.; Berendsen, H. J. C.; Herman, J.; Hol, W. G. J.; Postma, J. P. M. Proc. Narl. Acad. Sci. U.S.A. 1983,80,4315. (13) Levitt, M.; Sharon, R. Proc. Natl. Acad. Sci. U.S.A. 1988,85,7557. (14) Levitt, M. Chem. Scr. 1989,29A, 197. (1 5) Brunne, R.M.; Liepinsh, E.;Otting, G.; WUthrich, K.;van Gunsteren, W. F. J . Mol. Biol. 1993,211, 1040. (16) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swammathan, S.;Karplus, M. J . Comput. Chem. 1983,1, 187. (17) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D. J. Chem. Phys. 1983, 79,926. (18) Knapp, E. W.; Nilsson, L. Reaction Centers of Photosynthetic Bacteria; Michel-Beyerle, M.E., Ed.; Springer: Berlin, 1990; pp 437. (19) Brooks, C. L.,111; Karplus, M. J . Chem. Phys. 1983,79, 6312. (20) Briinger, A.; Brooks, C. L., 111; Karplus, M. Chem. Phys. Lett. 1984, 105,495. (21) Ryckaert, J.-P.: Ciccotti, G.;Berendsen, H. J. C. J. Comput. Phys. 1977,23, 327. (22) Mills, R. J . Phys. Chem. 1973, 77, 685. (23) Belch, A. C.; Berkowitz, M. Chem. Phys. Lett. 1985,113, 278. (24) King, G.; Warshel, A. J. Chem. Phys. 1989,91,3647. (25) Muegge, I.; Knapp, E. W. EBSA Workshop Palermo. Conference Proceedings of the Italian Physical Society; 1993, in prcss. (26) Kabsch, W. Acta Crystallogr. 1976,A32, 922. (27) Kraulis, P. J. J. Appl. Crystallogr. 1991, 24, 946. (28) Carslaw, H. S.; Jaeger, J. C. Conduction of Heat in Solids; Clarendon: Oxford, U.K., 1980.
