Variational Implicit-Solvent Modeling of Host-Guest Binding: A Case Study on Cucurbit[7]uril|

The synthetic host cucurbit[7]uril (CB[7]) binds aromatic guests or metal complexes with ultrahigh affinity compared with that typically displayed in protein-ligand binding. Due to its small size, CB[7] serves as an ideal receptor-ligand system for developing computational methods for molecular recognition. Here, we apply the recently developed variational implicit-solvent model (VISM), numerically evaluated by the level-set method, to study hydration effects in the high-affinity binding of the B2 bicyclo[2.2.2]octane derivative to CB[7]. For the unbound host, we find that the host cavity favors the hydrated state over the dry state due to electrostatic effects. For the guest binding, we find reasonable agreement to experimental binding affinities. Dissection of the individual VISM free-energy contributions shows that the major driving forces are water-mediated hydrophobic interactions and the intrinsic (vacuum) host-guest van der Waals interactions. These findings are in line with recent experiments and molecular dynamics simulations with explicit solvent. It is expected that the level-set VISM, with further refinement on the electrostatic descriptions, can efficiently predict molecular binding and recognition in a wide range of future applications.

Introduction

The synthetic cucurbit[7]uril (CB[7]) host molecule[1,2] has recently attracted experimental and theoretical attention due to its ultrahigh binding affinity of aromatic guests or metal complexes.[3−8] This host–guest system has various potential applications including catalysis, gas purification, crystal engineering, molecular machines, supramolecular polymers, drug transport, and gene transfection.[9−14] The detailed molecular reasons of the ultrahigh binding affinity are therefore of high interest and have been intensively explored.[15−17] Interestingly, CB[7] stands out in achieving high binding affinity among the cucurbit[n]uril hosts thanks to its good balance between the number of water molecules confined in the host cavity and the high energy for per water molecule induced by the incapability to form stable H-bonds network due to the limited confinement.[16] On the other hand, the host–guest system, cucurbit[7]uril (CB[7]) and B2 bicyclo[2.2.2]octane derivative (B2), is a highly suitable candidate for the development of computational approaches for biomolecular recognition because of its small size and relatively simple chemical complexity.

Theoretical studies of guest–host systems have been mainly based either on explicit-water molecular dynamics (MD) computer simulations,[6,18−20] or implicit-solvent models.[21,22] Traditional implicit-solvent models rely on the concept of solvent-accessible surfaces (SAS), solvent-excluded surfaces (SES), or van der Waals surfaces (vdWS), which define the interfaces between the solute and solvent regions in different ways.[23−27] Despite their wide and successful applications to many systems, the general accuracy and transferability are still questionable. Overall, the computational gain of eliminating explicit water is modest (a factor 3–5) compared to the sizable loss of information on hydration behavior in the solvation shells. In addition, a simple definition of SAS, SES, or vdWS often results in inaccurate estimation of the solvation free energy and fails in capturing hydrophobic interactions at molecular scales.[28]

To solve some of these problems, a new class of implicit-solvent models—the variational implicit-solvent model (VISM)[29,30]—has been developed in recent years. Coupled with the robust level-set numerical method, such a model allows a more physical description of the solute–solvent interface.[31−37] Central in the VISM is a free-energy functional of all possible solute–solvent interfaces, or dielectric boundaries, that separate the continuum solvent from all solute atoms. In a simple setting, such a free-energy functional consists of surface energy, solute–solvent van der Waals (vdW) interaction energy, and continuum electrostatic free energy, all depending solely on the geometry of a given solute–solvent interface. The minimization of the functional determines the solvation free energy and stable equilibrium solute–solvent interfaces. Our extensive numerical computations have shown that the level-set VISM can capture multiple local minima of hydrophobic hydration and is much more efficient than computer simulations (hours vs days).[31−37] Other related works based on variational formalisms of solvation have also recently shown advances in the description of solvation free energies.[38−40]

