Minimum Free Energy Path Research Articles

Efficient quantum mechanical minimum free energy path calculation by combining path integral hybrid Monte Carlo and climbing image nudged elastic band methods, and its application to the addition reaction of hydrogen isocyanide to formaldehyde.

We propose an efficient algorithm for a minimum free energy path calculation based on the path integral hybrid Monte Carlo (PIHMC) method by combining the climbing image-nudged elastic band (CI-NEB) and the thermodynamic integration (TI) methods. Here, the CI-NEB and the TI methods are used to find a transition state along the reaction path and evaluate the free energy path, respectively. Our algorithm is applied to the Walden inversion reaction of the hydronium ions (H3O+). The numerical results show that the computational effort by our algorithm is significantly reduced compared to that of the previously proposed algorithm combining PIHMC without losing accuracy. We also demonstrate the importance of temperature and isotope effects on the addition reaction of hydrogen isocyanide to formaldehyde. In this reaction, the nuclear quantum effect causes the structural change at the TS and decreases the energy barrier.

Theoretical Investigation into Polymorphic Transformation between β-HMX and δ-HMX by Finite Temperature String

Polymorphic transformation is important in chemical industries, in particular, in those involving explosive molecular crystals. However, due to simulating challenges in the rare event method and collective variables, understanding the transformation mechanism of molecular crystals with a complex structure at the molecular level is poor. In this work, with the constructed order parameters (OPs) and K-means clustering algorithm, the potential of mean force (PMF) along the minimum free-energy path connecting β-HMX and δ-HMX was calculated by the finite temperature string method in the collective variables (SMCV), the free-energy profile and nucleation kinetics were obtained by Markovian milestoning with Voronoi tessellations, and the temperature effect on nucleation was also clarified. The barriers of transformation were affected by the finite-size effects. The configuration with the lower potential barrier in the PMF corresponded to the critical nucleus. The time and free-energy barrier of the polymorphic transformation were reduced as the temperature increased, which was explained by the pre-exponential factor and nucleation rate. Thus, the polymorphic transformation of HMX could be controlled by the temperatures, as is consistent with previous experimental results. Finally, the HMX polymorph dependency of the impact sensitivity was discussed. This work provides an effective way to reveal the polymorphic transformation of the molecular crystal with a cyclic molecular structure, and further to prepare the desired explosive by controlling the transformation temperature.

Efficient acceleration of the convergence of the minimum free energy path via a path-planning generated initial guess.

We demonstrate that combining a shifted clustering algorithm with a fast-marching-based algorithm can generate accurate approximations of the minimum energy path (MEP) given a free energy landscape (FEL). Using this approximation as the initial guess for the MEP, followed by further refinement with the string method (referred to as the fast marching tree (FMT)-string combined approach), significantly reduces the number of iterations required for MEP convergence. This approach saves substantial time compared to using linear interpolation (LI) for the initial guess. Our method offers a viable solution for obtaining an effective initial guess of the MEP when an approximate or converged FEL is available. This work highlights the potential of applying FMT-based approaches to extract the MEP in chemical reactions.

DFT/MM Simulations for Cycloreversion Reaction of Cyclobutane Pyrimidine Dimer with Deprotonated and Protonated E283.

DNA photolyase targets the primary ultraviolet (UV)-induced DNA lesion─cyclobutane pyrimidine dimer (CPD), attaches to it, and catalyzes its dissociation. The catalytic mechanism of DNA photolyase and the role of the conserved residue E283 remain subjects of debate. This study employs two-dimensional potential energy surface maps and minimum free energy paths calculated at the ωB97XD/6-31G/MM level to elucidate these mechanisms. Results suggest that the catalytic process follows a sequential, stepwise reaction in which the C5-C5 and C6-C6 bonds are cleaved in order, facilitated by a protonated E283. Activation free energies for these cleavages are calculated at 4.4 and 4.2 kcal·mol-1, respectively. Protonation of E283 reduces electrostatic repulsion with CPD and forms dual hydrogen bonds with it and provides better solvation, stabilizing the CPD radical anion, particularly during intermediate state. This stabilization renders the initial splitting step exergonic, slows reverse reactions of the C5-C5 bond cleavage and electron transfer, and ensures a high quantum yield. Furthermore, the protonation state of E283 significantly affects the type of bond cleavage. Other residues in the active site were also investigated for their roles in the mechanism.

Surface-Accelerated String Method for Locating Minimum Free Energy Paths.

