20 years of polarizable force field development for biomolecular systems

(1)

2 0 Y E A R S O F P O L A R I Z A B L E

F O R C E F I E L D D E V E L O P M E N T

for biomolecular systems

thor van heesch

Supervised by Daan P. Geerke and Paola Gori-Giorg

(2)

contents 2

abstract

In this literature study we aimed to answer the following question: What has changed in the outlook on polarizable force field development during the last 20 years? The theory, history, methods, and applications of polarizable force fields have been discussed to address this question. This investigation showed that the quality of the force field potential is detrimental for any type of atomistic sampling technique. If the underlying energy function contains flaws, these flaws are one way or another embedded in the fast amount of data we seek to understand to repro-duce a wide range of quantifiable observables accurately. Subsequently, we identi-fied two key challenges associated with the development of polarizable force fields: First, how to determine transferable and accurate parameters sets, and second, how to advance the underlying physical model without computational overhead. Af-ter twenty year the classical avenue of force field development has the modeling community finally made the transition towards a general acceptance of the need to develop more physically sound models. During the last 5 years machine learning techniques emerged to provide new means to remove bottlenecks in the current process towards the development of an accurate force field potential.

1 introduction

Next-generation atomistic force fields include polarization effects for the simula-tion of biomolecular systems. How this apparent change happened is a different story, to explain this generational transition we begin our study in the late 1950s. Physicists Bernie Alder and Thomas Wainwright were the first to translate digital computation into the study of many particle systems.[1] Eventually, the offspring of their research brought the simulation method called molecular dynamics (MD) to reality. Currently, classical MD simulation methods are being applied to study a multitude of physical, chemical, and biological systems, ranging from pure liquids to large complex systems such as proteins and cell membranes. [2, 3] As a result atomistic simulations have become an important tool to understand fundamental processes of biological systems.

Since the pioneering work of Alder and Wainwright, computing performance has increased by more than trillion fold.[4] This rapid development of digital machines lead to the expansion of system sizes and increase of timescales. As this was not the only advancement, since technological advancement inspired in an equal man-ner the drive to search for faster, more efficient and accurate underlying physical models for our simulation methods. In response a diverse set of atomistic sim-ulations methods developed (co-)independently for the simulation of electrolytes, ionic liquids, metal organic frame works, biomolecular systems, and other types of nano-materials.

In this study we will keep our attention focused on the simulation of biomolecular systems. For this particular simulation field, the inclusion of explicit polarization effect has dominated the evolutionary process towards obtaining an improved de-scription of the systems under investigation. To understand the reasons why the biomolecular simulation community chose to include polarization effects into the

(4)

force fields: basics, caveats and extensions 4

atomistic model, we will let ourselves get inspired by the following remarkable ob-servation: The number of reviews that specifically address recent developments in polarizable force fields is numerous: [5, 6, 7, 8, 9, 10,11, 12, 13,14, 15,16]. This (incomplete) selection of publications account together for more than 2000 citations according to Google Scholar. The period over which these reviews span, from 2001 till the present, amounts to almost two decades. Therefore, it might be fruitful to ask ourselves the following question: what has changed in the outlook on polarizable force field development during the last 20 years?

First, this approach will allow us to assess the progress made and may show us the edge of where polarizable force fields are right now. Second, this method could expose whether there are persistent factors that hamper our progress towards a po-larizable force field for general use in biomolecular simulation. After the evaluation of this question we will discuss the open ends and look forward to how these new challenges can be solved. This will include the exploration of how Machine Learn-ing (ML) methods and novel samplLearn-ing techniques can aid towards the development of general use biomolecular force fields. Finally, we will summarise these efforts and provide an outlook on the current status of polarizable force field development for biomolecular simulations.

2 force fields: basics, caveats and extensions

Before we will proceed with an analysis of the aforementioned reviews we will provide the necessary background information to understand what is about to be discussed. The next sections will therefore explain the basics behind force fields together with the caveats of choosing a model that is based on fixed point charges, i.e., a non-polarizable force field. The first section discusses the components of the potential-energy function in terms of bonded and non-bonded interactions and seemingly continues with the caveats related to this classical approach. This next section will focus on the methods used to refine the model for its intent and purpose. Finally, the contributions of anisotropic charge distributions are discussed in the last section.

2.1 The Classical Approach

In molecular dynamics (MD) simulations are the atoms more often than not treated as point-like particles. Their reciprocal interactions in combination with the dy-namic equations of motion determine how the system will evolve over time. This simple approach is well suited to simulate the collective behavior of atoms in molec-ular structures and ensembles. Moreover, when a few assumptions are set aside, this level of theory can determine both the micro- and macroscopic properties of the respective system within the line of expectations. [17]

However, there is still plenty of room for uncertainty to develop during the course of the simulation. The biggest assumption that makes atomistic simulations to some extent a conjecture is the lack of an explicit expressions for electrons. In addition, the lack hereof assumes that the system is always in the electronic ground state.

(5)

As it where, these factors seem like a major shortcoming for creating a faithful molecular model. But, the simplicity of the atomistic model is also its greatest strength. Accounting for the electronic behavior implicitly allows for long simula-tion times up to the regime of µs, while maintaining relative low computasimula-tional cost.[18,19,20] Given the reason for treating atoms as simple particles, how do we resolve the need for inclusion of fundamental electrostatic effects without turning to on-the-fly quantum mechanical calculations?

Atomistic simulation methods have developed an implicit way to account for the existence and effects of electrons. The electronic energy is formulated as a para-metric function of the nuclear coordinates and corresponding parameters are sub-sequently fitted to experimental or higher level computational data. Such a param-eterised potential that describes all forces in the system is called a force field. This description allows us to think more pragmatic about what chemistry conveys in molecular simulations. Simply put, chemistry becomes the knowing of the energy as a function of nuclear coordinates and molecular properties become the knowing of how the energy changes upon adding a perturbation to the system. [17] Deter-mining the physics using this approach is justified in classical mechanics, because this is the basis of methods that give access to the Boltzmann weighted ensembles from which macroscopic properties of the system directly follow. [21] Still, there are numerous ways to construction of such energy functions. In the following section we will therefore discuss the most important components for the construction of a force field.

2.1.1 An atomistic-potential energy function

Figure 1: Schematic illustration of the terms in a classical fixed-charge force field, i.e. bond stretching (Ebond), bond-angle bending (Eangle), dihedral- angle torsion (Etor-sion), and improper dihedral-angle bending (Eimproper) as well as van der Waals (EvdW) and electrostatic (Eele) interactions. Figure courtesy to ref. [22]

