Recent Advances in Polarizable Force Fields for Macromolecules: Microsecond Simulations of Proteins Using the Classical Drude Oscillator Model.

In this Perspective, we summarize recent efforts to include the explicit treatment of induced electronic polarization in biomolecular force fields. Methods used to treat polarizability, including the induced dipole, fluctuating charge, and classical Drude oscillator models, are presented, including recent advances in force fields using those methods. This is followed by recent results obtained with the Drude model, including microsecond molecular dynamics (MD) simulations of multiple proteins in explicit solvent. Results show significant variability of backbone and side-chain dipole moments as a function of environment, including significant changes during individual simulations. Dipole moments of water in the vicinity of the proteins reveal small but systematic changes, with the direction of the changes dependent on the environment. Analyses of the full proteins show that the polarizable Drude model leads to larger values of the dielectric constant of the protein interior, especially in the case of hydrophobic regions. These results indicate that the inclusion of explicit electronic polarizability leads to significant differences in the physical forces affecting the structure and dynamics of proteins, which can be investigated in a computationally tractable fashion in the context of the Drude model.

Molecular modeling and simulations are important tools for exploring the structure, dynamics, and function of complex chemical and biochemical systems in condensed phases. Underlying these methods are force fields, which represent a computationally tractable approximation to the Born–Oppenheimer potential energy surface.[1] The functional form of the force fields (FFs) used to simulate macromolecular systems has remained largely unchanged during past decades, though the parameters in those FFs are continually being refined. In the widely used macromolecular FFs, the electrostatics term is described with Columbic interactions between fixed partial atomic point charges, referred to as an additive FF. This simplification is used largely to reduce the computational cost as well as to facilitate parametrization. In such additive fixed-charge FFs, electronic polarization effects are considered in a mean-field, average manner by empirically optimizing the partial atomic charges to overestimate dipole moments during the FF parametrization, thereby being representative of the condensed, typically aqueous, phase. While additive FFs have seen considerable use, the additive approximation significantly limits the accuracy of the method, for example, in treating highly polar versus hydrophobic environments. To achieve a more accurate description of the response of the charge distribution to variations in the surrounding electrostatic field, the explicit inclusion of electronic polarizability in the model is essential.

There are a variety of methods to include explicit electronic polarizability in FFs. An obvious approach is to treat the electronic response quantum mechanically. While this has been performed to treat local regions of macromolecular systems in the context of QM/MM methods[2−4] and preliminary FFs that include the treatment of electrostatics using semiempirical methods have been presented (e.g., XPOL),[5] these are computationally demanding approaches. The alternative is a more simple, computationally tractable extension of the potential energy function to mimic the quantum mechanical response of the electronic distribution to changes in the surrounding electric field. One way to perform this is to assign an induced point dipole at each atomic site.[2,6−11] Typically, the induced dipoles are determined by a self-consistent field (SCF) iterative procedure using isotropic atomic polarizabilities, following which the charge–charge, charge–dipole, and dipole–dipole interactions are computed to yield the electrostatic energy of the system. The computational bottleneck of induced dipole models is usually the SCF evaluation of mutual polarization. Recently, Wang et al. proposed an iAMEOBA model for water in which a single SCF step is performed with induced dipoles initially set to zero (ie. the induced dipoles are only determined by the electric field associated with the fixed charges or multipoles).[12] However, this approach neglects mutual polarization between the induced dipoles in the system, though it has been indicated that this contribution may be modeled in an average manner during the parametrization, an approach analogous to that performed with current additive FFs. Using this approach, the iAMEOBA water model was shown to outperform the original AMEOBA model in reproducing experimental measurements, although this is likely due to a parametrization algorithm that directly targets condensed-phase properties.[12,13] While the iAMEOBA approach is likely satisfactory for homogeneous systems, there are legitimate concerns that such an approximation will represent a significant compromise in accuracy in more polar, heterogeneous systems such as ion solvation.

Another approach to treat polarizability is the fluctuating charge model,[14−18] which treats partial atomic charges as dynamical variables based on the electronegativity at each atomic site. Charges are propagated according to the principle of electronegativity equalization,[19] and charge conservation constraints are assumed. Usually, charges are only allowed to redistribute within molecules to avoid unphysical intermolecular charge transfer, and an extended Lagrangian approach is adopted to propagate the charge variables during molecular dynamics (MD) simulations. One limitation of the model is the inability to describe out-of-plane polarization for a planar system such as benzene, although it is possible to add auxiliary out-of-plane charge sites to overcome this limitation. The CHeq FF is based on fluctuating charges, with parameters presented for proteins,[17,18] lipids,[20] and select carbohydrates.[21] In addition, CHeq has been applied to ligand binding to lysozyme,[22] ion solvation,[23] and lipid bilayer permeability.[24]