We present a surface-accelerated string method (SASM) to efficiently optimize low-dimensional reaction pathways from the sampling performed with expensive quantum mechanical/molecular mechanical (QM/MM) Hamiltonians. The SASM accelerates the convergence of the path using the aggregate sampling obtained from the current and previous string iterations, whereas approaches like the string method in collective variables (SMCV) or the modified string method in collective variables (MSMCV) update the path only from the sampling obtained from the current iteration. Furthermore, the SASM decouples the number of images used to perform sampling from the number of synthetic images used to represent the path. The path is optimized on the current best estimate of the free energy surface obtained from all available sampling, and the proposed set of new simulations is not restricted to being located along the optimized path. Instead, the umbrella potential placement is chosen to extend the range of the free energy surface and improve the quality of the free energy estimates near the path. In this manner, the SASM is shown to improve the exploration for a minimum free energy pathway in regions where the free energy surface is relatively flat. Furthermore, it improves the quality of the free energy profile when the string is discretized with too few images. We compare the SASM, SMCV, and MSMCV using 3 QM/MM applications: a ribozyme methyltransferase reaction using 2 reaction coordinates, the 2'-O-transphosphorylation reaction of Hammerhead ribozyme using 3 reaction coordinates, and a tautomeric reaction in B-DNA using 5 reaction coordinates. We show that SASM converges the paths using roughly 3 times less sampling than the SMCV and MSMCV methods. All three algorithms have been implemented in the FE-ToolKit package made freely available.

Fine-tuning activation specificity of G-protein-coupled receptors via automated path searching

Physics-based simulation methods can grant atomistic insights into the molecular origin of the function of biomolecules. However, the potential of such approaches has been hindered by their low efficiency, including in the design of selective agonists where simulations of myriad protein-ligand combinations are necessary. Here, we describe an automated input-free path searching protocol that offers (within 14 d using Graphics Processing Unit servers) a minimum free energy path (MFEP) defined in high-dimension configurational space for activating sphingosine-1-phosphate receptors (S1PRs) by arbitrary ligands. The free energy distributions along the MFEP for four distinct ligands and three S1PRs reached a remarkable agreement with Bioluminescence Resonance Energy Transfer (BRET) measurements of G-protein dissociation. In particular, the revealed transition state structures pointed out toward two S1PR3 residues F263/I284, that dictate the preference of existing agonists CBP307 and BAF312 on S1PR1/5. Swapping these residues between S1PR1 and S1PR3 reversed their response to the two agonists in BRET assays. These results inspired us to design improved agonists with both strong polar head and bulky hydrophobic tail for higher selectivity on S1PR1. Through merely three in silico iterations, our tool predicted a unique compound scaffold. BRET assays confirmed that both chiral forms activate S1PR1 at nanomolar concentration, 1 to 2 orders of magnitude less than those for S1PR3/5. Collectively, these results signify the promise of our approach in fine agonist design for G-protein-coupled receptors.

Improved Sampling of Adaptive Path Collective Variables by Stabilized Extended-System Dynamics.

Because of the complicated multistep nature of many biocatalytic reactions, an a priori definition of reaction coordinates is difficult. Therefore, we apply enhanced sampling algorithms along with adaptive path collective variables (PCVs), which converge to the minimum free energy path (MFEP) during the simulation. We show how PCVs can be combined with the highly efficient well-tempered metadynamics extended-system adaptive biasing force (WTM-eABF) hybrid sampling algorithm, offering dramatically increased sampling efficiency due to its fast adaptation to path updates. For this purpose, we address discontinuities of PCVs that can arise due to path shortcutting or path updates with a novel stabilization algorithm for extended-system methods. In addition, we show how the convergence of simulations can be further accelerated by utilizing the multistate Bennett's acceptance ratio (MBAR) estimator. These methods are applied to the first step of the enzymatic reaction mechanism of pseudouridine synthases, where the ability of path WTM-eABF to efficiently explore intricate molecular transitions is demonstrated.

Exhaustively Exploring the Prevalent Interaction Pathways of Ligands Targeting the Ligand-Binding Pocket of Farnesoid X Receptor via Combined Enhanced Sampling.