The majority of force fields describe the interactions between the atoms in the sys-tem via a potential energy function of the atomic coordinates. The total molecular potential energy is then composed of many terms and its exact form is unknown. Today classical biomolecular force fields use essentially the same gross approxima-tion for the potential energy as proposed by Levitt and Lifson in 1969: [23]

Epot(~r) = Espr(~r) + Eang(~r) + Edih(~r) + Eimp(~r) + EvdW(~r) + Eele(~r) (1) Here is Epot _{the potential-energy function of the system and ~r are the} coordi-nates of all atoms in the system. Espr, Eang, Edih, Eimp are the bonded terms

(6)

for bond stretching, bond-angle bending, dihedral-angle torsion, and improper dihedral-angle bending (or out-of-plane distortions) in the description of molecules. Together these terms account for all the energy contributions due to covalent inter-actions: Ecov _{= E}spr₍_{~r) + E}ang_{(~r) + E}dih₍_{~r) + E}imp₍_{~r). The interactions between} atoms that are not directly connected via covalent bonds and bond angles consist of the potential-energy term for all nonbonded interactions: Enb = EvdW + Eele. Where the van der Waals term, EvdW_{, accounts for dispersion and cocore} re-pulsion, and the latter term, Eele, for the electrostatic interactions of the system. Together with the potential-energy terms for the covalent interactions this again sums up to the full potential-energy expression of the force field:

Epot= Ecov+ Enb (2)

See figure 1 for a schematic illustration of the terms in a classical fixed-charged force field.

Without question, additional terms can be added to this general expression of the Epot_{, e.g., in the early days hydrogen-bonding interaction terms, or cross-terms} that describe coupling between the first three terms of the covalent interactions and other types of restraining potential-energy functions could be utilized to prevent the system from drifting away. [17] However, these extra descriptors are usually correction terms to account for omissions made in the other terms to still obtain realistic free energy profiles, instead of improving the description of the other terms. Despite the major importance of all the electronic effects, the partial charges of the atoms are fixed during simulations that are based upon classical force fields. Historically, this choice was made due to the high computational cost and added complexity associated with a more physics based force fields. [24]

Starting from the static decision the continuous development of fixed-charge force fields has actually shown to be surprisingly successful. A recent publication by S. Riniker illustrates the development of this class of force fields and provides an overview on the major force-field families, with detailed discussions about the dif-ferent covalent and nonbonded force-field terms, parametrization strategies and the historic background of these models. [22] Therefore, the aspects mentioned above revolving around non-additive force fields will only be touched upon in this litera-ture study. Instead, we will mainly focus on the effects and importance of extending the description of the classical electrostatics in biomolecular simulations.

2.1.2 Van der Waals potential

In (almost) all fixed-charged force-fields, the van der Waals interactions are de-scribed using a 12-6 Lennard-Jones (LJ) functional form

EvdW_ij (rij) =

C₁₂(i, j) r12_ij −

C₆(i, j)

r6_ij (3)

with C12(i, j) = 4ij(σij)12, C6(i, j) = 4ij(σij)6, where rij is the distance be-tween atoms i and j, and where ijdenotes the depth of the attractive well and σij the inter-particle distance where the potential changes sign.[25,26]

The function of van der Waals energy, Evdw_{, is to describe the repulsion or} at-traction between atoms that are not directly bonded. Often the Evdwis interpreted

(7)

as the non-polar part of the interaction not related to electrostatic energy due to charges.[17] Hence, the atoms in the system are grouped into atom types, which all have respective LJ parameters and typically represent atoms in a specific chemical environment, for example, a nitrogen in an amine group or oxygen in a hydroxyl-group, and so forth.

However, it should be noted that this widely used non-bonded term also has some (major) flaws embedded. Recently Frenkel et. al noted that over the course of time the LJ potential, which originally has been derived to describe the cohesive energy of crystals of noble gases, is used so often for systems where it is not expected to be particularly realistic that the disadvantages of the LJ potential have become very relevant these days.[27] The implications of the arguments made against the LJ 12-6 potential can be summarised as follows: the LJ 12-6 potential is anything but a well-defined standard and, in particular for proteins and nano-colloids it is not a good model. Mainly, because the range of attraction is too large compared to the effective diameter.

To top it off, Frenkel et al. argue there is not obvious indication of why a truncated LJ 12-6 potential should have any special merits that outweigh its disadvantages.[27] Instead, they argue that no such advantage exists and present a simple model pair potential to unify all the different unambiguous truncation methods used in the past. In contrast to all the shifted, truncated and interpolated LJ 12-6 models, once a cut-off distance, rc, is specified the model is uniquely defined.1 If required, their class of LJ-like potentials could be used for numerical studies of systems of particles with short-ranged attraction. For example, a set of thermodynamic and transport properties is reported for the cases rc=2.0 (“atomic liquid”) and rc=1.2 (“colloidal suspension”).[27] In short, the main message is that we should be careful with con-tinuing trends based old habits during the construction of next-generation force fields.

2.1.3 Electrostatic energy

The electrostatic interaction of the nonbonded interaction should account for the internal (re)distribution of the electrons, creating positive and negative parts of the molecule. [10] At the lowest order of approximation, this electronic behaviour can be modeled by assigning a bond dipole to the bond or, more conventionally, partial charges to each atom. For example, in classical fixed-charge force fields, this pairwise Coulomb interactions between the point charges of atoms i and j is considered as Eele_ij (r_ij) = qiqj 4πε0ε1 1 rij (4) where q is the partial charge, ε0 the electric constant, ε1the background dielectric permittivity, and rij the distance between atoms i and j.[22] The classical electro-static interactions can be viewed as the remaining long-range part, (1_r), of the quan-tum mechanical electrostatic interactions, after all short-range contributions like bonds, repulsion and dispersion are removed. Still, obtaining a good description of the electrostatic interaction between molecules (or between different parts of the same molecule) is one of the big challenges in force field developments. [10] The 1 The united model potential vanishes quadratically at a cut-off distance r_c and represents now specific substance. The authors stress that their model rather represent generic models. However, the authors mention, I cite, "very often the ubiquitious LJ 12-6 potential is used in exactly the same way."[27_]

(8)

predictive ability of MD simulations relies on the accuracy of the underlying force field.

2.2 The caveats of point-charge electrostatics

From equation 3 it is evident that not all types of electrostatic interactions are ac-counted for, and as such the deficiencies in the electrostatic potential that result from the fixed partial atomic charge approximation have been widely acknowl-edged over the past decade.[11, 10,28, 12, 29, 9, 13, 30] Both Jensen and Hagler independently made special efforts to summarize all the components that currently are unaccounted for in most “standard” force fields and are likely required to achieve experimental accuracy in (bio)molecular force fields.[17,31] The following 8 points summarize the combined efforts in addressing the caveats associated with fixed-point charge force fields:

i) Electronic polarization effects are not account for, as fixed-point charge

mod-els only two-body interactions are included. However, for polar bodies the three-body contribution is quite significant, perhaps 10–20% of the two-body term. [32] Moreover, the response of the molecular dipole moments to vari-ances in dielectric conditions are not taken into account by the fixed charge approximation models.[33,34]

ii) The geometric dependence of charge and higher order multi-poles is omitted. [35, 36] Studies have shown that both partial atomic charges and higher-order elec-tric moments have a significant dependence on the geometry. Consequently, these quantities cannot fulfill the requirement of transferability in parameter sets.

iii) The effects and energetics associated with charge transfer (CT) that occur

be-tween compounds, as for example in H-bonded complexes, are not accounted for in additive force fields. [37,38,39]

iv) The partial charge model cannot correctly model the electrostatics surrounding

a molecule, resulting in errors in the electrostatic potential ranging between 10to 20 kJ/mol.[17]

From section 2.1.2 we already concluded that the use of the standard LJ potential is not really justified. In addition, the following points will specifically address caveats hidden in the van der Waals potential and why these components should be accounted for or revisited.

v) The r−6 dispersion is an approximation omitting higher order terms, which

have been indicated to be important for close packed systems such as globular proteins. [40,41]

vi) Short range electrostatic penetration (CP) effects due to electron cloud overlap

are un accounted for. [42, 43, 44, 45] CP effects arises naturally when two atoms, at a typical van der Waals distance apart, penetrate each others electron clouds. The results is an attractive interaction, which cannot be described by a fixed partial charge model and the effect is instead accounted for implicitly in the van der Waals energy term.

(9)

vii) Another reason to revisit the van der Waals potential originates from the fact

that atoms are not always spherical, which results in anisotropy of van der Waals parameters. [46,47,48,45,49]

viii) The need to account for anisotropic interactions has been further

demon-strated by showing that inclusion of anisotropic exchange-repulsion, charge penetration, and dispersion effects, along with atomic multi-poles yielded substantive improvement in accounting for 91,000 dimer energies of 13 small organic compounds as well as several experimental systems. [50]

What should be noted from this exhaustive but most certainly incomplete list is that errors of a few kJ/mol are enough to draw the wrong qualitative conclusion from a simulation. In addition, as a more general observation about the points above, the origin of the errors made in standard force fields can be condensed to five topics that either require re-evaluation or an improved description. In summary, the 5 main caveats of point-charge electrostatics can be traced back to omissions and or lack thereof in polarization, dispersion, anisotropic, charge transfer and charge penetration effects, see figure2for an illustration of these 5 themes.

Figure 2: Electron cloud representations that illustrate 5 main caveats of point-charge electro-statics: electronic polarization, dispersion, anisotropy, charge transfer and charge penetration effects.

As mentioned in the first sentence of this study, the biomolecular simulation community chose to improve the classical potential by accounting for electronic polarization effects first. One of the main reasons is that we are interested in sim-ulations of solvated biological macro-molecules. In these types of simsim-ulations vary the electronic properties as a function of the environment, especially electrostatic interactions represent the dominant interaction in polar environments. The commu-nity therefore anticipated to yield a more physically realistic and consistent model if one would start including polarization effects.[51] Mainly because additive force fields do not provide information about the dependence of the charge distribution in the thermodynamic state of the system, nor do they resolve molecular motions causing fluctuations in the electric field.[52,53]

In addition, point charges are inadequate to reflect the electrostatic interactions in systems with varying polarity, such as proteins. Besides, omitting explicit po-larizability effects draw on when defining molecular interactions occurring in am-phiphilic environments,[54,55] as the contribution of electrostatic static interactions to the system will either be overestimated or underestimated.[56] However, when

(10)

electronic polarization effects are taken into account the dipoles moments of the molecular fragments in core of the system are actually susceptible to the total elec-tric field produced by the rest of the system as well as its direct molecular environ-ment. The inclusion of an extra set of polarizability parameters is therefore consid-ered to be the answer in treating hydrophobic versus highly polar environments in a transferable way.[54]

The perfect force field potential obviously does not only account for polarization effects, however, we hope that the message is clear: there is plenty of room to improve the classical force field. Still, it remains to be seen if a combination of tech-nological advancements and scientific goodwill will achieve a long overdue transition in atomistic simulations that will surpass the elementary description of electrostat-ics that atomic point-charges currently offer.[17, 31,30] In the remaining sections of this chapter we will therefore explore the physical phenomenon of polarizabil-ity in more depth and show how to upgrade the classical force field by explicitly accounting for polarization effects.

2.3 Physical phenomenon of polarizability

In general, the polarization of a body depends on the electric field strength. This relation can be explained as follows: All matter is built up of electrically charged particles. To the core they are either negative electrons or positive nuclei. Togther these particles combine to neutral atoms and molecules. But in other cases the combination results in charged particles, for example, ions either in solutions or crystals. [57] What these constituents have in common is that during the act of po-larization the electron density of a these particles is being reshuffled by an external electric field (usually generated by other molecules). [53]

Subsequently, the electric field polarizes the electron density away from the nuclei, creating a depletion of charge on one side of the molecule (δ+) and an increase in charge on the other side (δ−_{). The result is a more asymmetric distribution} across the electron density that causes a dipole along the field lines. The difference between the dipole moments before and after the application of the field is defined as the induced dipole moment. And the degree to which these induced dipoles align with (and augment) the electric field is referred to as polarizability. If a body shows an induced dipole moment differing from zero upon application of a uniform field, the body is said to be polarizable. [57]

What should be noticed is that the magnitude of the polarizability of a molecule correlates with the interaction of electrons and nuclei. The number of electrons affects how much influence the nuclei have on the overall charge distribution in the molecule. Less electrons means smaller and denser electron clouds, resulting in less pronounced shielding effects. For this reason are small atoms typically less prone to become polarized by an electric field, because their electrons are more tightly localized around their nucleus. In contrast, large negative ions are easily polarized, since diffuse electrons clouds and large atomic radii limit the interaction of the outer electrons with the nucleus.

In most cases the polarizability is well described as isotropic and scalar propor-tionality can be assumed. When the inducible dipole moments, ~µi, are assigned to (heavy) atoms i they can subsequently be determined by using equation 1, where

(11)

α_i is the polarizability of the species with units of Cm2V−1, and the electric field at that site, ~Ei. [53,52]

~

µi= αi(4πε0) ~Ei (5)

Eele= Eself+ ECoulomb

ECOS_self = Σ1 2kD,id

2 i