A third approach for the explicit treatment of polarization, which has been developed in our groups, is the classical Drude oscillator model.[25,26] The Drude model has also been referred to as the Shell[27] or the Charge-On-Spring model.[28] Briefly, a charged auxiliary particle (the Drude oscillator or particle) is attached to the atomic core of its parent atom via a harmonic spring with force constant kD. The displacement of the Drude particle relative to its parent atom in the presence of an electric field is given bywhere qD is the charge of the Drude particle. The induced dipole iswhich is equivalent to an atomic polarizability of α = qD2/kD.

In our Drude FF, the Drude particles are only associated with non-hydrogen atoms, which has been shown to reproduce molecular polarizabilities while maintaining computational efficiency.[29] For the sake of simplicity, a fixed value of 1000 kcal/mol/Å2 is used for the restoring force constant, kD, such that the charge qD is the parameter that governs the magnitude of α for a given non-hydrogen atom. While eq 2 describe isotropic atomic polarizability, anisotropy can be incorporated in the Drude model by considering (1/kD) as a tensor, that is, attributing an inhomogeneous spring constant based on a local reference frame for the virtual spring between the atom and its Drude particle. Such a treatment is implemented in the current Drude FF but only applied to hydrogen bond acceptors such as the peptide backbone carbonyl oxygen atoms. This and virtual particles representative of lone pairs are included to improve the treatment of nonbonded interactions as a function of orientation involving hydrogen bond acceptors.

An essential aspect of the Drude FF is that the dipole–dipole interactions between atoms involved in bond or valence angles (1–2 or 1–3 interactions) are explicitly included. However, for 1–2 and 1–3 interactions, Coulomb’s Law fails due to the short spatial separation. This problem may be overcome by applying an electrostatic shielding factor S to the electrostatic interactions as proposed by Thole[30]where r is the distance between the charge site i and j and αi, αj and ai, aj are the polarizabilities and Thole factors, respectively. The atom-based Thole factors introduced in the Drude model provide fine-tuning of near-field electrostatics, yielding improvements in the treatment of the orientation of molecular polarizabilities. Thus, atomic polarizabilities and Thole factors are FF parameters that, in addition to partial charges, need be optimized during the parametrization of the Drude polarizable FF. While the remainder of the potential energy function, including the bonded terms and the van der Waals (vdW) interaction, is identical to the CHARMM additive FF, the associated parameters also need to be optimized. This is due to the coupling of all terms in the energy function, such that changes in the electrostatic model require reoptimization of the vdW and internal bonded parameters to produce a carefully balanced empirical FF.

Parametrization started with the water model, initially the SWM4-DP model[25] followed by the SWM4-NDP model,[31] with the latter becoming hereafter the default model for all subsequent developments of the Drude FF. Parameters were then derived for model compounds that represent components of biomolecules. While the parametrization of the Drude FF shares similar target data with the CHARMM additive FF, additional target data such as the molecular polarizability were included to facilitate optimization of the electrostatic parameters. In general, the more sophisticated physical model and associated additional FF parameters in the Drude FF allow better reproduction of these target data including both gas-phase QM data and condensed-phase experimental data, an example being the treatment of polyalcohols.[32] Also, during the optimization, the atomic polarizabilities were scaled to reproduce pure solvent experimental dielectric constants,[33−35] yielding scaling factors ranging from 0.6 to 1.0, such that the Drude polarizabilities are equal to or lower than gas-phase experimental or QM values. Due to the nonadditive nature of the Drude FF, particular care was taken when transferring the parameters from model compounds to biomolecules, and often, target data for relatively large model compounds need to be included for further optimization. For example, the gas-phase relative energies between αR, PPII, and C5 conformations of the (Ala)5 peptide and interactions of water with the alanine dipeptide are included as additional target data for optimizing the parameters for the polypeptide backbone.

