Zhifeng F Jing

Chengwen Liu

Sara Y Cheng
email: sara.y.cheng@utexas.edu

Rui Qi

Brandon D Walker

Jean-Philip Piquemal

Pengyu Ren
email: pren@mail.utexas.edu

Polarizable force fields for biomolecular simulations: Recent advances and applications

Keywords: molecular dynamics simulation, polarizable force field, protein, nucleic acid, QM/MM

come   L'archive ouverte pluridisciplinaire

Introduction

Molecular dynamics (MD) simulation is indispensable tools for investigating physical properties of proteins, nucleic acids and designing new molecules and materials. [START_REF] Goh | Computational Methodologies for Real-Space Structural Refinement of Large Macromolecular Complexes[END_REF][START_REF] Mortier | The impact of molecular dynamics on drug design: applications for the characterization of ligandmacromolecule complexes[END_REF]128) Due to recent advances in computing hardware and improved simulation methods, the time and length scales of molecular dynamics simulations have been greatly extended. Noticeably, by combining GPU computing [START_REF] Carsten | Best bang for your buck: GPU nodes for GROMACS biomolecular simulations[END_REF][START_REF] Eastman | OpenMM 7: Rapid Development of High Performance Algorithms for Molecular Dynamics[END_REF]103) and enhanced sampling methods, [START_REF] Elber | Perspective: Computer simulations of long time dynamics[END_REF][START_REF] Faradjian | Computing time scales from reaction coordinates by milestoning[END_REF][START_REF] Husic | Markov State Models: From an Art to a Science[END_REF] molecular simulations are approaching the time scale of milliseconds and seconds, enabling the study of macromolecular interactions and folding with high fidelity. These advances not only lead to more reliable interpretation and predictions by computer simulations but also crucial for examining and improving the underlying physical models and simulation methods.

There has been much effort devoted to improving the potential-energy functions or force field (FF) used in MD simulations. It is believed in biology that amino acid sequences determine the structure, which then determines the function. The potential energy surface defines the physical driving forces underlying biomolecular structure and interactions. Force fields usually consist of several empirical energy terms including short-ranged bonded interactions and non-bonded interactions such as repulsion, dispersion and electrostatics. Electrostatics is both important and computationally expensive due to its long-range nature. To facilitate simulations of biomolecules with modest computational power, traditional force fields (FFs) use fixed point charge placed at atomic centers to represent the electrostatic interactions. The limitations of the fixed pointcharge force fields have been well recognized. [START_REF] Cisneros | Classical Electrostatics for Biomolecular Simulations[END_REF][START_REF] Lemkul | An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications[END_REF][START_REF] Bc | Polarizable force fields for molecular dynamics simulations of biomolecules[END_REF][START_REF] Marshall | Limiting assumptions in molecular modeling: electrostatics[END_REF], 89, 100) One significant approximation in traditional force fields is the omission of polarization, i.e. the response of the charge distribution to environment. This is problematic when applying the same set of charge parameters to different environments, such as aqueous solution, protein cavity, cell membrane and heterogeneous interfaces, where the charge distribution should change accordingly. Another approximation is the atom-centered point-charge model, whereas the realistic charge distribution should be smooth and anisotropic. To capture anisotropic features such as σ-holes, lone pairs and π-bonding, it is necessary to adopt higher-order multipolar electrostatics models [START_REF] Gresh | Anisotropic, Polarizable Molecular Mechanics Studies of Inter-and Intramolecular Interactions and Ligand-Macromolecule Complexes. A Bottom-Up Strategy[END_REF][START_REF] Jakobsen | Systematic Improvement of Potential-Derived Atomic Multipoles and Redundancy of the Electrostatic Parameter Space[END_REF][START_REF] Kolář | Computer Modeling of Halogen Bonds and Other σ-Hole Interactions[END_REF][START_REF] Kramer | Charge Anisotropy: Where Atomic Multipoles Matter Most[END_REF]122) and/or adding off-center sites. [START_REF] Gresh | Anisotropic, Polarizable Molecular Mechanics Studies of Inter-and Intramolecular Interactions and Ligand-Macromolecule Complexes. A Bottom-Up Strategy[END_REF][START_REF] Lemkul | An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications[END_REF] The effect of having atomic multipoles beyond fixed charges is of the same magnitude as the effect of polarization, suggesting that both should be included in force field development. [START_REF] Jakobsen | Systematic Improvement of Potential-Derived Atomic Multipoles and Redundancy of the Electrostatic Parameter Space[END_REF] Previously a number of reviews on polarizable force fields have been published. [START_REF] Cisneros | Classical Electrostatics for Biomolecular Simulations[END_REF][START_REF] Halgren | Polarizable force fields[END_REF][START_REF] Lemkul | An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications[END_REF][START_REF] Bc | Polarizable force fields for molecular dynamics simulations of biomolecules[END_REF][START_REF] Piotr | Polarization effects in molecular mechanical force fields[END_REF][START_REF] Ponder | Current Status of the AMOEBA Polarizable Force Field[END_REF]100,135) In this review, we will provide an update on recent progress in advanced electrostatic modeling, simulation algorithms, development and applications of accurate and efficient polarizable force fields for biomolecular simulations.

Electrostatic Models and Force Field Parameterization