Here, the atomic or molecular polarizability tensor can be represented as a tensor. In the case of isotropic polarizability, the tensor has three non- zero elements on the diagonal and αxx = αyy = αzz. However, we can also consider the case of anisotropic polarizability, the tensor is still symmetric (ααβ = αβα, α, β = x, y, z) with the diagonal elements generally not equal to each other and non-zero off-diagonal elements (in an arbitrary coordinate system). We can write the induced dipole moment as    µx µy µz   =    αxxαxyαxz αyxαyyαyz αzxαzyαzz       Ex Ey Ez    (6)

Often, the values are experimentally not known for the off-diagonal elements, although in the past various algorithms have been proposed to derive the atomic polarizability tensor in a molecule from ab initio calculations.[58,59] Yet, this may be not really necessary to capture the major polarization effects in case the polariz-ability is not very anisotropic, i.e., αxx, αyy, and αzz are not very different from each other. For example, for water the diagonal components of the molecular po-larizability are αxx = 1.415, αyy= 1.528, and αzz = 1.468 in units 4π0·10−3 nm3 [84], where the x-axis is perpendicular to the plane of the molecule and the y-axis is parallel to the line connecting the two hydrogen atoms. Since these values differ by less than 4% from their average, using an isotropically polarizable dipole in a water molecule is a simple, but yet accurate way to account for its polarizability. Furthermore, note that the polarizability of a single atom is just isotropic.

But if one considers the polarizability of a larger systems of covalently bonded atoms we have to consider the anisotropy. Meaning, the vector components of the induced dipole moments and the strength of the electric field may have different orientations. [52] This directional effect arises when a molecule or ion is not com-pletely spherical and thus the molecular polarizability cannot be described as a scalar component. Thus, for a complete description the use of a polarizability ten-sor is necessary (while still assuming that the effects remain linear). [57] When the anisotropic component the polarizability should be taken into account will be fur-ther discussed in the implementation methods of polarization methods in section 2.5.

2.4 Common implementation methods of electronic polarization

Polarization is a response property, in terms of molecular interactions polarization this means we are dealing with non-additivity, meaning that if molecule A is polar-ized by molecule B, A will interact differently with molecule C if A would not have

(12)

been polarized by B. We thus have to find a way such that our model can account for this non-additive response in its environment. Over the past decade this important question of how to explicitly treat polarization effects has been answered in various ways. Hence, multiple families of force fields have explored and implemented po-larization effects in order to model atom popo-larization during molecular simulations. However, the by far three most popular approaches are the induced point dipole (IPD), Drude oscillator (DO) and fluctuating charge (FQ) model. Figure3illustrates the electron clouds of the three different implementation methods compared to the fixed-charge model.

Figure 3: Electron cloud represented by the fixed-charge model and polarization models. Both the induced point dipole and the Drude oscillator models can represent the deformation of the electron cloud (dipole induction). The fluctuating charge model can represent the redistribution of charges within a molecule, while the induced dipole or the Drude oscillator models are required to represent the electron de-formation of monoatomic ions and out-of-plane polarization. Figure courtesy to ref. [13]

By accounting for the extra polarization term, the total electrostatic energy is now a sum of the Coulomb energy between all the charges and dipoles in the system and a self-energy term corresponding to the work needed to change the charge distribution: [13]

Eelec= Eself+ ECoulomb (6)

To account for the non-additive polarization effects we only have to define and solve the expressions for the self-energy corresponding to the polarization model. The following paragraphs will describe the three main polarizable models in more detail.

2.4.1 Fluctuating charge

In the fluctuating charge model is the partial charge of each atom placed at the site of the atomic nucleus. The basis of this model is the most similar to the classical non-polarizable force fields. However, the partial charges in the fluctuating charge (FQ) treatment are not static. Instead, they are free to change during the simulation. This dynamic flow of charges allows the model to account for polarization in the system. FQs are implemented by assigning fictitious masses to the charges that are evaluated as separate degrees of freedom in the dynamic equation. Charges keep flowing between atoms until instantaneous electronegativities on the atomic sites are equalized. [60] These instantaneous electronegativities depend on the

(13)

elec-force fields: basics, caveats and extensions 13

tronegativity of the atom, the hardness that represents the resistance of the electron to flow from or towards the atom and lastly the external electrostatic potential.

The main advantage of the FQ force field is that no extra interaction terms than in classical force fields are required, and the electrostatic interactions are calculated using a standard Coulomb potential. This cost effectively manner represents polar-ization that occurs in the direction of atomic bonds with a breeze. Through-bond and through-space polarization are actually treated equivalently in the FQ model. Use this same principle of treating polarization effects in such a way enables in principle the description of intermolecular charge transfer. [17] However, the major drawback of the FQ approach becomes apparent when polarization in any direction other than the bond should be accounted for. [11] Take for example an aromatic molecule such as benzene, the FQ cannot account for the charge polarization when the external field is perpendicular to the plane of the molecule. The lack of out-of-plane charge density representation can easily be solved by including additional point charges. [61,62] Up until now the charge equilibration (CHEQ) force field has put the most effort in applying the FQ model in biomolecular simulations. [63,64] However, the ABEEMσπ polarization force field that is also based upon the FQ principle has been made applicable for the simulation of base pairs with amino acid residue complexes.[65]

2.4.2 Charge on spring

The Drude oscillator (DO) or charge-on-spring (COS) works also on the basis of representing an inducible dipole moment at polarizable sites. [66,67,52] Each po-larizable atom is depicted by a pair of point charges. However, what is different compared to the ID model is that at the core of the atom one partial charge is in-stalled and the other charge is attached via a harmonic spring. Summing these two charges yields the total partial charge of the site and the magnitude of the total partial charge is kept constant during the whole of the simulation. This charge on a spring has no mass and is called the Drude particle and similar to the use of pseudo-atoms for modeling lone pairs. The Drude particle is free to move anywhere around the atomic center in response to the electrostatic environment. The resulting displacement gives rise to atomic dipole polarizability. [68] This approach is consid-ered to be the most intuitive from a chemical point of view, as the two particles carrying both partial charges can be considered as the nucleus and electron cloud of the atom that together mimic the polarization effects. Several groups have put considerable effort in the development of Drude oscillator-based polarizable force fields, but thus far only MacKerell et al. made their polarizable force field (based on the CHARMM package) suitable for the simulation of biomolecules. [51,12]

2.4.3 Polarizable point dipole

The next method adds a set of inducible point dipoles assigned to polarization sites, i.e., to an atomic center, lone-pair or interaction site between bonds, [69,70] and maintains the framework of fixed atomic charge force fields. To account for polarization the field due to the explicit charges is calculated first. Thereafter are the dipoles calculated as the field multiplied by a local polarizability tensor. [71] Induced dipoles themselves also create an electric field, in a response to the

(14)