An essential feature of an empirical FF is computational efficiency allowing for simulations of macromolecules with an explicit solvent representation on time scales of 100s of ns through microseconds. With the Drude model, this has been attained by maintaining the simplicity of the potential energy function described above as well as the implementation of an extended Lagrangian integrator with a dual thermostat, allowing for computationally efficient MD simulations while maintaining approximate SCF treatment of the polarizable degrees of freedom.[26] Unlike extended Lagrangians used in Carr–Parrinello (CP) MD or fluctuating charge model simulations, where fictitious masses are attributed to the electronic degrees of freedom, the Drude particle may be treated as a real particle that moves in space, allowing it to be assigned a mass, typically 0.4 amu, taken from its parent atom. In the dual-thermostat algorithm, the relative motion of each Drude–nucleus oscillator pair is coupled to a low-temperature thermostat (typically 1 K), such that the electronic degrees of freedom approach the adiabatic SCF limit during the MD simulation. The remainder of the system, including the center-of-mass motion of the Drude–nucleus pairs, is thermostated to the target temperature of the simulation (e.g., room temperature). To avoid polarization catastrophe, the distance d between the Drude particle and its parent atom is limited to typically 0.2 Å by imposing a hard wall constraint to d during the MD simulation.[36] This combination of a simple form of the potential energy function and an extended Lagrangian integrator implemented in an efficient parallelizable code, such as CHARMM, NAMD,[37] and ChemShell QM, has allowed for MD simulations on the order of 100s of ns on proteins, lipids, and DNA.[38,39] With NAMD, as compared to the additive FF, the computational overhead with the Drude model is approximately two-fold, which, when using a 1 versus 2 fs integration time step, yields an overall four-fold computational increase with the Drude model. In the remainder of this Perspective, we will present results from fully polarizable microsecond MD simulations of proteins using the Drude FF along with comparisons with the CHARMM36 (C36) additive protein FF.

Presented in Figure 1 are protein Cα (root-mean square deviations) RMSDs from 1 μs MD simulations of ubiquitin (1UBQ) and cold shock protein A (CspA, 1MJC) using the Drude and the C36 FFs. MD simulations were carried out in the NPT ensemble with NAMD, with simulation details provided in the Supporting Information. For the Drude simulations, a shorter time step of 1 fs was used versus 2 fs for the additive FF, which is due to the high-frequency motion of the Drude particles related to their small masses. Drude simulations were also performed on crambin (1EJG), the tight junction regulatory protein (3VQF), circular permutant of ribosomal proteion S6 (3ZZP), DNA methyltransferase associated protein (4IEJ), and protein G B3 domain (1P7E), with RMSD plotted in Figure S1 of the Supporting Information. The RMSDs of both ubiquitin and CspA were stable, and those from the Drude and C36 simulations were comparable. No systematic drift was observed, indicating the ability of Drude FF to maintain the protein folded structures on the microsecond time scale.

Figure 1

RMSD plots of 1 μs simulations of ubiquitin and cold-shock protein A. RMSDs were computed for Ca atoms in all residues, and results are presented as running 10 ns averages.

With polarizable FFs, the charge distribution of functional groups (i.e., dipole moments) is dominated by the electrostatic environment, which is fundamentally different from additive FFs where the dipoles only vary due to changes in the intramolecular geometry of the individual functional groups. This can be illustrated by examining dipole moments of the peptide backbone and protein side chains during the protein simulations. The backbone peptide bond group is defined as the C and O atoms of residue i and the N, H, Cα, Hα atoms of residue i + 1, yielding neutral charge in both FFs. The Drude simulations yield a mean value of 4.74 ± 0.31 D for peptide bonds in the helices and 5.14 ± 0.30 D for those in sheets, where the errors are the RMS fluctuations. The C36 simulations yield much narrower distributions, with mean values and RMS fluctuations of 3.82 ± 0.11 D for helical peptide bonds and 3.71 ± 0.09 D for sheet peptide bonds. The larger magnitude dipoles observed with the Drude model is notable because the parametrization of additive FFs typically involves systematic overestimation of dipole moments relative to that of the gas phase in an attempt to capture the effect of the environment in an average mean-field way.[40] While the enhanced favorable interactions associated with the larger dipoles are partly canceled by the positive self-energy for polarizing the atoms, the implication is that an explicit treatment of induced polarization captures local variations in the environment that are not realistically accounted for with an additive mean-field approximation. Similar results have been obtained for the dipole moments of nucleic acid bases in duplex DNA.[39] The ability of peptide bond dipoles to dynamically respond to the external electric field is likely to be important for the dynamics of protein and peptides. Recently, the Drude FF was found to reproduce the cooperativity of helix formation in the acetyl-(AAQAA)3-NH2 peptide, and such folding cooperativity was shown to be associated with enhanced dipole moments of the peptide backbone upon helix formation.[41]