Permanent Electrostatics Electrostatic interactions are essential in biomolecular systems. [START_REF] Cisneros | Classical Electrostatics for Biomolecular Simulations[END_REF]100) Atomic partial charges or multipoles (monopole, dipole and higher moments) are typically used to represent the "permanent" charge distributions in molecules and have significant impact on simulated structures and dynamics. In point-charge force fields, the partial charges are usually determined by fitting to quantum mechanical (QM) electrostatic potential (ESP) or interaction energy with water while implicitly accounting for polarization. [START_REF] Cerutti | Derivation of Fixed Partial Charges for Amino Acids Accommodating a Specific Water Model and Implicit Polarization[END_REF][START_REF] Debiec | Further along the Road Less Traveled: AMBER ff15ipq, an Original Protein Force Field Built on a Self-Consistent Physical Model[END_REF][START_REF] Mackerell | All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins[END_REF]119,131) For atomic multipoles, the parameters are typically derived from ab initio calculations using procedures such as Distributed Multiple Analysis (DMA), Atoms-in-Molecules (AIM) and Iterative Stockholder Analysis. [START_REF] Misquitta | Distributed Multipoles from a Robust Basis-Space Implementation of the Iterated Stockholder Atoms Procedure[END_REF]117) It is recognized that multipole parameters are redundant and subject to overfitting. [START_REF] Jakobsen | Searching the Force Field Electrostatic Multipole Parameter Space[END_REF] Jensen and coworkers showed by tensor decomposition that the number of multipole parameters that can be effectively determined from the electrostatic potential for peptide models is less than twice the number of atoms. [START_REF] Jakobsen | Searching the Force Field Electrostatic Multipole Parameter Space[END_REF] Accordingly, the number of multipole parameters can be significantly reduced without affecting the accuracy. In contrast to the standard charge fitting schemes, fitting points were placed on a single isodensity surface. [START_REF] Jakobsen | Systematic Improvement of Potential-Derived Atomic Multipoles and Redundancy of the Electrostatic Parameter Space[END_REF] In a similar fashion, Meuwly and coworkers developed a minimal distributed charge model (MDCM) based on offcentered point charges.(122) MDCM is capable of approximating the reference ab initio ESP with an accuracy as good as electric multipoles (EMPs).

Electronic Polarization

A significant advancement in modeling biomolecular electrostatics over the past decade is the explicit treatment of the polarization effect, allowing electrostatics to respond to chemical environments. Classical polarization models can be classified into two categories, one characterizing the charge redistribution within each atom, by either induced dipole [START_REF] Ponder | Current Status of the AMOEBA Polarizable Force Field[END_REF] or Drude oscillator [START_REF] Lemkul | An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications[END_REF] (also called charge-on-spring) and the other based on charge flow between atoms such as the fluctuation charge model (also known as charge equilibration, or chemical potential equilibration).(104)

In the induced dipole model, the induced dipole moment 𝜇 ind at each polarizable site, an atom or its lone-pair site, is proportional to the total electric field 𝐸 , 𝜇 ind = 𝛼𝐸 , where 𝛼 is the polarizability. [START_REF] Ponder | Current Status of the AMOEBA Polarizable Force Field[END_REF] The induced dipole also creates an electric field, and mutual polarization between induced dipoles exists. Therefore, polarization is non-additive and often solved iteratively via self-consistent field (SCF). The SCF typically requires several iterations for tight convergence and is thus computational expensive. A few new algorithms have been devised in the past few years to significantly accelerate the calculation of induced dipoles (see Polarization algorithms section below). To avoid the artificially strong interaction between point charges/dipoles at short range, the Thole model( 93) is often employed to screen the electrostatics interactions according to a smeared charge distribution.

In the Drude oscillator model,(53) a Drude particle carrying part of the atomic charge is attached to the core atom via a harmonic spring. The displacement of the Drude particle in response to an electric field will create a dipole moment. The positions of the Drude particles also need to be solved iteratively via SCF to ensure that the Drude particles are at ground state, although extended Lagrangian method has been employed in MD simulations to approximate the exact solution. Similar to the induced dipole model, the electrostatic interaction at short range is damped. [START_REF] Lemkul | An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications[END_REF] The Drude oscillator model can be considered as a finite-difference approximation to the induced dipole model. Recently, Huang et al. numerically established the equivalence of Drude oscillator and induced dipole models. [START_REF] Huang | Mapping the Drude polarizable force field onto a multipole and induced dipole model[END_REF] Due to the similarity between the two models, they can share many advanced polarization solvers.

The fluctuating charge (FQ) model is based on the electronegativity equalization principle. The atomic charges are redistributed to equalize the electronegativity/chemical potential at each site, which depends on the atomic electronegativity, hardness and the external electrostatic potential.(10, 104) FQ has been used to develop force fields for proteins(104) as well as inorganic materials. [START_REF] Rappe | UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations[END_REF]110) Theoretically both induced dipole and fluctuating charge can be included in the model.(116) However, there is no clear boundary between polarization and charge fluctuation. Mei et al. found that the contribution of fluctuating charge to polarization obtained from nine different population analysis schemes varies from 59.9% to 96.2%. [START_REF] Mei | Numerical Study on the Partitioning of the Molecular Polarizability into Fluctuating Charge and Induced Atomic Dipole Contributions[END_REF] In practice, either induced dipole or fluctuating charge model can reasonably reproduce molecular polarizability. [START_REF] Ponder | Current Status of the AMOEBA Polarizable Force Field[END_REF]104) Besides the explicit polarization models, methods for effective polarization in the framework of fixed-charge force fields have also been developed. Leontyev and Stuchebrukhov proved that if simulations involve only structurally similar configurations, inherently polarizable molecular systems can be described by equivalent non-polarizable fixed-charge models. [START_REF] Leontyev | Polarizable molecular interactions in condensed phase and their equivalent nonpolarizable models[END_REF] The fixedcharge models can be implemented by scaling the partial charges by a constant, which is equivalent to a uniform dielectric constant. [START_REF] Leontyev | Polarizable molecular interactions in condensed phase and their equivalent nonpolarizable models[END_REF] It is recognized that simple scaling is not a replacement of a well-built real polarizable force field. [START_REF] Leontyev | Polarizable molecular interactions in condensed phase and their equivalent nonpolarizable models[END_REF] Better and more rigorous ways to account for solvent polarization and reference self-energy in fixed-atomic charges have also been published. [START_REF] Cerutti | Derivation of Fixed Partial Charges for Amino Acids Accommodating a Specific Water Model and Implicit Polarization[END_REF][START_REF] Debiec | Further along the Road Less Traveled: AMBER ff15ipq, an Original Protein Force Field Built on a Self-Consistent Physical Model[END_REF]119) Recent studies have revealed the inadequacies in current polarizable force fields. [START_REF] Demerdash | Assessing many-body contributions to intermolecular interactions of the AMOEBA force field using energy decomposition analysis of electronic structure calculations[END_REF][START_REF] Mao | Assessing Ion-Water Interactions in the AMOEBA Force Field Using Energy Decomposition Analysis of Electronic Structure Calculations[END_REF] By studying water-water, water-ion, and water-trimers, the Thole damping (or the specific parameters used) was shown to produce incorrect distance dependence for polarization when compared with energy decomposition analysis (EDA).

