Heterogeneous Hydration of p53/MDM2 Complex.

Water-mediated interactions play critical roles in biomolecular recognition processes. Explicit solvent molecular dynamics (MD) simulations and the variational implicit-solvent model (VISM) are used to study those hydration properties during binding for the biologically important p53/MDM2 complex. Unlike simple model solutes, in such a realistic and heterogeneous solute-solvent system with both geometrical and chemical complexity, the local water distribution sensitively depends on nearby amino acid properties and the geometric shape of the protein. We show that the VISM can accurately describe the locations of high and low density solvation shells identified by the MD simulations and can explain them by a local coupling balance of solvent-solute interaction potentials and curvature. In particular, capillary transitions between local dry and wet hydration states in the binding pocket are captured for interdomain distance between 4 to 6 Å, right at the onset of binding. The underlying physical connection between geometry and polarity is illustrated and quantified. Our study offers a microscopic and physical insight into the heterogeneous hydration behavior of the biologically highly relevant p53/MDM2 system and demonstrates the fundamental importance of hydrophobic effects for biological binding processes. We hope our study can help to establish new design rules for drugs and medical substances.

Introduction

Water-mediated interactions play critical roles in biomolecular recognition processes, such as protein folding and protein–ligand binding.[1−6] In these processes, hydrophobic regions and apolar molecules are often driven together by the surrounding aqueous solution which has been noted since the famous publication by Kauzmann in 1959.[7,8] Usually, increasing the apolar binding surface area decreases its dissociation constant (increase the strength of binding).[9] During the process of optimizing lead molecules to drug, one of the important strategies for medicinal chemists is to utilize these hydration properties of ligands and receptors.[10] The importance of hydrophobic effects in quantitative structure–activity relationship (QSAR) has been documented since the first publication of QSAR in 1962.[11] Detailed understanding of the hydration properties for each important biological system is important for the rational design of better drugs.[12−15]

Theoretical approaches for modeling the hydrophobic effect, such as void volume theories, interpret that the hydrophobic effect is length-scale dependent.[16−23] At small length scales (e.g., solute radius smaller than roughly 0.5 nm), the solvation of solutes is mainly an entropically driven process in which the water hydrogen bonding network rearranges itself to adapt for the intercalation of solutes. On the other hand, at large length scales, a significant number of water hydrogen bonds around the solute are broken due to the restructuring at extended, relatively flat surfaces of the solutes, resulting into a domination of enthalpic contributions to the solvation process.[9]

Unlike generic theoretical models with uniform chemical detail, for proteins with both polar and apolar residues, the solute–solvent behavior can be quite complex.[24] In such a heterogeneous and dynamical environment, water and protein structures are sensitive to the local solute–solvent interface geometry, as well as van der Waals (vdW) and electrostatic interactions. The coupling between interface geometry and structural polarity is crucial to protein biological function. Characterizing the hydration properties surrounding the biomolecules and comparing them with those of bulk water have been a lasting scientific endeavor in the past decades. Significant progress has been made through the combination of various experiments, theoretical techniques, and models.[25]

Computer simulations, such as molecular dynamics (MD) simulation with explicit solvent models, can shed light on the delicate balance among these factors and illustrate the hydrophobic effects associated with ligand molecules[26,27] and induced polarization.[28] However, one drawback of molecular dynamics simulations is the relatively high computational demand for large systems such as biomolecules. Furthermore, the intrinsic fluctuations in all atom MD simulations can obscure the underlying intuitive physical pictures. In addition to the explicit solvent simulations, various implicit-solvent models are also developed as efficient alternatives to study the biomolecular hydration behaviors and the values of hydration free energy.

To date, most commonly used implicit-solvent models (e.g., Poisson–Boltzmann/Surface Area (PBSA)) are based on predefined solute–solvent interfaces (e.g., solvent accessible surface (SAS), solvent excluded surface (SES), or van der Waals (vdW) surface), in which polar and apolar interactions are assumed to be additive and decoupled.[29−34] Recently, a Gaussian-based approach was developed to obtain the distribution of the dielectric constant for a protein in solution.[35] By using the local dielectric constant values, the implicit-solvent models lead to better results than the ones using a uniform dielectric constant. There are also other microscopic theoretical treatments of water without predefined surface and based on the statistical mechanics theory of molecular liquid. The most notable one is the three-dimensional reference interaction site model (3D-RISM).[36,37] In this model, the solvent distribution function and correlation functions are self-consistently solved with a proper closure relation. It avoids the explicit molecular sampling of instantaneous solvent configurations because it works on average distribution function for the particular atoms. With recent advancement, 3D-RISM is able to calculate hydration free energy, equilibrium solvent distributions, and other thermodynamic quantities efficiently for biomolecular systems.[38]