It is well-known that the potency of a drug is heavily associated with its kinetic and thermodynamic properties with the target. Nuclear receptors (NRs), as an important target family, play important roles in regulating a variety of physiological processes in vivo. However, it is hard to understand the drug-NR interaction process because of the closed structure of the ligand-binding domain (LBD) of the NR proteins, which apparently hinders the rational design of drugs with controllable kinetic properties. Therefore, understanding the underlying mechanism of the ligand-NR interaction process seems necessary to help NR drug design. However, it is usually difficult for experimental approaches to interpret the kinetic process of drug-target interactions. Therefore, in silico methods were utilized to explore the optimal binding/dissociation pathways of the NR ligands. Specifically, farnesoid X receptor (FXR) is considered here as the target system since it has been an important target for the treatment of bile acid metabolism-associated diseases, and a series of structures cocrystallized with diverse scaffold ligands were resolved. By using random acceleration molecular dynamics (RAMD) simulation and umbrella sampling (US), 5 main dissociation pathways (pathways I-V) were identified in 11 representative FXR ligands, with most of them (9/11) preferring to go through Pathway III and the remaining two favoring escaping from Pathway I and IV. Furthermore, key residues functioning in the three main dissociation pathways were revealed by the kinetic residue energy analysis (KREA) based on the US trajectories, which may serve as road-marker residues for rapid identification of the (un)binding pathways of FXR ligands. Moreover, the preferred pathways explored by RAMD simulations are in good agreement with the minimum free energy path identified by the US simulations with the Pearson R = 0.76 between the predicted binding affinity and the experimental data, suggesting that RAMD is suitable for applying in large-scale (un)binding-pathway exploration in the case of ligands with obscure binding tunnels to the target.

Bridging semiempirical and abinitio QM/MM potentials by Gaussian process regression and its sparse variants for free energy simulation.

Free energy simulations that employ combined quantum mechanical and molecular mechanical (QM/MM) potentials at abinitio QM (AI) levels are computationally highly demanding. Here, we present a machine-learning-facilitated approach for obtaining AI/MM-quality free energy profiles at the cost of efficient semiempirical QM/MM (SE/MM) methods. Specifically, we use Gaussian process regression (GPR) to learn the potential energy corrections needed for an SE/MM level to match an AI/MM target along the minimum free energy path (MFEP). Force modification using gradients of the GPR potential allows us to improve configurational sampling and update the MFEP. To adaptively train our model, we further employ the sparse variational GP (SVGP) and streaming sparse GPR (SSGPR) methods, which efficiently incorporate previous sample information without significantly increasing the training data size. We applied the QM-(SS)GPR/MM method to the solution-phase SN2 Menshutkin reaction, NH3+CH3Cl→CH3NH3++Cl-, using AM1/MM and B3LYP/6-31+G(d,p)/MM as the base and target levels, respectively. For 4000 configurations sampled along the MFEP, the iteratively optimized AM1-SSGPR-4/MM model reduces the energy error in AM1/MM from 18.2 to 4.4 kcal/mol. Although not explicitly fitting forces, our method also reduces the key internal force errors from 25.5 to 11.1 kcal/mol/Å and from 30.2 to 10.3 kcal/mol/Å for the N-C and C-Cl bonds, respectively. Compared to the uncorrected simulations, the AM1-SSGPR-4/MM method lowers the predicted free energy barrier from 28.7 to 11.7 kcal/mol and decreases the reaction free energy from -12.4 to -41.9 kcal/mol, bringing these results into closer agreement with their AI/MM and experimental benchmarks.

Insights into generalization of the rate-limiting steps of the dehalogenation by LinB and DhaA: A computational approach

LinB and DhaA are well-known haloalkane dehalogenases (HLDs) capable of converting a plethora of halogenated alkanes, also those considered persistent pollutants. The dehalogenation reaction that these two enzymes catalyze has been studied to determine its rate-limiting step (rls) for the last two decades now. As a result, it has been determined that HLDs can show different rate-limiting steps for individual substrates, and at this point we do not have a basis for any generalization in this matter. Therefore, in this work we aimed at gaining insights into the enzymatic dehalogenation of selected dibromo- and bromochloro-ethanes and propanes by LinB and DhaA using computational approach to determine whether defined structural similarities of the substrates result in a unified mechanism and the same rls. By predicting halogen binding isotope effects (BIEs) as well as computing interaction energy for each HLD-ligand complex the nature of the protein-ligand interactions has been characterized. Furthermore, C and Br kinetic isotope effects (KIEs) as well as the minimum free energy paths (MFEPs) were computed to investigate the chemical reaction for the selected systems. Accuracy of the approach and robustness of the computational predictions were validated by measuring KIEs on the selected reactions. Overall results strongly indicate that any generalization with respect to the enzymatic process involving various ligands in the case of DhaA is impossible, even if the considered ligands are structurally very similar as those analyzed in the present study. Moreover, even small structural differences such as changing of one of the (non-leaving) halogen substituents may lead to significant changes in the enzymatic process and result in a different rls in the case of LinB. It has also been demonstrated that KIEs themselves cannot be used as rls indicators in the reactions catalyzed by the studied HLDs.