The dipole moments of amino acid side chains are also typically larger in the Drude model as compared to those for the additive FF, although this varies based on residue type, as illustrated in Figures S2 and S3 of the Supporting Information. Exceptions occur with nonpolar residues such as Phe, Val, and Ile, with their side-chain dipole magnitudes being similar between the two models. Notably, the polarizable FF shows wide variability between residues of a given type during the MD simulations as well as in individual residues themselves, as previously observed for Trp residues in lysozyme.[38] An example is shown in Figure 2, where dipole moments for four Gln residues in ubiquitin are plotted as a function of time. Transitions in dipole moments on a time scale of tens of nanoseconds are observed for the buried Gln41 residue, correlated with the side chain rotating about its χ1 dihedral angle (Figure 2C). Similar rotation also occurs in the additive C36 simulation (Figure 2D); however, due to the fixed charge nature, no significant changes in the dipole moment are observed. Images of the two environments of Gln41 sampled in the Drude simulation (Figure 2E and F) show that its side chain carries a smaller dipole moment when pointing toward a helix and larger dipole moment when it forms a hydrogen bond with the carbonyl oxygen of Pro38 in a loop region. In contrast, the dipole moment of surface residue Gln2 is stable throughout both the Drude and C36 simulations. However, the Drude model yields a much larger absolute value (5.8 D) than C36 (4.2 D).

Figure 2

Side-chain dipole moments of glutamine residues in ubiquitin during 1 μs MD simulations with the Drude (A) and the C36 (B) FFs. The χ1 dihedral angles for Gln41 are plotted for the Drude (C) and the C36 (D) simulations, and images of Gln41 (Licorice) in ubiquitin (cartoon) are shown from the 400 (E) and 500 ns (F) time points. Residues containing atoms within 3 Å of the Gln41 side chain are shown in lines. Dipole moments are running averages over 10 ns.

The distribution of water dipole moments during the Drude simulations was also analyzed. For the Drude simulation, the SWM4-NDP water model has a dipole moment of 1.85 D in the gas phase and a much larger average value of 2.46 D in the bulk phase (black dashed lines in Figure 3) due to the mutual polarization. The additive TIP3P water model has a fixed dipole moment of 2.35 D. RMS fluctuation of the bulk water dipole moments in the polarizable FF is 0.16 D, similar to the value of 0.22 D[42] obtained from CPMD simulations using the Bader approach.[43] The distribution of water dipole moments in bulk is compared with those in the first solvation shell of selected residues in Figure 3, as defined as any water within 2.4 Å of any non-hydrogen atom in the protein. Larger dipole moments are observed for water in the proximity of negatively charged Asp and Glu side chains, while smaller values are obtained near positively charged Arg and Lys. Water close to the peptide backbone groups and the hydrophobic side chains also carries slightly smaller dipole moments compared to those of bulk water. Such a difference in the dipole distributions, though small as compared to water in the bulk phase, indicates that water can probe and respond to the complex electric environment of the protein in the polarizable FF. Ab initio MD simulations and QM calculations of ion solvation in water found that the average dipole moment of water in the first hydration shell of K+ is 0.22 D smaller than the average bulk value, and those around Cl– are 0.11 D smaller.[44,45] The importance of variations in water dipoles were shown in a recent study on DNA base flipping using the Drude DNA FF, where changes in water dipole moments in the first solvation shell around the flipping bases occur during the transition from the Watson–Crick base paired to flipped conformations.[46]

Figure 3

Probability density distribution of water dipole moments in bulk (black lines) and in the first solvation shell in the vicinity of (A) arginine and lysine side chain groups, (B) aspartic acid and glutamic acid side chain groups, (C) peptide backbone groups and polar (Gln, Thr, Asn, and Ser) side-chain groups, and (D) hydrophobic (Val, Ala, Phe, Trp, Tyr, Met, Ile, and Leu) side-chain groups.

We also examined the dipole moments and the molecular polarizabilities of the entire proteins, which are closely related to their dielectric properties (Table 1). The traces of protein polarizability tensors during the 1 μs Drude simulations for ubiquitin and CspA are shown in Figure S4 of the Supporting Information, showing them to be stable along the 1 μs MD trajectories, with fluctuation occurring on the nanosecond time scale. Dielectric constants obtained by modeling the proteins as spheres with radii based on the radius of gyration[47] are summarized in Table 1. The average εinf over the six proteins considered is 2.0 in the Drude model, which equals the commonly assumed value.[47,48] In contrast, some studies have suggested that a higher value should be used due to the aromatic groups in proteins or the higher density of proteins.[49,50] However, the present results indicate that this is not necessary. While εinf among the different proteins is similar, their variation will be important when studying properties dominated by the protein dielectric relaxation, for example, as when computing the reorganization energy for electron transfer in proteins using Marcus theory.

