# Fast Atomic Charge Calculation for Implementation into a Polarizable Force Field and Application to an Ion Channel Protein

Authors: Raiker Witter; Margit Möllhoff; Frank-Thomas Koch; Ulrich Sternberg
Journal: Journal of Chemistry
 (2015)

---

## Abstract

Polarization of atoms plays a substantial role in molecular interactions. Class I and II force fields mostly calculate with fixed atomic charges which can cause inadequate descriptions for highly charged molecules, for example, ion channels or metalloproteins. Changes in charge distributions can be included into molecular mechanics calculations by various methods. Here, we present a very fast computational quantum mechanical method, the Bond Polarization Theory (BPT). Atomic charges are obtained via a charge calculation method that depend on the 3D structure of the system in a similar way as atomic charges ofab initio calculations. Different methods of population analysis and charge calculation methods and their dependence on the basis set were investigated. A refined parameterization yielded excellent correlation of R
=
0.9967. The method was implemented in the force field COSMOS-NMR and applied to the histidine-tryptophan-complex of the transmembrane domain of the M2 protein channel of influenza A virus. Our calculations show that moderate changes of side chain torsion angle χ
1 and small variations of χ
2 of Trp-41 are necessary to switch from the inactivated into the activated state; and a rough two-side jump model of His-37 is supported for proton gating in accordance with a flipping mechanism.


---

## Body

## 1. Introduction

Substantial effort has been put in the development of polarizable force fields to understand the function of biomolecules in their realistic environment, that is, biomembranes or aqueous solutions. Standard Class I and II molecular mechanics force fields, used for chemical and biological applications, describe electrostatic interactions in terms of fixed, in most cases atom-centred, partial charges. However, polarization plays a substantial role in shifting charges within real molecular systems. The interaction energies of molecular systems are considerably influenced by polarization due to molecular rearrangements that cause charges to move in the region of interaction. The neglect of polarizations is one of the most severe drawbacks of force fields using fixed atomic charges; and many attempts have been made to develop polarizable force fields (for a review see, e.g., Cieplak et al. [1] and references therein).Partial atomic charges in molecules are not directly measurable and they are not observable in the context of quantum theory. They give a quantified description of the electron distribution in molecules in a simplified view. Atomic charges can only be determined indirectly by measuring strongly charge dependent observables, for example, the dipole moment, the electric field gradient, or the magnetic shielding. Charges induce correct tendencies and magnitudes of interactions utilized for the calculation of electrostatic interactions by applying Coulomb’s law. Therefore, the concept of partial atomic charges is extremely helpful in many fields of biology, chemistry, and physics.Widely used force fields (and collections) like AMBER [2, 3], CHARMM [4–6], GROMACS [7–10], GROMOS [11–13], and MMX [14] use fixed atomic charges which result in inadequate descriptions for some applications, especially in case of highly charged molecules. Efforts were made to improve existing force fields with applications tested by the developers, for example, including fluctuating charges in OPLS/PROSA [15], applying induced dipole models in AMBER [16], and introducing multipole expansions as well as additional potentials in AMOEBA [17] or utilizing Drude oscillator models in CHARMM [18].Methods, based on the combination of quantum mechanics (QM) and molecular mechanics (MM), received more attention due to the development of increasingly powerful computation options and fast data processing. Although QM/MM methods always bear the problem of truncation effects caused by the discontinuity at the borderline between the QM and MM treated section within a molecular system, QM/MM methods have been a valuable option to model chemical reactions (breaking and forming of chemical bonds as well as states far away from equilibrium) [19]. The results allow valuable insights into, for example, proton transfer reactions and transition states in biological interesting molecular systems due to the consideration of improved coupling potentials and a critical data analysis. However, quantitative comparisons to experimental and calculated data require careful analysis; see [20]. An excellent article regarding the capabilities and limitations to treat complex biomolecular processes was published by Riccardi et al. [21]. Emphasis was given to the prevention of misinterpretations when only a minimum QM region was defined; moreover, consistent treatment of the electrostatics and sufficient sampling are crucial for the success of the method.In our approach, partial atomic charges are calculated using the Bond Polarization Theory (BPT) [22–25]. Two parameters per bond type are required to calculate partial charges from bond orbitals. These two parameters linearly enter into the BPT formula and therefore can be determined from a least square solution of a system of equations, provided partial charges are known fromab initio calculations. These two parameters have a physical meaning: One parameter is the bond partial charge of a nonpolarized bond and the other parameter is the change of the charge due to polarization. The polarization energy for each bond is calculated within the BPT framework. BPT charge calculations are easily integrated into a force field by constructing a bond orbital for each bond defined in a molecule. The capability to describe polarization within a force field model potentially serves as key to a better understanding of structural based function of biological systems [23, 26]. The concept of polarization has been also applied to fast calculation of chemical shift tensors, for nuclear magnetic resonance spectroscopic investigations, and for direct structure elucidation using them [25, 27–33].The semiempirical treatment of charge distributions integrated into a force field improves the description of short- and long-range electrostatics within a large molecular system and this should bridge the gap between QM/MM methods and force field calculations with fixed charges. Previously, we developed a force field (COSMOS-NMR) [34] with polarisable BPT-derived charges to describe intermolecular interactions for most organic compounds: Molecules consisting of C, O, N, and H atoms [26]. Although the BPT is not able to describe chemical reactions, the COSMOS-NMR force field provides a novel approach to treat biomolecular systems inclusive polarization effects with consistent electrostatics and allows sufficient sampling due to the low computational effort; for benchmark simulations see Schneider et al. [35].The interaction of Valinomycin with potassium ions was studied earlier [26] and offers an explanation regarding the selectivity of the ion transporter protein for K+ versus Na+. The calculations yielded interaction energies of the cyclic depsipeptide complex that are in the order of the hydration energy; that is, Valinomycin is able to strip off the hydration sphere of a K+ ion by forming an interaction complex lower in energy compared to the ion in its hydration shell. The K+ selectivity of Valinomycin in nature provides qualitative proof of our quantitative model. The calculations pointed to polarization effects that are mainly attributed to carbonyls of the backbone of Valinomycin, which matches experimental results in related systems [36].In this work, one of our objectives is the extension of our semiempirical charge calculation method to F, Si, P, S, and Cl containing compounds by a reparameterization. Also, some criticism of the BPT charge calculation was voiced due to the choice of a minimum basis set (STO-3G) and population analysis (Mulliken) to obtain reference atomic charges for the parameterization. Thus, several alternative parameterizations of the BPT charge calculation were tested for best performance based onab initio calculations using different definitions of partial charges and different sized basis sets as reference data.The new method was applied to the molecular model of a proton conducting ion channel where charges and polarizations are of concern for the function of the gating mechanism. The influenza A virus contains a tetrameric bundle of proteins in its envelope that form a membrane spanning proton channel, the M2 channel. The channel conducts protons essential for viral replication [37–40]. After the virus enters a cell and becomes embedded in an acidic compartment (pH 5-6), protons are imported into the virion balancing the pH-gradient across the membrane between lumen and cytoplasm. The resulting change in protein-protein and protein-lipid interaction triggers the uncoating of the virus. The channel formation can be inhibited by the antiviral drugs Amantadine [41] and Rimantadine [42]. Nevertheless, the virus becomes increasingly drug resistant [43] due to mutations. Therefore there exists a priority to clearly understand the structural basis for the opening and closing mechanism of the channel in order to find new drugs against the virus and at the same time addressing drug resistant viral strains [44].A wealth of information is available on the structure and structural details concerning the function of the M2 channel: structures from experiments have been discussed [44–58]; several solid state NMR structures of the M2 transmembrane domain (TMD) [44, 50, 52–54, 59, 60] and solution state NMR structures of M2 in micelles [46, 61, 62] as well as X-ray crystal structures of M2 TMD in octyl-glucoside [47, 63] have been published, and there are also a broad variety of calculated structural models available [49, 57, 58]. These proteins highlight different features in structure and function which seem to be strongly dependent on the membrane environment, pH adjustment, and temperature [56]. Also, a heterogeneous model was reported showing a helix-kink angle distribution [50].One feature these structures have in common is the core of the channel: It comprises four Trp-41 (indole rings) and four His-37 (imidazoles) amino acids. The latter groups are either protonated or deprotonated at physiologic range between pH 5 and pH 8 [64] which is essential for the channel function/gating mechanism. His-37 initiates the opening of the channel [65, 66] when it becomes positively charged through protonation. His-37 residues appear to be physically involved in proton conducting, while Trp-41 groups mechanically open and close the channel [67].The time-limiting step of the channel opening is the protonation/deprotonation of His-37 (imidazole side chain) proceeding on a ms-to-sec time-scale [68, 69]. Khurana et al. [49] performed MD simulations on four structures/channel models in different protonation states to show the relaxation of the His-37 groups. The authors propose the channel action as a peristaltic water pump with a “sphincter” formed by the Val-27 residues. There was no explanation how the protonation state of His-37 is “recognized” by hydrophobic Val-27 residues.Interestingly, the M2 channel has a surprising low conductivity. Two different mechanisms of proton conduction through the channel have been suggested: (1) A continuous “water-wire” enabling proton “hopping” [70] which is very fast [57, 71] or (2) a shuttling/sphincter mechanism that directly involves the protonation/deprotonation of His-37. The latter mechanism is more consistent with the relative slow conductance rate of 104 ions per second [49, 69, 72].There is evidence for a gating function of tryptophan rings [32, 73] from NMR investigations on membrane oriented channels with 15N and 19F labeled Trp-41 groups. The NMR investigations provided fluor-fluor distances for the opened and closed state. Cation-π interactions exist between the protonated imidazole ring of His-37 and the indole ring of Trp-41 [66]. Fluorescence investigations on M2 also demonstrated a pH dependence of the Trp-41 and His-37 interactions [74].A two-site Trp-41 exchange model has been introduced to explain the dynamics of the system. TheJ
N
C
γ
3 coupling of Trp-41 has been measured to be 2.6 Hz. However, Pérez et al. reported an error on the Karplus relation for calculating J
N
C
γ
3 of 0.8 Hz [75], which renders any conclusion concerning the conformational state uncertain. The relative large error could be an indication for the simultaneous occurrence of mixed protonation states within a sample, causing different conformational states of the side chain residues His-37 and Trp-41 at the active site, the “inner” gate. Therefore, the precise conformational states of the channel core are a key to understand the functional mechanics of the M2 channel.In this paper, we provide evidence for a plausible opening mechanism for the M2 channel: The tryptophan gate of the channel is pushed open by the histidines as consequence of their conformational change triggered by protonation. A conformational search of the gating model at different states offers the energetic rationale for a detailed opening mechanism by including all mutual polarizations in the calculation of the electrostatic energy.


