Article pubs.acs.org/JPCB
Parametrization of 2,2,2-Trifluoroethanol Based on the Generalized Amber Force Field Provides Realistic Agreement between Experimental and Calculated Properties of Pure Liquid as Well as Water-Mixed Solutions Jiří Vymětal and Jiří Vondrásě k* Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic (AS CR), Flemingovo nám. 2, 166 10 Praha 6, Czech Republic S Supporting Information *
ABSTRACT: We present a novel force field model of 2,2,2-trifluoroethanol (TFE) based on the generalized AMBER force field. The model was exhaustively parametrized to reproduce liquid-state properties of pure TFE, namely, density, enthalpy of vaporization, self-diffusion coefficient, and population of trans and gauche conformers. The model predicts excellently other liquid-state properties such as shear viscosity, thermal expansion coefficient, and isotropic compressibility. The resulting model describes unexpectedly well the state equation of the liquid region in the range of 100 K and 10 MPa. More importantly, the proposed TFE model was optimized for use in combination with the TIP4P/Ew and TIP4P/2005 water models. It does not manifest excessive aggregation, which is known for other models, and therefore, it is supposed to more realistically describe the behavior of TFE/water mixtures. This was demonstrated by means of the Kirkwood−Buff theory of solutions and reasonable agreement with experimental data. We explored a considerable part of the parameter space and systematically tested individual combinations of parameters for performance in combination with the TIP4P/Ew and TIP4P/2005 water models. We observed ambiguity in parameters describing pure liquid TFE; however, most of them failed for TFE/water mixtures. We clearly demonstrated the necessity for balanced TFE−TFE, TFE−water, and water−water interactions which can be acquired only by employing implicit polarization correction in the course of parametrization.
■
proteins.15−21 However, the explanatory power of MD directly reflects the quality of the force field parameters utilized. Chitra and Smith tested several force field models of TFE15,22,23 for performance in simulations of pure liquid TFE as well as TFE/ water mixtures.24 Depending on the particular force field, TFE models varied in a broad range of parameters, specifically Lennard-Jones parameters for fluorine and partial atomic charges on the hydroxyl and fluoromethyl groups. None of their TFE models simultaneously and correctly reproduced all properties of pure liquid. More importantly, their models of TFE/water mixtures indicated overexaggerated self-aggregation of TFE molecules. The commonly used family of AMBER force fields offers an automated protocol for parametrization of small organic compounds with the general AMBER force field (GAFF).25 A TFE model represented by GAFF parameters reproduced pure liquid properties in moderate agreement with experimental values,26 leaving room for further refinement. The obvious shortcomings are underestimated density (or overestimated molar volume) and a low self-diffusion coefficient. Although successful simulation of a peptide in mixed solvent modeled by GAFF parameters has been reported,20 characterization of
INTRODUCTION 2,2,2-Trifluoroethanol (TFE) ranks among the most wellestablished cosolvents used in experimental characterization of peptides and proteins. The effects of TFE involve stabilization of the peptide conformation and solubilization of peptides and proteins.1 TFE can induce helix formation in peptides that are unstructured in standard buffers,2 facilitate β-sheet-to-helix transitions3,4 and allow measurement of the helical propensities of different peptide sequences in a qualitative manner.5 However, TFE can also stabilize other structural elements such as β-hairpins6 or intermolecular strands.7 Several mechanisms explaining how TFE affects the behavior of polypeptides in solution have been proposed. One class of hypotheses assume direct interaction and binding of TFE molecules to polypeptides and direct competition with water molecules. The TFE−polypeptide interactions themselves or the stability of the internal hydrogen bonds that they promote due to the local hydrophobic environment is usually attributed to structural and thermodynamic changes.8−12 Self-aggregation of TFE molecules in TFE/water mixtures can also contribute to polypeptide behavior in solution.13 Another class of hypotheses suggest an indirect effect via the thermodynamics of folding; i.e., the chaotropic effect of TFE on the water solvation layer is different in the folded and unfolded states.14 Molecular dynamics (MD) can be a valuable tool to elucidate the atomistic mechanisms of TFE’s effects on peptides and © XXXX American Chemical Society
Received: June 12, 2014 Revised: August 6, 2014
A
dx.doi.org/10.1021/jp505861b | J. Phys. Chem. B XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry B
Article
Lorentz−Berthelot combination rules, which use AMBER force fields for calculation of heterogeneous Lennard-Jones interactions. However, the performance of classic combination rules has been questioned for fluorinated compounds.42,43 Finally, the following parameters, which to a large extent determine the liquid-state properties, were included for optimization: 1. Lennard-Jones parameter σ of fluorine (σF), which drives the density of the liquid phase. 2. Lennard-Jones parameter ε of fluorine (εF), which has a major effect on the heat of vaporization. 3. Polarity of O−H bonds (μOH), which significantly influences the self-diffusion coefficient and heat of vaporization. 4. Additional torsional term for the C1−C2−O−H dihedral angle (τ), which rectifies the ratio of gauche and trans conformers of the TFE molecule. 5. Polarity of C−F bonds (μCF), which was added after we identified a need for more polar TFE molecules. Parameters representing the polarity of covalent bonds were applied by bond-contribution-based decomposition of partial charges. First, the initial partial charges qi on each atom were distributed into the dipole contribution μij:
TFE/water mixtures themselves revealed that TFE selfaggregation closely depends on the water model used.27 Performance with a TIP4P/Ew water model28 was poor, leading to complete separation of both liquids despite their full miscibility observed in experimental studies.29−32 Even if compatibility with TIP4P/Ew as a promising water model for further development of AMBER force fields33,34 is favored, the current GAFF model of TFE must be optimized. The Kirkwood−Buff theory of solution provides a suitable connection between computer simulations and experimental thermodynamic data. Moreover, Kirkwood−Buff integrals (KBIs) shed light on the microscopic structure of the mixtures and allow us to estimate the affinities of individual compounds and quantify their aggregation tendencies.35−37 A correct and balanced description of TFE−water interactions in mixed solvent should be the first step in the process of modeling more complex interactions with biomolecules. The TFE solutions are suspected to form transient clusters,13,38 which could influence interpretation of experiments and simulations. Therefore, TFE parameters able to correctly describe the structural, thermodynamic, and kinetic properties of TFE/water mixtures are desired. In this study, we addressed the refinement of the GAFF TFE model by optimizing force field parameters to establish better agreement with the experimental properties of pure TFE in its liquid phase. Moreover, we propose to enhance the TFE model by implicit polarization correction similarly to the successful parametrization of the TIP4P/Ew and TIP4P/2005.39 This step allows us to reasonably fit the experimental properties of the liquid phase, which can be reproduced by many combinations of optimized force field parameters. Nevertheless, only a few of them were found compatible with TIP4P/Ew and TIP4P/2005 in simulation of TFE/water mixtures, as assessed by application of the Kirkwood−Buff theory. Finally, we describe the best TFE model, which reproduces the liquid-state properties as well as the behavior of TFE in binary mixtures.
qi =
∑ μij j
μij = −μji
(1) (2)
The summation is evaluated through all atoms j bonded directly to the atom i. Consequently, the bond dipoles μCF and μOH were varied in the optimization procedure. The signs of μCF and μOH were adjusted to simultaneously represent the partial atomic charges on the fluorine and hydrogen atoms, respectively. Optimized Properties. The pure-liquid-state properties used for parameter refinement included density, enthalpy of vaporization, self-diffusion coefficient, and ratio of trans and gauche conformers at a temperature of 298.15 K and pressure of 1 bar. All these properties could be easily calculated with reasonable precision from a single molecular dynamics simulation of liquid accompanied by a single simulation of one TFE molecule in the gas phase. Consequently, they were evaluated using objective functions as described comprehensively in the Appendix. A detailed description of how the individual properties were calculated is provided, as well as the necessary details of workflow for optimization. Setup of ab Initio Calculations. All ab initio quantum electronic calculations were performed in Gaussian 09.44 The polarizability tensor and dipole moment were analyzed on the density functional theory (DFT) level using the B3LYP functional45 and diffuse aug-cc-pVTZ basis sets commonly used for prediction of electronic properties. Full or constrained geometrical optimization of the molecular structure was conducted before each analysis, utilizing the same method and basis set. Setup of MD Simulations. All MD simulations were performed using the Gromacs simulation package, version 4.5.5.46 All liquid-phase simulations used the following setup: Equations of motion were integrated by the leapfrog method with 2 fs time steps with constraints on all bond lengths enforced by a fourth-order LINCS algorithm with one iteration. Periodic boundary conditions were employed together with a
■
METHODS Parametrization Strategy. Our model was intended to reproduce the pure liquid properties of TFE and the properties of a TFE/water mixture over the entire dilution range. The parametrization process was divided into three stages. First, the force field parameters for optimization were varied and selected on the basis of the sensitivity of the liquid-state properties. Second, the liquid-state properties were optimized by an exhaustive search in parameter space using objective functions. Finally, the best-performing sets of parameters were critically tested in simulations of TFE/water solutions, and the quality of optimized parameters was assessed by properties not involved in the optimization procedure. Selection of Parameters for Optimization. All relevant parameters (bond lengths, bond angles, corresponding force constants, dihedral angles, and Lennard-Jones parameters) were adopted from GAFF.25 Partial charge atom distribution was obtained from the RED server database40,41 as a high-quality multiconformational electrostatic fit. Starting with these parameters, a preliminary set of simulations of liquid TFE at standard conditions were conducted to assess the influence of individual parameters on the targeted properties. Additional criteria were taken into account for the optimized set, namely, keeping the set reasonably small, no modification of abundant atom types, and preservation of compatibility with GAFF and other AMBER force fields. This also involves accepting B
dx.doi.org/10.1021/jp505861b | J. Phys. Chem. B XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry B
Article
grid search algorithm for neighbor searching updated every five steps. The electrostatic interactions and forces were evaluated by the particle mesh Ewald (PME) method with fourth-order interpolation, a Fourier grid spacing of 1.0 nm, and a real-space cutoff of 1.0 nm. The Lennard-Jones potential was treated with a cutoff of 1.4 nm, compensated by a long-range correction to the energy and pressure tensor. The simulation temperature was maintained by a V-rescale thermostat with a time constant of 1 ps−1. The pressure was controlled by isotropic Berendsen coupling with a time constant of 1 ps−1 and compressibility factor of 4.5 × 10−5 bar−1. The total linear momentum of the system was removed every 10 steps. Gas-phase simulations were performed with a stochastic dynamic integrator (sd) with a time step of 2 fs, an inverse friction constant of 1 ps, and constraints on all bond lengths. Electrostatic and Lennard-Jones interactions were evaluated by direct summation without any cutoffs. The total angular and linear momentums of the molecule were removed every 10 steps. Simulations of pure liquid TFE intended for parametrization and rapid routine evaluation of liquid-state properties involved 500 molecules of TFE in a cubic box with approximate edges of 3.9 nm. The simulations started from the equilibrated liquid structure obtained by initial GAFF parameters. The first 2 ns was considered as the equilibration phase, and the remaining 10 ns was used for analysis. In cases in which temperatures and pressures were different from the standard conditions (298.15 K and 1 bar) used in parametrization, the lengths of the equilibration and production phases were chosen to be 5 and 15 ns, respectively. The trajectories were stored every 4000 steps; other observables (energy and its components, volume of the box, and pressure tensor) were logged at 50-step intervals. Final validation of the liquid-state properties of the proposed TFE model was conducted in cubic boxes with approximate edges of 6.7 nm containing 2500 TFE molecules. Of a total simulation time of 60 ns, the last 50 ns was used for analysis. The storage frequency remained the same as for simulations in smaller boxes. Simulations of single-molecule TFE in the gas phase necessary for evaluation of the enthalpy of vaporization were conducted for 40 ns. The first 2 ns was discarded; the rest was used for analysis. Both the trajectory and energies were stored every 50 steps. Simulations of TFE/water mixtures were performed in two different ways. First, rapid scanning for aggregation utilized smaller boxes (approximately 125 nm3). Second, for reliable calculation of KBIs, larger boxes (approximately 500 nm3) were used. The counts of TFE and water molecules for the examined concentrations in both variants are listed in Table 1. Random uniform distribution of TFE and water molecules served as the starting structure for all simulations. Therefore, the first 50 ns was discarded for equilibration, and the analyses were performed on the rest of the trajectory up to 200 ns of total simulation time. In both cases, the trajectories and other observables were stored every 5000 and 100 steps, respectively. Unlike the other simulations, calculation of shear viscosity was conducted by Gromacs, version 4.6.3, in which bug no. 1244 concerning viscosity evaluation has been fixed. These simulations were conducted with boxes containing 2500 TFE molecules for 5 ns, starting from a well-equilibrated liquid structure. Soft-core interactions were used in solvation free energy simulation to remove singularities in the interaction potential
Table 1. Content of Boxes with TFE/Water Mixtures x(TFE)
N(TFE), small
N(water), small
N(TFE), large
N(water), large
0.0365 0.0661 0.113 0.13 0.188 0.271 0.354 0.519 0.645 0.773
137 230 352 390 501 623 715 846 915 970
3616 3249 2763 2608 2163 1675 1305 784 503 285
561 942 1440 1596 2051 2551 2930 3464 3751 3975
14821 13301 11306 10682 8861 6862 5347 3210 2064 1167
during interpolation between states. The soft-core interactions employed the following parameters: sc-alpha = 0.5, sc-r-power = 6, sc-power = 1, and sc-sigma = 0.3. The decoupling of TFE− water interactions was scheduled in 23 steps (λ values), first turning off Coulombic terms followed by Lennard-Jones terms. For the production run, 10 ns was used for each λ value. The data for analysis were collected in 0.2 ps intervals. The system for solvation free energy calculation was composed of 1 TFE molecule and 814 water molecules.
■
RESULTS Ab Initio Survey of the TFE Molecule in the Gas Phase. Table 2 shows the magnitude of the dipole moment, isotropic polarizability, and individual components of the polarizability tensor corresponding to the optimized geometry of the two main conformers of TFE in the gas phase at the B3LYP/aug-cc-pVTZ level. The dipole moments of both conformers differ significantly. The trans conformer possesses a total dipole moment almost 2-fold higher than that of the gauche conformer. However, the calculated polarizability of both conformers differs negligibly. The isotropic portion of the polarizability tensor (35.2 and 35.0 au for the trans and gauche conformers, respectively) and all its components showed remarkably similar values for the trans and gauche conformers (see Table 2). The diagonal components of the polarizability tensors have almost the same magnitude, and the values of offdiagonal components are negligible, suggesting isotropic polarizable character of the TFE molecule. This finding justified usage of the simple term for polarization costs in the gas phase in our study (eqs 6 and 7). The energetic difference found between gauche and trans conformers (8 kJ mol−1) agrees well with the previously reported values.47−51 Parametrization of the Pure Liquid Phase. An initial exhaustive 5-dimensional scan of the parameters involved in optimization was performed for the values listed in Table 3. All combinations of the initial parameters were examined in 1080 simulations. The corresponding densities, self-diffusion coefficients, enthalpies of vaporization, and populations of the trans conformer were analyzed, and their dependencies on the examined force field parameters were fitted in analytic forms using eq 4 in the Appendix. Then five selected force field parameters were optimized to reproduce the fitted properties incorporated into the objective function. The obtained parameter values were repeatedly tested in simulations and used to improve the fits. In total, an additional 1416 points from simulations enhanced the initial regular grid. The selected properties were reproduced with a polynomial formula (eq 4) in the whole parameter space with root-mean-square deviations C
dx.doi.org/10.1021/jp505861b | J. Phys. Chem. B XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry B
Article
Table 2. Electronic Properties of trans and gauche Conformers of TFEa conformer
|μ⃖ | (D)
trans gauche
3.39 1.83
α̅ (au) 35.2 35.0
αxx
αxy
αyy
αxz
αyz
αzz
38.1 36.8
−0.7 −1.3
34.7 34.9
0.0 −0.1
0.0 0.2
33.0 33.3
The magnitude of the dipole moment |μ⃖|, isotropic polarizability α̅ , and individual components of the dipole−dipole polarizability tensor αxx, αxy, αyy, αxz, αyz, and αzz as obtained at the B3LYP/aug-cc-pVTZ level are shown.
a
according to the positions on the μOH−μCF plot of the objective function (Figure 1B). Similarly, the area in parameter space involving the aforementioned minimum on the noncorrected plot (Figure 1A) will be further referenced as the noncorrected (N) region. To confirm our assumption of strong coupling between optimized parameters, we analyzed correlations between them in simulations reproducing the liquid properties of TFE (i.e., values of the objective function of