Atomic Polarizability

The atomic polarizabilities used in polarizable force fields are empirically chosen to reproduce the molecular polarizabilities of model molecules. Verstraelen et al. developed the Atom-Condensed Kohn-Sham DFT model approximated to Second order (ACKS2) to compute the polarization parameters directly computed as expectation values of an electronic wave function.(125) ACKS2 includes an electronic kinetic energy term to overcome the limitation of Electronegativity Equalization Method (EEM) (and other fluctuation charge models) that the distance of charge transfer for dielectric systems needs to be restrained.(125) Wang et al. developed a method to calculate the atomic polarizabilities by fitting to the electrostatic potentials (ESPs) under external electric field obtained from quantum mechanical (QM) calculations within the linear response theory.(129) The molecular polarizability obtained from this method is comparable to those from directly fitting to molecular polarizabilities. For polar molecules such as water, polarization is highly anisotropic. To account for the anisotropy of polarizability, atomic polarizability tensors are used.(117, 123) The use of gas-phase clusters for deriving atomic polarizabilities can overestimate polarization in condensed phase.(126) Vosmeer et al. used a combined QM/MM approach to estimate condensed-phase polarization.(127) The obtained polarizabilities for water and ethanol were found to be close to those used in previous water and methanol models.(127)

Parametrization of Polarizable Force Field

In principle, polarizable force fields are more accurate and transferable than non-polarizable force fields when applied to contrasting dielectric environments, as evidenced by the better agreement with QM on gas-phase interaction energies and successful applications in various problems such as ion solvation, [START_REF] Grossfield | Ion Solvation Thermodynamics from Simulation with a Polarizable Force Field[END_REF] protein-ligand binding [START_REF] Jiao | Calculation of protein-ligand binding free energy by using a polarizable potential[END_REF] and pKa prediction.(73) However, there are cases where polarizable force fields are comparable to or even worse than nonpolarizable force fields, due to the poor quality of parameters in the polarizable force fields. [START_REF] Bc | Polarizable force fields for molecular dynamics simulations of biomolecules[END_REF]125) In general, deriving accurate force field parameters is challenging because of the large parameter space, non-linear interdependencies of parameters and limitation in the amount and quality of experimental and ab initio reference data.(132) The parameterization could be even more difficult for polarizable force fields because of the additional parameters. On the other hand, by improving the physics and utilizing high-level QM data, the parameterization process can be made more robust,(120, 125) leading to more accurate, transferable and reproducible force fields. Nonetheless, due to the limitation of ab initio methods, recently biological force fields have heavily relied on experimental data such as NMR data to refine their parameters for proteins and nucleic acids. [START_REF] Lemkul | Polarizable Force Field for DNA Based on the Classical Drude Oscillator: I. Refinement Using Quantum Mechanical Base Stacking and Conformational Energetics[END_REF]111,146) Significant efforts have been made to design systematic and automatic approaches for the parameterization of force fields. Typically, the development of a force field consists of determining reference data (QM and experimental properties), defining an objective function to measure the quality of force field, and optimizing a large set of parameters to improve the objective function. When the nonlinear interdependency between parameters is nontrivial, sophisticated optimization methods can be utilized. Methods in artificial intelligence, such as evolutionary algorithms, have been applied to optimize force field parameters.(58, 130) Wang et al.(132) developed the ForceBalance software to tackle several problems in force field development. Specially tuned objective function, regularization and gradient-based optimization algorithm were used to improve the optimization results.(132) To choose subset of parameters to optimize, one can utilize sensitivity analysis or test optimization using cheaper objective functions.(133, 143) These automated algorithms can save substantial human efforts. However, the quality of resulting force fields depends critically on the reference data set and various "weights" assigned to different reference properties, due to the imperfection in the underlying models. Additionally, overfitting can potentially lead to problems in transferability. 27) independently developed a PME algorithm based on spherical tensors. The algorithm when applied to quadrupoles only slows down the calculation by 1.5 to 2 times compared to a charge-only model, in part because a shorter real-space cutoff is possible with fast-decaying higher multipole moments. The result is quite encouraging considering that the charge-dipole-traceless quadrupole model has nine degrees of freedom for each atom. [START_REF] Giese | Multipolar Ewald Methods, 1: Theory, Accuracy, and Performance[END_REF] This manifests the advantage of point multipole over an equivalent representation by a set of point charges. Lin derived a general formula for EMP based on Cartesian tensors. [START_REF] Lin | Generalized and efficient algorithm for computing multipole energies and gradients based on Cartesian tensors[END_REF] It was argued that this algorithm is more efficient than the spherical tensor formalism in terms of implementation, since the latter needs coordinate transformation at every time step. [START_REF] Lin | Generalized and efficient algorithm for computing multipole energies and gradients based on Cartesian tensors[END_REF] The isotropic periodic sum (IPS) method uses the so-called isotropic periodic images to represent the remote structure, so that the sum of interactions with periodic images can be solved analytically for most potential functions.(139) Compared to PME, IPS achieves 2-to 3-fold increase in efficiency for a dihydrofolate reductase (DHFR) system. IPS also has a better scaling O(N) compared to O(N log N) for PME.(139)

Efficient methods for polarizable MD simulations