perma-force fields: basics, caveats and extensions 14

nent charges, and mutual polarization between induced dipoles is present. [53] Elec-trostatic interactions can therefore no longer be approximated by just charge-charge interactions described by a standard Coulomb potential. An iterative procedure is usually employed to solve the non-additive electrostatic energy term that should also includes the charge-dipole and dipole-dipole interactions. [71, 72] The most frequently used polarizable force field based on inducible dipoles is implemented in the AMBER simulation software package as the AMBER ff02 force field. [73,72] The other major induced dipole model is the AMOEBA force field.[8]

A brief background on the use of Thole screening factors. At the start of the devel-opment of the induced-dipole models it was discovered that the model is suscepti-ble to a so-called polarization catastrophe, although this phenomenon is not exclusive to this type of model. At short range point dipoles could cause artificially strong in-teractions between one another. To alleviate this problem, Thole proposed to screen the polarization interactions if two atoms cross this distance threshold with a damp-ening function, which later became known as the Thole model. [8,74]. Even though the Thole models has been considered effective for many years, however, recent studies revealed much of the inadequacies related to the use of the Thole damping function.[75,76]

One alternative strategy to prevent the “polarization catastrophe” from happen-ing that has been used in the "early days" for solvent models is done by damphappen-ing the linear dependence between the polarization and the electrostatic field. In this case the polarizabilities are damped when the electrostatic field has higher magnitude than a predefined threshold to prevent the unlimited polarization.[77,78] More re-cently, simple modifications ought to remedy some of the inadequacies of the Thole model, which resulted in improvements by separating the parameters of the nonlin-ear effects and optimizing the exponent in the many-body energy function.[79]

2.5 Accounting for anisotropic interactions

Above, we have discussed that one of the major short-comings of most force fields is the inability to model anisotropic charge distributions. Point charges are not able to reproduce relevant interaction features, or they are unable to properly estimate the directionality of the interactions.[80] Having the proper knowledge of about this spatial charge effect embedded in the force field is critical for determining the equilibrium geometry and energy of molecular complexes. Anisotropy becomes especially important for charge distribution such as σ-holes, lone pairs, and aro-matic systems. A few examples of planar molecules that have obvious polarizability anisotropy are shown in figure4. Considering these type and conformational spe-cific charge effects, it is no easy task to develop a polarizable force field that properly models all anisotropic effects. During the quest of how to account for anisotropy it is therefore important to ask the questions of what is the right balance for includ-ing anisotropic effects and how much rigor should be applied to account for these spatial charge effects.

Anisotropy can be accounted for in both polarizable and non-polarizable force fields, however, the difference is made by to what extend the force field can describe the dynamic character of the charge. In principle an anisotropic object would be stretched or deformed by a force to varying degrees depend on its orientation. In a

(15)

Figure 4: Representative molecules that have substantial polarization anisotropies. Figure courtesy to ref. [81]

static manner one could account for this deformation before the actual deformation occurs by introducing extra point charges that are usually fixed in magnitude and position based on the type of interaction they should represent.

For example, anisotropy can be included by the addition of off-centered point charges to represent lone pairs. Early on the addition of off-centered charges was already found to be both effective in improving structures and interaction energies as well as an easy to implement method.[82] This approach of adding virtual charge interaction sites has subsequently been applied to model electron lone pairs.[12] This simple solution is able to mimic the complicated directionality of H-bonding interactions reasonably well. [11] In addition, the same remedy can also be ap-plied to describe σ-holes. The reason being that the σ-hole is actually a region with positive electrostatic potential on halogen atoms, which has the possibility to interact with a lone pair on a hetero-atom, thereby forming the so-called halogen bond. In the classical fixed-point charge model, the halogen atoms have a spheri-cal negative electrostatic potential, but by simply attaching an off-centered positive charge to the halogen atom an extra interaction site is created.[83, 84] Similarly, off-centered charges can also be used to account for π-bonding by attaching two negative charges to each heavy atom, although this approach has found to be com-putationally inefficient.[84]

Still, one might wonder, with how much confidence is the position of such a vir-tual site determined? Some methods used ab initio derived dipole and quadruple moments to fit optimal positions of the charges.[85] However, even after refine-ments this method yields a large number of possible locations for the placement of these charges. Other methods are exclusively applicable for specific types of bonding, such as the σ-hole and halogen-bonds, limiting the transferability of the efforts.[86,83] Or in a like-wise manner the approaches are tailored towards indi-vidual elements, such as sulfur.[87] Another path is to use experimental data in conjunction to parameterize the positions of the off-centered charges.[87]

Adding off-centered charges has shown to be an effective, albeit ad hoc, solution for describing different interaction types, including H-bonding. Still, even if the QM/DFT community has not come to a consensus about the complex nature of hy-drogen bonds,[88] how can off-centered charges properly describe its mechanism? One aspect where it becomes evident that adding off-centered point charges are

(16)

no panacea is illustrated by how poorly this method reproduces the electrostatic po-tential (ESP) of a molecule.[80] By this reasoning, a more recent avenue aimed to determine off center point charges by optimizing the recreation of the QM ESP with additional point charges.[89]2 3 Their approach shows promise, as this method was able to reduce for a set of small molecules the atom’s mean ESP error by 65.8%. Despite this, the authors note that they cannot say with certainty that their method is appropriate for all possible organic molecules, even though it has already been applied successfully to over 100 small organic molecules. On the other hand, a minimal distributed charge model has also been developed to determine a minimal number of off-centered point charges to approach the reference ESP.[92] Overall, there is no clear answer to whether more or less off-centered charges improve the anisotropic charge description.

Nevertheless, if the static approach is not sufficient enough to model anisotropy, how do we account for the anisotropic charge deformation? Polarizable force fields are more flexible than their non-polarizable counter-parts and are therefore able to address the anisotropy problem dynamically. For example, in the Drude Oscillator model the atomic polarizability, α, of a given atom is given by equation 2.5.[68] Here is qDthe charge on a Drude particle and kD the force constant for the spring that connects the Drude particle.

α = q 2 D kD

(6) In this model could, kD, for atoms that act as H bond acceptors actually be treated as a vector rather than a scalar, allowing for an anisotropic representation of po-larizabilities. [11] However, a disadvantage of the Drude particle method is the introduction of extra charges, which means more interactions to evaluate. And thus in practice, a constant value of kDis used for all atoms that results in qD determin-ing the polarizability of an individual atom.