The variational implicit-solvent model (VISM) developed by Dzubiella et al.[33,34,39,40] provides a self-consistent description of the molecular solvation with the contributions from the solute–solvent interface geometry and van der Waals and electrostatic interactions coupled together in a physical based free-energy functional within the purely surface-based implicit-solvent model framework. An equilibrium solute–solvent interface is determined by minimizing the solvation free-energy functional that balances the interface local geometry with various hydration free energy contributions.[33,39,41] A similar approach has also been put forward by Wei et al.[42,43] In their model, the expression of free energy and the optimization algorithm are different from VISM.

Both theoretical and experimental studies over the past decades showed a complicated picture of hydration behaviors around biomolecules. In this study, we apply a joint MD and VISM analysis to gain certain quantitative and qualitative insights into the heterogeneous hydration, the solute–solvent interface, and individual water molecule behavior around proteins. Most previous studies of hydration behavior focused on either solute geometrical complexity with uniform chemical details[41,44,45] or heterogeneous chemical details with a simple domain interface.[27] Hua et al. modeled parallel plates consisting of both hydrophobic and hydrophilic particles.[46] They showed that the behavior of water between the two plates strongly depends on the distribution of the hydrophobic and hydrophilic particles on the plates. In the studies presented here, we choose a biologically important and realistic system, the p53/MDM2 protein complex with both geometrical and chemical detail complexities, as an example to investigate hydration properties near complex biomolecular binding sites.

The receptor protein MDM2 acts as a negative regulator of the tumor suppressor protein p53, and loss of the p53 function is observed in ∼50% of human malignancy.[47] To understand the hydration properties of the p53/MDM2 complex during the process of association is very important for the understanding of the p53/MDM2 binding mechanism. Furthermore, inhibitor design to reactivate p53 as tumor suppressor protein will also benefit from this study.

This system exhibits a strong hydrophobic character at the p53/MDM2 binding interface (70% of the residues at the binding interface are apolar).[48−50] However, the edge of the binding pocket is decorated by polar hydrophilic residues. Therefore, it provides an excellent example to study the heterogeneous hydration properties in a protein–protein binding process at both atomic and continuum levels. In Figure 1, we show the molecular surface of this system. The red color represents charged hydrophilic residues, the pink color represents neutral hydrophilic residues, and the blue color represents hydrophobic residues. In Figure 1A, it is apparent that the hydrophobic areas of the two domains are in contact with each other in the bound state. In Figure 1B, we separate the two domains and expose the interfacial area. One can clearly see that the binding interfaces are essentially hydrophobic (blue).

Figure 1

p53/MDM2 complex (PDB code: 1ycr.pdb). (A) Molecular surface of the complex in the bound state. The domain interface is buried with hydrophobic residues in blue, charged residues in red, and neutral hydrophilic residues in pink. (B) Two domains of the complex are separated and rotated in such a way that the binding interfaces face the reader.

In this work, MD simulations and VISM calculations are used to study the heterogeneous hydration behaviors around the protein during the p53/MDM2 binding. Through MD simulations, we construct water density profiles from which we observe high water density solvation shells near the hydrophilic residues and depleted water density regions near the hydrophobic binding pocket of MDM2. Furthermore, we identify a bimodal state for the water occupancy behavior in the binding pocket when the interdomain distance is between 4 and 6 Å. A solute–solvent hydration energy density map is used to illustrate the relation between the water density and the various solute–solvent interactions. As of now, six clinical trials are underway to evaluate the clinical benefits for inhibiting p53/MDM2 binding.[51] We hope our study can further the understanding of the binding affinity and kinetics between p53 and MDM2, which eventually could help to design next generation p53/MDM2 inhibitors.

Theory and Methods

System Preparation

The initial structure of the p53/MDM2 complex is taken from the Protein Data Bank (PDB code: 1ycr.pdb). Hydrogen atoms are added at pH = 7.0 using the Protein Preparation workflow in Maestro.[52] Acetyl (ACE) and N-methyl amide (NMA) are incorporated to cap the N- and C-termini, respectively. A series of configurations for the protein complex are produced by separating the two domains along the axis through their geometrical centers. The interdomain separation d is chosen as the reaction coordinate. The increment for interdomain distance is 2 Å in this study. The value of d = 0 Å corresponds to the native crystal structure of the complex.

For MD simulations, each configuration is then placed in an orthorhombic box with a minimum distance of 10.0 Å to the boundary of box and hydrated with a pre-equilibrated box of TIP3P water using the System Builder module of the Desmond package.[53−55] All overlapping solvent molecules are removed, and counterions are added to maintain charge neutrality.

