Computing the Free Energy Barriers for Less by Sampling with a Coarse Reference Potential while Retaining Accuracy of the Target Fine Model

Abstract

Proposed in this contribution is a protocol for calculating fine-physics (e.g., ab initio QM/MM) free-energy surfaces at a high level of accuracy locally (e.g., only at reactants and at the transition state for computing the activation barrier) from targeted fine-physics sampling and extensive exploratory coarse-physics sampling. The full free-energy surface is still computed but at a lower level of accuracy from coarse-physics sampling. The method is analytically derived in terms of the umbrella sampling and the free-energy perturbation methods which are combined with the thermodynamic cycle and the targeted sampling strategy of the paradynamics approach. The algorithm starts by computing low-accuracy fine-physics free-energy surfaces from the coarse-physics sampling in order to identify the reaction path and to select regions for targeted sampling. Thus, the algorithm does not rely on the coarse-physics minimum free-energy reaction path. Next, segments of high-accuracy free-energy surface are computed locally at selected regions from the targeted fine-physics sampling and are positioned relative to the coarse-physics free-energy shifts. The positioning is done by averaging the free-energy perturbations computed with multistep linear response approximation method. This method is analytically shown to provide results of the thermodynamic integration and the free-energy interpolation methods, while being extremely simple in implementation. Incorporating the metadynamics sampling to the algorithm is also briefly outlined. The application is demonstrated by calculating the B3LYP//6-31G*/MM free-energy barrier for an enzymatic reaction using a semiempirical PM6/MM reference potential. These modifications allow computing the activation free energies at a significantly reduced computational cost but at the same level of accuracy compared to computing full potential of mean force.

Figure 1

Calculating the activation free-energy by positioning local PMF regions relative to the coarse-physics free-energy shifts instead of computing full PMF from fine-physics sampling. PMF calculations involve sampling with a series of biased fine-physics potentials T₁..T_M.. along the reaction coordinate (cyan circles). These mapping points are sequentially connected using the multistep free-energy perturbation ΔF(T₁ → T_M), what is depicted by cyan arrows. An alternative approach involves computing PMF locally (blue dots) using fewer fine-physics potentials (blue stars). The free-energy difference ΔF(T₁ → T_M) is instead computed from the thermodynamic cycle, shown by blue and red arrows. Red arrows show connecting a sequence of biased coarse-physics potentials R₁..R_M.. using the multistep free-energy perturbation ΔF(R₁ → R_M). Blue arrows show the free-energy perturbation of switching from the coarse-physics potential to the fine-physics potential (both having the same bias) at reactants and at the transition state, ΔF(R_i → T_i). This allows for the positioning of two local PMF regions (blue dots) and to determine the activation free-energy (shown with a black arrow).

Figure 2

Iterative accuracy improvement of the fine-physics free-energy surface computed from coarse-physics sampling. The first step involves computing a low-accuracy free-energy surface from the coarse-physics sampling by reweighting and by free-energy perturbation approaches (red upper block). All ensemble averages are computed with the coarse-physics reference potential, <···>R. Once the surface is constructed, the regions of interest (reactants and the transition state) are identified. In the next step, a limited fine-physics sampling is performed at those regions. This allows computing the ensemble averages with fine-physics target potential <···>T, which are used to calculate a high-accuracy free-energy surface at those regions (blue medium block). After the second step, the accuracy becomes limited by positioning the local regions. In the third step, a limited sampling is performed in the same regions with a linear combination of the two potentials. With the use of the generated data, the free-energy perturbation is further improved by including the corresponding average for the intermediate ensemble.

Figure 3

Active center of haloalkane dehalogenase. The protein backbone is shown in green; the catalytic residue, aspartate, and the substrate are shown as sticks. The dotted red line depicts the direction of the nucleophilic attack on the substrate, 1, 2-dichloroethane.

Figure 4

Reaction scheme for the studied reaction.

Figure 5

Selecting regions for targeted fine-physics sampling from low-accuracy B3LYP/MM free-energy surfaces computed from sampling with PM6/MM reference potential for the benchmark enzymatic reaction. These surfaces are used for detecting reactants and the transition state (shown by shaded gray areas). The first surface (shown in green) is computed by reweighting the reaction coordinate distribution from PM6/MM sampling using the umbrella sampling method. Another low-accuracy surface (shown in black) is calculated by adding the linear expansion of the free-energy perturbation to the PM6/MM high-accuracy PMF.

Figure 6

Improving the accuracy of probability distributions of the reaction coordinate at selected regions of the reaction path by performing local sampling with B3LYP/MM potential for the benchmark system. The upper plot shows local regions of the free-energy surface computed from the targeted sampling: green and magenta triangles show free-energy penalties of incrementally moving along the reaction coordinate relative to the first (leftmost) bias at each region. They are generated with multistep linear response approximation (M-LRA) and with WHAM, respectively. Red and blue lines are the local PMF computed using the weight-averaged umbrella sampling approach (red) and WHAM (blue), respectively. The lower plot shows biased distributions of the reaction coordinate for the data points used in calculations.

Figure 7

Improving accuracy of the free-energy shifts associated with changing the center of bias ξ_m⁰ while moving along the reaction coordinate. All changes are computed relative to the leftmost (first) bias. Red stars correspond to free energies computed from sampling with biased PM6/MM reference potential (R) using a multistep free-energy perturbation. The cyan circles represent the exact searched solution and correspond to free energies computed from sampling with a biased B3LYP/MM target potential (T) using a multistep free-energy perturbation. Green stars correspond to the PM6/MM free-energy shifts to which the linear expansion of the free-energy perturbation to B3LYP/MM potential is added (at corresponding biases). This perturbation is computed as PM6 ensemble average of the energy gap between B3LYP and PM6, ΔE. These shifts are used in constructing the low-accuracy B3LYP free-energy surface. After performing a limited fine-physics sampling with B3LYP target potential (T), a more-accuracte estimate for these shifts is computed (purple ▲) using two-step linear response approximation. These estimates are used for positioning the local PMF regions (see the main text). Next, the estimates are improved further with a three-step perturbation (blue ▲), which involves computing the average energy gap with intermediate potentials M = 0.5R + 0.5T.

Figure 8

Increasing accuracy of the B3LYP activation free-energy barrier computed by positioning the local PMF regions from the targeted sampling and its comparison to the full PMF. Full B3LYP/MM PMF (dashed gray line) is computed from 64 molecular dynamics B3LYP trajectories. Black line is a low-accuracy PMF computed from 64 PM6 trajectories using the averaged energy gap with the B3LYP potential. It is used to identify regions for targeted sampling with B3LYP. Local B3LYP PMFs are computed from this targeted sampling, which includes 5 trajectories at reactants, 5 trajectories at the transition state, and 3 trajectories at products (blue line), but the relative positions of these regions is not known. They are found by computing the free-energy perturbations of switching from PM6 reference potential (R) to B3LYP target potential (T), shown by black arrows. The local PMF regions, positioned using a two-step linear response approximation, are shown in red. Three-step free-energy perturbation slightly further improves the positioning, and the activation free-energy estimate (shown in green). The three-step estimate involves additional targeted sampling at the same regions with 0.5R + 0.5T potential.