Estimation of frequency factors for the calculation of kinetic isotope effects from classical and path integral free energy simulations.

We use the modified Bigeleisen-Mayer equation to compute kinetic isotope effect values for non-enzymatic phosphoryl transfer reactions from classical and path integral molecular dynamics umbrella sampling. The modified form of the Bigeleisen-Mayer equation consists of a ratio of imaginary mode vibrational frequencies and a contribution arising from the isotopic substitution's effect on the activation free energy, which can be computed from path integral simulation. In the present study, we describe a practical method for estimating the frequency ratio correction directly from umbrella sampling in a manner that does not require normal mode analysis of many geometry optimized structures. Instead, the method relates the frequency ratio to the change in the mass weighted coordinate representation of the minimum free energy path at the transition state induced by isotopic substitution. The method is applied to the calculation of 16/18O and 32/34S primary kinetic isotope effect values for six non-enzymatic phosphoryl transfer reactions. We demonstrate that the results are consistent with the analysis of geometry optimized transition state ensembles using the traditional Bigeleisen-Mayer equation. The method thus presents a new practical tool to enable facile calculation of kinetic isotope effect values for complex chemical reactions in the condensed phase.

Automated Path Searching Reveals the Mechanism of Hydrolysis Enhancement by T4 Lysozyme Mutants

Bacteriophage T4 lysozyme (T4L) is a glycosidase that is widely applied as a natural antimicrobial agent in the food industry. Due to its wide applications and small size, T4L has been regarded as a model system for understanding protein dynamics and for large-scale protein engineering. Through structural insights from the single conformation of T4L, a series of mutations (L99A,G113A,R119P) have been introduced, which have successfully raised the fractional population of its only hydrolysis-competent excited state to 96%. However, the actual impact of these substitutions on its dynamics remains unclear, largely due to the lack of highly efficient sampling algorithms. Here, using our recently developed travelling-salesman-based automated path searching (TAPS), we located the minimum-free-energy path (MFEP) for the transition of three T4L mutants from their ground states to their excited states. All three mutants share a three-step transition: the flipping of F114, the rearrangement of α0/α1 helices, and final refinement. Remarkably, the MFEP revealed that the effects of the mutations are drastically beyond the expectations of their original design: (a) the G113A substitution not only enhances helicity but also fills the hydrophobic Cavity I and reduces the free energy barrier for flipping F114; (b) R119P barely changes the stability of the ground state but stabilizes the excited state through rarely reported polar contacts S117OG:N132ND2, E11OE1:R145NH1, and E11OE2:Q105NE2; (c) the residue W138 flips into Cavity I and further stabilizes the excited state for the triple mutant L99A,G113A,R119P. These novel insights that were unexpected in the original mutant design indicated the necessity of incorporating path searching into the workflow of rational protein engineering.

Multiscale Modeling of Grain-Boundary Motion in Cylinder-Forming Block Copolymers

Using the combination of a soft, coarse-grained, particle-based model, a free-energy functional that depends on the local composition, and a lattice model of local, metastable states, we study the structure and motion of a grain boundary between two orthogonal grains of cylindrical domains in asymmetric block copolymers. The particle-based model provides direct insights into the elementary class of transitions of the self-assembled morphology in the course of grain-boundary translation. These processes are correlated in space and time. We identify a minimal set of transitions, whose free-energy changes and barriers are obtained by describing the system by a free-energy functional of the local composition and calculating the minimum free-energy path (MFEP). The spatiotemporal correlation arises from the dependence of the free-energy characteristics on the local environment. We use this information to parametrize a lattice model of the correlated processes in the course of grain-boundary motion. This allows us to investigate the grain-boundary motion by kinetic Monte Carlo (kMC) simulation and determine its free-energy landscape. Grain-boundary motion proceeds by nucleating a two-dimensional, anisotropic cluster inside the plane of the grain boundary.

Identifying a Feasible Transition Pathway between Two Conformational States for a Protein.