Another approach to improve the dynamic description of anisotropic charge dis-tribution is the use of atomic multipoles, the model that is utilized by the polariz-able point dipole force field approaches.[8, 9, 10,13] Multipoles are series expan-sions that can represent arbitrary angular distributions. Atomic multipoles trun-cated at quadrupole have shown to be sufficient to model common chemical inter-actions, such as the previously mentioned σ-holes, lone pairs, and π-bonding.[30] However these higher-order moments also come at a certain cost. Furthermore it has also been shown, in contrast to common assumptions, that by itself the multi-pole approach cannot completely account for all energetically important effects of atomic-level anisotropy.[50] See figure 5 for a comparison of the ESP based upon point-charges, multipoles and QM approaches.

A more recent approach that aimed for the reproduction of ab initio polarizabil-ity anisotropy is the work of Wang et al. in which a set of atomic polarizabilpolarizabil-ity parameters for a new polarizable Gaussian model (pGM) has been developed.[81] Instead of fitting the ESP, this method takes advantage of the ability of the Gaussian model to screen all short-range electrostatic interactions and subsequently fit the ab initio molecular polarizability tensors directly. The authors suggest that by using 2 This QM ESP fit is very similar to an slightly older approaches, [90,91] although difference is made

in how both methods treat symmetry. The new method maintains the symmetry around the atom and bonding environment when placing the off center point charges.

(17)

knitting the reviews and perspectives together 17

Figure 5: Comparison of electrostatic anisotropy illustrated by the ESP of bromobenzene with (i) off-centered charges and (ii) atomic multipoles, and (iii) the reference QM potential. Figure courtesy to ref. [13]

Gaussian functions to seamlessly treat multipoles and electron penetration effects, a promising polarization framework has been developed. The main reason being that the pGM will most likely improve the stability in charge fitting, energy, and force calculations and the accuracy of multi-body polarization.

3 knitting the reviews and perspectives together

In the following chapter we will discuss the content and outlooks of 12 reviews and perspectives that addressed the developments of polarizable force fields over the past 20 years.

3.1 Polarizability: A smoking gun?

The first review in our list already concluded that without a doubt polarizable force fields provide a far more superior physical description compared to their fixed-charged counterparts.[5] The future seemed bright as the authors, Halgren and Damm, noticed that after twenty years the field finally started to gain traction within the (bio)molecular simulations community. Despite the technical limitations of the 2000’s, there was still no real proof or reason to believe that including polariz-ability would improve in accuracy, for example, in the calculation of ligand-receptor binding affinities, i.e., there was as of yet no "smoking gun".[5]

Six years later in 2007 Warshel et al. made special efforts to review a wide range of applications where the inclusion of polarizability is important.[6] Unfortunately, the trial’s long-awaited smoking gun still failed to surface, as it turned out that the inclusion of induced dipoles made no major difference in binding calculations of

(18)

neutral molecules.[93] At the time, the convergence problems were still more prob-lematic than the errors associated with the implicit inclusion of the induced dipoles in the parametrization procedure. Another highlight in the review described the procedure of a systematic (polarizable and non-polarizable) force field calibration by calculations on cation solvation energies with water and cation binding sites of proteins.[6] Interestingly, the polarizable models performed better only when the system was moved from water to other environments and even then only in the case of divalent ions (when dealing with ions that are in contact with water).

What should have been noticed by now, is that a general biomolecular force field was still far far away. For example, the study above about the cation solvation en-ergies had to carry out a complete empirical parameterization procedure on the system of interest. Nevertheless, theses early (not especially convincing) studies seem like a small steps, but the subtle differences between polarizable and non-polarizable force fields stimulated a more widespread realization of the importance to account for polarization effects. To wrap up this review, Warshel et al. men-tioned in their outlook that consistent quantum mechanical studies with QM/MM inclusion of the rest of the environment should be extremely useful in further force field development, as such, the effect of the induced dipoles can be separated from charge-transfer effects.[6]

3.2 New branches of electronic polarization

In 2009 already Cieplak et al. touched exactly on this latter point, as the review focused for a major part on the progress of how electronic polarization effects are in-corporated into force fields.[7] At this point in time, roughly five active branches of polarizable force fields methods were being in development: the fluctuating charge (FQ), Drude oscillator, induced point dipoles, electronic polarization via quantum mechanical treatment (QM) or mixed QM/MM,[94,95,96] and polarization treat-ment using a continuum solvent.[97] Most efforts had been related to explicitly induced dipoles and fewer to Drude oscillator and FQ approaches. In general, the parameterization of all these aforementioned approaches were considered to be a much slower process compared with additive counterparts. Therefore, the applica-tion of polarizable force fields was sparse and became even less due to the extra associated computational costs.[7] The reason for the upstart of these many differ-ent branches was attributed by Cieplak et al. to the fact that the force fields of ten years ago were "overdue for much-needed advances". However, the question of what approach to include polarizability will become the method of choice or is most suited for biomolecular simulations remained to be seen.

From 2000 till 2010 the much needed development occurred and molecular force fields were about to approach a generational transition, moving away from well-established and well-tuned, but physically less sound, non-additive point charge models towards more realistic and expensive polarizable models. Ponder et al. wondered if the new polarizable force field parameterizations actually reached a new level of predictive power over their non-polarizable predecessors.[8] In their review they addressed the validation of the AMOEBA force field, which was a leading example of the next-gen force fields at at the moment of the publication.[8, 98,99,100,101]

(19)

AMOEBA focused on the description of intermolecular interactions, which is es-pecially important for protein-ligand binding predictions. The result was a much improved reproduction of structural and thermodynamic properties. For example, in study that investigated the binding free energies of a series of benzamidine-like inhibitors to trypsin with both a polarizable and non-polarizable potential.[102,103] The study revealed that the polarization between water and benzamidine was re-sponsible for 4.5 kcal/mol out of the total 45.8 kcal/mol hydration free energy. On the other hand, the polarization between trypsin-in-water and benzamidine lowered the binding energy by 22.4 kcal/mol. These diverging results can be explained as follows: when polarization effects are included in the description of the medium, the protein becomes responsive and can screen the permanent electrostatic interac-tions, thereby weakening the attraction between benzamidine and trypsin.

The results from the protein-ligand binding studies with inclusion of polarizabil-ity looked very promising, although the hydration free energies (HFE) of the drug-like compounds in the SAMPL 2009 test revealed much larger errors.[104] The main source of these errors in the HFE were attributed to uncertainty in the experimental data as well as to problems with halogenated molecules and nitro compounds.4

Point dipole and Drude oscillator, both use point charge models to account for the permanent electrostatics, which in essence limits the physical character of the method. Real atoms are anisotropic while the point charge models are intrinsically isotropic.[30] Lone pairs, π-clouds and σ-holes are some examples of anisotropic in-teraction sites that are caused by specific electron distributions in molecules. These effects can be incorporated within point charge models by adding extra point charges, but require extra computation or extra attention during the development and pa-rameterization process.[105] In comparison„ multipole models naturally capture any non-spherical contribution of the atomic charge density, due to the inclusion of monopoles (charges) and higher order terms such as dipoles and quadrupoles. Still, Ponder et al. acknowledged that further fine-tuning would be necessary to accu-rately describe dynamical properties that are not sampled at ambient conditions as well as in the description of aromatic interactions (paradoxically interactions involv-ing anisotropy) as well as for solvent models.