In this work, we apply the level-set VISM to the guest B2 binding to the host CB[7]; see Figure 1 for an illustration of the system. We use the efficient Coulomb-field approximation (CFA) with an empirical definition of the dielectric boundary to approximate the electrostatic part of the solvation free energy.[36,37] We calculate surface hydration maps, describe the curvature properties, and investigate the binding free energy along the reaction coordinate of the symmetry axis of the system. Our VISM estimates of the total free energy of binding are in good agreement with experimental and previous computational results. We further separate the individual contributions in the potential of mean force and show that hydrophobic water effects and host–guest intrinsic (vacuum) van der Waals (vdW) interactions are the two dominant driving forces for binding. A key feature of the host CB[7] is its toroidal shape. This leads to complex hydration patterns that are hardly describable by traditional, fixed-surface type implicit-solvent models. Our level-set VISM can, however, capture details of such complex patterns systematically. The extensive numerical computations reported here demonstrate particularly that the level-set VISM is promising in the efficient study of noncovalent molecular binding. With further refinement, it may thus constitute a powerful tool for calculating insightful hydration patterns and accurate binding free energies in a wide variety of aqueous molecular systems.

Figure 1

Bicyclo[2.2.2]octane derivative (B2) and cucurbit[7]uril (CB[7]). Left: Side view of the bound state. Middle: Top view of the bound state. Right: The distance d (Å) between the geometric centers of CB[7] and B2.

The rest of this paper is organized as follows: In section II, we briefly recall the VISM, the method of shifting a dielectric boundary, and the VISM level-set implementation. We also describe the reaction coordinates and the potential of mean force for the system CB[7]–B2. In section III, we report our computational results of applying our theory and method to the host–guest system CB[7]–B2. Finally, in section IV, we draw our conclusions.

Theory and Methods

Variational Implicit-Solvent Model

We divide geometrically the underlying system into three parts: the solute region Ωm, the solvent region Ωw, and the corresponding solute–solvent interface Γ, cf. Figure 2. We use this interface as the dielectric boundary. We assume that there are N atoms constituting the solute molecules that are located at positions 1,..., inside Ωm and carry point charges Q1,...,Q, respectively.

Figure 2

Schematic view of a solvation system with an implicit solvent. The solute region Ωm and the solvent region Ωw are separated by a solute–solvent interface (i.e., the dielectric boundary) denoted by Γ. The solute atoms are located at x1,...,x, carrying partial charges Q1,...,Q, respectively.

In the VISM, the equilibrium interface Γmin is obtained by minimizing the solvation free-energy functional of all possible solute–solvent interfaces ΓHere, the term P vol (Ωm) is the energy of creating the solute cavity, with P the difference of pressures inside and outside the solutes. This term can be neglected for systems at nanometer scale at a normal pressure condition (P = 1 bar). The surface integral term is the surface energy, where γ is the liquid–vapor surface tension. For small solutes with large curvature, a curvature correction typically applies, which we assume is linear in mean curvature.[41,42] The surface tension is therefore given bywhere γ0 is the surface tension for a planar interface, τ is the curvature correction coefficient or the Tolman coefficient,[41] and H is the mean curvature defined as the average of the two principal curvatures. We denote by Ggeom[Γ] the sum of the first two terms in eq 2.1, and call it the geometrical part of the free energy.

For each i, U(|x – x|) in the ith volume integral in eq 2.1 is the potential of vdW interaction between the solute particle at x and a solvent molecule at x. As we model the solvent molecules as a continuum material with a homogeneous bulk density ρw, the interaction between the ith solute particle and the solvent is then described by a volume integral over the solvent region Ωw. As usual, each U is taken to be a pairwise Lennard-Jones (LJ) potentialThe parameters ε of energy and σ of length can vary with different solute atoms as in conventional force fields. We shall denote by GvdW[Γ] the third term (i.e., the summation term) in eq 2.1, and call it the vdW part of the free energy.

The last term Gelec[Γ] in eq 2.1 is the electrostatic part of the solvation free energy, which is defined by the Born cycle[43] as the difference of the electrostatic energies of the vacuum state and the solvated state. Our approach to electrostatic energies is based on the Coulomb-field approximation (CFA) of the electric displacement field, where solvent polarization effects are neglected.[36,37] With such an approximation, we can express the electrostatic part of the solvation free energy as a simple volume integral over the solvent region that can be arbitrarily shaped:where ε0 is the vacuum permittivity and εw and εm are the relative permittivities of the solvent and solutes, respectively. This is similar to the generalized Born model.[44,45] However, in our approach, we do not need to compute generalized Born radii, introducing therefore no additional parameters in the calculation. For methodological details and performance of the CFA we refer to our previous work.[36,37]