## 2. Materials and Methods

### 2.1. Bond Polarization Theory for Partial Atomic Charges

The BPT (Bond Polarization Theory) is the starting point for the calculation of atomic charges, describing the dependence of expectation values from polarization [76]. The definition of a local charge operator q
^
A (site A) is a prerequisite to calculate atomic charges q
A:(1a)
q
A
=
Ψ
q
^
A
Ψ
=
Ψ
e
Z
A
-
∑
j
A
j
j
Ψ
.The expectation value is expressed as [23, 24, 26](1b)
q
A
=
∑
i
A
-
B
n
i
q
i
+
n
i
2
A
i
∑
x
≠
A
x
∈
C
χ
A
i
q
x
R
⃑
x
-
r
⃑
χ
A
i
-
χ
B
i
q
x
R
⃑
x
-
r
⃑
χ
B
i
.The general atomic charge operator is defined by the difference between the nuclear charge and the electron population. The latter depends strongly on the kind of atomic orbitals, j, associated with atom A. The molecular wave function Ψ is build up as a configuration interaction (CI) series of Slater determinants of bond orbitals (BOs). The BOs are linear combinations of hybrid atomic orbitals (AOs) χ
A
i and χ
B
i. The first sum in (1b) counts all bonds A
-
B of atom A denoted by i. The variable n
i refers to the occupation number of ith bonds; and q
i is the corresponding bond partial charge. The parameter A
i quantifies the change of the charge due to bond polarization. The outer sum adds the contribution of each bond of A to the total charge, consisting of the bond partial charge and the correction due to polarization caused by the surrounding atoms. Therefore, the inner sum, accounting for all polarizing atoms (denoted as system C), has to exclude atom A itself. The bond partial charges q
i and their variation caused by bond polarization A
i are the parameters that have to be determined by a parametrization procedure. This is achieved with known charges q
A. Two parameters need to be determined for each type of bond/respective combination of elements A
-
B. A sufficient large set of charges is necessary to obtain an overdetermined system of equations (see (1b)). The least square solution of this set of equations will provide the parameters q
i and A
i.The occupation number of bondsn
i is introduced to treat molecular systems with conjugated bonds between two atoms A and B. An empirical relationship between the contraction of an ideal single bond R
A
B and the valence of the bond n
i
=
n
i
(
A
-
B
) was applied to estimate the bond occupation numbers n
i (see (1b)) between two atoms A and B via the corrected/real bond length (r
A
-
B):(2)
r
A
-
B
=
R
A
-
B
-
0.37
ln
⁡
n
i
A
-
B
.The ideal single bond lengths are obtained from a parameterization of O’Keeffe and Brese [77].BPT chargesq
A are calculated from all other atomic charges q
X of the molecular system (see (1b)). The resulting system of equations has the dimension of the number of atoms in the molecular system.