For VISM calculations, the same input solute structures with partial charges as in the MD simulations are used, so that comparisons between explicit and implicit modeling can be made.

Molecular Dynamics Simulations

All molecular dynamics (MD) simulations are performed using the Desmond package.[53−55] The OPLS 2005 force field[56,57] is used to model the protein interactions, and the TIP3P model[58] is used for the water. Particle-mesh Ewald method[59] (PME) is used to calculate the long-range electrostatic interactions with grid spacing of 0.8 Å. van der Waals and short-range electrostatic interactions are smoothly truncated at 9.0 Å. The Nose-Hoover thermostat[60] is used to maintain the constant simulation temperature and the Martina-Tobias-Klein method[61] is used to control the pressure. The equations of motion are integrated using the multistep RESPA integrator[62] with an inner time step of 2.0 fs for bonded interactions and nonbonded interactions within the short-range cutoff. An outer time step of 6.0 fs is used for nonbonded interactions beyond the cutoff. Periodic boundary conditions are applied.

To focus on the water behavior around the p53/MDM2 surface, we constrain the protein molecules so that their conformational fluctuations are not convoluted with the process. These constraints will also facilitate the grid-based water density mapping around the protein and inside of the binding pocket and the direct comparison with the VISM calculations. The systems are equilibrated with the default protocol provided in Desmond. After the equilibration, a 40 ns NPT production simulation is performed for each configuration at temperature 300 K and pressure 1.01 bar with a 20 kcal/(mol·Å2) harmonic potential restraint on solute atoms, and the simulation trajectories are saved in 4-ps intervals for analysis.

Level-Set Variational Implicit-Solvent Method (VISM)

The details of VISM with the Coulomb-field approximation (CFA) are described extensively in the previous publications.[33,34,39,40] In short, we optimize the free energy of the solvation system as a functional of all possible solute–solvent interfaces Γ.where N solute atoms are located at x1, ..., x inside Ω and with point charges Q1, ..., Q, respectively. (In our system of p53/MDM2, N = 1688.) The first term Pvol(Ω) is the volumetric part of the energy for creating the solute cavity Ω with P being the pressure difference between the solvent liquid and solute vapor. The second term is the surface energy, where γ(x) is the surface tension given by γ(x) = γ0(1 – 2τH(x)), where γ0 is the constant macroscopic surface tension for a planar solvent liquid–vapor interface which is set as 0.127 kBT/Å2 at 300 K according to the TIP3P water simulation.[63] τ is the first order correction coefficient often termed as the Tolman coefficient[64] which is set as 1.0 Å,[41,65] and H(x) is the mean curvature defined as the average of the two principal curvatures. The third term is the energy of the van der Waals interaction between the solute atoms and the continuum solvent. The parameter ρ is the solvent density which is set as bulk water density ρ0 = 0.0333/Å3. The fourth term is the electrostatic contribution to the solvation free energy. It is defined by the Born cycle[66] as the difference of the energies of two states, where ε0 is the vacuum permittivity, ε is the relative permittivity of the solute molecule, and ε is the relative permittivity of the solvent which are defined by the VISM solute–solvent interface Γ by ε(x) = ε if x ∈ Ω and ε(x) = ε if x ∈ Ω.

Now the free energy G[Γ] determines the effective boundary force, −δΓG[Γ], acting on the VISM solute–solvent surface Γ, where δΓ is the variational derivative with respect to the location change of Γ. It is only the normal component of this force that can affect the motion of such a solute–solvent surface. We denote by n = n(x) the unit normal vector at a point x on the solute–solvent surface Γ, pointing from the solute region Ω to the solvent region Ω. Then the normal component of the effective boundary force is given by[67−69]where K = K(x) is the Gaussian curvature, defined as the product of the two principal curvatures, at a point x on Γ. This force will be used as the “normal velocity” in our level-set numerical calculations.

To minimize the free-energy functional (1), we choose two different types of initial surfaces, a loose and a tight initial surface.[39,40] Both of them enclose all the solute atoms located at X1, ..., X. The tight initial surface is defined by the van der Waals (vdW) surface. The loose initial surface is often set to be a large sphere enclosing all the solute atoms. In this study, it is chosen so that the closest solute atom (from the vdW sphere edge) is at least 1.5 water diameters away form the surface. The initial interface can have a very large value of the free energy. It is subsequently moved in the direction of steepest descent of the free energy by the level-set method until a minimum is reached. The starting point of the level-set method is the representation of a surface Γ using the (zero) level-set of a function ϕ = ϕ(x): Γ = {x: ϕ(x) = 0}.[70−72] The motion of a moving surface Γ = Γ(t) with t denoting the time is then tracked by the evolution of the level-set function ϕ = ϕ(x, t) whose zero level-set is Γ(t) at each t. Such evolution is determined by the level-set equation