Boateng developed two Cartesian tree algorithms which employ Taylor approximations and hierarchical clustering. [START_REF] Boateng | Mesh-free hierarchical clustering methods for fast evaluation of electrostatic interactions of point multipoles[END_REF] The algorithms are suitable for simulations with free-space boundary conditions, and speed up the evaluation of point multipole interactions by an order of magnitude compared to direct sum. [START_REF] Boateng | Mesh-free hierarchical clustering methods for fast evaluation of electrostatic interactions of point multipoles[END_REF] Algorithms for Evaluating Polarization Self-consistent field (SCF) iterations provides rigorous solution to polarization energy and gradient at the group state, which is needed for structure optimization, QM/MM application and energy conservation in MD simulations etc. However, full SCF calculation is computationally demanding. To speed up the SCF calculation, Brooks and coworkers developed an empirical extrapolation scheme based on perturbation theory.(113) They showed that the fourth order perturbation method (OPT4) achieves the best compromise between accuracy and efficiency, with a cost similar to that of three SCF iterations. Truncated conjugate gradient (TCG) by Aviat et al. is another recent development to accelerate the SCF calculation, which provides a solution at user-chosen cost and accuracy.(4) This method can be combined with preconditioned CG, or one additional Picard fixed point iteration after the last step. Tests on various systems demonstrate that three to four iterations provide excellent accuracy. Noticeably, TCG produces analytical forces of the corresponding energy, and thus avoid energy drift and permits large time steps in MD.(4) Beran and coworkers proposed divide-and-conquer JI (DC-JI) with overlapping blocks to accelerate the polarization solution for use with PME and multipoles in spherical harmonics. [START_REF] Nocito | Massively Parallel Implementation of Divide-and-Conquer Jacobi Iterations Using Particle-Mesh Ewald for Force Field Polarization[END_REF] The algorithm showed 20-30% improvement in speed compared to PCG or JI/DIIS and Cartesian based multipoles. [START_REF] Nocito | Massively Parallel Implementation of Divide-and-Conquer Jacobi Iterations Using Particle-Mesh Ewald for Force Field Polarization[END_REF] Extended Lagrangian is an alternative method where an additional set of electronic degrees of freedoms is propagated to approximate the SCF solution. Due to stability issues, early extended Lagrangian methods only permit small integration time step of 1 fs or less. Recently, Albaugh et al. introduced an iteration-free method, inertial extended Lagrangian with 0 SCF (iEL-0SCF). [START_REF] Albaugh | A New Method for Treating Drude Polarization in Classical Molecular Simulation[END_REF][START_REF] Albaugh | Accurate Classical Polarization Solution with No Self-Consistent Field Iterations[END_REF] In this scheme, the auxiliary dipoles drive the time evolution of real dipoles that stays close to the true SCF solution. This method allows for 6 fs time step for single-point polarizable water.(2) When used to simulate the flexible AMOEBA water model with the same 1.0 fs time step, iEL-0SCF is twice as fast as standard SCF algorithm. Future work will combine iEL-0SCF with RESPA. [START_REF] Albaugh | Accurate Classical Polarization Solution with No Self-Consistent Field Iterations[END_REF] Enhanced Molecular Dynamics and Sampling Algorithms Enhanced sampling techniques are needed for achieving the necessary sampling efficiency of biomolecular systems regardless of how force fields are used. [START_REF] Elber | Perspective: Computer simulations of long time dynamics[END_REF] Several efforts have been made in recent years to accelerate the simulations of polarizable force fields. Multiple time step algorithms have been developed to allow for very large time steps in molecular dynamics simulations. (59, 78, 83) In the extreme case, the computation speed can be accelerated by 10 to 20 times. [START_REF] Margul | A Stochastic, Resonance-Free Multiple Time-Step Algorithm for Polarizable Models That Permits Very Large Time Steps[END_REF] Dual force field approach introduced by Schnieders and coworkers, (86) takes advantage of the sampling efficiency of the fixed-point charge model (OPLS-AA) and accuracy of polarizable force fields (AMOEBA) to compute the absolute crystal decomposition thermodynamics. A similar procedure was used by Shirts and coworkers [START_REF] Dybeck | Effects of a More Accurate Polarizable Hamiltonian on Polymorph Free Energies Computed Efficiently by Reweighting Point-Charge Potentials[END_REF] to indirectly calculate the free energy of three benzene polymorphs by AMOEBA. (19) There have also been significant advances in thermodynamic and kinetic reweighting methods, [START_REF] Chodera | Dynamical reweighting: Improved estimates of dynamical properties from simulations at multiple temperatures[END_REF]136) which can in principle be combined with the dual-force field methods. Orthogonal space random walk (OSRW) and orthogonal space tempering (OST) by Yang and coworkers(71) allows more effective sampling of conformational transitions in aqueous solution, and has been utilized on crystal(108) and host-guest (5) systems with AMOEBA force field.

Recent Development of Polarizable Force Fields for Biomolecules

Over the past decades, several polarizable force fields have been developed for biological systems, including AMBER,(12, 130) AMOEBA, [START_REF] Peng | Accurate Evaluation of Ion Conductivity of the Gramicidin A Channel Using a Polarizable Force Field without Any Corrections[END_REF][START_REF] Ponder | Current Status of the AMOEBA Polarizable Force Field[END_REF] CHARMM Drude,(53) CHARMM fluctuating charge,(14, 104) SIBFA, GEM, [START_REF] Gresh | Anisotropic, Polarizable Molecular Mechanics Studies of Inter-and Intramolecular Interactions and Ligand-Macromolecule Complexes. A Bottom-Up Strategy[END_REF] and ABEEMσπ. [START_REF] Liu | Development of the ABEEMσπ Polarization Force Field for Base Pairs with Amino Acid Residue Complexes[END_REF]142) Their coverage and software implementation are summarized in Table 1. Most of the force fields are supported on GPU platforms, [START_REF] Carsten | Best bang for your buck: GPU nodes for GROMACS biomolecular simulations[END_REF][START_REF] Eastman | OpenMM 7: Rapid Development of High Performance Algorithms for Molecular Dynamics[END_REF][START_REF] Harger | Tinker-OpenMM: Absolute and relative alchemical free energies using AMOEBA on GPUs[END_REF][START_REF] Huang | Molecular dynamics simulations using the drude polarizable force field on GPUs with OpenMM: Implementation, validation, and benchmarks[END_REF]103) which provides two orders of magnitude acceleration compared to CPU and permits routine access to the microsecond time scale. Tinker-HP is a massively parallel package for polarizable MD simulations of large systems on supercomputers.(50) Below we will only overview some of recent developments, and the readers are referred to the respective literature for more details. AMBER AMBER ff02pol( 12) is one of the earliest polarizable force fields for proteins and nucleic acids. Point charge and simple induced dipole model with no damping were employed in ff02pol. Later ff12pol with Thole-style damping functions was developed to improve the accuracy of intermolecular interaction energies.(130)

AMOEBA