Proteins usually need to transit between different conformational states to fulfill their biological functions. In the mechanistic study of such transition processes by molecular dynamics simulations, identification of the minimum free energy path (MFEP) can substantially reduce the sampling space, thus enabling rigorous thermodynamic evaluation of the process. Conventionally, the MFEP is derived by iterative local optimization from an initial path, which is typically generated by simple brute force techniques like the targeted molecular dynamics (tMD). Therefore, the quality of the initial path determines the successfulness of MFEP estimation. In this work, we propose a method to improve derivation of the initial path. Through iterative relaxation-biasing simulations in a bidirectional manner, this method can construct a feasible transition pathway connecting two known states for a protein. Evaluation on small, fast-folding proteins against long equilibrium trajectories supports the good sampling efficiency of our method. When applied to larger proteins including the catalytic domain of human c-Src kinase as well as the converter domain of myosin VI, the paths generated by our method deviate significantly from those computed with the generic tMD approach. More importantly, free energy profiles and intermediate states obtained from our paths exhibit remarkable improvements over those from tMD paths with respect to both physical rationality and consistency with a priori knowledge.

Cover Image

Protein ScienceVolume 31, Issue 5 Cover ImageFree Access Cover Image First published: 04 April 2022 https://doi.org/10.1002/pro.4302AboutPDF ToolsRequest permissionExport citationAdd to favoritesTrack citation ShareShare Give accessShare full text accessShare full-text accessPlease review our Terms and Conditions of Use and check box below to share full-text version of article.I have read and accept the Wiley Online Library Terms and Conditions of UseShareable LinkUse the link below to share a full-text version of this article with your friends and colleagues. Learn more.Copy URL Share a linkShare onFacebookTwitterLinked InRedditWechat Abstract Minimum free-energy path for the early stages of domain dissociation in the orange carotenoid protein (OCP). See Sharawy M, Pigni NB, May ER, Gascón JA. A favorable path to domain separation in the orange carotenoid protein. Protein Science. 2022;31:850–63. https://doi.org/10.1002/pro.4273 Volume31, Issue5May 2022 RelatedInformation

An investigation of the kinetics and thermodynamics of NaCl nucleation through composite clusters.

Having a good understanding of nucleation is critical for the control of many important processes, such as polymorph selection during crystallization. However, a complete picture of the molecular-level mechanisms of nucleation remains elusive. In this work, we take an in-depth look at the NaCl homogeneous nucleation mechanism through thermodynamics. Distinguished from the classical nucleation theory, we calculate the free energy of nucleation as a function of two nucleus size coordinates: crystalline and amorphous cluster sizes. The free energy surface reveals a thermodynamic preference for a nonclassical mechanism of nucleation through a composite cluster, where the crystalline nucleus is surrounded by an amorphous layer. The thickness of the amorphous layer increases with an increase in supersaturation. The computed free energy landscape agrees well with the composite cluster-free energy model, through which phase specific thermodynamic properties are evaluated. As the supersaturation increases, there is a change in stability of the amorphous phase relative to the solution phase, resulting in a change from one-step to two-step mechanism, seen clearly from the free energy profile along the minimum free energy path crossing the transition curve. By obtaining phase-specific diffusion coefficients, we construct the full mesoscopic model and present a clear roadmap for NaCl nucleation.

Molecular Mechanism of Dual Intercalation in Sac7d-DNA Complexation.

Sac7d belongs to the hyperthermophilc chromosomal protein family, and it is very stable with regard to heat and acidic environments. Unlike many other DNA-protein complexes, the present one is a nonspecific complexation where two amino acids (AA), VAL26 and MET29, are found to intercalate into the same base pair of DNA. Here, we have carried out multiple short molecular dynamic simulations to calculate the distribution of nonspecific protein-DNA aggregates to find the most probable state, which was subsequently used to construct the free energy landscape of protein intercalation into DNA. Analysis of trajectories along the minimum free energy path revealed mechanistic details such as rotation of the protein, simultaneous intercalation of two amino acids, and bending/kinking of the DNA. Moreover, the results indicate a strong interdependency between the intercalating amino acids such that the deintercalation of one AA leads to a spontaneous deintercalation of the other.

Self-Assembly of Model Triblock Janus Colloidal Particles in Two Dimensions.