Methods for Data Analysis

To study the water density profiles around the p53/MDM2 complex, a lattice with grid size of 0.8 Å is constructed and TIP3P water oxygen atoms from MD simulations are assigned to the nearest grid points.[73,74]

To understand the origin of the inhomogeneous water distribution from the solute–solvent interaction perspective, a solute–solvent hydration energy density map is constructed. From this energy density map, we directly identify the hydrophobic and hydrophilic regions around the protein. In this study, we use the VISM free energy functional to fulfill this purpose. In principle, the solute–solvent hydration energy density map can also be obtained from MD results. However, it is much more desirable to calculate the solute–solvent hydration energy density from efficient methods than long time MD simulations.

From the VISM solvation free energy functional (eq 1), the van der Waals and electrostatic solvation free energy contributions can be considered as the integration of solute–solvent interactions over the entire solvent region. Here, we define a solute–solvent hydration energy density aswhere

The vdW solute–solvent hydration energy density φvdw(x) and electrostatic solute–solvent hydration energy densities φelec(x) are obtained by using the same formula as the individual component of the solvation free energy functional (i.e., eq 1) in VISM without the volume integration.

In order to study the aqueous behavior only in the hydrophobic binding pocket of MDM2, we used the differences between the contracted level-set VISM surface with a loose initial surface and the molecular surface to define the binding pocket (Figure 2).

Figure 2

Binding pocket region is defined by the differences between molecular surface (green) and a contracted VISM surface with loose initial surface (red).

Results and Discussion

Water Density, van der Waals (vdW), and Electrostatic Contributions in the Hydration Process of p53/MDM2

The water density profile is constructed from a 40 ns MD simulation with TIP3P water. Figure 3 shows the cross-section of water density profile across the target protein MDM2 binding pocket when interdomain distance d = 12 Å. In this figure, the protein complex is represented by its molecular surface with the same color code as Figure 1. The value of the local water density is represented by colors in the legend. White color indicates zero density, blue represents half of the bulk water density, yellow represents the bulk water density, and red represents twice of the bulk water density.

Figure 3

(A) Cross-section of water density profile through the binding pocket based on the MD simulation at interdomain distance d = 12 Å. The molecular surface is represented with the same color code as Figure 1 (Blue for hydrophobic residues, red for charged residues, pink for neutral hydrophilic residues). The colors in the legend represent the average water density (range from 0 to 2 and in units of bulk density ρ0). (B) The same cross-section water density profile as (A) without showing the protein.

It is observed that high density hydration layers are formed around most parts of the solute molecules as shown by the red color. The water density of the hydration layer close to hydrophilic residues (red molecular surface) is much higher and more discrete than that near hydrophobic residues (blue molecular surface). When near the convex protein surface, even around the hydrophobic residues, high density and continued hydration layers are observed. When the protein surface is concave and residues are hydrophobic near the binding pocket at the interdomain region, the hydration layers form at the entrance of the pocket away from the molecular surface and water density drops to zero at the bottom of the hydrophobic pocket. These hydration shells behave differently according to local protein surface geometry and chemical details.

Such a phenomenon as high density layers formed far away from the molecular surface is also observed in several simpler model systems previously[41,45] and is attributed to the liquid–vapor interface in the hydrophobic binding pocket, and the water density fluctuating between higher than bulk to zero. These behaviors are explained by capillary evaporation and capillary condensation driven by solute–solvent interactions.[75]

To understand quantitatively the relations between the hydration water density distribution around the protein and the solute–solvent interactions, we analyze and compare the water density profile ρ(x) and the solute–solvent hydration energy density φvdw+elec(x). The solute–solvent interactions contributed hydration energy density consists of vdW and electrostatic components between solute and solvent. Details can be found in the Theory and Methods part.

In Figure 4A, it shows the cross-section of water density profile ρ(x) across the binding pocket with the location of molecular surface (black line inside). Figure 3B shows the solute–solvent hydration energy density map φvdw+elec(x). In this map, the red color represents φvdw+elec(x) larger than +0.1kBT/Å3. The blue color represents φvdw+elec(x) smaller than −0.1 kBT/Å3. By comparing Figure 3A and Figure 4B, we find that the high density water layers form at the region with the solute–solvent hydration energy density less than −0.1 kBT/Å3 shown as deep blue in Figure 4B. The individual components of this solute–solvent hydration energy density (vdW and electrostatic contributions) are shown in Figure 4C,D, respectively.

