Local Ion Densities can Influence Transition Paths of Molecular Binding

Abstract

Improper reaction coordinates can pose significant problems for path-based binding free energy calculations. Particularly, omission of long timescale motions can lead to over-estimation of the energetic barriers between the bound and unbound states. Many methods exist to construct the optimal reaction coordinate using a pre-defined basis set of features. Although simulations are typically conducted in explicit solvent, the solvent atoms are often excluded by these feature sets-resulting in little being known about their role in reaction coordinates, and ultimately, their role in determining (un)binding rates and free energies. In this work, analysis is done on an extensive set of host-guest unbinding trajectories, working to characterize differences between high and low probability unbinding trajectories with a focus on solvent-based features, including host-ion interactions, guest-ion interactions and location-dependent ion densities. We find that differences in ion densities as well as guest-ion interactions strongly correlate with differences in the probabilities of reactive paths that are used to determine free energies of (un)binding and play a significant role in the unbinding process.

Keywords: SAMPL system; binding affinity; free energy; ligand unbinding; mechanisms; molecular dynamics; weighted ensemble.

FIGURE 1

The OA-G3 and OA-G6 systems. The OA host molecule (left). The G3 (top right) and G6 (bottom right) guest molecules.

FIGURE 2

General WE Framework. Every circle represents a trajectory in the ensemble. Colors represent conformations and circle size represents probability, with all trajectories beginning with the same conformation and probability. Trajectories are run for a predetermined number of steps (dynamics), followed by a resampling step containing merging and cloning procedures. This cycle repeats until the end of the simulation.

FIGURE 3

Analysis of t ₀ poses. (A) The OA host molecule with the G6 ligand in the starting pose (multi-color) and example t ₀ poses (pastels). Some atoms from the host have been removed in A for clarity. (B) Average probabilities from -30 cycles to the final unbinding event organized by unbinding probability for the 2020 OA-G6 data set. (C) The average number of cycles between t ₀ and the unbinding event for OA-G3 (blue) and the 2020 OA-G6 data set (gray) organized by unbinding probability.

FIGURE 4

Feature Analysis. (A) A visualization of the region of space considered for the guest-ion features using the G6 ligand. The maximum distance for the 5 Å scale is in gray and the 3 Å scale is in blue (top). The two logistic functions used to determine the molecule counts (bottom). (B–E) Molecule counts for Na⁺ ions with results organized by both time and exit point probability. The legend in C applies to all four plots. The average total ion count (5 Å scale) around the upper negative charges of the host for (B) OA-G6 and (C) OA-G3. The average total ion count (5 Å scale) around the guest for (D) OA-G6 and (E) OA-G3. OA-G6 results correspond to the 2020 data set.

FIGURE 5

Ion Density Analysis. (A) A diagram showing the simulation box and the cylindrical space above the host where the number of ions (η(t)) is determined. The equation for calculating the autocorrelation of this quantity (C(τ)) is shown. (B) An autocorrelation plot of the cylindrical ion density (η(t)) is calculated using all reactive and non-reactive trajectory data. (C) The average number of ions in the cylinder space above the host for OA-G3 (dark blue), OA-G6 (2020, gray), and OA-G6 (2018, cyan). The average from 4500 random simulation cycles is shown in red.

FIGURE 6

Exit Point Analysis. (A) Unbinding event locations for exit points with probabilities 10^–7 (red) and 10^–12 (gray VolMap) for OA-G6. (B) Unbinding event locations for exit points with probabilities 10^–6 (red) and 10^–12 (gray VolMap) for OA-G3. The surfaces show a density contour (Isoval) of 0.0001 in both panels.