A Kirkwood-Buff Force Field for the Aromatic Amino Acids
In a continuation of our efforts to develop a united atom non-polarizable protein force field based upon the solution theory of Kirkwood and Buff i.e., the Kirkwood-Buff Force Field (KBFF) approach, we present KBFF models for the side chains of phenylalanine, tyrosine, tryptophan, and histidine, including both tautomers of neutral histidine and doubly-protonated histidine. The force fields were specifically designed to reproduce the thermodynamic properties of mixtures over the full composition range in an attempt to provide an improved description of intermolecular interactions. The models were developed by careful parameterization of the solution phase partial charges to reproduce the experimental Kirkwood-Buff integrals for mixtures of solutes representative of the amino acid sidechains in solution. The KBFF parameters and simulated thermodynamic and structural properties are presented for the following eleven binary mixtures: benzene + methanol, benzene + toluene, toluene + methanol, toluene + phenol, toluene + p-cresol, pyrrole + methanol, indole + methanol, pyridine + methanol, pyridine + water, histidine + water, and histidine hydrochloride + water. It is argued that the present approach and models provide a reasonable description of intermolecular interactions which ensures that the required balance between solute-solute, solute-solvent, and solvent-solvent distributions is obtained.
Introduction<!><!>Introduction<!>Kirkwood-Buff Analysis of Experimental Data<!>Experimental Sources for Composition and Activity Data<!>Parameter Development<!>Molecular Dynamics Simulations<!>Analysis of Other Simulated Solution Properties<!>Results and Discussion<!>Benzene<!>Solution Mixtures<!>Benzene (1) + Methanol (2)<!>Toluene (1) + Methanol (2)<!>Benzene (1) + Toluene (2)<!>Phenol (1) + Toluene (2) / p-Cresol (1) + Toluene (2)<!>Pyrrole (1) + Methanol (2) / Indole (1) + Methanol (2)<!>Pyridine (1) + MOH (2) / Pyridine (1) + Water (2)<!>Histidine (2) + Water (1) / HisHCl (2) + Water (1)<!>Thermodynamic properties<!>Pure liquid properties<!>Solution properties<!>Conclusions
<p>Essentially all biochemical processes, and likewise many of industrial interest, take place in solution mixtures. Attempts to understand solutions have led researchers to develop, and often later disband, many theories and models over the last two centuries.1 Despite much advancement, our present understanding of solution mixtures is far from satisfactory. One way to measure our progress in this area is by the repeated evaluation of models developed to mimic the behaviour of solution mixtures and, in particular, their experimental properties. Additionally, accurate modelling allows for the future prediction of behaviour in situations where the experimental work is prohibitive on some grounds. One of the ultimate goals is to create solution models that would only require the composition and state of the system as input, yet the models would be able to reliably output that system's thermodynamic properties. The focus of the present work lies in the development and application of such models to eventually obtain thermodynamic and structural information for peptides and small proteins in various physiological environments using classical simulation techniques.</p><p>Several force fields (FFs) for the simulation of peptides and proteins have been proposed.2–6 However, in recent comparisons of established protein FFs a significant degree of reduced solvation, i.e. overestimation of solute-solute and solvent-solvent interactions, was frequently observed.7–15 Sometimes this leads to obviously erroneous results, such as the spontaneous formation of aggregates or phase separation in systems which should be fully miscible.7–9 Unfortunately, in many studies the errors are not immediately apparent from a visual inspection of the trajectory, but can often be inferred from, for example, the rotational diffusion of proteins.15 This imbalance in interaction strength is especially problematic when preferential interactions are of interest,16–21 or when considering simulations of systems for which little quantitative experimental data is available for comparison.</p><p>The origin of such misbalanced interactions is often related to the approximate description of electrostatic interactions provided by non-polarizable FF electrostatic parameters.7, 9 Consequently, polarizable FFs have been developed which, in principle, should perform better than non-polarizable FFs. However, they are computationally expensive and are not without their own set of difficulties.10 Thus, non-polarizable FFs still remain in widespread use when simulating large systems and/or attempting to study long time events. Popular methods to develop nonpolarizable partial charges for use in condensed phases have included the use of scaled gas phase quantum mechanical partial charges based upon calculations of the restrained-electrostatic potential, the use of scaled gas phase charges from ab initio calculations of the minimum interaction energies and geometries between dimers, and the parameterization of partial charges to reproduce thermodynamic and structural properties of pure liquids and the free enthalpies of solvation in water and other solvents.11 The primary disadvantage of each of these approaches is that the molecule is being studied in either the gas phase, under infinitely dilute liquid solution conditions, or as a pure liquid. Hence, because the molecule being parameterized is generally not considered in a variety of environments, the chosen charges may not optimally describe intermolecular interactions over a range of compositions, i.e. where solute-solute interactions become more important.</p><p>It is widely acknowledged that biomolecular force fields are imperfect and parameters are constantly being refined.11 However, the presence of several existing force fields suggests that the biological simulation community does not need another FF, unless it has been developed with a markedly different approach. One possibility involves a focus on appropriately balancing solute-solute, solute-solvent, and solvent-solvent interactions by studying these interactions far beyond the infinitely dilute range. Naturally, an important benchmark for such a FF would be its ability to reproduce solution activities. This is the role the Kirkwood-Buff Force Field (KBFF) seeks to fill. Since the solvation process distorts the electronic structure of the solute, with different parts of a molecule becoming polarized to different extents, we believe it is extremely difficult to develop accurate FF parameters, specifically effective condensed phase charge distributions, by either simply scaling gas phase charges or just reproducing pure liquid or infinitely dilute solution properties.22 Instead, we have developed a series of models where the best effective solution phase charges are derived empirically in order to fit experimental liquid solution data for mixtures of model compounds over a range of liquid phase compositions.7, 8, 23–29 Specifically, our approach attempts to obtain an effective charge distribution that best reproduces the experimentally derived Kirkwood-Buff integrals (KBIs) for representative solution mixtures.</p><!><p>The integrals provide additional data for use during parameterization.</p><p>One has easy access to changes in the solute activity and thereby quantitative information concerning intermolecular distributions, which are a sensitive property of solution mixtures.1, 30</p><p>Typically one does not have to sacrifice agreement for other properties of interest in order to reproduce the integrals.7, 8, 23–29</p><!><p>A variety of solute models have been developed using this type of approach.7, 8, 23–29 In our opinion it is important for a FF to be able to reproduce the features displayed by the experimental KBIs (see Results), because this assures that the FF is accurately balancing intermolecular interactions due to the KBIs' dependence on the solution activity. Furthermore, FFs that can reproduce the experimental KBIs could be used to gain insights into the very underlying solution microstructure that determined the KBIs. This is advantageous since the microstructure of solutions is generally difficult to determine experimentally. Thus, these FFs have the potential to significantly enhance the predictive power of solution molecular dynamics simulations.</p><p>We should note that outside of the biological FF community there has been significant work to develop FFs that reproduce solution activities, a property not normally used in the development of biological force fields. Notably the TraPPE FFs developed by the Siepmann group are specifically designed to model vapor-liquid coexistence curves.12–14 However, while they are attempting to reproduce the activities of each component in solution, only one composition is considered for a given T and p. For biological simulations it is more important to reproduce activities over a range of compositions for a relatively narrow window of (near-physiological) temperatures and pressures. Thus while the TraPPE approach is well suited for numerous thermodynamic and engineering purposes, it is of limited use in the typical biological context.</p><p>As we continue to develop a FF specifically designed to reproduce the experimental KBIs for condensed phase solution mixtures applicable to the simulation of peptides and small proteins,7, 8, 23–25, 27–29 models for the aromatic amino acid sidechains are required. In this work, toluene (Tol) was used as a model for the phenylalanine (Phe) sidechain, p-cresol (pCr) as a model for the tyrosine (Tyr) sidechain, and pyrrole (Pyrr) and indole (Ind) as models of the tryptophan (Trp) sidechain. Additionally, KBFF models for neutral histidine (His) and doubly-protonated histidine (HisH) were developed. First however, benzene (Ben) was studied as a model system to develop parameters for a fully conjugated carbon atom (Car). Then the new Car atom type was used in the parameterization of the partial charges for the other aromatics.</p><p>Beyond the mere necessity of developing these aromatic parameters to achieve completion of the KBFF for proteins, the systems studied are interesting in their own right. Aromatic molecules have unique features that make them challenging to model accurately, such as their ability to exhibit π-effects (cation-π, polar-π, and π donor-acceptor). This issue receives considerable attention due to the role π-effects frequently play in molecular recognition.31, 32 The negative electrostatic potential on the surface of the aromatic ring and positive electrostatic potential along the ring edges allows for a favourable interaction between cations or other polar group with the face of the ring. These effects are not explicitly taken into account in any classical FFs due to a lack of orbital descriptions, but the effects undoubtedly play an important role in many of the systems the community attempts to study using these simple FFs.</p><p>Previously, the Smith group has developed models for the sidechains of other amino acids by studying their analogues solvated in water. Unfortunately, many of the solute models that are representative of the aromatic sidechains are only sparingly soluble in water, making the study of aqueous solutions troublesome. To avoid this problem we have switched to solutions of aromatics in a solvent of methanol (MOH) or Tol where necessary. A KBFF model for MOH was previously developed by Weerasinghe and Smith from a study of aqueous solutions employing the SPC/E water model.28, 33 Hence, in this work we have made the usual assumption that our FF was fully transferable, i.e. a good description of Ben + MOH solutions coupled with a good description of MOH + water solutions provides for a reasonable description of Ben + water solutions.</p><!><p>Kirkwood-Buff (KB) theory was used to extract information concerning the microstructure of binary solutions from experimental thermodynamic data available in the literature. The process by which this was achieved has been thoroughly detailed elsewhere.7, 8, 25, 27–29, 34 Briefly, the KBIs (Gii, Gij = Gji, and Gjj) were obtained from the composition dependent chemical potential derivatives (µij), partial molar volumes (V̄i), and isothermal compressibility (κT) of binary mixtures at constant pressure (p) and temperature (T) according to1
(1)Gij=kTκT−δijρi+ρkT(1−ρiV¯i)(1−ρjV¯j)ρiρjμij where δij is the Kroenecker delta function, ρi = Ni/V is the number density of i molecules, ρ = ρi + ρj, µij = (δµi/δNj)p,T,Ni, k is the Boltzmann constant, T is the temperature of the liquid solution, Ni is the number of i molecules in the solution, and V is the total volume of the solution. Experimental activity data was obtained from literature values for the activity coefficients or excess molar Gibbs energies, GE, and were fit to an analytical function to facilitate differentiation (if the fit was not already performed in the original study). Partial molar volumes were calculated from experimental data for the solution's mass density or volume as a function of composition according to23
(2)V¯i=Vm+xjdVmdxi where Vm is the molar volume of the solution, and xi is the mole fraction of component i in the solution. For systems in which there was no experimental data available for the density or partial molar volumes as a function of composition it was assumed, based on precedence,34 that V¯i=Vi°, leading to Vm+xiVi°+xjVj0 and VE = 0. Here Vi° is the molar volume of pure component i and VE is the excess volume. It has also been shown that the value of κT has an insignificant effect on the calculation of the KBIs at ambient temperatures and pressures.34 Thus, for all systems the solution κT was calculated according to (3)κT=ρiV¯iκT,i+ρjV¯jκT,j where κT,i is the isothermal compressibility of pure component i. The KBIs extracted in the above manner are most sensitive to the experimental activity or Gibbs excess data.34, 35 Unfortunately, it is not always possible to accurately determine the appropriate errors associated with this data. In general, the errors associated with the Gij values tend to increase as the concentration of i or j decreases. To a large degree the errors in both the experimental and simulated KBIs can be minimized by the use of excess coordination numbers, Nij = ρjGij (see later discussion).</p><!><p>Benzene (1) + Methanol (2): Partial molar volumes were set equal to the molar volumes of the pure components.36 The Wilson equation was used to fit the GE data at 308 K.37 Toluene (1) + Methanol (2): Partial molar volumes were set equal to the molar volumes of the pure components.36 Activity coefficients were reported directly in an experimental isothermal vapor liquid equilibria (VLE) study at 318 K and 1 bar.38 From the activity data, GE was calculated as a function of composition and fit using the Redlich – Kister equation with a nonlinear least squares method involving three parameters. Toluene (1) + Benzene (2): Partial molar volumes were set equal to the molar volumes of the pure components.39 Activity data was obtained from experimental isothermal VLE data at 313 K using the Wilson equation.40 Toluene (1) + Phenol (2) / Toluene (1) + p-Cresol (2): Isothermal VLE data provided a fit for VE using the Redlich – Kister equation.41 The NRTL GE model was used to calculate the activity coefficients.41 Pyrrole (1) + Methanol (2): Partial molar volumes were set equal to the molar volumes of the pure components.39 Activity data was taken from experimental liquid-liquid equilibria (LLE) data at 298 K and fit using the NRTL equation for the ternary system methanol (1) + hexadecane (2) + pyrrole (3), where x2 was simply set equal to zero.39 Indole (1) + Methanol (2): Partial molar volumes were set equal to the molar volumes of the pure components.39 Activity data at 298 K was taken from the same source as for the pyrrole + MOH system in which activity data were provided for solutions up to values of xInd = 0.74. Pyridine (1) + Methanol (2) / Pyridine (1) + Water (2): The VE and GE data were obtained from measurements of vapor pressures and densities and modelled using the Redlich – Kister power series equation.42 Histidine (2) + Water (1) / Histidine Hydrochloride (2) + Water (1): Experimental data for His and HisH mixtures is sparse, presumably because histidine is only sparingly soluble in water. Volume data for His + water was obtained from a density versus composition study at 298 K and 1 bar.43 Activity coefficients were obtained from a vapor pressure osmometry study at 298 K and 1 bar with a maximum concentration of mHis = 0.25.44 We studied only the mHis = 0.25 composition with a tautomeric ratio of 15% HisA to 85% HisB.45–49 Partial molar volumes were set equal to the molar volumes of the pure components. Activity coefficients were obtained from the same study as for His + water at 298 K and 1 bar with a maximum mHisHCl = 0.63.44 We studied two compositions, mHisHCl = 0.3 and 0.6 due to the low solubility of HisHCl and the absence of a complete set of experimental data for imidazole. Analysis of the experimental and simulated data for salts is complicated by the presence of multiple ionic species in solution. We adopted the indistinguishable ion approach, which has been outlined in detail elsewhere.8, 28, 50</p><!><p>The KBFF approach uses a simple Lennard-Jones 6–12 plus Coulomb potential for all nonbonded interactions. Bond, angle, and torsion parameters were adopted (with permission) from Gromos with the exception of the H-Car-Car-Car improper dihedrals for six-membered rings, which were removed in order to improve the translational temperature of the aromatics (see Results for further discussion).51 The CH, CH2, and CH3 united atom parameters were also adopted from the Gromos force field,51 while polar hydrogen and heteroatom parameters were taken from our previous models. The Npy atom type was developed following the same scheme as described in detail previously,7, 27 in which the σ and ε parameters were determined by using the correlation between atomic size and atomic hybrid components of molecular polarizabilities.52</p><p>Initial studies with our existing sp2 carbon parameters resulted in poor agreement with experiment for the liquid densities. Hence, a new aromatic carbon type was developed to reproduce several properties of pure benzene (density, crystal structure, heat of vaporization). Traditionally, we have not used heat of vaporization data for pure liquids during our parameter development as we consider the KBFF models to be polarization corrected models in a similar fashion to the SPC/E water model,33 which is the water model of choice. However, the polarization corrections for benzene will be small and therefore the comparison with experiment should be reliable. The effective solution phase partial atomic charges used in the Coulomb potential were parameterized to reproduce the KBIs of binary mixtures. In our development of the partial charges for the aromatics, the charge group concept was adopted in an attempt to decrease the number of variables used in the parameterization process. Starting with an initial guess for the partial charges of the molecule of interest, the simulated KBIs were compared with those extracted from the experimental data for the same system for a range of solution compositions. Simulated KBIs were calculated by integration of the center of mass (com) based radial distribution functions (rdfs), which were obtained from the simulations in the usual manner. The required KBIs are given by30
(4)Gij=4π∫0∞[gijμVT(r)−1]r2dr≈∫0Rcut[gijNpT(r)−1]r2dr</p><p>The simulated KBI values were obtained by carefully averaging over a range of distances in which the com rdfs were essentially unity and the KBIs had reached a reasonably constant value. This small range of distances, typically corresponding to one solvation shell, was not necessarily the same for each system or for each of the three KBIs within a specific system. Some of the systems displayed significant aggregation at some compositions, often in agreement with experiment. However, the simulation results in which the KBIs were not converged using the present system sizes were not reported because they were considered statistically unreliable. Error estimates were obtained from the averages of multiple 5 ns runs.</p><p>Once both simulated and experimental KBIs were obtained, they were compared. If the solution KBIs agreed closely with experimental values, then the parameterization process was complete for that molecule. If not, a new set of partial charges was chosen and the process was repeated. Thus, the charges on the atoms were adjusted in an iterative manner to best reproduce the KBIs for solution mixtures over a range of compositions. Once the parameterization process was complete, other thermodynamic properties of the new models were assessed. The specific parameter development for each system is discussed in the Results section. The final atom type LJ parameters are shown in Table 1 and the final partial charges are provided in Figure 1.</p><!><p>All binary mixtures were simulated with classical molecular dynamics techniques using the Gromacs program in the isothermal isobaric (NpT) ensemble at a pressure of 1 bar.55 The simulation temperature for each system was chosen to match the temperature of the experimental GE data, which corresponded to 298 K for most systems except Ben + MOH (308 K), Tol + MOH (318 K), Ben + Tol (313 K), Tol + PhOH (333 K), and Tol + pCr (333 K). All systems were simulated at mole fractions of x1 = 0.00, 0.25, 0.50, 0.75, and 1.00. Additionally Ben + MOH, Tol + MOH, Tol + PhOH, Tol + pCr were simulated at x1 = 0.90, and Py + water was simulated at x1 = 0.10, as these regions displayed interesting behaviour. His + water was simulated at mHis = 0.25 and the HisHCl + water salt system at mHisHCl = 0.30 and 0.60. All systems were simulated in a 6 nm cubic box except for His and HisHCl, which were simulated in 8 nm cubic boxes due to the low molality of the solute. The Berendsen weak coupling technique was used to modulate the temperature and pressure with relaxation times of 0.1 ps and 5.0 ps, respectively, with a 4.5×10−5 bar−1 compressibility.56 All bonds were constrained using the Settle and LINCS algorithms for water and non-water molecules, respectively.57, 58 The use of bond constraints allowed for a 2 fs time step to be used for the integration of the equations of motion, which was performed using the Leap-Frog algorithm.59 The particle-mesh-Ewald technique was used to calculate the electrostatic interactions using a convergence parameter of 2.603 nm−1, cubic interpolation, a maximum fast fourier transform grid spacing of 0.12 nm for the reciprocal space sum, and tinfoil boundary conditions.60 The simulations employed a twin range cut-off of 1.2 and 1.5 nm for the real space electrostatic and van der Waals interactions, respectively.60 The neighbour list was updated every 10 steps. Each system was defined by a cubic simulation box into which molecules were randomly placed to generate the initial configuration for each mixture. The steepest descent method was then used to perform ≥1000 steps of energy minimization. This was followed by extensive equilibration that was continued until the intermolecular potential energy contributions and com rdfs displayed no systematic drift with time. Configurations were saved every 1.0 ps for analysis. All mixture simulations were run for at least 20 ns to ensure reasonable precision for the KBIs. We found that relatively long simulation times were necessary to yield reasonably precise data, presumably due to the large size of aromatics (slow diffusion) when compared to our previously studied molecules. Non-bonded interactions were excluded for neighbouring atoms up to six bonds away for benzene, toluene, phenol, and p-cresol, up to three bonds away for indole, pyrrole, and pyridine, and up to two bonds away for His and HisH. For phenol, a force constant of 7.24 kJ/mol with a multiplicity of two was used to constrain the angle between the CE-CZ-OH and CZ-OH-HH planes near 180° in the CE-CZ-OH-HH proper dihedral angle. The force constant was chosen based on the experimentally obtained values for this rotational barrier in phenol.61–63</p><!><p>The enthalpy of vaporization (ΔvapH) was calculated for the pure liquids according to (5)ΔvapH=−(Etot−Eintra)N+RT assuming that (1) the intermolecular interactions in the gas phase are negligible and (2) the intramolecular interactions are equal in the gas and liquid phases. Here, Etot is the total potential energy of the molecules in the liquid phase, Eintra is the intramolecular potential energy of the molecules in the liquid phase, N is the number of molecules in the liquid phase, and RT is the pV work contribution from the gas phase (ideal gas approximation). The pV work term is assumed to be zero for the liquid phase (incompressible liquid approximation).64, 65 The dielectric constant was determined from the dipole moment fluctuations using a reaction field permittivity of infinity as implemented by the g_dipoles Gromacs code.29, 66 For liquids, which are nearly incompressible, ΔmixV ≈ 0 thus ΔmixH ≈ ΔmixE, where the subscript 'mix' refers to mixing throughout.65 The enthalpy of mixing (ΔmixH) for each system was thus calculated according to64
(6)ΔmixH=EsinNsin−∑ixiEi°Ni° where Ei° is the total potential energy of pure i and Ni° is the number of i molecules in the pure i simulation and the subscript 'sln' refers to the solution properties. Finite difference κT values were obtained by performing additional 5 ns long simulations at pressures of 250 and 500 bar.64 Finite difference isobaric thermal expansion (αp) values were obtained from five additional 5 ns long simulations at 5 K temperature intervals encompassing the desired αp temperature. Translational self-diffusion coefficient (Dt) values were computed using the Einstein relation from the calculation of the mean square displacement (msd) of atoms from their initial positions, as implemented by the Gromacs g_msd code.</p><!><p>In this section we present and describe the comparison with the available experimental data. The primary goal was a consistent set of nonbonded parameters, developed to reproduce the experimental KBIs, for simple models that may be used for future biomolecular simulations. Unfortunately, while perfect agreement with experiment was the goal, the use of relatively simple (nonpolarizable) potentials necessarily introduced a degree of approximation and consequently some systems displayed some deviation from experiment. In particular, it would be expected to be difficult to fully reproduce all of the subtle polarization effects between the extremes in composition for binary mixtures of polar and nonpolar compounds.</p><!><p>Initially, we intended to use the same LJ parameters, previously developed for amides,25 for the carbon in Ben and to simply adjust the polarity of the Car-H bonds in Ben to reproduce the experimental KBIs. We found that the KBIs were indeed very sensitive to the charge distribution. For example, in a binary mixture of Ben + MOH, the degree of MOH-MOH (GMOH-MOH) aggregation noticeably differed depending upon the charge distribution across the Car-H bond. The difference was even observable from a simple visual inspection of the trajectories. However, the density and ΔvapH of pure Ben were far from the experimental values regardless of the choice of charge. Thus, we found that it was necessary to develop a new atom type for a fully conjugated carbon, Car, to produce better agreement for these properties as well.</p><p>Starting with an initial guess for the charge distribution (Car=- 0.14 |e|, taken from GROMOS) the Ben LJ σ and ε of the sp2 amide carbon type were adjusted to best reproduce the experimental liquid density of 0.8629 g/cm3 at 308 K and ΔvapH of 33.85 kJ/mol at 298 K.36, 67 The experimental KBIs were then extracted from simulations of Ben + MOH mixtures at 308 K, analysed, and the charge distribution was parameterized.</p><p>Although Ben is not considered a polar molecule, large differences were observed in the simulated KBIs as changes were made to the charge distribution (Figure 2). We believe this was related to small changes in the point charges representing a cumulatively large change in the overall charge distribution of the molecule due to the six Car-H bonds per molecule. Additionally, the change in charge parameters would have an even more pronounced effect on Ben–Ben interactions. As can be seen in Figure 2 there was a systematic reduction in the aggregation of MOH-MOH, as was evident from the decreasing GMOH-MOH as the charge distribution was increased. Starting with a carbon charge of q = −0.10 |e|, an incremental reduction in the KBIs occurred on scaling to q = −0.14 |e|, yielding a systematic improvement across the full composition range for each of the three KBIs. This was also observed from a visual inspection of the snapshots of these mixtures in which the number of MOH that were free in solution increased as the carbon charge was increased (results not shown). It should be noted that the simulated KBIs displayed for the q = −0.13 |e| distribution in Figure 2 are not in perfect agreement with those of the final KBFF due to the iterative procedure by which the σ and ε were later adjusted.</p><p>As we have developed the KBFF models, we have repeatedly seen a similar sensitivity of the KBIs on the charge distribution across many other systems.7, 8, 23–25, 27–29 We believe that the responsiveness of the KBIs to changes in the charge distributions is a major advantage of the current approach. The charge differences between different FFs for Ben and Phe are very small (Table 2), although an exact comparison across FFs is not straightforward due to the effects of other variables including different LJ parameters and different van der Waals and electrostatic cut-off distances.</p><p>Figure 3 indicates that Ben + MOH is a highly non-ideal system in which the large negative GBen-MOH at high xBen (middle panel) and the corresponding large positive GMOH-MOH in the same region of the composition range (bottom panel), quantify the reduced mixing between Ben and MOH and the corresponding microscopic self-association of MOH molecules, respectively.</p><p>The final Car parameters were found to be σ = 0.381 nm, ε = 0.33 kJ/mol, qcar=−0.13, and qH=+0.13. Fortunately, this was a q very near to that used in the parameterization of σ and ε and our final simulated Ben density and ΔvapH remained in good agreement with experiment despite iterating on q after choosing σ and ε. The final simulated values for the density and ΔvapH for benzene were 0.862 g/cm3 and 32.88 kJ/mol, respectively. The benzene crystal dimensions were also well reproduced with experimental (x, y, z) values of (37.049 Å, 48.664 Å, 35.264 Å) compared to our simulated (x, y, z) values of (37.271 Å, 48.955 Å, 35.475 Å), and an experimental density of 1.2119 g/cm3 compared to the simulated density of 1.1904 g/cm3. These results provided confidence that the parameters describing the size of an aromatic carbon, σ, and the atoms energy of most favourable interaction, ε, were appropriately chosen.72</p><p>In summary, the major features of this solution were accurately reproduced by MD simulations with the above parameters, but only after the partial charges were carefully parameterized to appropriately balance the solute-solute, solute-solvent, and solvent-solvent interactions.</p><!><p>The simulated com rdfs for all the systems studied here are shown in Figure 4. It is noticeable that the major changes in the rdfs with composition occurred for the alcohol (major hydrogen bonding) component. The rdfs for each mixture will be discussed in more detail in the sections below. Figure 5 compares the calculated and simulated KBIs after a multiplication by the number density. The quantity ρjGij is often referred to as an excess coordination number and is denoted as Nij. As mentioned previously, the use of Nij values provides a dimensionless quantity which diminishes the noise in the KBIs by effectively removing the uncertainty that is inherent at low concentrations of j in both the experimental and simulated KBI data. This uncertainty is attributed to sampling errors related to very few j molecules at high i concentration (and vice versa), and the relative difficulty of large aromatic molecules to move past each other at high aromatic concentrations. In addition, the Nij values have a well-defined physical interpretation, i.e. the excess or deficit in the number of j molecules within a certain volume of space away from a central i molecule when compared to the number that would have been found had the central i not been present.</p><!><p>The major feature in the KBIs for Ben + MOH, as provided in Figures 3 & 5, was the increase in GMOH-MOH (and NMOH-MOH) as xMOH decreases. The simulations using the new models provided excellent agreement with the experimental data. This feature was also mirrored in the gMOH-MOH rdf displayed in Figure 4 where the first peak dramatically increased in magnitude. Noticeably, gBen-Ben and gBen-MOH remained comparatively constant over the full concentration range.</p><p>A large positive Nij, such as that observed from the experimental peak of almost 10 for NMOH-MOH at xBen ≈ 0.7, can be interpreted as an indication of excessive clustering of MOH molecules at this composition. Physically, this corresponds to larger and larger deviations from a random distribution of methanol molecules as the number of possible hydrogen bonding partners decreases. The subsequent final decrease in NMOH-MOH presumably indicating an increasing entropic cost of maintaining these enthalpically favourable hydrogen bonded structures. If Ben + MOH were a symmetric ideal (SI) system, the GMOH-MOH (NMOH-MOH) would have been 2 cm3mol−1 (0) according to1
(7)GijSI=RTκT−Vi−Vj+ρ1V12+ρ2V22</p><p>The differences from the SI values were large and serve to quantify the deviation of this mixture from an ideal (essentially random) solution mixture.</p><p>A visual inspection of the xBen = 0.90 system revealed that MOH predominantly formed MOH-MOH chain structures over MOH-MOH cyclic structures. This was an observation that agreed with other simulations and experiments.73–76 However, a rigorous analysis of the percentages or lifetimes of the chain and cyclic structures was not performed. A more thorough discussion of this microstructure can be found in one of our recent publications.26</p><!><p>Although there is only a small difference between the atomic structure of Ben and Tol, the corresponding experimental NMOH-MOH values displayed a significantly greater peak, over twice the magnitude of that in the Ben + MOH system, at xTol ≈ 0.7. The major peak was in the same position, but much narrower than the peak in Ben + MOH. The simulations reproduced this increase in the NMOH-MOH maximum for the Tol + MOH, although the full increase was not reproduced quantitatively.</p><p>The charge distribution for the Car-CH3 bond in Tol was tested to determine if the best model should include polarity across the C-CH3 bond, or if the KBIs would be better reproduced with zero charges across this bond. From fundamental chemical principles it seems reasonable that one could argue either way. The electronegativity difference between sp2 and sp3 carbons is minimal, supporting charges of 0.00 and 0.00 on both the Car and CH3 sites. On the contrary, our simple model for aromaticity on the Ben ring would be disturbed under such an arrangement. We tested Tol + MOH (and Ben + Tol) with the Car-CH3 charges set equal to zero, against Tol + MOH (and Ben + Tol) using the same charges for Car-CH3 as those across the Car-H of the benzene ring. The results showed a clear improvement when the polarity across the bond was retained. Thus, the only way in which the final KBFF Tol and Ben models differ from one another is in their LJ parameters as they have identical coulombic parameters.</p><!><p>Ben + Tol is one of the quintessential ideal systems. The Ben + Tol com rdfs and Nij were not included in Figures 4 and 5 because we have discussed the excellent simulation of this ideality previously.77 The Tol + Ben com rdfs displayed a first solvation shell at ~0.6 nm, however the gBen-Ben and gTol-Tol distributions had slightly different magnitudes with a small increase in the probability for finding a Ben ~0.6 nm away from a central Ben than a Tol ~0.6 nm away from a central Tol. This was likely due to the tighter packing that is possible between Ben molecules due to the absence of a methyl group when compared to Tol. The simulated Gij and Nij curves were also in agreement with the values predicted by Equation 7 and further supported the ideality of the system.77</p><!><p>We first attempted to use the same MOH partial charges for the Car-O3-H charge group in phenol and cresol; however, this resulted in relatively poor agreement with the experimental KBIs. After several attempts the charges which appeared to best represent the experimental KBIs corresponded to a simple scaling of the MOH charges by a factor of 0.80. The difference in the partial charge distribution between aromatic and aliphatic alcohols may be attributed to both the aromaticity of the PhOH and pCr molecules and the large size of PhOH and pCr when compared to MOH, which allows for fewer hydrogen bonding partners to be formed in these solutions. An increase in the first peak of gTol-Tol was observed as xTol approaches one. Cyclic clusters of PhOH or pCr molecules were each observed in an arrangement vaguely reminiscent of reverse micelles. A more thorough discussion of this microstructure can be found in our recent publication.26</p><!><p>As indicated in Figure 1 we were not forced to adjust the charges around the ring to obtain good agreement with the experimental KBIs. Instead it was possible to retain the simple group concept and merely adjust the dipole across the N-H bond of indole. The same N-H dipole provided satisfactory results for both pyrrole and indole molecules. The com rdfs for these systems showed structured secondary solvation shells and a significant increase in the first solvation shell in all three com rdfs as the concentration of Pyrr or Ind was increased. The scale of these changes was much smaller than that observed for the Ben + MOH and Tol + MOH systems, so this did not necessarily suggest these systems were highly aggregating. Indeed, both the experimental and simulated excess coordination numbers did not indicate an excessive degree of association between species in these solutions.</p><!><p>In the comparison of the Py + MOH and Py + water systems, the first solvation shell of gwater-water exhibited a higher peak than gMOH-MOH at a slightly smaller radius due to the smaller excluded volume of water when compared to MOH. The Py + water Nij were reproduced quite well, except at xPy = 0.1. Experimentally, there was a high degree of self-aggregation at this composition. However, our simulations did not appear to fully capture this aggregation. In particular, the self-association of Py (NPy-Py) was underestimated, while the solvation of Py by water (NPy-water) was overestimated. Further attempts to improve the Py + water KBIs at low Py concentrations were unsuccessful and typically resulted in a corresponding reduction in the agreement with the Py + MOH KBIs and the Py + water ΔmixH (see later).</p><!><p>As mentioned previously, the experimental data for His and HisH mixtures is sparse, presumably because histidine is only sparingly soluble in water. Thus, the charge distributions were developed based upon the previous studies of Pyrr + MOH and additional studies of Py + MOH and Py + water over the full composition range, and then tested using simulations of His + water at 0.25 mHis and HisHCl + water at 0.3 and 0.6 mHisHCl. The N-H Pyrr charge distribution obtained for use in Trp was used as an initial guess of the charges on the protonated sp2 nitrogen (N2) present in His and HisH. The charges on HisH were then adjusted to achieve better reproduction of the experimental KBIs. Py was studied to model the charges of neutral histidine's unprotonated sp2 nitrogen (Npy) and the charges on the two Car-H that flank this nitrogen. The Pyrr and Py results were then combined to develop a crude approximation for the charges in His. His + water was then simulated to ensure that the amino acids did not aggregate excessively – a minimum requirement.</p><p>The His + water and HisHCl + water systems both exhibited similar gwater-water with a sharp first solvation shell and well-structured second and third solvation shells. gHis-water did not display any strong peaks, while gHisHCl-water had a sharp first solvation shell at a distance of ~0.3 nm with a peak slightly less than 2.25. gHis-His and gHisHCl-HisHCl both showed first solvation shells at a distance of ~0.7 nm with broad shoulders, however the gHisHCl-HisHCl peak at mHisHCl = 0.3 was ~4 and the gHis-His peak was ~2. Even in a relatively large simulation box with a length of ~8 nm, the gHisHCl-HisHCl rdf did not reach unity. Thus the errors in NHisHCl-HisHCl were significant and they were not reported due to their unreliability.</p><p>The Nij were well reproduced in the His + water system. Although we only considered one composition, we believe that this demonstrated a modest description of neutral histidine in our simulations. The HisHCl + water Nij results were also quite good, especially when considering the errors associated with the NHisHCl-water values. Although we did not gain meaningful information on the NHisHCl-HisHCl values, due to an unconverged gHisHCl-HisHCl, a visual inspection of the mixture did not reveal any significant self-aggregation of the amino acids with the final set of partial charges. It should be noted that ion pairing and histidine clustering were observed with other sets of trial charges, but that this clustering was not present using final charge distributions. Our analysis of the site-site coordination numbers for the His + water and HisHCl + water systems did not exhibit any interesting features (data not shown). They displayed small coordination numbers with values significantly less than one, further indicating weak interactions between His-His and HisHCl-HisHCl in the systems studied.</p><!><p>These KBFF models were designed to reproduce solution properties, with a particular focus on the experimental KBIs. However, the KBFF models developed for previous solutes have been shown to be competitive with existing FFs in their ability to reproduce other properties of solution mixtures and the properties of pure liquids.7, 8, 23–27 Traditionally, solution properties are not necessarily well reproduced with other FFs unless, of course, a certain property was the focus of the parameterization effort.9, 27 Here we explored the ability of the KBFF to reproduce (where data is available) the experimental ΔvapH, αp, κT, and dielectric constants (ε) of the pure liquids and the translational self-diffusion coefficient (Dt) and ΔmixH and ΔmixVE of the solutions. No values were explicitly corrected for polarization or vibrational frequency changes.</p><!><p>As mentioned previously, the iterative nature of the approach used for obtaining the parameters of Car, in which the charge distribution was adjusted after the choice of the LJ σ and ε was made, resulted in small deviations from experiment for ΔvapH. Additional errors in our results may have been due, at least in part, to the omission of polarizability and a lack of quantum corrections for vibrations, since not all the necessary data was available for calculation of these contributions. However, the SPC/E water model and the KBFF MOH model have been corrected for polarizability and have been shown to accurately reproduce the ΔvapH data.29, 96 Thus, polarizability was at least implicitly accounted for in systems where one species was either water or MOH. Nevertheless, the enthalpy of vaporization data was in good agreement with experiment considering only the Ben values were used during the parameterization procedure. The values of the αp (Table 3) were also in good agreement with experiment.</p><p>The calculated compressibilities of the pure liquids (Table 3) were consistently lower than the experimental κT values. The cause of this deviation was unclear. During an early iteration of the parameters for the Car of Ben, the compressibility of pure Ben was observed to be severely low (Ben κT,Exp=9.67×10−5 bar−1, Ben κT,Sim=4.39×10−5 bar−1), which is actually lower than the experimental and simulated values of water (water κT,Exp=4.5×10−5 bar−1). This was related to an inability to equilibrate and/or maintain the translational temperature (Ttrans) of the aromatic molecules. The value of Ttrans was consistently too low for the pure aromatics (output Ttrans=273 K for pure Ben, with TSim=308 K). Many attempts were made to elucidate the main cause for the low Ttrans of the pure aromatics, e.g. testing different temperature coupling algorithms, removing bond constraints while using a timestep of 0.5 fs, running canonical ensemble (NVT) simulations to eliminate the possibility of incorrect coupling between T and p induced by the timestep, etc. The only improvement (increase) in the values of Ttrans and κT was obtained on removal of the H-Car-Car-Car improper torsional potentials in Ben and the other aromatic molecules. It should be noted that the improper torsional potentials that were removed are not required to maintain planarity in other force fields.3 This led to an output Ttrans=299 K (TSim=308 K) and an improved κT,Sim=7.5×10−5 bar−1 for Ben. Of the properties we investigated, the removal of the improper potentials only affected Ttrans and κT for the pure liquids, with the properties of mixtures being essentially unaffected.</p><p>Values for the dielectric constants of the pure liquids are provided in Table 3. These values agreed moderately well with experiment, with the exclusion of pCr in which the appropriate trend in the values on going from PhOH to pCr was not observed. Simulated dielectric constants as a function of composition are shown in Figure 6 and exhibited the expected trends with composition. Unfortunately, we could not find any appropriate experimental data for comparison.</p><!><p>The ΔmixH results are shown in Figure 7 for the systems in which experimental data was available. The simulated results underestimated the mixing enthalpy for Ben + MOH solutions, although the Tol + MOH values were reasonable. Unfortunately, the simulated values were not sensitive to changes in the charge distribution. Recent investigations seem to indicate that explicit polarization may be required for obtaining accurate enthalpies of mixing, and that even reasonable models can provide incorrect signs.97, 98 The ΔmixH for Py + MOH indicated that mixing of the solute and solvent was too unfavourable, while the ΔmixH for Py + water was well reproduced across the full composition range. We are unsure of the reason for this disagreement; however we must point out that the ΔmixH is a property that is nontrivial to model correctly, and similar deviations from experiment are common in the literature.97, 99 Through iterative tests on the aromatic carbon's σ and ε, and through iterative tests on the partial charges of Ben, the ΔmixH changes for Ben + MOH mixtures were always small. Hence, it was not clear which parameters should be modified in order for the model to reach agreement with experiment, and it appeared doubtful that a single Py model could reproduce both the water and methanol data.</p><p>In Figure 8 we provide a comparison of the excess volume of mixing for several of the solutions investigated here. The results appear to be very reasonable, especially in comparison to other FFs.100 Previous studies have also indicated systematic deviations from experiment for the solution volume (or density) data.100 There appears to be no such trend in the models developed here.</p><p>The experimental translational self-diffusion coefficients for Ben + MOH and Ben + Tol were available and were compared to the simulation results obtained using the KBFF models (Figure 9).104 For Ben + MOH, Dt,Ben was in close agreement with the experimentally observed trend that as the concentration of Ben increases Dt,Ben decreases, except at xBen = 0.90 where we saw a small increase in Dt,Ben.104 Experimentally, Dt,MOH decreases as xBen increases from xBen = 0.00 to xBen ~ 0.50, and then the Dt,MOH increases as xBen approaches 1.0.104 Our model for Ben + MOH does not capture the increase in Dt,MOH at low xMOH, but the low simulated values seemed to agree with the presence of MOH-MOH clusters observed in our simulations and by others,73–75, 105, 106 which one would expect to diminish the rate of diffusion of MOH. The experimental work referenced for the Ben + Tol system describes the diffusion in this system using a single linear fit in which the values of Dt,Ben were assumed equal to the values of Dt,Tol over the entire composition range. Simulated values for Ben + Tol spanned the experimental values across the entire range of compositions. We note that our results were not corrected for finite size effects, which typically result in 5–10% larger diffusion coefficients.107</p><!><p>We have presented the Kirkwood Buff Force Field parameters for the side chains of the aromatic amino acids and have shown how simulations of these models compare with the available experimental thermodynamic and structural data in solutions over a range of compositions. The ability of our simple nonpolarizable force fields to reproduce the KBIs over the full composition range is a non-trivial result.8, 19, 21, 25, 28, 108 For example, the ability to accurately model benzene + methanol and toluene + methanol activities over the range from pure methanol to pure benzene or pure toluene suggests that we will have a good model for phenylalanine in a variety of environments, i.e. as a free amino acid, as a solvent exposed residue, or as a residue buried deeply within a hydrophobic protein core. Furthermore, the ability of the classical models used in the KBFF to mimic the experimentally observed aromatic π-effects can now be tested and compared with other FFs. Finally, reasonable results were obtained for other solution properties not included in the parameter development.</p><p>The choice of a force field should depend upon which properties of a system a particular researcher wishes to probe. The overall message of this work is not meant to be that traditional force fields are fundamentally inferior to the KBFF models. We would be greatly amiss not to recognize that in many areas biological simulations using traditional force fields have certainly been successful at reproducing experimental properties, and even providing new hypotheses for experimentalists to subsequently test. Where we believe we have made an important improvement over traditional force fields is in quantifying the balance of intermolecular interactions. Before adopting the KBFF models it is always prudent to compare results from other force fields and experiment with results obtained using the KBFF for the properties of interest.76</p>