The AMOEBA polarizable force fields employ atomic induced dipole to model polarization and atomic multipoles up to quadrupole to represent the permanent electrostatics. AMOEBA force fields have been widely used to simulate water, ions, organic molecules and proteins. (35, 77, An automatic and systematic approach for the parameterization of AMOEBA using the ForceBalance package(132) has also been explored. Overall the AMOEBA water model reparameterized (AMOEBA14)(51) using ForceBalance(132) better reproduces high-level quantum mechanical (QM) data and experimental condensed-phase properties compared to the original AMOEBA03. Variations of the functional form were devised to improve the computational speed, including the direct polarization (iAMOEBA)(132) and united atom models (uAMOEBA).(95) Both iAMOEBA and uAMOEBA, parameterized using ForceBalance, have comparable accuracy to AMOEBA03 for predicting gas-phase and liquid properties. As an example, the liquid densities over a wide temperature range predicted by different water models are compared in Figure 1.

To improve the accuracy and transferability and mitigate the tedious parameterization process, the next-generation AMOEBA force field focuses on calibrating each energy component to highlevel QM energy decomposition such as Symmetry-Adapted Perturbation Theory (SAPT), and using automated optimization methods(132) for parameterization at large scales. For electrostatic interactions, the point charge or multipole model fails at close distances where electron clouds overlap. In this situation "charge penetration" (CP) effect must be considered. By utilizing empirical smearing functions either for charge-charge interactions only(134) or higher order multipoles,(98) the charge-penetration correction can be accurately captured. For polarization, the Thole damping function used in AMOEBA(93) was improved to better capture the explicit many-body interactions for a range of molecules at different intermolecular distances. [START_REF] Liu | Capturing Many-body Interactions with Classical Dipole Induction Models[END_REF] The polarization model also offers a way to separate the polarization energy from the charge-transfer energy in a physically consistent way. For vdW interactions, the buffered-14-7 potential used in AMOEBA is re-parametrized by targeting the SAPT exchange-repulsion and dispersion energy.

SIBFA and GEM

SIBFA(31) (Sum of Interaction Between Fragments Ab Initio

) is an ab initio polarizable force field formulated as a sum of electrostatic multipole, short-range repulsion, polarization, charge transfer and dispersion contributions, each of which is designed to reproduce its QM counterpart. It was first developed to deal with divalent cations metalloproteins(32) but extended halogen compounds [START_REF] El | Could an anisotropic molecular mechanics/dynamics potential account for sigma hole effects in the complexes of halogenated compounds?[END_REF] and nucleic acids. [START_REF] Gresh | Stacked and H-Bonded Cytosine Dimers. Analysis of the Intermolecular Interaction Energies by Parallel Quantum Chemistry and Polarizable Molecular Mechanics[END_REF] SIBFA is implemented into Tinker-HP(50) for massively parallel MD simulations. Recent developments include the Gaussian electrostatic model (GEM), which provides a more faithful representation of ab initio electron density.( 31) GEM has been incorporated into SIBFA [START_REF] Gresh | Anisotropic, Polarizable Molecular Mechanics Studies of Inter-and Intramolecular Interactions and Ligand-Macromolecule Complexes. A Bottom-Up Strategy[END_REF] and AMOEBA. [START_REF] Duke | GEM*: A Molecular Electronic Density-Based Force Field for Molecular Dynamics Simulations[END_REF] Recent Application of Polarizable Force Field Small Molecules Ionic liquid systems have received much attention because of their excellent thermal and electrochemical stability and good solvation properties. [START_REF] Kan | Polarization Effects on the Cellulose Dissolution in Ionic Liquids: Molecular Dynamics Simulations with Polarization Model and Integrated Tempering Enhanced Sampling Method[END_REF]115) Due to their charged nature, they are studied with MD simulations employing polarizable force fields. The study of dissolution of cellulose in ionic liquids shows that the conformational changes with polarizable model are broader than those with non-polarizable models. 119 Comparison to other fixed-point charge models, AMOEBA for ionic liquid resulted a better agreement with experiments on ionic liquid densities, enthalpies of vaporization, and diffusion coefficients. 120 Busch et al. studied a highly concentrated aqueous solution of proline using neutron diffraction experiments and MD simulations employing AMOEBA and CHARMM force fields. Detailed structural analysis revealed the existence of proline-proline dimers, which explains well the experimental observation. Compared to non-polarizable CHARMM force field, the polarizable AMOEBA simulation gives better agreement with the EPSR fits to the diffraction data, which is similar to ab initio (CPAIMD) methods.( 7)

Ions

The association of Mg 2+ and H2PO4 -in water may give insights into our understanding of Mg and phosphate-containing biomolecules, e.g. DNA, RNA and ATP. A recent simulation study shows that the binding free energy between Mg 2+ and H2PO4 -determined by AMOEBA simulations (-2.23 kcal/mol) closely match the experimental value (-1.7 kcal mol -1 ).( 76) Another recent quantum calculation which used a mixed explicit/continuum solvent model gave a value of -3.3 kcal mol -1 , while non-polarizable force field over-predicted the binding free energy by a factor of ten.(105) These results again emphasized the importance of polarization in highly charged systems.

Solvation of salt ions in the non-aqueous solvent has significant implications for understanding ion transport in cellulose. Noskov and coworkers performed combined experimental and computational analysis of the solvation of LiCl salt in N-methyl-acetamide (NMA). They found that polarizable Drude oscillator model was capable of reproducing energetics and geometries of the gas-phase clusters, and yielded qualitative agreement with experimental data on the concentration-dependence of solvation enthalpies. Polarization also has a dramatic impact on the computed potential of mean force (PMF) for ion permeation. [START_REF] Manin | LiCl solvation in N-methylacetamide (NMA) as a model for understanding Li(+) binding to an amide plane[END_REF] Protein and Peptide Structures Polarization plays a role in stabilizing helices by enhancing the dipole moment of peptide bonds. [START_REF] Lemkul | An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications[END_REF][START_REF] Ouyang | When are Many-Body Effects Significant?[END_REF] The enhanced unfolding of amyloid A-beta peptide is shown to be due to the mutated side chains altering the local peptide-bond dipole moments leading to local destabilization of the alpha-helix. [START_REF] Lemkul | Induced Dipole-Dipole Interactions Influence the Unfolding Pathways of Wild-Type and Mutant Amyloid β-Peptides[END_REF] Cui and coworkers studied the titration response of buried residues in staphylococcal nuclease mutants by MD simulations and Adaptive Poisson-Boltzmann Solver (APBS) calculation.