Figure 4

(A) Cross-section of water density profile ρ(x) (range from 0 to 2 and in units of bulk density ρ0) through the MDM2 binding pocket, black line inside is the location of molecular surface (MS) from MD simulations, (B) solute–solvent hydration energy density map (φvdw+elec(x)) and (C) vdW (φvdw(x)), (D) electrostatic (φelec(x)) (range from −0.1 to 0.1 and in units of kBT/Å3) individual parts derived from VISM free energy functional. Figures are for the p53/MDM2 at interdomain distance of 12 Å.

Most hydrophilic areas shown in red in Figure 4A originate from the favorable electrostatic solute–solvent interactions. This is demonstrated in Figure 4D where the blue color depicts the strongly attractive solute–solvent electrostatic interaction only. For this reason, more water molecules accumulate in those regions until strong solute–solvent vdW repulsion (shown red in Figure 4C) counter-balances the condensation. Thus water densities in these solvation shells are much higher near the polar hydrophilic residues. In contrast, the white color at the binding pocket region indicates that the electrostatic contribution is nearly zero, i.e., an apolar region. The vdW force is the dominant interaction between solute and solvent at this region. The lack of strong competition between attractive and repulsive solute–solvent interactions leads to a lower water density solvation shell minimizing the surface energy and the high density hydration layers form far from the molecular surface at this concave region. The comparison between water density distribution from extensive MD simulations and the hydration energy density from VISM free energy functional shows that we can qualitatively identify the hydrophobic and hydrophilic regions from the solute–solvent hydration energy density map without extra-simulations and it can be achieved in minutes just based on the information of 3D protein structure.[33,34,39,40]

One of the VISM’s advantages is that it can track the solute–solvent interfacial geometry through the coupling of the various solute–solvent interactions in a free energy functional and is solved self-consistently.[33,34,39,40] In the following study, we investigate the interplay between solute–solvent interactions and the solute–solvent interface geometry by numeric minimization of the VISM free energy functional (i.e., eq 1) with respect to the solute–solvent interface.

Figure 5 shows the superposition of equilibrium VISM surfaces from both loose and tight initial surfaces[33,34,39,40] with the water density profile and the solute–solvent hydration energy density map at the interdomain distance of 12 Å. The two equilibrium VISM surfaces correspond to two stable states of solvation and are the results of the complex relations between polar, apolar, and surface energy contributions to the hydration process. When superimposing the equilibrium VISM surfaces on the MD water density profile ρ(x) in Figure 5A, we find that most part of the VISM surface corresponds to the high density hydration shell regardless of the initial conditions except for the interdomain region. In Figure 5B, we also superimpose the equilibrium VISM surfaces on the solute–solvent hydration energy density map φvdw+elec(x) to show the relation between water density and solute–solvent hydration energy density.

Figure 5

Water density profile (range from 0 to 2 and in units of bulk density ρ0) (A) from MD simulations and solute–solvent hydration energy density map (range from −0.1 to 0.1 and in units of kBT/Å3) (B) from VISM free energy functional are superimposing with the equilibrium VISM surfaces (depicted by the thick orange–white–blue lines) which are obtained from loose and tight initial surfaces, respectively. Figures are for the p53/MDM2 at interdomain distance 12 Å.

Furthermore, in Figure 6A, it shows the average water density (in units of average bulk water density ρ0) vs the signed distance to the p53/MDM2 equilibrium VISM surface (a tight initial surface was used here). A negative value represents the distance to the equilibrium VISM surface from inside and vice versa. The water density peak at zero suggests that the equilibrium VISM surface is largely located at the first hydration shell and the water density vanishes to nearly zero inside the VISM surface (<−1.4 Å). However, in the MDM2 binding pocket, the VISM surface captures the high energy density shell much further away from molecular surface (Figure 5B). Comparing MD results with VISM surfaces, we can identify the delicate hydration balance for hydrophobic and hydrophilic effects near protein surfaces with both geometrical and chemical complexities. We find VISM is able to properly describe hydrophobic binding pockets which are of ultimate interest in most cases. In contrast, traditional implicit-solvent models based on predetermined surfaces, such as PBSA, are only able to account for the geometrical shape defined by the vdW radii.

Figure 6

(A) Water density profile ρ(x) (in units of average bulk water density ρ0) in 1D obtained from the MD simulation vs the distance to the VISM surface of a tight initial. (B) Solute–solvent hydration energy density and individual components from the VISM free energy functional vs the distance to the VISM surfaces obtained by a tight initial surface. A negative distance represents the distance from the VISM surface to the solute direction, and a positive value represents the other way.