A simplified two-dimensional effective-solvent model of triblock Janus particles, consisting of three interaction sites in a linear configuration, a core particle, and two particles modeling the attractive patches at the poles, is developed to study the mechanism of nucleation and self-assembly in triblock Janus particles. The potential energy parameters are tuned against phase transition temperatures and free energy barriers to the nucleation of crystalline phases, calculated from our previous detailed model of Janus particles. Vapor-liquid equilibria and critical temperatures are calculated by grand-canonical molecular dynamics simulations for particles of different patch widths. With metadynamics, phase equilibria, mechanism of nucleation, and free energy barriers to nucleation are investigated. The minimum free energy path to nucleation indicates two steps. The first step, with a higher free energy increase, consists of the densification of the fluid into a disordered cluster. In the second step, of a lower free energy barrier, the inner particles of the disordered cluster reorient to form a crystalline nucleus. This two-step mechanism of nucleation of a kagome lattice is in complete agreement with the experiment and with our previous simulations using a detailed model of Janus particles. Large systems at a slight supersaturation generate multiple crystalline domains, which are misaligned at the grain boundaries. In complete agreement with the experiment and with previous simulation results, we observe a two-step mechanism for crystal growth: melting of the smaller (less stable) crystallites to a fluid followed by recrystallization at the surface of neighboring bigger (more stable) crystallites. A comparison of the present softer modeling of a Janus particle with harder models in the literature for self-assembly of Janus particles indicates that softer potentials stabilize open lattices (e.g., kagome) more than dense lattices (e.g., hexagonal). Also, experimental locations of phase transition points and barrier heights to nucleation are better reproduced by the present model than by the existing simple models.

Semiflexible oligomers crystallize via a cooperative phase transition.

Semicrystalline polymers are ubiquitous, yet despite their fundamental and industrial importance, the theory of homogeneous nucleation from a melt remains a subject of debate. A key component of the controversy is that polymer crystallization is a non-equilibrium process, making it difficult to distinguish between effects that are purely kinetic and those that arise from the underlying thermodynamics. Due to computational cost constraints, simulations of polymer crystallization typically employ non-equilibrium molecular dynamics techniques with large degrees of undercooling that further exacerbate the coupling between thermodynamics and kinetics. In a departure from this approach, in this study, we isolate the near-equilibrium nucleation behavior of a simple model of a melt of short, semiflexible oligomers. We employ several Monte Carlo methods and compute a phase diagram in the temperature-density plane along with two-dimensional free energy landscapes (FELs) that characterize the nucleation behavior. The phase diagram shows the existence of ordered nematic and crystalline phases in addition to the disordered melt phase. The minimum free energy path in the FEL for the melt-crystal transition shows a cooperative transition, where nematic order and monomer positional order move in tandem as the system crystallizes. This near-equilibrium phase transition mechanism broadly agrees with recent evidence that polymer stiffness plays an important role in crystallization but differs in the specifics of the mechanism from several recent theories. We conclude that the computation of multidimensional FELs for models that are larger and more fine-grained will be important for evaluating and refining theories of homogeneous nucleation for polymer crystallization.

Doubly Polarized QM/MM with Machine Learning Chaperone Polarizability.

A major shortcoming of semiempirical (SE) molecular orbital methods is their severe underestimation of molecular polarizability compared with experimental and ab initio (AI) benchmark data. In a combined quantum mechanical and molecular mechanical (QM/MM) treatment of solution-phase reactions, solute described by SE methods therefore tends to generate inadequate electronic polarization response to solvent electric fields, which often leads to large errors in free energy profiles. To address this problem, here we present a hybrid framework that improves the response property of SE/MM methods through high-level molecular-polarizability fitting. Specifically, we place on QM atoms a set of corrective polarizabilities (referred to as chaperone polarizabilities), whose magnitudes are determined from machine learning (ML) to reproduce the condensed-phase AI molecular polarizability along the minimum free energy path. These chaperone polarizabilities are then used in a machinery similar to a polarizable force field calculation to compensate for the missing polarization energy in the conventional SE/MM simulations. Because QM atoms in this treatment host SE wave functions as well as classical polarizabilities, both polarized by MM electric fields, we name this method doubly polarized QM/MM (dp-QM/MM). We demonstrate the new method on the free energy simulations of the Menshutkin reaction in water. Using AM1/MM as a base method, we show that ML chaperones greatly reduce the error in the solute molecular polarizability from 6.78 to 0.03 Å3 with respect to the density functional theory benchmark. The chaperone correction leads to ∼10 kcal/mol of additional polarization energy in the product region, bringing the simulated free energy profiles to closer agreement with the experimental results. Furthermore, the solute-solvent radial distribution functions show that the chaperone polarizabilities modify the free energy profiles through enhanced solvation corrections when the system evolves from the charge-neutral reactant state to the charge-separated transition and product states. These results suggest that the dp-QM/MM method, enabled by ML chaperone polarizabilities, provides a very physical remedy for the underpolarization problem in SE/MM-based free energy simulations.