Noticeably larger structural disruption is observed upon ionization of some mutants when the Drude force field is used, compared with the non-polarizable force fields. However, due to the limited amount of experimental data for comparison, it is difficult to tell which force field is closer to reality.(144) Polarizability of nonpolar solvent also affects the protein stability. With polarization, the alphahelix is stabilized compared to beta-hairpin by about 1 kJ mol -1 per residue for methanol and chloroform and by about 2 kJ mol -1 per residue for carbon tetrachloride. This highlights that inclusion of polarizability in models for less polar and nonpolar solvents or protein environments is as important as including polarizability in models for liquid water. [START_REF] Lin | Effects of Polarizable Solvent Models upon the Relative Stability of an α-Helical and a β-Hairpin Structure of an Alanine Decapeptide[END_REF] Additionally, an MD study of alkali-acetate solutions at various concentrations indicated that polarizable force fields may be needed to accurately capture behavior of protein in electrolyte solutions using MD simulations. [START_REF] Ahlstrand | Computer simulations of alkali-acetate solutions: Accuracy of the forcefields in difference concentrations[END_REF] Protein-Ion and Protein-Ligand Binding Several recent studies have been focused on capturing the interactions of ions with proteins and nucleic acids. Using the AMOEBA polarizable-force field, many-body effects were shown to be important for ion-selectivity in Mg and Ca protein complexes. [START_REF] Huang | Mapping the Drude polarizable force field onto a multipole and induced dipole model[END_REF] With the Drude polarizable force field, the secondary coordination shells of proteins were shown to be perturbed in cationdependent manner, with significant delocalization and long-range effects of charge transfer and polarization on Ca 2+ binding.(87) Mehandzhiyski et al. showed that near equilibrium, charge transfer between metal ions and deprotonated carboxylic acids are significant. [START_REF] Mehandzhiyski | Density Functional Theory Study on the Interactions of Metal Ions with Long Chain Deprotonated Carboxylic Acids[END_REF] 2D Free energy profiles for Zn-binding to a voltage-gated proton channel (Hv1) calculated with the Drude force field were consistent with the voltage clamp fluorometry data, supporting the existence of two Zn 2+ -binding sites and the involvement of different amino acid residues in the two binding sites. [START_REF] Qiu | Molecular mechanism of Zn<sup>2+</sup> inhibition of a voltage-gated proton channel[END_REF] Several apical iodide translocation pathways have been proposed for iodide efflux out of thyroid follicular cells, including a pathway mediated by the sodium-coupled monocarboxylate transporter 1 (SMCT1), which remains controversial. Vergara-Jaque et al. evaluated the structural and functional similarities between SMCT1 and the well-studied sodium-iodide symporter (NIS) that mediates the first step of iodide entry into the thyroid. These results suggest that wild-type hSMCT1 in the inward-facing conformation may bind iodide only very weakly, which may have implications for its ability to transport iodide.(124) Qi et al. showed that the AMOEBA force field could accurately predict the binding free energy between phosphate and the phosphate binding protein (PBP). [START_REF] Qi | Elucidating the Phosphate Binding Mode of PBP: The Critical Effect of Buffer Solution[END_REF] By considering the interaction between phosphate and the buffer ligands, the thermodynamic stability between two protonation states of phosphate in the binding pocket was established, which has been ambiguous from analyses of experimental measurements.

Diffusion and Permeation of Small Molecules