3.3 Descriptions of electrostatics

Regarding the latter obstacle involving solvation free energies, the description of solvation and electrostatics go hand in hand when studying biomolecular folding, binding, enzyme catalysis and dynamics. Consequently, understanding the rela-tion between these two phenomena is ought to be important for the development of polarizable force fields. The in 2012 published review by Ren and co-workers discussed the advances made in the solvation of biomolecules with a sharp aim at computational biophysics,[9] and thus including those with interest in the ad-vancement of polarizable force fields. By now several groups have performed MD simulations using polarizable force fields to study ion behavior or to determine ion solvation energies employing the three different approaches for including polariz-4 The atomic polarizability values in the AMOEBA force field were taken from original work done by Thole, which did not include halogen atoms.[74] The nitro compounds were challenging because of large

bond length changes between the gas and liquid phase, as well as intricate “push-pull” polarization that was not captured in their “simple” polarization model.

(20)

ability effects mentioned above.[100,106,107,108] Despite all the efforts, modeling explicit ions remained a considerable challenge due to the subtle nature and com-plex dynamics of ions.[9]

Without regard to the abstinence5

in the development in classical polarizable force field methods for ions, the development of advanced classical electrostatic model beyond simple polarization was a great area of interest in the early 2010’s. In addition to polarization effect more quantum effect became the topic of interest, for example, local charge-transfer (CT) and penetration effects displayed an important role for short-range molecular interactions in water, aromatics and high-valence ions.[109,45,110,111] Inclusion of empirical, additive terms for CT and penetration effects have been shown to be rather effective, due to their short-range nature.[112, 109, 110, 113] Mainly, because these interactions can be treated with a local cut-offs, resulting in negligible additional computational cost relative to polarizable electrostatics computed with particle-mesh Ewald summation. In light of these improvements, Ren et al. (almost just as how Warshel et al. did) made the prognosis that the future direction would start to focus more on applying ab initio treatments and polarizable force fields or hybrid QM/MM approaches.[9]

3.4 Solvation and polarization

One year later in 2013 Cisneros et al. reviewed methods that accurately describe electrostatics in classical biomolecular simulations in explicit solvents. [10] Apart from the computational methods that have been developed to deal with long-range nature of the electrostatic interaction, the ways of representing the molecular elec-tronic charge cloud beyond the fixed point-charge representation were discussed extensively. Note, I focus solely on the aspect of including polarizability, but Cis-neros et al. also showed examples where the molecular electronic clouds can be accurately described by the inclusion of higher-order multipoles and or continuous electrostatic functions.[10,114] Again, the authors state that in spite of the presence of various models for treating polarization effects, the progress toward the develop-ment and applications of general-purpose polarizable force fields was still limited. The limitations were once more attributed to concerns related to computational speed and the lack of understanding of the importance of polarization effects. An-other facet blocking the way towards a widely applicable polarizable force fields are found in the challenges associated to other representations in the potential energy function, such as the van der Waals interactions and the short-range valence term.

At the time, most polarizable force fields were derived from both gas-phase ab ini-tio data and experimental properties for parametrizaini-tion, although to a differing ex-tent. The basis of deriving more accurate representations of the molecular electronic clouds depend (often) upon the fitting the electrostatic property of interest (usually the atomic charge) to a molecular electrostatic potential (MEP) obtained from ab initio, density functional, or semi-empirical wave function methods. However, the individual contribution of atoms in a molecule to the polarizability is not physically observable; it is only the molecular polarizability that can be measured.[30] Still,

po-5 Judging from publication dates of the articles mentioned above (ranging between 2003 and 2007) and the year of publication of the review article itself, the progress in describing ions seemed to be relatively slow.

(21)

larization may be described entirely in terms of classical electrostatics, and as such, it should be distinguished from the dispersion interaction that arises from instanta-neous fluctuations of molecular charge distributions, which is entirely quantum in origin.[10] Even though the overall molecular polarizabilities are recovered exactly in the ab initio approach, determining such atomic polarizabilities from quantum mechanics is not trivial. On the other hand fixed charge force field simulation tech-nology has become rather mature over the years, therefore leaving much room for development in the parameterization process of polarizable force fields. Cisneros et al. expected therefore to see growth of polarizable force fields along with their in-creasing application to unconventional molecular systems, mainly where traditional force fields have been challenged.

3.5 The rise of new challenges

In 2015 Baker made an overview on the systems to which polarizable force fields have been applied as well as the pro’s and cons that arise due to the inclusion of polarizability effects.[11] By the time of this publication a shift in the application and coverage of polarizable force fields has occurred, as parameters for fully opti-mized polarizable force fields were being published along with associated code im-plemented in standard simulation software. Together cover AMBERff02, AMOEBA, CHARMM Drude and CHEQ all major classes of biomolecules6

(either with com-plete or partial coverage),[72, 99, 8, 51,115,116] but Baker notes not all teething problems have been eliminated and the relative slow sampling rate employing po-larizable force fields remained an issue. For future directions in the development of polarizable force fields, Baker categorized issues that need more attention into three groups: parameterization, sampling and protein-ligand binding.

One of the challenges for (both polarizable and non-polarizable) force field devel-opers should be the identification of any weaknesses in parameter sets. Automated schemes for the optimization of force field parameters circumvent the need of man-ual involvement, thereby reducing human errors and resulting in systematic and reproducible parameters.[11, 117, 118] The following two studies show successful applications of such automated parameterization protocols. First, the force balance method for obtaining a single parameter set has shown to be able to optimize pa-rameters for a simplified AMOEBA water models.[119,120] Secondly, QM target data has been used to parameterize small molecules in an automated manner for both polarizable and non polarizable force fields.[121]

The second challenge addresses the slow reputation of polarizable force fields. In absolute terms this is neither true nor false, but compared to the non-polarizable counterparts this statement is true. Of course this issue also depends on the research question, as to whether a more accurate description is more important or adequate sampling. Enhanced sampling methods and Hamiltonian replica exchange between simulations of two different force fields are proposed by Baker as a possible resolve for the sampling problem.[11,122,123]

The last challenge of protein-ligand binding is actually a combination of both the previous issues, as the challenge come from both the accuracy of the force fields and the quantity of sampling. The free energy estimates will be poor if the two 6 _{Major classes: Proteins, Nucleic Acids, Lipids, Carbohydrates}

