Efficient Sampling of Cavity Hydration in Proteins with Nonequilibrium Grand Canonical Monte Carlo and Polarizable Force Fields

Abstract

Prediction of the hydration levels of protein cavities and active sites is important to both mechanistic analysis and ligand design. Due to the unique microscopic environment of these buried water molecules, a polarizable model is expected to be crucial for an accurate treatment of protein internal hydration in simulations. Here we adapt a nonequilibrium candidate Monte Carlo approach for conducting grand canonical Monte Carlo simulations with the Drude polarizable force field. The GPU implementation enables the efficient sampling of internal cavity hydration levels in biomolecular systems. We also develop an enhanced sampling approach referred to as B-walking, which satisfies detailed balance and readily combines with grand canonical integration to efficiently calculate quantitative binding free energies of water to protein cavities. Applications of these developments are illustrated in a solvent box and the polar ligand binding site in trypsin. Our simulation results show that including electronic polarization leads to a modest but clear improvement in the description of water position and occupancy compared to the crystal structure. The B-walking approach enhances the range of water sampling in different chemical potential windows and thus improves the accuracy of water binding free energy calculations.

Figure 1:

A thermodynamic cycle used to evaluate the binding free energy of $N_f$ water molecules to a GCMC region in a protein (indicated by the dashed circle), $\mathrm{\Delta }{F}_{\mathit{bind}}\left(N_f\right)$.

Figure 2:

The benzamidine binding site of trypsin. (a) Location of the benzamidine binding site in trypsin (PDB ID 5MOQ). The dashed circle indicates the GCMC region, which is a spherical cavity of 6 Å radius centered at the average position of $C\alpha$ atoms of Ala221 and Gly226. (b) The binding mode of benzamidine in trypsin and nearby water molecules resolved in the combined X-ray/neutron crystal structure. Hydrogen bonding interactions are indicated by dashed lines with the distances (in Å) between the hydrogen and the hydrogen bonding acceptor atoms.

Figure 3:

Results of the hydration level in the benzamidine binding site of trypsin from ~560 ns equilibrium molecular dynamics (MD) simulations using the CHARMM-Drude 2019 force field. (a) The average radial distribution function (RDF), $g(r)$, of the water oxygen around the two benzamidine nitrogen atoms. The dashed curve indicates the integrated RDF. (b) The RDF of water oxygen around the average position of $C\alpha$ atoms of Ala221 and Gly226 (cavity center). (c) The instantaneous number of water within 6 Å from the cavity center. (d) The distribution of the number of water within 6 Å from the cavity center. (e-f) Comparison of water positions from the crystal structure with the clustering analysis results based on MD simulations. The crystal waters are shown as sticks, and the spheres represent the locations of the clusters with occupancies >10% (see Table 3). The occupancy of the cluster increases as the color changes from blue to red. Different numbers of frames from the MD trajectories are used in the analysis, skipping every $n_{\mathit{skip}}=100$ frames in (e), and every $n_{\mathit{skip}}=10$ frames in (f).

Figure 4:

Results of the hydration level in the benzamidine binding site in trypsin from GCNCMC simulations at the $B_{eq}$ value using the CHARMM-Drude 2019 force field. (a-d) The same set of radial distribution functions as in Fig. 3, but from GCNCMC simulations. (a-b): results for GCNCMC simulations with different $n_{\mathit{pert}}$ and $n_{\mathit{prop}}$ combinations; each simulation protocol is carried out with 3 replicas and each replica involves 1000 GCNCMC iterations. (c-d): results for the protocol $n_{\mathit{pert}}=149$ and $n_{\mathit{prop}}=100$ with 10 replicas, and each replica includes 1000 GCNCMC iterations; the shaded region illustrates the standard deviation of the mean. (e-f) The distribution of the number of water molecules within the GCMC region. (e): results from different $n_{\mathit{pert}}$ and $n_{\mathit{prop}}$ combinations with 3 replicas each; (f): results for $n_{\mathit{pert}}=149$ and $n_{\mathit{prop}}=100$ with 10 replicas.

Figure 5:

Comparison of water positions at the benzamidine binding site of trypsin from the crystal structure with the clustering analyses results based on different GCNCMC simulations. The crystal waters are shown as sticks, and the spheres represent the locations of the clusters with occupancies >20%. The occupancy of the cluster increases as the color changes from blue to purple to red. (a) $n_{\mathit{pert}}=99$ and $n_{\mathit{prop}}=150$; (b) $n_{\mathit{pert}}=149$ and $n_{\mathit{prop}}=100$; (c) $n_{\mathit{pert}}=199$ and $n_{\mathit{prop}}=50$.

Figure 6:

Results of grand canonical integration for water binding at the benzamidine binding site of trypsin. (a-b) represent the results from regular GCNCMC simulations, while (c-d) are the results with $B$-walking with Adam’s parameters ($B$) in the range of [−8.05, −3.55]. For the binding free energy curves (b,d), the solid line and the shaded area represent the mean and the 95% confidence interval, respectively.