Polarization effects are essential for capturing ion transport as shown by calculating the potential of mean force of Li transport through a narrow ion channel [START_REF] Manin | LiCl solvation in N-methylacetamide (NMA) as a model for understanding Li(+) binding to an amide plane[END_REF]. The free energy path for an oxygen molecule to travel along E. Coli AlkB tunnels has been determined with AMBER and AMOEBA. Both PMFs indicate passive transport of O-2 from the surface of the protein. However, the inclusion of explicit polarization shows a very large barrier for diffusion of the co-substrate out of the active site, compared with the non-polarizable potential. Also, the results suggest that the mutation of a conserved residue along the tunnel, Y178, has dramatic effects on the dynamics of AlkB and the transport of O-2 along the tunnel.(121) Zhu et al. studied the permeation behavior of 2-aminoethoxydiphenyl borate (2-APB), a broadspectrum modulator for some membrane proteins. They showed that the protonation state and therefore the polarity of the drug is critical for its partition, and that the drug is likely to switch between different protonation states along its permeation pathway. By changing the degrees of freedom, protonation further affects the thermodynamic of the permeation pathway of 2-APB, leading to different entropic contributions. A survey of 54 analog structures with the similar backbone to 2-APB showed that delicate balance between entropy and polarity plays an important role in the potency of drugs.(147) Ion Channels By explicitly introducing the multipole terms and polarization into the electrostatic potentials, the permeation free energy barrier of K+ through the gA channel is considerably reduced compared to the overestimated results obtained from the fixed-charge model. Moreover; the estimated maximum conductance, without any corrections, for both K + and Na + passing through the gA channel is much closer to the experimental results than any classical MD simulations, demonstrating the power of AMOEBA in investigating the membrane proteins. [START_REF] Peng | Accurate Evaluation of Ion Conductivity of the Gramicidin A Channel Using a Polarizable Force Field without Any Corrections[END_REF] Voltage-gated sodium (Na-v) channels play vital roles in the signal transduction of excitable cells. Upon activation of a Nav channel, the change of transmembrane voltage triggers conformational change of the voltage sensing domain, which then elicits opening of the pore domain and thus allows an influx of Na+ ions. Description of this process with atomistic details is in urgent demand. In this work, the partial activation process of the voltage sensing domain of a prokaryotic Nav channel using a polarizable force field was simulated. It was not only observed the conformational change of the voltage sensing domain from resting to preactive state, but also rigorously estimated the free energy profile along the identified reaction pathway. Comparison with the control simulation using an additive force field indicates that voltage-gating thermodynamics of Na-v channels may be inaccurately described without considering the electrostatic polarization effect.(118) Roux and coworkers investigated the properties of an ion channel from the Gram-positive bacterium Tsukamurella paurometabola with a selectivity filter formed by an uncommon prolinerich sequence. Electrophysiological recordings show that it is a non-selective cation channel and that its activity depends on Ca 2+ concentration. In the crystal structure, the selectivity filter adopts a novel conformation with Ca 2+ ions bound within the filter near the pore helix where they are coordinated by backbone oxygen atoms, a recurrent motif found in multiple proteins. The binding of Ca 2+ ion in the selectivity filter controls the widening of the pore as shown in crystal structures and in molecular dynamics simulations. The structural, functional and computational data provide a characterization of this calcium-gated cationic channel. [START_REF] Dhakshnamoorthy | Structural and functional characterization of a calcium-activated cation channel from Tsukamurella paurometabola[END_REF] Nucleic Acids Studies have shown the importance of polarization for capturing the flexibility and stacking behavior of nucleic acids with ions.(34, 106, 146) Gresh et. al. showed that SIBFA, QC multipoles, and explicit representation of lone pairs is essential to account for coulomb anisotropies and exchange repulsion when studying the stacking of cytosine dimers and a doubly H-bonded dimer. [START_REF] Gresh | Stacked and H-Bonded Cytosine Dimers. Analysis of the Intermolecular Interaction Energies by Parallel Quantum Chemistry and Polarizable Molecular Mechanics[END_REF] For DFT-D3 optimization RNA shows higher flexibility as compared with the MM showing that DFT-D3 methods complement MD results and can be used to benchmark faster computational methods. [START_REF] Kruse | QM Computations on Complete Nucleic Acids Building Blocks: Analysis of the Sarcin-Ricin RNA Motif Using DFT-D3, HF-3c, PM6-D3H, and MM Approaches[END_REF] Polarization effects significant when introducing a second ion to a G-stem quadruplex indicating a delicate balance between electrostatic and induction.(28) Gao et al. showed that orbital overlap is vital for capturing short hydrogen bonding [START_REF] Gao | Improved Polarizable Dipole-Dipole Interaction Model for Hydrogen Bonding, Stacking, T-Shaped, and X-H...pi Interactions[END_REF] Simulations with Drude polarizable force field yielded near-quantitative agreement with experimental measurements of the equilibrium between the base-paired and flipped states. Free energy barriers to base flipping are reduced by changes in dipole moments of both the flipped bases that favor solvation of the bases in the open state and water molecules adjacent to the flipping base. [START_REF] Lemkul | Induced Polarization Influences the Fundamental Forces in DNA Base Flipping[END_REF] Polarizable force field using the classical Drude oscillator better reproduced experimental solution X-ray scattering for DNA compared to non-polarizable AMBER parmbsc0 and CHARMM36 force fields. The simulations also indicate that the conformational properties of DNA in solution are sensitive to the type of monovalent ion. The primary conformational mode associated with the variations is a contraction of the DNA minor groove width with decreasing cation size.(107) Song et al. showed that polarization of the nucleobases by K + enhanced electrostatic attraction between the base and ions.(114) This increased attractive interaction is critical to stabilizing the stem-loop junction ions in G-DNA. With non-polarizable force fields, the top and bottom cations would be released into the solvent within just a few nanoseconds, and an incorrect bifurcated bonding geometry of G-DNA will be adopted which is not observed in experiments.

Kinetics

There have been limited kinetics studies by polarizable force fields. Lin et al. showed that the relaxation rate of proteins was overestimated by one order of magnitude by fixed-charge force fields; while relaxation can be slowed down by using polarizable FF (AMBER12pol), it cannot make up for the gap in timescales between experiments and simulations. [START_REF] Lin | Implementing electrostatic polarization cannot fill the gap between experimental and theoretical measurements for the ultrafast fluorescence decay of myoglobin[END_REF] These results indicate the certain areas for improvement of force fields, such as consideration of conformational transition rates.(44) Interaction with Electric Field Electronic polarization is essential for modeling the interaction with electric field, such as in the simulations of THz spectra. AMOEBA force field was used to simulate the THz spectra of two zwitterionic amino acids (glycine and valine) in aqueous solution. After detailed check of the THz spectral assignments, the mode-specific spectral decomposition into intramolecular solute motions, and solute-water cross-correlation modes, the authors found promising agreement of AMOEBA and ab initial molecular dynamics (AIMD) data for both systems.( 23)

QM/MM

Polarizable force fields (e.g. AMOEBA, GEM and CHARMM-Drude) have been applied to the hybrid QM/MM method to better describe the environment of the QM region. [START_REF] Ke | LICHEM: A QM/MM program for simulations with multipolar and polarizable force fields[END_REF][START_REF] Gökcan | QM/MM Simulations with the Gaussian Electrostatic Model: A Density-based Polarizable Potential[END_REF][START_REF] Lipparini | Polarizable Molecular Dynamics in a Polarizable Continuum Solvent[END_REF][START_REF] Loco | Hybrid QM/MM Molecular Dynamics with AMOEBA Polarizable Embedding[END_REF][START_REF] Loco | A QM/MM Approach Using the AMOEBA Polarizable Embedding: From Ground State Energies to Electronic Excitations[END_REF]145) The methods have been implemented in software interfaces, such as Gaussian/TINKER,(69) Psi4/TINKER(30), Q-Chem/CHARMM.(145) The use of polarizable force field improves both ground-state energy and structure [START_REF] Gökcan | QM/MM Simulations with the Gaussian Electrostatic Model: A Density-based Polarizable Potential[END_REF][START_REF] Loco | A QM/MM Approach Using the AMOEBA Polarizable Embedding: From Ground State Energies to Electronic Excitations[END_REF] and excited-state spectral properties.(67, 69, 145) Loco et al. used QM/MM simulations with B3LYP and AMOEBA to study the color tuning in Carotenoid pigment-crustacyanin complexes.(67) It was found that polarizable force field and MD simulations are necessary to obtain quantitative predictions of the spectrum. The high color tunability of the pigment-protein complex was explained by the bond length alternation in the long-chain carotenoids modulated by the dynamical protein environment.

Summary and Outlook

Polarizable force fields have grown steadily during the past few years in terms of computational efficiency, model accuracy and applications to biomolecular systems. The AMEOBA force field that has recently been extended to DNA and RNA shows an improved description of the conformational ensemble in different environments. The CHARMM Drude force field has recently been refined for DNA and extended to carbohydrates and halogenated molecules. Advances in GPU computing, polarization and simulation algorithms have provided access to the microsecond time scale with polarizable force fields, and the computational overhead compared to fixedcharge force field has been significantly reduced.