In summary, our VISM free-energy functional with the CFA of electrostatics is given by

Effective Dielectric Boundary

It is well-known that in implicit-solvent descriptions the accuracy of continuum electrostatics depends sensitively on the choice of a dielectric boundary.[21,30,37,46−48] Such a boundary does not necessarily overlap with a solvent-accessible or solvent-excluded boundary. It is believed that the effective position of a dielectric boundary is related to the first peak of the charge distribution of the solvent surrounding a charged molecule.[46,49,50] However, the charge distribution depends on solute charge. Consider, for example, two oppositely charged ions of the same size with an equal amount of charges. Due to the different signs of charge, the first peak of charge distribution of the asymmetric and dipolar water molecule surrounding one ion is different from that of the other ion. This charge asymmetry implies that the dielectric boundary for the solvation of one of the ions in water should be different from that for the other ion, possibly located closer to the anion for which water hydrogen atoms point inward. To address this issue, we employ an empirical parallel shift of the VISM surface toward the solute region by a fitting parameter ξ, when we calculate the electrostatic part of the free energy. The amount ξ should be comparable to the size of the solvent molecule. The assumption underlying such an adjustment is that the effective dielectric boundary is closely related to the VISM surface. In our recent work,[37] we find by fitting the solvation free energy of a large set of small organic molecules of mixed anionic and cationic nature that an average value of ξ close to 1.4 Å works quantitatively well. In this work, we consistently find that ξ = 1.35 Å is well suited for a reasonable description of the electrostatics of the host–guest system, based on comparisons to previous implicit-water calculations.[5] This indicates that the value of ξ is transferable and context-independent.

Level-Set Relaxation

To numerically minimize the free-energy functional eq 2.4, we employ the level-set method that we developed previously.[31−37] We begin with an initial surface that encloses all the solute atoms located at x1,...,x. We then move the surface in the direction of steepest descent of the free energy by solving the level-set equation until a steady state 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}.[51−53] With this representation of the surface, the unit normal n = n(x), the mean curvature H = H(x), and the Gaussian curvature K = K(x) of a point x at the surface are then given by n = ∇ϕ/|∇ϕ|, H = (1/2)∇·n, and K = n·adj (∇2ϕ)n, respectively, where ∇2ϕ is the 3 × 3 Hessian matrix of the function ϕ whose entries are all the second order partial derivatives ∂2ϕ of the level-set function ϕ, and adj (∇2ϕ) is the adjoint matrix of the Hessian ∇2ϕ. 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 equationwhere v = v(x,t) is the normal velocity (i.e., the normal component of velocity) of a point x on the surface at time t.

When we apply the level-set method to relax the VISM free-energy functional G[Γ], we choose the normal velocity v to be proportional to the normal component of the effective boundary force that is given by – δΓG[Γ], the negative functional derivative of the free-energy functional G[Γ] with respect to the location change of the boundary Γ. We set the constant of proportionality, which is the mobility, to be the unity. This normal component of the boundary force is given by[31,32,35,36,54,55]

As the VISM free-energy functional is quite nonconvex due to the possible energy barriers existing in an underlying molecular system,[56] different initial surfaces can lead to different local minimizers that are of practical interest. In order to capture the latter, we design two types of initial solute–solvent interfaces. The first one is a tight wrap: a surface that is close to the vdW surface of the solute atoms. The second one is a loose wrap: a surface that loosely encloses all the solute atoms. Typical choices of the loose and tight initial surfaces for the host only and for the host–guest bound state, respectively, are displayed in Figure 3.

Figure 3

Initial surfaces in the level-set VISM computations. (a) A loose initial surface for the host. (b) A tight initial surface for the host. (c) A loose initial surface for a host–guest bound state. (d) A tight initial surface for a host–guest bound state.

Host–Guest System and Force-Field Parameters