Figure 6B shows the solute–solvent hydration energy density map φvdw+elec(x) (black line) and the individual components φvdwc(x) (green line) and φelec(x) (red line) as a function of signed distance to the equilibrium VISM surface. By comparison of these two figures, one can see that the solute–solvent interactions are essentially the main driving forces for the formation of high density hydration layer.

There are 1688 atoms in the p53/MDM2 complex. In this study, it takes about a week on average for 40 ns MD simulation with 8 CPUs in parallel for each configuration. In contrast, it takes about 1 h for tight initial surface and 3 to 4 h for loose initial surface with single CPU processor for the same configuration in MD simulations. The calculation speed of VISM strongly depends on the initial surfaces and convergence threshold.[39,40] We are working on speeding up the VISM calculation through numerical algorithm improvement, better initial conditions, and parallelization.

Hydration Properties inside of Hydrophobic MDM2 Binding Pocket

Interesting hydration behaviors happen inside of the hydrophobic MDM2 binding pocket. In the following study, we define the region of MDM2 binding pocket as the area between the molecular surface and a contracted VISM surface from a loose initial surface of MDM2 (details in methods of analysis section).

In Figure 7A, it shows the water density distribution inside of the binding pocket of MDM2. In this figure, the blue color represents density lower than half of the average bulk density ρ0 and the red color represents regions with water density 1.5 times higher than the average bulk density ρ0. When the transactivation domain of p53 binds to MDM2, the key residue L26 occupies the left region in the figure, W23 occupies the middle, and F19 occupies the right. The left pocket is formed by hydrophobic residues Ile99, Leu57, and Leu54. The middle one is formed by Val75, Leu57, Ile61, Val93, and Ile99, and the right pocket is formed by Val93, Val75, and Il61. In Figure 7A, we find three isolated high water density regions in the p53-L26 binding pocket, four high water density regions in the W23 binding pocket, and two high water density regions in the F19 binding pocket. These high density regions will be occupied by the three p53 residues when they bind. In Figure 7B, we show φvdw+elec(x) inside of the binding pocket, the blue color represents hydration energy density lower than −0.1 kBT/Å3 indicating hydrophilic regions, the red color represents energy density higher than +0.1 kBT/Å3 indicating hydrophobic regions, and the white color represents the zero hydration energy density level. When p53 binds to this pocket and water molecules are replaced by ligand atoms, a substantial amount of entropy can be gained by occupying the red high water density with ordered water regions (red in Figure 7A). Enthalpic gain is more pronounced at the blue region where water molecules are relatively close to the surface of the protein. In Figure 7B, this figure illustrates that the low φvdw+elec(x) regions locate near the critical binding regions for the three p53 key residues.

Figure 7

(A) Water density profile inside of the MDM2 binding pocket from MD simulations with p53 peptide absence. The entrance of the binding pocket is at the front and the bottom is at the back. In part A, the blue color represents density lower than 0.5ρ0 and the red color represents density higher than 1.5ρ0. (B) Solute–solvent hydration energy density map φvdw+elec(x) from VISM free energy functional inside of the binding pocket. In part B, the blue color represents the value lower than −0.1 kBT/Å3, the red color represents values higher than +0.1 kBT/Å3, and the white color represents the zero energy level.

In Figure 8A, we show the distribution of water densities inside of the binding pocket in a range of interdomain distances from 4 to 12 Å. In this figure, when the interdomain distance is larger than 6 Å, the water density inside of the cavity has a one-state distribution, and the peaks are located around 0.5–0.6ρ0. This is consistent with the overall hydrophobic nature of the binding pocket. When the interdomain distance decreases to 5 Å, an interesting two-state water density profile emerges with peaks at 0.2ρ0 and 0.5ρ0. They correspond to partially dry and wet states for this pocket. As the interdomain distance decreases further to 4 Å, the water density profile peaks at zero, corresponding to a completely dry state.

Figure 8

(A) Water density distribution inside of the binding pocket obtained from MD simulations when different interdomain distances range from 4 to 12 Å. Blue line represents the bulk water distribution in the same volume as d = 12 Å. (B) Normalized autocorrelation function (ACF) C(t) of occupancy fluctuations in the pocket with the p53/MDM2 interdomain distance increasing from 4 to 12 Å. The blue line represents the bulk water autocorrelation in the same volume as d = 12 Å.

