Improving sampling of binding free energy differences between covalently bound ligands in alternate binding pockets using MT-REXEE

Abstract

The primary limitation for the application of alchemical free energy methods to a wider variety of complex molecular systems is achieving reasonable sampling. Flexible binding complexes often have high free energy barriers, which require prohibitively long simulations to sample sufficiently to obtain reliable free energy estimates. An example of such a system is the complex formed between FabB, an elongating β-ketoacyl-acyl carrier protein (ACP) synthase (KS) from Escherichia coli and ACP, which carries acyl chains of varying lengths. Previous experimental evidence suggests that growing acyl chains can bind to at least two pockets. With the multiple topology replica exchange of expanded ensemble (MT-REXEE) enhanced sampling approach, we can obtain highly efficient sampling of both pockets by adaptively growing and shrinking the chains in the simulation ensemble, allowing each simulation to visit chain lengths where transitions between the pockets do occur. This enables unbiased sampling of alternate configurational states for large complex systems. Using the new swapping approach gives significantly enhanced sampling even for this simpler problem, as demonstrated by faster convergence of free energy estimates. This case study demonstrates the utility of MT-REXEE and its open-source implementation for systems that feature high free energy barriers for a subset of ligands of interest, demonstrating a valuable addition to the existing stable of enhanced sampling methods.

Figures (16)

• Figure 1: (A) A set of MT-REXEE simulations that perform transformations of the acyl substrate from C2 to C4, C4 to C6, and C6 to C8. Each simulation, which visits a range of alchemical states $\lambda$ in [0,1], has configurations that should (in the limit of sufficient sampling) visit both pocket A (orange) and pocket B (brown). The dashed boxes indicate states in each individual simulated transformation, and the dark blue arrows show where configurational swaps can occur, at the alchemical end states. (B) When we separate the transformations by end state and configurational state, we can estimate each of the FE differences we compute from a MT-REXEE simulation. We compute the alchemical FE difference for each transformation independently in each pocket (light blue) as well as the configurational FE difference at at least one end state (pink). This enables us to compute the relative FE difference between any combination of alchemical end state and configurational state.
• Figure 2: Random range swapping yields statistically equivalent free energies while significantly decreasing the round-trip time through alchemical states. (A) A diagram of the sequence of ligand transformations used in binding free energy estimates with the protein MUP1 to validate the random range swapping scheme. The curved arrows show where MT-REXEE swaps occur. (B) The structure of MUP1 in complex with ligand A. (C) We compare the binding free energy estimates using EE, MT-REXEE with the exhaustive swapping scheme, and MT-REXEE with the RR swapping scheme. We observe no significant differences between the free energy estimates for all three methods as the values are all within statistical error of one another. (D) We observe a significant (p$<$0.001) decrease in the round-trip times from 10-20 ns to 3-5 ns for both solvent and MUP1 complex simulations using the RR swapping scheme compared to the previous exhaustive scheme.
• Figure 3: (A) The binding pocket of the FabB--acyl covalent intermediate is shown with a 12 carbon acyl substrate in both pocket A (orange) and pocket B (brown) in the context of the overall protein. The structures shown are the centroids from MT-REXEE simulations of the C10-12 transformation for $\lambda=1$. (B) Occupancy of pocket B as a function of substrate chain length using both standard MD simulations and MT-REXEE enhanced sampling as described in section \ref{['sect:methods_config']}. Standard MD is not able to access pocket B within 500 ns for any substrate longer than 4 carbons, but we readily see sampling in pocket B for all chain lengths within only 100 ns of MT-REXEE simulation. (C) A close-up image of the binding pocket in (A) highlighting several residues on both the catalytically active monomer and adjacent monomer that restrict the acyl chain's conformational flexibility and prevent interconversion during MD simulations of at least 500 ns.
• Figure 4: The convergence of binding free energy simulations initiated in pockets A and B demonstrates MT-REXEE can increase sampling of the conformational space. (A) Comparison of estimates from individual independent EE calculations with estimates of the same transformations in MT-REXEE simulations. MT-REXEE's significantly improved sampling leads to moderate differences between the two methods, with MT-REXEE expected to be more accurate. Estimates for binding free energy using MT-REXEE with the initial acyl chain configurations in either pocket A (brown) or pocket B (orange) for the calculation of the RFE in pocket A (B) and pocket B (C), showing independence of starting configuration. (D) Cumulative relative binding free energy can be obtained for all chain lengths relative to a C2 chain for both pockets using simulations initiated in both pockets A and B (Figure \ref{['SI:total_FE']}). Thus, from a single set of MT-REXEE simulations, we can obtain the full chain length and pocket dependence on binding FE. For all figures, we plot the mean of n=3 (n=6 for D) replicas and error bars represent the standard error of the mean.
• Figure 5: (A) The G107M mutation (red) blocks access to a region of pocket A which destabilizes chain binding to this pocket, but has little direct affect on pocket B. We show centroids for the C10 to C12 transformation for the chain in pocket A (brown) and pocket B (orange). To assess the effect of this mutation on the binding complex, we compare the binding free energy differences by chain length for the WT and G107M variant in pocket A (B) and pocket B (C). (D) To highlight the difference in stability between the two pockets we then compare the difference in individual relative binding free energies between pocket A and B for all chain lengths for the WT and G107M variant. There is a clear bias in the mutant FabB variant against the growth of the acyl chain from 6 to 8 and from 8 to 10 carbons.