We model the CB[7] host and B2 guest as rigid entities along the distance reaction coordinate d as shown in Figure 1. Equilibrium coordinates of the host–guest bound state are taken from a recent MD simulation study.[57] Notice that the relative orientation between the host and guest, based on such equilibrium coordinates, is used along the reaction coordinate in our VISM calculations. Hence, the only force-field parameters that have to be defined are the atomic partial charges and the solute–water LJ interaction parameters. The partial charges of CB[7] and B2 are calculated using the RESP program.[58,59] To ensure symmetry, charge equivalence is enforced on each one of the seven units of the cucurbituril host. The molecular electrostatic potential is calculated at the HF/6-31G* level. In order to investigate the effects of charge on hydration and free energy, we will study both charged and uncharged cases. In the uncharged case, we set the values for all the charges of the host and guest to be zero.

The LJ parameters for CB[7] and B2 are assigned using the academic CHARMM force field.[60,61] We use the TIP4P water model to determine the water parameters and employ the Lorentz–Berthelot mixing rules to define the LJ interaction between the water and the individual solute atoms in our VISM functional (eq 2.1).[42] Notice that we assign only one vdW interaction to a single water molecule and use LJ parameters of the oxygen atom to represent water molecule. Specifically, we use T = 300 K, P = 1 bar, γ0 = 0.1315 kBT/Å2, ρw = 0.0331 Å–3, εm = 1, and εw = 71. The unit of energy is kBT. We choose the Tolman coefficient τ = 0.76 Å appropriate for TIP4P water.[42] With this set of parameters, the VISM theory is complete and used essentially as a fit-parameter free water model in this work.

Potential of Mean Force

We study the effective interaction or potential of mean force (PMF) between the two binding partners along a simple reaction coordinate d. We shall define the PMF to include the Coulomb and vdW interactions between two binding partners in vacuum. We define this reaction coordinate to be the distance between the geometrical centers of the host and guest molecules, respectively, cf. Figure 1. The reference state is conveniently chosen with dref = ∞, i.e., the atoms in the host molecule are at infinite separation from those in the guest molecule. The relative positions of all atoms in the host or the guest are all fixed for every d. We consider the position of the guest molecule along the symmetry axis that passes through the center of the host and is perpendicular to the nearly circular host molecule. We allow the guest to move along this axis of reaction coordinate to the entrance of the circular host in both negative and positive directions, as shown in Figure 1. This pathway may not be the most probable pathway in the real system. However, close to a binding state (small d), this geometry is meaningful. The analysis of such a pathway can shed light on hydration mechanisms before binding, such as imposed energetic desolvation barriers.

Consider a VISM optimal (i.e., free-energy minimizing) surface Γ that corresponds to a reaction coordinate d. We denote by Γξ the dielectric boundary that is obtained by shifting the VISM surface Γ by ξ toward to the solutes CB[7], or B2, or the combined system. As in our previous work,[36] we define the (total) PMF as the sum of its separate contributions[36]Different parts, Ggeom, GvdW, and Gelec, of the free-energy functional G are defined in eqs 2.1 and 2.3. Effectively, we evaluate Ggeom[Γ∞] as follows. We fix the locations and partial charges of all CB[7] atoms and calculate the optimal VISM surface ΓCB[7] by minimizing the corresponding free-energy functional GCB[7] which only has the CB[7] atoms. Similarly, we obtain the optimal VISM surface ΓB2 by minimizing the corresponding free-energy functional GB2, which only has the B2 atoms. We then calculateThe other terms related to Γ∞ or Γ∞ξ can be calculated similarly. The notation x ∈ CB[7] and x ∈ B2 mean that x is a CB[7] atom and x is a B2 atom, respectively. The double-sum terms in GvdWpmf(d) and Gelecpmf(d) account for the vacuum vdW interactions and charge–charge interactions of atoms in CB[7] and those in B2, respectively.

The VISM binding free energy is given by Gtotpmf(0). Note that the contributions of configurational entropy are not included in our model. We use fixed coordinates of the host–guest system and thus neglect the translational and configurational fluctuations which are significantly restricted upon binding. However, an accurate estimate of these contributions exists[61] and will be added to the VISM binding energy to obtain the total binding energy for comparison to experimental values.

