QM/MM Minimum Free Energy Path: Methodology and Application to Triosephosphate Isomerase

Abstract

Structural and energetic changes are two important characteristic properties of a chemical reaction process. In the condensed phase, studying these two properties is very challenging because of the great computational cost associated with the quantum mechanical calculations and phase space sampling. Although the combined quantum mechanics/molecular mechanics (QM/MM) approach significantly reduces the amount of the quantum mechanical calculations and facilitates the simulation of solution phase and enzyme catalyzed reactions, the required quantum mechanical calculations remain quite expensive and extensive sampling can be achieved routinely only with semiempirical quantum mechanical methods. QM/MM simulations with ab initio QM methods, therefore, are often restricted to narrow regions of the potential energy surface such as the reactant, product and transition state, or the minimum energy path. Such ab initio QM/MM calculations have previously been performed with the QM/MM-Free Energy (QM/MM-FE) method of Zhang et al.1 to generate the free energy profile along the reaction coordinate using free energy perturbation calculations at fixed structures of the QM subsystems. Results obtained with the QM/MM-FE method depend on the determination of the minimum energy reaction path, which is based on local conformations of the protein/solvent environment and can be difficult to obtain in practice. To overcome the difficulties associated with the QM/MM-FE method and to further enhance the sampling of the MM environment conformations, we develop here a new method to determine the QM/MM minimum free energy path (QM/MM-MFEP) for chemical reaction processes in solution and in enzymes. Within the QM/MM framework, we express the free energy of the system as a function of the QM conformation, thus leading to a simplified potential of mean force (PMF) description for the thermodynamics of the system. The free energy difference between two QM conformations is evaluated by the QM/MM free energy perturbation method. The free energy gradients with respect to the QM degrees of freedom are calculated from molecular dynamics simulations at given QM conformations. With the free energy and free energy gradients in hand, we further implement chain-of-conformation optimization algorithms in the search for the reaction path on the free energy surface without specifying a reaction coordinate. This method thus efficiently provides a unique minimum free energy path for solution and enzyme reactions, with structural and energetic properties being determined simultaneously. To further incorporate the dynamic contributions of the QM subsystem into the simulations, we develop the reaction path potential of Lu, et al.2 for the minimum free energy path. The combination of the methods developed here presents a comprehensive and accurate treatment for the simulation of reaction processes in solution and in enzymes with ab initio QM/MM methods. The method has been demonstrated on the first step of the reaction of the enzyme triosephosphate isomerase with good agreement with previous studies.

Figure 1

Active site and the reaction catalyzed by the enzyme triosephosphate isomerase (TIM).

Figure 2

Convergences of the free energy gradients of all QM atoms in the TIM system. The X-axis is the simulation time, while the Y-axis is the free energy gradient defined in Eq. 9.

Figure 3

Free energy profile from NEB optimization by the QM/MM-MFEP method for the proton transfer reaction in TIM. To make the results comparable to previous studies, the X-axis R is defined as the distance differences R_C–H – R_O–H, instead of the indices of the conformations used in the simulations.

Figure 4

Free energy profile from Ayala-Schlegel MEP optimizations by the QM/MM-MFEP method for the proton transfer reaction in TIM. The X-axis is defined in the same way as in Figure 2.

Figure 5

Comparison between the energies calculated for different MM environments with the same QM conformation. The X-axis is the energy calculated by direct QM/MM method in which the exact QM electron density is solved for each QM/MM conformation; the Y-axis is the energy calculated by the RPP method with polarizable QM charges. The RPP is constructed using the mean-field of all MM charges as the external electrostatic potential for the QM part. For better display, the energies on both axes have been shifted by the same amount.