(22)

previous issues are not accounted for. At the time of publication, polarizable force fields were not being applied large scale free energy calculations due to the lack of parameters for a diverse range of small molecules. However, the AMOEBA method (where a set for small molecules was already available) displayed the capability to produce reliable estimates of experimental binding free energies.[124] As a final remark Baker states that if automated methods for generating force field parame-ters become successfully applicable to polarizable force fields, it is likely that the accuracy of calculating binding free energies will be further enhanced.

3.6 Parameterization or polarization?

In 2016 Lemkul et al. described the latest developments in the Drude force field parametrization and application.[12] The main focus is placed on the Drude-2013 polarizable force field for proteins, DNA, lipids and carbohydrates. A set of small organic compounds have been used as target data for the parameterization process. QM calculations and condensed-phase experimental data were targeted during the iterative refinement protocol that previously was developed for the parameteriza-tion of the additive CHARMM force field, which yielded a transferable parameter set across molecules.[125] To achieve such a robust model and transferability, the parameterization protocol involved flexible tuning of Thole screening factors on a per-atom basis, the use of atom-pair-specific LJ parameters, and scaling of gas-phase polarizabilities based on the nature of the model compounds.[12]

Shortly going back to Baker, as he noted that compared to the non-additive counter parts polarizable force fields tend to include more individual parameters per atom, which is of course inherent to the more complex additive potential func-tion. [11] Nevertheless, granted that polarizable force fields represent a more phys-ically sound potential, atoms should be more adaptive to a diverse set of surround-ings. The actual implication should therefore be that less distinct atom types are necessary for a polarizable force field than the non-additive counterparts. As of 2015the CHARMM Drude polarizable force field required as many distinct atom types as the non-polarizable analogue, which indicates towards a point of future improvements.[126]

By 2011 the parameterization of small molecules and ions for the CHARMM Drude model seemed as good as finished.[127] However, the full release of the Drude-2013 force field, was delayed (as the name speaks for itself) by 2 years. The extra time required for the fully applicable polarizable protein force field to become available was caused by challenges associated with non-additive effects that arose when going from small molecules to larger, polymeric macro-molecules.[12] For instance, the interactions between positive and negative charged residues showed small imbalances, as these were not treated at the small molecule level. Conse-quently, this lack in description of multiple charged residues caused nonphysical interactions that could result into a polarization catastrophe.

(23)

3.7 Enough response: how far away?

The amount of reviews specifically aimed at the development of polarizable force fields remained relatively sparse, until a burst of three consecutive publications appeared in 2019. Before I will discuss these three recent literature studies in further detail, it might be a good idea to wrap up the conclusions, challenges and issues made thus far. After this pre-evaluation we hopefully find breakthroughs to the problems in the publications of 2019.

In the early 2000 a general-purpose polarizable force fields was still far from re-ality and the development was rather limited compared to the fixed-counterparts. Mainly, because the parameterization of polarizable force fields was much slower with considerable more computational costs. Another issue associated with the pa-rameterization process is the non trivial way atomic polarizabilities are obtained from quantum mechanical data with the addition of the need for less distinct atom types than the non-additive counterparts. As the development of methods to ac-count for polarizability effects started to move on apace many other issues arose, such as how to handle anisotropic interactions, charge penetration effects, dynamic ionic systems and the limitations of other representations in the potential energy function. Finally, challenges associated with non-additive effects arise when go-ing from small molecules to larger, polymeric macro-molecules and/with multiple charged residues.

Apart from the often bleak but hopeful outlooks, considerable progress was made towards a generally applicable polarizable force field for biomolecular simulations. As automated schemes for the optimization of force field parameters as well as enhanced sampling methods and other algorithmic improvements began to soften the computational burden of development and application of the more expensive polarizable models. And at last, complete polarizable models appeared of which the CHARMM Drude and the AMOEBA force fields are two examples that cover a large array of system types and also have been through various iterations of testing and improvements.

3.8 The last perspectives

By 2016, however, considerable efforts are still required to address all problems associated with the inclusion of polarizability in one package and the amount of up-to-date software is still limited compared to the number of non-polarizable force fields. What has become clear, at this point, is that the fixed-point charge treatment has come out of due and polarizable force fields start to ascend in the mainstream of the biomolecular simulations. Obviously, there are still pieces of the polarizable force field puzzle that have been omitted so far. Therefore, we will quickly continue with our final course of reviews and afterwards I hope to construct a clear basis of how (some of) the current issues have been resolved and what challenges still need to be tackled in the coming years.

3.8.1 Jing et al.

The first review of 2019 was published by Jing et al. and as always, the review provides an update on the recent advances and applications, as well as future

(24)

di-knitting the reviews and perspectives together 24

rections of polarizable force fields in biomolecular simulations.[13] Advances in the utilization of computational resources, especially the application of the graphical processing unit (GPU), have made microsecond MD simulation using polarizable force fields accessible.[128,129,130] The accuracy and coverage of polarizable force fields have also improved in recent years, see figure6for the polarizable force fields easily accessible and applicable for biomolecular simulations. The AMOEBA force field has recently been extended to DNA and RNA and has also become more flexi-ble with sampling in different environments.[131] The CHARMM Drude force field has been refined for DNA and recently also includes carbohydrates and halogenated molecule.[84,132]

Figure 6: Schematic of polarizable force fields for biomolecules and available software. Fig-ure courtesy to ref. [13]

With regard to the development in systematic and automatic parameterization approaches, Jing et al. expect that these protocols will profit from the continuous advancement in QM methods and machine learning (ML) approaches that become more frequently adopted in chemistry.[13] According to Jing et al. the combina-tion of polarizable force fields with enhanced sampling methods has not yet been explored fully, but starts to gain more interest. This field of methods that poten-tially can extend the simulation timescales of large biomolecular complexes, is il-lustrated with examples of orthogonal space sampling, Markov state models, and mile-stoning.[133,134,135]7Still, the authors note that further improvements of the underlying physical models are necessary, particularly for short-range interactions such as charge penetration (CP) and charge transfer (CT).[44,42]

3.8.2 Bedrov et al.

In May of 2019 Bedrov et al. published a review (with over 435 references herein) in which they discuss the different polarization models in MD simulations of ionic materials, a subject that is important for many applications in chemistry, biology, as well as energy storage and conversion.[14] Subsequently, the review provides the pros and cons of the different polarization treatments for highly ionic materials. Furthermore, the authors address and compare the methods and strategies for the 7 Further strengthened by the recent publication of Celerse et al. in which the steered molecular dynamics methodology has been implemented in the framework of the massively parallel Tinker-HP software allowing for both long polarizable and non-polarizable MD simulations of large protein.[136_]