We remark that for a given reaction coordinate d there can be multiple stable equilibrium interfaces Γ that are local minimizers of the VISM free-energy functional. Different local minimizers Γ for the same coordinate d define multiple values G[Γ] of VISM local minimum free energies. Therefore, the PMF can have multiple branches along the reaction coordinate d and, hence, can lead to hysteresis. Once these local minimizers have been obtained, an ensemble averaging can be calculated to estimate the free energy of the host–guest system. Suppose there are M equilibrium hydration surfaces, {Γ}, at distance d. The actual free energy of the system can then be expressed as[33,34]where G0 is a constant such that G(+∞) = 0.

Results and Discussion

Equilibrium Interface of the Isolated Host and Bound Host–Guest

In Figure 4, we show the final equilibrium interfaces of the isolated host obtained by our level-set VISM with both tight and loose initial surfaces. We observe that the loose initial surface leads to a final equilibrium dry state where no water molecules are in the host cavity and that the tight initial surface leads to a wet state where the cavity is hydrated. Such a bimodal hydration behavior has been previously also observed in concave hydrophobic pockets due to capillary evaporation.[34] Due to the toroidal shape and narrow confinement featured by the host such a capillary mechanism may play a role in this system. Note, however, that those mechanisms are highly sensitive to local vdW or electrostatic potentials.[29,30] In Figure 5a the VISM surface is shown for the bound host–guest system. There, for both loose and tight initial surfaces, we obtain the same equilibrium surface. This independence of the initial surface is reflected also in the same value of total VISM solvation free energy listed in Table 1 and 2.

Figure 4

Equilibrium surfaces of the host with electrostatics starting from (a) a loose initial and (b) a tight initial. Color represents mean curvature, where red renders convex curvature and blue renders concave curvature. Green interpolates the red and blue and essentially represents flat patches.

Figure 5

Equilibrium surfaces of the host–guest system CB[7]–B2 in the bound state for the case of (a) including electrostatics and (b) no electrostatics. In both cases, both choices of the initial surface (loose or tight) reach the same equilibrium state. The color code is the same as in Figure 4.

Table 1

Individual Solvation Contributions to the Host–Guest Solvation Free Energy with Loose Initial Surfaces (units kBT)

systems	Ggeom[Γ]	GvdW[Γ]	Gelec[Γ]	total solvation free energy
noncharged CB[7]	91.8	–88.0	0.0	3.8
charged CB[7]	91.2	–86.3	–142.6	–137.7
noncharged CB[7]–B2	91.4	–93.2	0.0	–1.8
charged CB[7]–B2	90.5	–91.0	–156.4	–156.9
charged B2	28.0	–18.3	–26.8	–17.1

Table 2

Individual Solvation Contributions to the Host–Guest Solvation Free Energy with Tight Initial Surfaces (units kBT)

systems	Ggeom[Γ]	GvdW[Γ]	Gelec[Γ]	total solvation free energy
noncharged CB[7]	104.9	–97.2	0.00	7.7
charged CB[7]	103.9	–95.2	–164.5	–155.8
noncharged CB[7]–B2	91.4	–93.2	0.00	–1.8
charged CB[7]–B2	90.5	–91.0	–156.4	–156.9
charged B2	28.0	–18.3	–26.8	–17.1

To evaluate charge effects on hydration, we also present the interfaces with electrostatic solute–solvent interactions switched-off in the numerical computation. In Figure 5b and Figure 6, we show the final equilibrium interfaces of the bound host–guest system and the isolated host without electrostatics. Comparing the equilibrium surfaces obtained with and without electrostatics, we observe that the electrostatic interaction attracts polar water molecules to the charged solutes, resulting much tighter final surfaces.

Figure 6

Equilibrium surfaces of the host without electrostatics. (a) Calculations start from a loose initial. (b) Calculations start from a tight initial. The color code is the same as in Figure 4.