In Figure 8B, we show the water occupancy autocorrelation function (ACF), C(t) = ⟨δN(t)δN(0)⟩/⟨δN2⟩. Bulk water ACF is shown with the blue line as a comparison.[45] In the recent work of Setny et al.,[45] a more than 1 order of magnitude shift in time scale for the ACF in a purely hydrophobic pocket was reported. Here, we find that the ACF of occupancy fluctuation in the MDM2 binding pocket depends on the interdomain distances. A very long correlation time exists for a particular interdomain distance.[76] When the separation deviates from there, the correlation time reduces quickly. In the case of p53/MDM2, the distance is around 5 Å. When the interdomain distance is larger than that value, the binding pocket begins to get flooded by bulk water giving shorter correlation times. When the interdomain distance is smaller than 5 Å, the water also fluctuates faster. Under this condition, the waters inside the pocket prefer to evacuate the binding site quickly. Interestingly but not surprisingly, when around d = 5 Å, the drying and wetting processes compete and result in long correlation times (∼400 ps). As a result, a bimodal distribution with two peaks at 0.2ρ0 and 0.5ρ0 arises at this distance as shown in Figure 7A with population ratio of 0.41/0.59 which is obtained by Gaussian distribution fitting of the bimodal density probability.

In order to interpret the population of dry and wet states in 40 ns MD simulation at certain distance, we artificially define the dry state as the pocket water density smaller than 0.2ρ0 and the wet state as pocket water density larger than 0.5ρ0. Figure 9 shows the water density profile by averaging MD frames which are selected by the dry (A) and wet (B) state criteria. In Figure 8A, a water depleted dry region is observed in the interdomain region which is flooded in Figure 8B. During the 40 ns simulation, we find 19% frames in which the pocket water density is smaller than 0.2ρ0 and 24% frames with larger than 0.5ρ0 pocket water density. In this interdomain distance, the dry and wet states show comparable probability.

Figure 9

(A) Water density profile by averaging MD trajectories in which the binding pocket water density is less than 0.2ρ0. (B) Water density profile by averaging MD trajectories in which the binding pocket water density is larger than 0.5ρ0. The black line inside represents the location of the molecular surface. Figures are for the p53/MDM2 at an interdomain distance of 5 Å.

In Figure 10, we also calculated the relative population of dry and wet states for different interdomain distances ranging from 4 to 12 Å. From these p53/MDM2MD simulations, we find that the binding pocket will prefer a flooded wet state when the interdomain distance is larger than 6 Å and a depleted dry state when the interdomain distance is smaller than 4 Å. The interesting dry–wet transition likely takes place when the p53/MDM2 interdomain distance is between 4 and 6 Å. The dehydration at the hydrophobic interface has profound implications to the kinetics of p53/MDM2 binding.[45] MDM2 has positively charged residues (K51 and H73) at the edge of the binding pocket, and p53 transactivation domain (TAD) has two corresponding negatively charged partners (E17 and E28). The electrostatic interactions help to orientate ligand and receptor during binding. However, Schon et al.,[77] showed that the association rate of p53/MDM2 was independent of the ionic strength and the truncation of E17 and E28 did not reduce the binding affinity. On the basis of these results, they hypothesized that the binding of p53/MDM2 is dominated by the dehydration of the hydrophobic interface, hydrophobic interactions, and interface rearrangement while electrostatic contribution is less pronounced.[77]

Figure 10

Probability of having a dry and wet pocket when interdomain distances ranging from 4 to 12 Å in MD simulations. Water density smaller than 0.2ρ0 for dry pocket (blue line) and larger than 0.5ρ0 for wet pocket (red line).

In the current study, the binding process of the p53 transactivation domain (TAD) and MDM2 is modeled by conformational constrained protein domains. Although the p53 TAD has been shown as predominantly disordered before binding in solution,[78] a high percentage of α-helical secondary structure is required for stable binding with MDM2. The high affinity of stapled p53 peptide also demonstrates that preorganization of TAD is beneficial to the binding.[79,80] Therefore, the current frozen conformational modeling system can provide important insights to the binding process from the hydration perspective. In addition, the constrained MD facilitates grid based water density analysis. It should be noted that the static protein simulations and analysis cannot capture every aspect of the actual binding process, such as the long relaxation time for the pocket water due to large scale protein conformational dynamics. The realistic p53-MDM2 binding involves global conformational changes.[81] The coupling of protein dynamics to the surrounding water dynamics is of immense interest, and a completed picture will be undoubtedly much more complex than that presented in this study.

Conclusions

In summary, we have studied the aqueous behavior around the p53/MDM2 complex in detail using MD simulations and VISM calculations. We show that water density is significantly lower in regions near the concave hydrophobic binding pocket and the hydrophobic ligand residues. Near the hydrophobic cavity, the first solvation layer forms at the entrance of the pocket. The water density drops from higher than bulk to zero at the pocket bottom. The water in the hydrophobic binding pocket behaves like the liquid–vapor interface and the water density fluctuates between high and low values.