The applications of polarizable force fields have provided many new insights. Recent studies using polarizable force fields have demonstrated the critical role of polarization for the stability of nucleic acids and proteins, base-pair flipping, ion distribution around DNA, diffusion and permeation of small molecules. In general, simulations with polarization force fields agree better with experiments.

There are several future directions for polarizable force fields. The underlying physical models, particularly for short-range interactions such as charge penetration [START_REF] Rackers | An optimized charge penetration model for use with the AMOEBA force field[END_REF]134) and charge transfer, can be further improved and incorporated to achieve significantly better accuracy and transferability with little additional computational cost. Such models will allow the utilization of ab initio EDA to systematically calibrate the energy components and robust parameterization, and reduce the reliance on error cancellation. Systematic and automatic methods for assigning atom types, parameterization using (ab initio and experimental) reference data set are crucial for reducing human efforts and errors and improving reproducibility. An exciting new direction is to combine polarizable force fields with enhanced sampling methods such as orthogonal space sampling (OSS),(71) Markov state models (MSMs) and Milestoning, [START_REF] Elber | Perspective: Computer simulations of long time dynamics[END_REF][START_REF] Husic | Markov State Models: From an Art to a Science[END_REF] which will significantly extend the time and length scales of polarization force fields simulations to areas such as protein and nucleic acids conformational dynamics. These studies would provide crucial feedback to the force field development and insights into our understanding of the intermolecular forces and how they affect the structure and properties of biomolecular systems.

  Algorithms for Computing Long Range Electrostatics Particle Mesh Ewald (PME) is an efficient algorithm for calculating electrostatics interactions under periodic conditions. Recently, several groups have developed generalized and efficient PME algorithms for electric multipoles (EMP) of arbitrary order.(27, 63, 112) Simmonett et al.(112) and Giese et al.(

  101, 102, 111, 137) Mu et al. showed that the σ-hole effect can be captured by AMOEBA.(85) Recently, Zhang et al. developed the AMOEBA force field for DNA and RNA.(146) The force field was extensively validated through 35 microseconds of MD simulations. The simulated solution and crystal structures of DNA duplexes, RNA duplexes and hairpins agree with NMR structures with RMSDs < 2.0 Å. Notably, the interconversion between A-and B-form DNAs was observed in ethanol-water mixtures, (see Figure 2) indicating a balanced description of the stabilities of different forms. Clavaguéra and coworkers developed the AMOEBA force field for Fe(II) and the heme cofactor in ferrous and ferric form.(109, 138) The parameters were validated for energy calculation of larger clusters and MD simulations of cytochromes, showing good agreement with DFT and NMR data. To match the energy components from ab initio calculations, Xia et al. incorporated an explicit charge-transfer term into the AMOEBA force field for Fe(III).(140) For the transition metal ions Cu 2+ and Zn 2+, AMOEBA-VB model was derived.(141) This model generates correct ion-ligand geometry and energetics for both QM gas-phase clusters and the coordination of first solvation shell structure of their aqueous solutions. To better model the water ligand exchange rate around Mg 2+ , Kurnikov and Kurnikova (49) treated the polarizability of AMOEBA water as variables according to the distance between water and Mg 2+ .

( 96 )

 96 CHARMM Drude force fieldThe CHARMM Drude force field utilize the Drude oscillator model for polarization and off-center charges to represent anisotropic charge distributions.[START_REF] Lemkul | An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications[END_REF] The Drude force field covers the majority of molecules commonly used in molecular simulations, including small organic molecules, protein/peptide, DNA, and lipid. (See Ref[START_REF] Lemkul | An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications[END_REF] and references therein) Lin et al. improved the Drude force field for both aliphatic and aromatic halogenated molecules by including off-site charges, anisotropic polarizability on halogen and vdW parameter on the Drude particle.[START_REF] Lin | Polarizable Empirical Force Field for Halogen-Containing Compounds Based on the Classical Drude Oscillator[END_REF] The Drude model for DNA has been refined to resolve problems of the previous version Drude-2013, such as the weak base stacking in A-and B-DNA, the unwinding of Z-DNA.[START_REF] Lemkul | Polarizable Force Field for DNA Based on the Classical Drude Oscillator: I. Refinement Using Quantum Mechanical Base Stacking and Conformational Energetics[END_REF] Ions and water models have been adjusted accordingly to obtain better compatibility with DNA model by fitting to QM energy profiles and aqueous solution properties.(53) Similar strategies were used to develop the ion-protein model.[START_REF] Lopes | Polarizable Force Field for Peptides and Proteins Based on the Classical Drude Oscillator[END_REF] It is notable that the Drude model for lipids has also appeared very recently.[START_REF] Li | Drude Polarizable Force Field for Molecular Dynamics Simulations of Saturated and Unsaturated Zwitterionic Lipids[END_REF]

Figure 1 .

 1 Figure 1. Density of liquid water over the temperature range of ~250-370 K at atmospheric pressure. The data were reproduced from the original papers by using WebPlotDigitizer (https://automeris.io/WebPlotDigitizer).

Figure 2 .

 2 Figure 2. Transition from A-DNA to B-DNA in ethanol/water solution as captured by AMOEBA simulations. Figure taken from reference(146).

Table 1 .

 1 Polarizable force fields for biomolecules and available software

	Force field	Systems	Software
	AMBER	Protein, nucleic acids AMBER
	ff02pol/ff12pol
	AMOEBA	Protein, DNA, RNA,	Tinker, Tinker-HP,
		lipids	AMBER, OpenMM
	CHARMM Drude	Protein, DNA, lipid,	CHARMM, NAMD,
		carbohydrates	GROMACS, OpenMM
	CHARMM FQ	Protein, lipids,	CHARMM
		carbohydrates
	SIBFA	Protein, nucleic acids Tinker-HP
	ABEEMσπ	Protein, nucleic acids Tinker (modified)

Acknowledgement

We are grateful for support by the National Institutes of Health (Nos. R01GM106137 and R01GM114237), the Robert A. Welch Founda-tion (No. F-1691) and CPRIT(RP160657).

Literature Cited