The hydration behavior described by VISM is consistent with recent explicit solvent MD simulations of the identical apolar and polar systems.[57] Results for the MD-derived hydration maps for the unbound host with and without electrostatic interactions are shown in Figure 7. Clearly, in the nonpolar case, the hydration in and around the host is much weaker than in the electrostatic case. The MD simulations indeed found enhanced fluctuations in the toroidal confinement of the host,[57] and the host was found mostly in the wet state. These behaviors are very well predicted by our level-set VISM computations. For bound system, we also observe that, analogous to the unbound case, the electrostatic interaction attracts water molecules closer to the solute atoms. As shown in Figure 7, with the electrostatics, the water density is much higher in the vicinity of the solute and the water molecules distribute closer to solute atoms, in both the bound and unbound systems.

Figure 7

Effect of electrostatics on hydration: Top view and side view of average water density maps from MD simulations of bound (left) and unbound (right) host–guest systems, with the partial charges of the host included (blue) and removed (red).

In Tables 1 and 2, we list values of different parts of the solvation free energy for different systems, charged or uncharged, and isolated or bound, corresponding to loose and tight initial surfaces, respectively. For the host only system, we observe that the solvation free energy of the wet state is much lower than that of the dehydrated state, with a difference of 18.1 kBT. This implies that the wet state is more favorable for the host cavity. Comparing the energy differences of each contribution of the solvation free energy, we also see that the electrostatics plays a more important role than the nonpolar parts. Moreover, the host only system switches its favorable hydration state from the dry state to the wet state owing to the switched-on electrostatics.

PMF and Binding Affinity

In Figures 8 and 9 we show the equilibrium surfaces for the system with varying distance d between CB[7] and B2 along the reaction coordinate that we defined in II E, using loose initial surfaces and tight initial surfaces, respectively. When the guest is at or close to the center of the host, our level-set relaxation reaches the same final equilibrium surface, with a loose or tight initial surface, cf. Figure 8a and Figure 9a. With a tight initial surface, the VISM equilibrium surface of the host–guest system CB[7]–B2 gradually breaks up when the center-to-center distance between CB[7] and B2 becomes greater than 10 Å. In contrast, with a loose initial surface, the final VISM equilibrium surface of the host–guest system still sticks together until their center-to-center distance exceeds 11 Å. This is due to the bimodal hydration behavior of the host cavity. Note that for small distances the part of surface at the intersection of host and guest has a high concave curvature, showing squeezed and probably unfavorable water hydration.

Figure 8

Equilibrium surfaces of the CB[7]–B2 with loose initial surfaces, with varying center-to-center distance d (in Å) along the reaction coordinate. The color code is the same as in Figure 6.

Figure 9

Equilibrium surfaces of the CB[7]–B2 with tight initial surfaces, with varying center-to-center distance d (in Å) along the reaction coordinate. The color code is the same as in Figure 6.

In Figure 10, we plot the total VISM solvation free energy and its components along the reaction coordinate d. At each distance d, we show both the VISM values corresponding to a loose initial and tight initial surfaces. A negative value of d means that B2 is placed from the side of CB[7] different from that with a positive value of d. Note that these energy curves are symmetric in d → −d, resulting from the symmetry of the host–guest system CB[7]–B2. Moreover, we observe clearly differences in these energy curves between loose and tight initials. Such differences originate from the bimodal hydration behavior that also found in the isolated (unbound) host alone. For the unbound host, the wet state is always favored over the dry one when a tight initial surface is used. The difference of the solvation free energy of bound (d = 0) and unbound cases (|d| = 14) is about 16.3 kBT. Thus, the solvation part of the total PMF strongly disfavors binding. This is consistent with the previous prediction[61] of a 20.1 kBT difference.

Figure 10

Plots along the reaction coordinate of (a) the minimum VISM solvation free energy G[Γ], (b) the geometric part Ggeom[Γ], (c) the solute–solvent vdW interaction part GvdW[Γ], and (d) the electrostatic part Gelec[Γ] of the total solvation free energy G[Γ].