Table 1

Protein Dielectric Constants Computed from MD Simulationsa

	εinf	ε
		entire protein		protein interior		hydrophobic core
	Drude	Drude	C36	Drude	C36	Drude	C36
1UBQ	2.10 ± 0.00	10.3 ± 1.1	11.3 ± 1.1	2.6 ± 0.1	2.6 ± 0.0	2.1 ± 0.1	1.4 ± 0.0
1MJC	2.00 ± 0.00	10.8 ± 0.8	9.0 ± 0.5	3.6 ± 0.4	2.5 ± 0.1	4.6 ± 0.5	2.8 ± 0.2
1EJG	1.76 ± 0.00	2.5 ± 0.1	2.1 ± 0.0	2.7 ± 0.0	1.8 ± 0.1	2.2 ± 0.0	1.6 ± 0.1
3VQF	2.10 ± 0.00	7.0 ± 0.5	5.7 ± 0.4	3.1 ± 0.2	3.2 ± 0.2	3.7 ± 0.3	2.1 ± 0.1
3ZZP	1.91 ± 0.00	8.3 ± 0.3	5.8 ± 0.3	5.0 ± 0.2	3.1 ± 0.2	2.4 ± 0.2	1.6 ± 0.1
1P7E	1.91 ± 0.00	8.7 ± 0.4	5.2 ± 0.4	4.3 ± 0.3	1.8 ± 0.0	1.2 ± 0.0	1.0 ± 0.0

Optical dielectric constants εinf were computed from the protein polarizability according to the Clausius–Mosotti equation, and the molecular polarizability tensor for a given conformation was computed by applying an external electric field along the x, y, and z directions, relaxing the Drude particles and then evaluating the changes in dipole moments. Static dielectric constants of entire proteins, the protein interior (residues with a relative SASA less than 50%), and the hydrophobic core (residues with a relative SASA less than 20%) were computed from the Drude and the C36 simulations. Proteins were modeled as spheres whose radius r is given by the radius of gyration Rg multiplied by (5/3)1/2.[47]

Additional analysis involved estimates of the static dielectric constant ε of different regions of the proteins from Kirkwood g factors based on the fluctuation of protein dipole moments as presented by Simonson and Brooks.[47,51] The protein interior is defined using the relative solvent-accessible surface area (SASA) as a criterion, and the computation of ε’s is detailed in the Supporting Information. As shown in Table 1, dielectric constants vary across the different proteins. The small value of 2.5 for crambin is consistent with the fact that this protein is not soluble in water. For the other proteins, ε ranges from 7 to 12, much larger than the value of 2–5 from the measurements on dry protein powders[52] and consistent with previous calculations from MD simulations using additive FFs.[47] Notably, ε values of the protein interiors are much smaller than the ε of the entire protein, and the Drude model yields dielectric constants that are larger or similar to the C36 model. For the protein hydrophobic cores, defined as residues with a relative SASA less than 20%, the Drude dielectrics are systematically larger than those from C36. These results are consistent with ε values for the entire proteins being dominated by charged surface residues,[15] while for the buried residues in the core, additional electronic degrees of freedom of the charge distribution become important. Finally, we note that the dielectric constant is a macroscopic quantity, and estimating its value in microscopic regions of proteins, which are anisotropic and inhomogeneous, remains challenging.[53,54] Furthermore, ε is dependent on the actual protein site being considered and how the site is defined.[53,55] This analysis highlights the importance of using an explicit, fully polarizable model to obtain a detailed and accurate description of the response of proteins to their own electric fields as well as those of the surrounding environment.

In this Perspective, we summarized results on protein simulations using the Drude-2013 polarizable FF. The Drude model incorporates a more physically realistic treatment of electronic polarizability versus additive FFs and yet maintains a computational efficiency that is crucial for a FF, as illustrated here by microsecond simulations of globular proteins. The Drude polarizable protein model, along with currently available parameters for DNA,[39,56] dipalmitoylphosphatidylcholine,[57] and hexapyranoses,[58] as well as ongoing efforts to parametrize RNA and more lipid and carbohydrate molecules, will allow for studies of heterogeneous biomolecular systems using a fully polarizable FF.