Machine-Learning-Assisted Free Energy Simulation of Solution-Phase and Enzyme Reactions

Abstract

Despite recent advances in the development of machine learning potentials (MLPs) for biomolecular simulations, there has been limited effort on developing stable and accurate MLPs for enzymatic reactions. Here we report a protocol for performing machine-learning-assisted free energy simulation of solution-phase and enzyme reactions at the ab initio quantum-mechanical/molecular-mechanical (ai-QM/MM) level of accuracy. Within our protocol, the MLP is built to reproduce the ai-QM/MM energy and forces on both QM (reactive) and MM (solvent/enzyme) atoms. As an alternative strategy, a delta machine learning potential (ΔMLP) is trained to reproduce the differences between the ai-QM/MM and semiempirical (se) QM/MM energies and forces. To account for the effect of the condensed-phase environment in both MLP and ΔMLP, the DeePMD representation of a molecular system is extended to incorporate the external electrostatic potential and field on each QM atom. Using the Menshutkin and chorismate mutase reactions as examples, we show that the developed MLP and ΔMLP reproduce the ai-QM/MM energy and forces with errors that on average are less than 1.0 kcal/mol and 1.0 kcal mol^-1 Å^-1, respectively, for representative configurations along the reaction pathway. For both reactions, MLP/ΔMLP-based simulations yielded free energy profiles that differed by less than 1.0 kcal/mol from the reference ai-QM/MM results at only a fraction of the computational cost.

Figure 1:

Workflow for the training of MLP and its use to generate energy and forces for MD simulations. Q_i is the ESP charge on QM atom i, which is fitted using the electrostatic potential ϕ_B on inner MM atom positions.

Figure 2:

Schemes for (a) Menshutkin and (b) chorismate mutase reactions.

Figure 3:

Conservation of the total energy in 100ps NVE simulations of the chorismate mutase reaction using A) PM3*, B) MLP, and C) PM3*+ΔMLP models. In each figure, the line shown in orange indicates the drift of energy (see the values mentioned in the main text for each model).

Figure 4:

Accuracy of MLP (top), PM3*+ΔMLP (middle), PM3 and PM3* (bottom) energy, forces, electrostatic potential (ϕ) and electric field ($\mathcal{E}$) for the 2,000 testing configurations for aqueous Menshutkin reaction. In each figure, the reference values are obtained from the B3LYP/6–31G*/MM calculations. The root-mean-square error (RMSE) value is also shown for each method.

Figure 5:

Distribution of high-level (B3LYP/6–31G*) and low-level (PM3*, MLP, and PM3*+ΔMLP energy differences for configuration collected from B3LYP/MM MD trajectories (blue) or low-level MD trajectories (orange) of the aqueous Menshutkin reaction.

Figure 6:

(A) Sampled pathway and (B) potential of mean force for the aqueous Menshutkin reaction based on umbrella sampling using PM3*+ΔMLP and MLP potentials in comparison to PM3* and B3LYP/6–31G* results. The pathways are represented by the average bond lengths for each of the 80 windows in the umbrella sampling simulations (shown in Figure S2). The stars show the locations of the transition states on the pathways.

Figure 7:

Accuracy of MLP (top), PM3*+ΔMLP (middle), PM3 and PM3* (bottom) energy and forces for the 2,000 testing configurations for the chorismate mutase reaction.

Figure 8:

Distribution of high-level (B3LYP/6–31G*) and low-level (PM3*, MLP, and PM3*+ΔMLP energy differences for configuration collected from B3LYP/MM MD trajectories (blue) or low-level MD trajectories (orange) of chorismate mutase reaction.

Figure 9:

(A) Sampled pathway and (B) potential of mean force for the chorismate mutase reaction based on umbrella sampling using PM3*, PM3*+ΔMLP and MLP (the 2nd iteration) potentials in comparison to B3LYP/6–31G* results. The pathways are represented by the average bond lengths for each of the 80 windows in the umbrella sampling simulations (in Figure S3). The stars show the locations of the transition states on the pathways.