The individual contributions to the solvation free energy are shown in Figure 10b, c, and d. We see that the surface energy strongly drives the binding, due to the desolvation of the unfavorable solute upon binding. The difference in end-point energies, that is, between bound cases (d = 0) and unbound cases (|d| = 14), is large and favorable: −28.9 kBT for the loose initial. The corresponding difference for the tight initial is larger, about −41.2 kBT, because unfavorably bound water is freed upon binding. This is in contrast to the loose initial case where the host was already void of water before binding. The solute–water vdW interaction of the solvation free energy, shown in Figure 10c, disfavors binding by roughly 12.9 kBT for a loose initial surface. The reason for this energy penalty is that favorable solute–water vdW interactions are reduced upon binding. Again the energy value depends on loose and tight initial surfaces. With a tight initial surface, VISM predicts that water molecules are released upon binding. Therefore, the energy penalty is larger, up to 21.8 kBT. This is in line with the conclusion drawn from MD simulations that the release of the unfavorable water inside the host cavity plays an important role upon binding.[16,17] The nonpolar part (i.e., the sum of surface contribution and the solute–water vdW interaction) of the end-point free energy is −15.7 kBT for the loose case and −19.4 kBT for the tight case. This predicts that the nonpolar part of solvation favors the host–guest binding, agreeing qualitatively with the value of −4.4 kBT as in previous implicit water calculations.[61]

We now consider the effect of electrostatics on the host–guest binding. Figure 10d shows the electrostatic part of the solvation free energy. Again, the energy value depends on loose and tight initial surfaces. We observe that the electrostatics disfavors binding, with 13.6 kBT energy penalty for the loose initial and 35.7 kBT energy penalty for the tight initial. Notice that our electrostatic part of the PMF includes the Coulombic interaction between the host and guest in vacuum. Simple calculations indicate that such interaction has a favorable contribution of a value up to −12.2 kBT to the binding free energy. Figure 11 displays the electrostatic part of the PMF. It shows an electrostatic penalty of 1.4 kBT on binding for the loose case and 23.6 kBT for the tight case, in contrast to 10.6 kBT reported in the previous work.[61] All of these data imply that the penalty that stems from the electrostatic part of the solvation free energy is partially balanced by the attractive Coulombic interaction in vacuum. Such balance leads to a weaker, unfavorable electrostatic contribution to the binding free energy.[4,5,61]

Figure 11

Electrostatic part of the PMF.

In Figure 12, we show the host–guest vdW interaction in vacuum. It naturally vanishes for large distances while it displays a highly favorable vdW attraction of about −38.7 kBT in the bound state. Moghaddam et al.[61] calculated −57.9 kBT. The discrepancy may be due to differences in the positional coordinates of the bound host–guest complex. In our d-resolved energy a high barrier is found at |d| ≃ 3 Å, stemming from too close atomic overlaps at these distances. This is a consequence of our chosen reaction pathway and must be considered artificial. While future studies should consider relaxation and fluctuations along the pathway, we believe that the solvation free energy discussed above does not suffer from these intrasolute overlaps as they are restricted to the internal cavity regions not in touch with water.

Figure 12

Vacuum host–guest vdW interaction energy.

Figure 13 shows the total PMF profiles including the vacuum vdW part, obtained with loose initials, tight initials, as well as their ensemble averages. The averaged profile almost overlaps with the tight case, indicating that the results with tight initials are much more favorable. The end-point free-energy difference for the ensemble averaged profile is about −35.3 kBT, comparing with −49.6 kBT reported in the work of Moghaddam et al.[61] The desolvation barrier around d = 10 Å is physical while the large barriers at d ≃ 4 Å are likely to be artificial, as discussed above.

Figure 13

Potential of mean force.

For the comparison to experiments, the configurational entropy penalty upon binding to the end-point free-energy differences has to be added to our total PMF to consider the restriction of configurational degrees of freedom upon binding. Moghaddam et al.[61] estimated that these contribute to a large free-energy penalty, about 29.4 kBT, including the correction to account for the experimental standard free energy. Thus, the total binding free energy of about −5.9 kBT predicted by our level-set VISM has to be compared to the experimental value from titration measurements of about −22.6 kBT.[61] The gap between two predictions can be attributed to the LJ force field parameters or the approximations in the CFA. Overall, our VISM results agree qualitatively with experiments and MD simulations.