In addition, we show that the local water density is correlated with the local solute–solvent hydration energy density. The interplay among polar, apolar, and geometric factors forms the energy minimum at the solute–solvent interface. Inside of the hydrophobic binding cavity, the averaged water density is lower than the bulk, and water molecules are heterogeneously distributed inside the pocket. The water occupancy in the concave binding pocket fluctuates between dry and wet states as a function of interdomain distance. This fluctuation has important implications to the kinetics of the p53/MDM2 binding process.

As one of the most studied tumor proliferation pathways, detailed molecular binding mechanism and kinetics are of ultimate interest for developing therapeutics for activating p53. Although it is generally known for the hydrophobic nature of binding, we illustrate that the balance of hydrophobicity, shape, and polarity plays a key role in hydration. Numerous successful drugs have also demonstrated that balanced hydrophilic and hydrophobic properties are critical for drug development.[82] Our results hopefully can provide a rational strategy to the design of new drugs with balanced hydrophobic and polar properties and taking advantage of the bimodal solvation process in the p53/MDM2 complex. For instance, one can optimize molecules that minimize the probability of bimodal solvation to gain fast onset kinetics.

This study also illustrates the physical connection between the newly developed VISM model and explicit water simulations. The VISM surface mostly corresponds to the first high water density shell around the protein. In a protein binding site with complex geometry and charge distributions, it often defines the interface where bulk waters are greatly perturbed by the protein residues. In the case of the MDM2 binding site, it encloses the region that has 40% less water than the surrounding bulk.

There have been many attempts to locate hot spots for water on the protein surface to find the potential binding sites, from both explicit and implicit-solvent aspects.[83−87] The kinetics of hydration sites (both drying and wetting processes) strongly depends on the protein surface geometrical and physiochemical properties. For instance, the water molecules can be trapped and take a longer time to escape in deep and polar cavities than the ones in shallow hydrophobic areas. The variety of protein surfaces leads to relatively broad distributions of thermodynamics properties for hydration sites. In explicit solvent MD simulation studies, pure water solvents as well as mixed solvents are used to investigate the solvation of protein–ligand binding sites.[84] Correct hydration site predictions are only possible if the bound solvent molecules can be exchanged at a reasonably faster speed relative to the simulation time scale. Seco et al. have studied the exchange frequency by counting all the occurrences of solute–solvent contacts for eight proteins including p53 peptide and MDM2.[84] They have shown that solvent interaction withthe p53 binding site is mostly short-lived (up to 9 ns at most).[84] Beuming et al. have produced consistent convergence in the hydration site-free energies across a broad range of targets in their internal studies with 5.0 kcal/mol·Å2 harmonic restraint and 2 ns production simulations.[86] Of course, sufficient sampling is not always warranted in every system by MD simulation. For example, water molecule permeation through the Na/glucose cotransporter (SGLT1) under equilibrium require hundreds of nanosecond simulations due to cotransporting partners, long transporting channel, and protein conformational change.[88] In p53/MDM2, sampling is less problematic as the pocket is relatively shallow so that the time scale to reach equilibrium hydration is shorter than SGLT1. Although an extremely long kinetic trap for water cannot be completely ruled out, we believe that 40 ns production MD simulation after equilibration provides us sufficient hydration information for this particular system.

All our calculations are based on the constrained protein surfaces to emphasize the role of water fluctuations during the hydration processes. It is similar to the early Bagchi-Zewail model.[89−91] On the basis of their theoretical models, they revealed that the hydration dynamics directly relate to the water residence time near this protein surface. The hydration dynamics at the protein surface was characterized by two time scales, an ultrafast bulk like time scale followed by a slow one 10 times longer. The fast time scale arises mainly due to the reorientation and translation of hydrating water.[91] It has been argued that protein fluctuations contribute mainly to the slow component. Li et al.[92] showed that the protein–water relaxation time had nearly the same amplitude with the frozen protein. It indicates that the water response is not qualitatively modified by a constrained protein. Although the constrained protein surface limits rearrangements of local water networks, it does not alter the hydration dynamics qualitatively.[92] Of course, protein flexibility is still required to observe the slow component. In the current static model, important hydration effects, such as the different relaxation time scales for the pocket water due to the protein chemical and geometrical environment, the fluctuations between bound and bulk water, can be captured. It has also been shown that constrained protein simulations facilitate the grid mapping of water density.[85] Incorporating protein full flexibility will undoubtedly offer a richer and more complex picture of actual water behaviors around protein. We are actively pursuing this direction and hope to report our results in a future publication.