### 2.2. Parameterization

A set of 163 molecules consisting of H, C, N, O, F, Si, P, S, and Cl atoms (with 25 bond types: H-C, H-N, H-O, C-F, C-N, C=N, C-C, C=C, C-O, C=O, Cl-C, P-O, P=O, Si-H, Si-C, Si-O, Si-Cl, S-H, S-C, S=C, S-O, S=O, S-S, Zn-O, and Zn-N) was selected to calibrate the BPT charge parametersq
i (22 numbers) and A
i (25 parameters) with 25 valences n
i and global scaling factor for the electrostatic energy.Based on experimentally available structural data [78], the initial molecular geometries were optimized with GAUSSIAN 98 [79] using hybrid functional option (DFT/B3LYP) for a 6-31G(d,p) basis set. The distances were checked against the original experimental data from literature [79, 80]. No major deviations in bond distances (≥0.02 Å) and bond angles (≥2°) were observed.The BPT method was parameterized using three different atomic charge models: Mulliken population analysis (MPA), charges derived from electrostatic potentials (ESP), and charges from a natural population analysis (NPA). The ESP method does not fit into the context of methods for population analysis since these charges are fitted to represent the electrostatic potential at the molecular surface. It is possible to derive the partial charge as expectation value of an operator (see (1a)) mathematically correct only for the MPA and NPA method. The basis set dependence of the charges was investigated using 11 different basis sets for each charge model (MPA, ESP, and NPA) (see Supplementary Material available online at http://dx.doi.org/10.1155/2015/908204). Table 1 shows the correlation coefficients of the parameterization of the BPT charge calculation of six selected charge calculation models; only results of the best and worst performing basis sets are listed. The overall best performance has been obtained with NPA/6-31G(d,p).Table 1
Correlation coefficients (R) and standard deviations (SD) of the parameterization of the BPT charge calculation based on selected different charge models MPA (Mulliken population analysis), ESP (charges derived from electrostatic potentials), and NPA (natural population analysis) of the best and worst performing basis sets for each model. The overall best performance, NPA/6-31G(d,p), is marked in bold. (Δ
Q is the total range of the charge distribution.)
Charge calculation model
Basis set
Correlation coefficient
R
Standard deviationSD [e]
Relative standard deviationSD/Δ
Q [%]
MPA
STO-6G
0.9944
0.0244
10.50
6-311++G(d,p)
0.8182
0.1463
47.05
ESP
3-21G
0.9768
0.0818
20.90
6-31+G(d,p)
0.9414
0.1155
31.75
NPA
6-31G(d,p)
0.9967
0.0437
8.10
6-311+G(d,p)
0.9795
0.0959
19.69For a special application, zinc was introduced as new element to prove the general applicability, the quality of the correlation, and the stability of the results. The new training set for the parameterization consisted of the same 163 molecules that was used earlier plus 12 additional zinc compounds. Again, the NPA/6-31G(d,p) basis set yielded the highest correlation coefficient compared to different charge models (MPA/STO-6G and ESP/3-21G) that showed also high correlation coefficients (see Table1) for the 175 compounds. Figure 1 depicts the correlation of the BPT charge calculation against NPA charges from DFT calculations with the 6-31G(d,p) basis set for 175 compounds of the calibration set. The ESP method did not perform such well as the initial parameterization.Figure 1
Correlation of NPA atomic charges calculated with BPT and DFT employing a 6-31G(d,p) basis set. The calibration set included 175 compounds containing the elements H, C, N, O, F, Si, P, S, Cl, and Zn. The correlation coefficient isR
=
0
.
9961 and the standard deviation SD = 0.05 e.A correlation coefficient of 0.9961 was obtained between BPT/NPA and DFT/NPA, and the standard deviation was 0.05 e. The inclusion of rather exotic zinc compounds did not drastically change correlation coefficient and kept the BPT parameters nearly unchanged. This demonstrates that (1) the BPT is able to reproduce atomic charges fromab initio calculations in a very satisfactory manner and (2) the interpretation of the BPT parameters regarding the polarization is physically correct.One major point to be addressed is the basis set dependence of the calculated atomic charges. The correlations of BPT/NPA parameters for different compact basis sets (3-21G and 6-31G(d,p)) are shown in Figure2. The correlation of the bond partial charges q
i yields a correlation coefficient of R
=
0.9973 (SD = 0.03 e), while the polarization parameters A
i show a R-value of 0.9969 and a standard deviation of 0.26 e/H. This indicates that sufficient parameter stability is maintained regarding the usage of different compact basis sets (no polarization functions are used); no basis set dependence of the calculated atomic charges was found in our tests.Figure 2
Correlation of BPT/NPA parameters using different compact basis sets. (a) Correlation ofA
i/3-21G and A
i/6-31G(d,p) with R
=
0.9969 and SD = 0.26 e/H. (b) Correlation of q
i/3-21G and q
i/6-31G(d,p) with R
=
0.9973 and SD = 0.03 e.
(a)
(b)In a next step, the BPT/NPA method was compared with DFT/NPA results of a pseudopeptide zinc complex that was not included in the training set, consisting of 64 atoms (H, C, N, O, and Zn). The correlation of the calculated atomic charges for the zinc complex (correlation coefficientR
=
0.9882, standard deviation SD = 0.07 e) is shown in Figure 3.Figure 3
Correlation between BPT atomic charges and DFT charges calculated employing the NPA method and a 6-31G(d,p) basis set for a pseudopeptide zinc complex not included in the training set. The correlation coefficient isR
=
0.9882 and the standard deviation SD = 0.07 e.A detailed comparison of all correlation coefficients regarding different charge models, basis set sizes, and numbers of the corresponding BPT parameters is given in Supplementary Material.One major focus of this work was to improve the description of intermolecular interactions for biomolecular applications. Intermolecular interactions are mainly described by the sum of Coulomb and Van der Waals energy terms within the COSMOS-NMR force field with fluctuating charges (for details see [26]). Since interacting electron distributions are replaced by potential functions, this simple approach will lead to a systematic bias in interaction energies.A scaling factor was introduced into the potential energy term in order to fine-tune the interaction energies. This scaling factor was determined iteratively by comparing interaction energies calculated with the COSMOS-NMR force field with interaction energies from the web service based force field evaluation suite by Halgren [81–83]. Among other types of data, this data base containsab initio calculated intermolecular interaction energies of small dimeric molecules. 66 hydrogen-bonded dimers were optimized with the COSMOS-NMR force field [25]. The intermolecular interaction energies were calculated as the difference between the energy of the dimer and the sum of the energies of the two monomers. The optimized structures were compared to the quantum mechanically optimized geometries from Halgren [81–83] by calculating the RMS deviations for the whole structures as well as for the hydrogen bond structures (distances and angles). A scaling factor was calculated as arithmetic mean of the scaling factors of each interaction pair. The whole procedure was repeated until the difference of the scaling factor per iteration was smaller than 10−5, which is 103 times below the error limit of the calculated BPT charges. The calculated scaling factor for the Coulomb energy of the COSMOS-NMR force field is 0.87245. A comparison of the mean deviations for the interaction energies and hydrogen bonding geometries from several force fields calculations including our own developed COSMOS force field is given in Table 2.Table 2
Comparison of hydrogen-bonded dimer geometries and interaction energies to HF/6-31G∗ data  (this table is an extension of a web-based force field evaluation suite by Halgren [81–83] and compares the COSMOS-NMR force field with data provided by the web page http://server.ccl.net/cca/data/ff_evaluation_suite/).
Force field
Interaction energy [kcal/mol]
X⋯Z distance [Å]
X-H⋯Z [°]
Mean deviations
COSMOS-NMR
0.41
−0.04
10.6
MMFF94
−0.25
0.02
14.6
MMFF94s
−0.24
0.03
18.1
CFF95
1.10
0.21
27.0
CVFF∗
2.24
0.23
26.1
MSI CHARMM
0.60
0.08
34.0
CHARMM 22∗
0.00
0.05
15.0
OPLS
−0.49
0.02
22.9
AMBER
0.08
0.04
26.9
MM2
1.92
0.06
39.5
MM3∗
2.12
0.22
34.0
∗Each of the force fields handles all 66 hydrogen-bonded dimers except for (i) CVFF, for which the Cerius2 implementation could not perform the atom typing for imidazole, (ii) MM3, which lacks parameters for 12 of the dimers, and (iii) CHARMM 22, which lacked parameters for about half the systems even with extensions made locally. Some problems occurred with the interaction energies of two cyclic dimers that were not included in the comparison (oxalic acid/water, cyclic ammonia dimer) and three charged dimers (formamidinium/water, formaldehyde iminium/water, and bidentate/water) that show larger deviations of the interaction energies (≥2 kcal/mol).The COSMOS-NMR force field shows slight overestimation of the interaction energies compared to the quantum mechanical data, whereas the hydrogen bond geometries show relative small deviations from the quantum mechanical calculated geometries; theX
⋯
Z distances are too short (0.04 Å) and the X
-
H
⋯
Z angles are larger (10.6°) on average than the HF/6
-
31
G
∗ calculated data.


### 2.3. Program Implementation

The routines for the BPT charge calculations were integrated into the COSMOS-NMR force field to enable the recalculation of the charges at every step of a molecular mechanics calculation. A graphical user interface for Windows (Win32) was implemented that integrates COSMOS-NMR into the general modelling and graphics program COSMOS [34].All calculations can be performed using the COSMOS-backend program (C++) that was compiled for several operating systems as, for instance, UNIX, LINUX, AIX, and Win32. The COSMOS-backend is controlled by command line parameters and a project file (
.
c
o
s
∗
). The project file (ACII) can be edited or alternatively generated using the GUI program COSMOS. A MPI parallelized version is available [35] for MD calculations.The most time consuming part of the charge calculation is the evaluation of the integrals over bond orbitals describing the bond polarization energies (see (1a) and (1b)). The expressions for the integrals that are needed as coefficients of the system of equations are given in the supporting material. These integrals are used to build up a matrix of the dimension equal to the number of atoms N. The computational time for setting up all N charge equations is proportional to the number of bonds multiplied by N
2. Charge calculation within the BPT framework means solving this set of linear equations for which the number of floating point operations is proportional to N
3. A time profile for the routines was generated and showed that the solution of the equation is about ten times faster than the calculation of all integrals. Highly optimized standard packages (LAPACK, see, e.g., [35]) are used to speed up the matrix inversion.


## 2.1. Bond Polarization Theory for Partial Atomic Charges

The BPT (Bond Polarization Theory) is the starting point for the calculation of atomic charges, describing the dependence of expectation values from polarization [76]. The definition of a local charge operator q
^
A (site A) is a prerequisite to calculate atomic charges q
A:(1a)
q
A
=
Ψ
q
^
A
Ψ
=
Ψ
e
Z
A
-
∑
j
A
j
j
Ψ
.The expectation value is expressed as [23, 24, 26](1b)
q
A
=
∑
i
A
-
B
n
i
q
i
+
n
i
2
A
i
∑
x
≠
A
x
∈
C
χ
A
i
q
x
R
⃑
x
-
r
⃑
χ
A
i
-
χ
B
i
q
x
R
⃑
x
-
r
⃑
χ
B
i
.The general atomic charge operator is defined by the difference between the nuclear charge and the electron population. The latter depends strongly on the kind of atomic orbitals, j, associated with atom A. The molecular wave function Ψ is build up as a configuration interaction (CI) series of Slater determinants of bond orbitals (BOs). The BOs are linear combinations of hybrid atomic orbitals (AOs) χ
A
i and χ
B
i. The first sum in (1b) counts all bonds A
-
B of atom A denoted by i. The variable n
i refers to the occupation number of ith bonds; and q
i is the corresponding bond partial charge. The parameter A
i quantifies the change of the charge due to bond polarization. The outer sum adds the contribution of each bond of A to the total charge, consisting of the bond partial charge and the correction due to polarization caused by the surrounding atoms. Therefore, the inner sum, accounting for all polarizing atoms (denoted as system C), has to exclude atom A itself. The bond partial charges q
i and their variation caused by bond polarization A
i are the parameters that have to be determined by a parametrization procedure. This is achieved with known charges q
A. Two parameters need to be determined for each type of bond/respective combination of elements A
-
B. A sufficient large set of charges is necessary to obtain an overdetermined system of equations (see (1b)). The least square solution of this set of equations will provide the parameters q
i and A
i.The occupation number of bondsn
i is introduced to treat molecular systems with conjugated bonds between two atoms A and B. An empirical relationship between the contraction of an ideal single bond R
A
B and the valence of the bond n
i
=
n
i
(
A
-
B
) was applied to estimate the bond occupation numbers n
i (see (1b)) between two atoms A and B via the corrected/real bond length (r
A
-
B):(2)
r
A
-
B
=
R
A
-
B
-
0.37
ln
⁡
n
i
A
-
B
.The ideal single bond lengths are obtained from a parameterization of O’Keeffe and Brese [77].BPT chargesq
A are calculated from all other atomic charges q
X of the molecular system (see (1b)). The resulting system of equations has the dimension of the number of atoms in the molecular system.


## 2.2. Parameterization

A set of 163 molecules consisting of H, C, N, O, F, Si, P, S, and Cl atoms (with 25 bond types: H-C, H-N, H-O, C-F, C-N, C=N, C-C, C=C, C-O, C=O, Cl-C, P-O, P=O, Si-H, Si-C, Si-O, Si-Cl, S-H, S-C, S=C, S-O, S=O, S-S, Zn-O, and Zn-N) was selected to calibrate the BPT charge parametersq
i (22 numbers) and A
i (25 parameters) with 25 valences n
i and global scaling factor for the electrostatic energy.Based on experimentally available structural data [78], the initial molecular geometries were optimized with GAUSSIAN 98 [79] using hybrid functional option (DFT/B3LYP) for a 6-31G(d,p) basis set. The distances were checked against the original experimental data from literature [79, 80]. No major deviations in bond distances (≥0.02 Å) and bond angles (≥2°) were observed.The BPT method was parameterized using three different atomic charge models: Mulliken population analysis (MPA), charges derived from electrostatic potentials (ESP), and charges from a natural population analysis (NPA). The ESP method does not fit into the context of methods for population analysis since these charges are fitted to represent the electrostatic potential at the molecular surface. It is possible to derive the partial charge as expectation value of an operator (see (1a)) mathematically correct only for the MPA and NPA method. The basis set dependence of the charges was investigated using 11 different basis sets for each charge model (MPA, ESP, and NPA) (see Supplementary Material available online at http://dx.doi.org/10.1155/2015/908204). Table 1 shows the correlation coefficients of the parameterization of the BPT charge calculation of six selected charge calculation models; only results of the best and worst performing basis sets are listed. The overall best performance has been obtained with NPA/6-31G(d,p).Table 1
Correlation coefficients (R) and standard deviations (SD) of the parameterization of the BPT charge calculation based on selected different charge models MPA (Mulliken population analysis), ESP (charges derived from electrostatic potentials), and NPA (natural population analysis) of the best and worst performing basis sets for each model. The overall best performance, NPA/6-31G(d,p), is marked in bold. (Δ
Q is the total range of the charge distribution.)
Charge calculation model
Basis set
Correlation coefficient
R
Standard deviationSD [e]
Relative standard deviationSD/Δ
Q [%]
MPA
STO-6G
0.9944
0.0244
10.50
6-311++G(d,p)
0.8182
0.1463
47.05
ESP
3-21G
0.9768
0.0818
20.90
6-31+G(d,p)
0.9414
0.1155
31.75
NPA
6-31G(d,p)
0.9967
0.0437
8.10
6-311+G(d,p)
0.9795
0.0959
19.69For a special application, zinc was introduced as new element to prove the general applicability, the quality of the correlation, and the stability of the results. The new training set for the parameterization consisted of the same 163 molecules that was used earlier plus 12 additional zinc compounds. Again, the NPA/6-31G(d,p) basis set yielded the highest correlation coefficient compared to different charge models (MPA/STO-6G and ESP/3-21G) that showed also high correlation coefficients (see Table1) for the 175 compounds. Figure 1 depicts the correlation of the BPT charge calculation against NPA charges from DFT calculations with the 6-31G(d,p) basis set for 175 compounds of the calibration set. The ESP method did not perform such well as the initial parameterization.Figure 1
Correlation of NPA atomic charges calculated with BPT and DFT employing a 6-31G(d,p) basis set. The calibration set included 175 compounds containing the elements H, C, N, O, F, Si, P, S, Cl, and Zn. The correlation coefficient isR
=
0
.
9961 and the standard deviation SD = 0.05 e.A correlation coefficient of 0.9961 was obtained between BPT/NPA and DFT/NPA, and the standard deviation was 0.05 e. The inclusion of rather exotic zinc compounds did not drastically change correlation coefficient and kept the BPT parameters nearly unchanged. This demonstrates that (1) the BPT is able to reproduce atomic charges fromab initio calculations in a very satisfactory manner and (2) the interpretation of the BPT parameters regarding the polarization is physically correct.One major point to be addressed is the basis set dependence of the calculated atomic charges. The correlations of BPT/NPA parameters for different compact basis sets (3-21G and 6-31G(d,p)) are shown in Figure2. The correlation of the bond partial charges q
i yields a correlation coefficient of R
=
0.9973 (SD = 0.03 e), while the polarization parameters A
i show a R-value of 0.9969 and a standard deviation of 0.26 e/H. This indicates that sufficient parameter stability is maintained regarding the usage of different compact basis sets (no polarization functions are used); no basis set dependence of the calculated atomic charges was found in our tests.Figure 2
Correlation of BPT/NPA parameters using different compact basis sets. (a) Correlation ofA
i/3-21G and A
i/6-31G(d,p) with R
=
0.9969 and SD = 0.26 e/H. (b) Correlation of q
i/3-21G and q
i/6-31G(d,p) with R
=
0.9973 and SD = 0.03 e.
(a)
(b)In a next step, the BPT/NPA method was compared with DFT/NPA results of a pseudopeptide zinc complex that was not included in the training set, consisting of 64 atoms (H, C, N, O, and Zn). The correlation of the calculated atomic charges for the zinc complex (correlation coefficientR
=
0.9882, standard deviation SD = 0.07 e) is shown in Figure 3.Figure 3
Correlation between BPT atomic charges and DFT charges calculated employing the NPA method and a 6-31G(d,p) basis set for a pseudopeptide zinc complex not included in the training set. The correlation coefficient isR
=
0.9882 and the standard deviation SD = 0.07 e.A detailed comparison of all correlation coefficients regarding different charge models, basis set sizes, and numbers of the corresponding BPT parameters is given in Supplementary Material.One major focus of this work was to improve the description of intermolecular interactions for biomolecular applications. Intermolecular interactions are mainly described by the sum of Coulomb and Van der Waals energy terms within the COSMOS-NMR force field with fluctuating charges (for details see [26]). Since interacting electron distributions are replaced by potential functions, this simple approach will lead to a systematic bias in interaction energies.A scaling factor was introduced into the potential energy term in order to fine-tune the interaction energies. This scaling factor was determined iteratively by comparing interaction energies calculated with the COSMOS-NMR force field with interaction energies from the web service based force field evaluation suite by Halgren [81–83]. Among other types of data, this data base containsab initio calculated intermolecular interaction energies of small dimeric molecules. 66 hydrogen-bonded dimers were optimized with the COSMOS-NMR force field [25]. The intermolecular interaction energies were calculated as the difference between the energy of the dimer and the sum of the energies of the two monomers. The optimized structures were compared to the quantum mechanically optimized geometries from Halgren [81–83] by calculating the RMS deviations for the whole structures as well as for the hydrogen bond structures (distances and angles). A scaling factor was calculated as arithmetic mean of the scaling factors of each interaction pair. The whole procedure was repeated until the difference of the scaling factor per iteration was smaller than 10−5, which is 103 times below the error limit of the calculated BPT charges. The calculated scaling factor for the Coulomb energy of the COSMOS-NMR force field is 0.87245. A comparison of the mean deviations for the interaction energies and hydrogen bonding geometries from several force fields calculations including our own developed COSMOS force field is given in Table 2.Table 2
Comparison of hydrogen-bonded dimer geometries and interaction energies to HF/6-31G∗ data  (this table is an extension of a web-based force field evaluation suite by Halgren [81–83] and compares the COSMOS-NMR force field with data provided by the web page http://server.ccl.net/cca/data/ff_evaluation_suite/).
Force field
Interaction energy [kcal/mol]
X⋯Z distance [Å]
X-H⋯Z [°]
Mean deviations
COSMOS-NMR
0.41
−0.04
10.6
MMFF94
−0.25
0.02
14.6
MMFF94s
−0.24
0.03
18.1
CFF95
1.10
0.21
27.0
CVFF∗
2.24
0.23
26.1
MSI CHARMM
0.60
0.08
34.0
CHARMM 22∗
0.00
0.05
15.0
OPLS
−0.49
0.02
22.9
AMBER
0.08
0.04
26.9
MM2
1.92
0.06
39.5
MM3∗
2.12
0.22
34.0
∗Each of the force fields handles all 66 hydrogen-bonded dimers except for (i) CVFF, for which the Cerius2 implementation could not perform the atom typing for imidazole, (ii) MM3, which lacks parameters for 12 of the dimers, and (iii) CHARMM 22, which lacked parameters for about half the systems even with extensions made locally. Some problems occurred with the interaction energies of two cyclic dimers that were not included in the comparison (oxalic acid/water, cyclic ammonia dimer) and three charged dimers (formamidinium/water, formaldehyde iminium/water, and bidentate/water) that show larger deviations of the interaction energies (≥2 kcal/mol).The COSMOS-NMR force field shows slight overestimation of the interaction energies compared to the quantum mechanical data, whereas the hydrogen bond geometries show relative small deviations from the quantum mechanical calculated geometries; theX
⋯
Z distances are too short (0.04 Å) and the X
-
H
⋯
Z angles are larger (10.6°) on average than the HF/6
-
31
G
∗ calculated data.


## 2.3. Program Implementation

The routines for the BPT charge calculations were integrated into the COSMOS-NMR force field to enable the recalculation of the charges at every step of a molecular mechanics calculation. A graphical user interface for Windows (Win32) was implemented that integrates COSMOS-NMR into the general modelling and graphics program COSMOS [34].All calculations can be performed using the COSMOS-backend program (C++) that was compiled for several operating systems as, for instance, UNIX, LINUX, AIX, and Win32. The COSMOS-backend is controlled by command line parameters and a project file (
.
c
o
s
∗
). The project file (ACII) can be edited or alternatively generated using the GUI program COSMOS. A MPI parallelized version is available [35] for MD calculations.The most time consuming part of the charge calculation is the evaluation of the integrals over bond orbitals describing the bond polarization energies (see (1a) and (1b)). The expressions for the integrals that are needed as coefficients of the system of equations are given in the supporting material. These integrals are used to build up a matrix of the dimension equal to the number of atoms N. The computational time for setting up all N charge equations is proportional to the number of bonds multiplied by N
2. Charge calculation within the BPT framework means solving this set of linear equations for which the number of floating point operations is proportional to N
3. A time profile for the routines was generated and showed that the solution of the equation is about ten times faster than the calculation of all integrals. Highly optimized standard packages (LAPACK, see, e.g., [35]) are used to speed up the matrix inversion.


## 3. Application and Results

The membrane protein M2 of influenza A virus consists of 97 residues including a 24-residue N-terminal and a 54-residue C-terminal segment. The transmembrane domain (TMD) consists of 19 residues that form anα-helical secondary structure. The α-helices associate with a tetrameric channel forming a pore in a lipid bilayer as illustrated in Figures 4 and 6 [74, 84, 85].Figure 4
Transmembrane part of the M2 proton channel; illustrated are the His-37 (green) and Trp-41 (blue) interactions that control the gating mechanism of the four helices in a bundle. The conformations of the Nishimura model at pH 8 of His-37 (χ
1, χ
2) = (−177°, +172°) and Trp-41 (χ
1, χ
2) = (−177°, −105°) are used (PDB-code: 1NYJ). The His-37 and Trp-41 conformations in this model cause strong Van der Waals contacts. (a) Bottom view. (b) Side view.
(a)
(b)Investigations under different conditions indicate that the tilt angle of the helix backbone orientation varies between 15° and 45° [46, 47, 86–93]. Channels without Amantadine show much less variation of the tilt angle. A tilt angle between 32° and 38° was detected for the channel in lipid DMPC. The tertiary structure seems not to change significantly during the switch between an open and closed state. Thus, we use the same backbone structure for our calculations (DMPC/water box) at different pH. A 25-residue peptide (SSDPLVVAASIIGILH37LILW41ILDRL) containing the hydrophobic TMD (P25-L46) was investigated. This TMD model and its equivalent mutants are widely used for structural investigations [52]. Here, we focused on the pH-dependent conformational change of the His-37 and Trp-41 side chains inside the channel pore.The COSMOS-NMR force field is especially suited for calculations of charged systems since all mutual polarizations are included into the electrostatic energy. The PDB structure 1NYJ [52] was used as initial model of M2 TMD. The protonation states determined by Hu et al. [64] were applied for the His-37-tetramer. The charged groups were modeled according to the pH values for the opened and closed channel configuration. Figure 5 illustrates Trp-41; the side chain torsion angles of His-37 are named accordingly.Figure 5
The side chain torsion angles of TRP-41:χ
1
(
N
,
C
α
,
C
β
,
C
γ
) and χ
2
(
C
α
,
C
β
,
C
γ
,
C
δ
1
).Figure 6
(a) The initial structure of M2 TMD from the Protein Data Base (1NYJ). Trp-41 are shown in blue, His-37 in green, hydrophilic/charged amino acid side chains in light blue, all other side chains in yellow. (b) View along the channel with a network colored by the electrostatic potential (from blue to red corresponding to negative to positive). (c) Equilibrated box with a size of 603 Ǻ3, containing 128 DMPC lipid, 3000 water molecules, and the M2 TMD structure in the centre.
(a)
(b)
(c)COSMOS-NMR force field energies were calculated for different His-37 and Trp-41 side chain conformations at pH 5 and mapped onto the total force field energy of the channel.Four His-37 have an average total charge of +1.3 at pH 8.0 and +3.0 at pH 5.0. The charged tetrameric M2 TMD structure was placed into an equilibrated box with a size of 603 Ǻ3, containing 128 DMPC lipid and 3000 water molecules; see Figure 6(c). All charges and energies, including the coordinate-dependent electrostatic interactions, were calculated for conformations of the four His-37 and Trp-41 residues under these conditions.Figures7 and 8 show conformational energy maps at pH 8 (closed channel) and at pH 5 (open channel), respectively. Panels (a), (d), and (g) picture the calculated total energies as a function of the His-37 (a) and Trp-41 ((d) and (g)) side chain torsion angles χ
1
(
N
,
C
α
,
C
β
,
C
γ
) and χ
2
(
C
α
,
C
β
,
C
γ
,
N
π
), as well as χ
1
(
N
,
C
α
,
C
β
,
C
γ
) and χ
2
C
α
,
C
β
,
C
γ
,
C
δ
1 in different configurations. The dark areas represent allowed rotamers within the channel in a DMPC/water bilayer environment. Figures 7(b), 7(e), 7(h), 8(b), 8(e) and 8(h) depict the calculated Van der Waals and (c), (f), and (i) the electrostatic energy contour plots for the different combinations of side chain torsion angles, χ
1 and χ
2.Figure 7
Torsion angle: energy maps of His-37 with fixed Trp-41 at (χ
1, χ
2) = (−100°, +110°) in (a), (b), and (c) and maps of Trp-41 with fixed His-37 at (χ
1, χ
2) = (+120°, +110°) in (d), (e), and (f) and (−180°, −180°) in (g), (h), and (i) for pH 8.0/closed state. Panels (a), (d), and (g) represent total energies, (b), (e), and (h) Van der Waals energies, and (c), (f), and (i) electrostatic energies. All energy maps are shown within a 5 MJ/mol range from minimum (black) to higher energies (dark grey, grey, light grey, and white).
(a)
His-37
(b)
His-37
(c)
His-37
(d)
Trp-41
(e)
Trp-41
(f)
Trp-41
(g)
Trp-41
(h)
Trp-41
(i)
Trp-41Figure 8
Torsion angle: energy maps of His-37 with fixed Trp-41 (χ
1, χ
2) = (−50°, +115°) in (a), (b), and (c) and maps of Trp-41 with fixed His-37 (χ
1, χ
2) = (+95°, +130°) in (d), (e), and (f) and (+100°, −50°) in (g), (h), and (i) at pH 5.0/open state. Panels (a), (d), and (g) represent total energies, (b), (e), and (h) Van der Waals energies, and (c), (f), and (i) electrostatic energies. All energy maps are shown within a 5 MJ/mol range from minimum (black) to higher energies (dark grey, grey, light grey, and white).
(a)
His-37
(b)
His-37
(c)
His-37
(d)
Trp-41
(e)
Trp-41
(f)
Trp-41
(g)
Trp-41
(h)
Trp-41
(i)
Trp-41The calculations allow only an estimate of forbidden conformational areas for the side chains of the gating residues His-37 and Trp-41 and provide reasons for the conformational restrictions that are concluded from the energy decomposition analysis, since the structures are not geometry optimized.An analysis of the calculations for the inactivated/closed state at pH 8 yields the following results: (1) Dihedral angels of Trp-41 were fixed according to experimental results (χ
1, χ
2) = (−100 ± 10°, +110 ± 10°) by Witter et al. [32]: (i) There is a broad forbidden area (between −110° and 60°) for the χ
1 angles of His-37 due to Coulomb interaction (see Figure 7(c)). (ii) Van der Waals interactions cause two main allowed conformational areas (see Figure 7(b)). (iii) From these two conformational areas each can be split into two rough torsion angles combinations of His-37: (χ
1, χ
2) = (−175°, +50°)/(+120°, +110°) and (−175°, −150°)/(+110°, −65°) which can be considered as rough energetically allowed/preferred conformations. The angle distribution around these values is between 20° and 50°. Transitions between these two regions possibly take place by 180° flips of χ
2 which correspond to flipping mechanism [72]. (2) Two representing angles: (χ
1, χ
2) = (+120°, +110°) and (−180°, −180°) of His-37 were used for the analysis of the torsion angles of the mechanically gating Trp-41 side chains. The resulting angle distributions for Trp-41 are shown in Figures 7(d)–7(i). (i) No major restrictions of the allowed conformational areas are attributed to the Coulomb interactions of Trp-41 (see widespread dark patches, Figures 7(f) and 7(i)). (ii) Van der Waals interactions cause major forbidden areas in the conformational space due to the bulky side chains of Trp-41 residues. (iii) Torsion angles of (χ
1, χ
2) = (−115° ± 50°, +130° ± 30°) are identified as energetically lowest conformations, which are close to the initial experimental values for the Trp-41 torsion angle calculation; see above: (χ
1, χ
2) = (−100° ± 10°, +110° ± 10°) [32]. Nevertheless view values with higher energies might be possible, for example, (−100°, −100°), depending on His-37 conformation. The resulting structure of a closed channel configuration is pictured in Figure 9.Figure 9
Side chain conformations of Trp-41 and His-37 in the transmembrane channel structure of the homotetrameric M2 protein. Representatives of (a) closed state at pH 8 for Trp-41: (χ
1, χ
2) = (−115° ± 50°, +130° ± 30°), His-37: (χ
1, χ
2) = (−180° ± 20°, −140° ± 20°) and (b) open state at pH 5 for Trp-41: (χ
1, χ
2) = (−60° ± 20°, +110° ± 20°) and His-37: (χ
1, χ
2) = (+95° ± 20°, +130° ± 20°) are shown.
(a)
(b)A similar analysis of the calculations for the activated/open state at pH 5 yields the following results: The panels in Figure8 show similar energy maps as panels in Figure 7 that lead to corresponding conclusions regarding forbidden conformational areas for the side chains of the gating residues His-37 and Trp-41. 1. His-37 was analyzed for fixed dihedral angels of Trp-41 according to experimental results (χ
1, χ
2) = (−50° ± 10°, +115° ± 10°) [32]; see Figures 8(a), 8(b), and 8(c). (i) Again, the Coulomb interactions cause a large forbidden conformational area for the χ
1 angles of His-37 (panel (c)). (ii) There are mostly only two areas for the torsion angles of His-37, (χ
1, χ
2) = (+95°, +130°) and (+100°, −50°), that keep the side chain locked in the open position of the channel due to Van der Waals interaction of the protonated His-37 side chains. The angle distribution is roughly 20°. The conformational difference is just about a 180° flip. These two dihedral angles of His-37, (+95°, +130°) and (+100°, −50°), were used for the conformational analysis of the side chain angles of Trp-41; see Figures 8(d)–8(i). (i) Again, restrictions of the conformational space are mainly caused by the bulky side chains of Trp-41 due to Van der Waals interaction that lock the side chain. The position with the lowest energy can be found around the initial experimental open state position (χ
1, χ
2) = (−70° ± 20, +110° ± 20). Still other conformations at higher energies are valid. During proton gating molecules of water (not considered in these calculations) interacting with His-37 side chains will restrict Trp-41 to such open configuration, presented in Figure 9.Figures7(a), 7(d), 7(g), 8(a), 8(d), and 8(g) show the sum of the electrostatic and Van der Waals energies. The total energy is dominated by Van der Waals interactions that are altered by the electrostatic contributions.Four Trp-41 and four His-37 residues participate actively in the gating of the channel structure. The function of the Trp-41 residues is a mechanical closing and opening mechanism of the channel. A change in pH causes the protonation/deprotonation of His-37 and causes a structural response of His-37 and Trp-41 driven by electrostatic and Van der Waals interactions between these residues. We conclude from our calculations that the bulky Trp-41 side chains enhance the steric effect of the His-37 residues that are directly involved in the H+ gating. In the open state the conformational space of His-37 is more restricted such that they are locked for gating. It also implies the possibility of a two-side jump model which supports the function of the imidazole rings flipping around χ
2 by almost 180°, thereby able to perform the “turning step” of passing the protons through the channel. At the closed state the His-37 side chains are locked due to hydrogen bridges and are not accessible for H+ gating. When the channel opens, the Trp-41 (folded back into the sides of the channel) are exposing the His side chains that only now are able to bind the hydrogen ions for passing through the open channel. That is in accordance with the flipping mechanism [72].


## 4. Discussion

BPT procedures for atomic charge calculation were integrated into traditional molecular mechanic force fields. Thus, a more detailed description of Coulomb interactions becomes feasible within the COSMOS-NMR force field without defining special QM regions. A careful parameterization of the BPT procedure allows the calculation of atomic charges virtually onab initio level for very large molecular systems. The charges depend on all mutual polarizations and provide with a correlation coefficient of 0.9967 improved accesses to electrostatic interaction energies.The COSMOS-NMR force field was applied to the proton channel of the influenza A virus. A detailed structure-function relation of the channel and its gating residues was obtained and attributed to Van der Waals and Coulomb interactions of the gating residues His-37 and Trp-41. Experimentally derived torsion angles of Trp-41 at pH 8 and pH 5 were theoretically approved. The presented states of the virus channel are not single valued or very fixed but represent preferred conformational distributions. His-37 exists in mainly two different conformational distributions which can be transformed by 180° flips aroundχ
2. His-37 are pointed to the centre of the channel at the open state. At the closed state they are locked via hydrogen bridges to the channel’s side. To change Trp-41 from the inactivated into the activated state, only moderate changes of torsion angle are necessary. With our analysis we confirm former results obtained with 15N-, 19F-NMR [32, 73] and the gating function of His-37 for proton conduction related to a His-37 flipping [72].


## Supporting Information

All coordinates of used structures within this publication are available from the authors upon request. The COSMOS-backend can be also obtained from the authors. The program is available without charge for scientific applications (request to Ulrich Sternberg: postmaster@COSMOS-software.de).



---

*Source: 908204-2015-09-30.xml*