It has been reported extensively in literature that the hydrophobic effect, size complementary, and host–guest direct interaction can all account for the high binding affinity for Cucurbit[n]urils.[4,5,16,61] Moreover, it has been shown that the expulsion of perturbed water from the host cavity is one of the major driving forces for high affinity binding.[16,61] According to our numerical computations, the surface (water induced) part and direct host–guest interaction in vacuum, including vdW contributions and electrostatic contributions, play important roles in binding.

Conclusions

We have applied our level-set VISM to the host–guest system CB[7]–B2. The toroidal shape of the host gives rise to complex hydration patterns within the host cavity. Some of such patterns can be hardly captured consistently by traditional, fixed-surface type, implicit-solvent models. Our VISM descriptions of the host hydration is in line with recent computer simulations which point to highly perturbed water molecules in the host binding pocket.[16,17] The dry state and wet state predicted by our level-set VISM with loose initials and tight initials correspond to different local minima of the free-energy functional. These explain the enhanced fluctuations of water molecules inside the host cavity in MD simulations. A similar behavior has been recently found in the simulation of a ligand binding to a purely nonpolar hemispherical pocket[34] with a remarkable enthalpic thermodynamic signature.[62]

To show the influence of electrostatics on the hydration, we have performed numerical computations on the host–guest system with partial charges switched-off and switched-on, respectively. Our computational results have shown that the electrostatic interaction enhances wetting by pushing the water molecules to charged atoms, leading to a much tighter hydration surface. These are consistent with the MD simulations and our previous level-set VISM results.

The VISM has been in general successful in qualitative predictions of the binding free energies. We have found that the surface energy and the direct host–guest interactions in vacuum account most for the host–guest binding, while the electrostatic part of the solvation and solute–water vdW interactions disfavor the binding. The ensemble averaged total PMF obtained by the VISM predicts a favorable free-energy binding affinity. All of these agree qualitatively with existing studies in the literature.[4,5,16,17,61]

We discuss now several issues related to our approach and applications to the system CB[7]–B2. First, we have used the CFA for the electrostatic part of the free energy. Although such an approach is quite accurate qualitatively as we have found in this work, it does not include the effects of water polarization by dielectric voids nor those of ionic screening. Therefore, the description of the electrostatic contribution needs to be improved. A natural way to do so will be to incorporate the Poisson–Boltzmann (PB) equation in our level-set VISM. Note that we have used a boundary shift as a parameter to deal with the issue of a proper definition of the dielectric boundary. However, this issue remains also in the PB approach. More systematic studies are needed to resolve this issue in our level-set VISM as well as in other implicit-solvent models.

Second, a reasonable pathway (i.e., reaction coordinate) for the guest binding to the host should be introduced. In this work we have chosen the symmetry axis perpendicular to the circular molecules of the fixed host as coordinate for binding a rigid guest. We have found a higher barrier in the profile of the host–guest vdW interactions at d ≃ 3 Å, cf. Figure 12, probably due to artificially overlapping atoms in the approaching host and guest. This indicates that our binding pathway is not most reasonable. It will be for future studies to identify a probable binding pathway.

Third, the configurational change of the solutes, such as translation, vibration, and rotation, should be included in the VISM framework. It is reported that there is a large configurational entropy penalty upon binding. In our current model, however, the host and guest are treated as static rigid bodies. We should consider these motions of the solutes in our future work.

Finally, the occurrence of hydration fluctuations and polymodal hydration behavior requires a treatment of surface fluctuations in the level-set VISM model. Only with that a proper equilibrium (Boltzmann-weighted) sampling of all solvation states is possible. This remains one of the grand challenges for implicit-solvent modeling.

On the computational side, we note that each of our level-set VISM relaxation usually takes minutes to hours. The computational time can very much depend on the initial surfaces and numerical resolutions in the level-set relaxation. It is clear that our level-set VISM approach has not yet been as efficient as generating a SAS, SES, or vdWS of a charged molecule. However, the VISM calculations are clearly much more efficient than MD simulations: hours to days. Moreover, when compared to MD simulation results, our level-set VISM can predict qualitatively well the polymodal hydration free energies, dry and wet states, and binding affinities. These are rather difficult to be described by SAS, SES, or vdWS approaches in systems with complex hydration patterns. We are currently working to improve our computational efficiency by parallelizing our codes.