AMBER and CHARMM Force Fields Inconsistently ... - ACS Publications

Article Cite This: J. Chem. Theory Comput. 2019, 15, 665−679

pubs.acs.org/JCTC

AMBER and CHARMM Force Fields Inconsistently Portray the Microscopic Details of Phosphorylation Jirí̌ Vymeť al, Veronika Juraś kova,́ † and Jirí̌ Vondraś ě k* Institute of Organic Chemistry and Biochemistry of the Czech Academy of Sciences, Flemingovo náměstí 542/2, 166 10 Praha 6, Czech Republic

Downloaded via BUFFALO STATE on August 2, 2019 at 09:36:51 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Phosphorylation of serine, threonine, and tyrosine is one of the most frequently occurring and crucial post-translational modifications of proteins often associated with important structural and functional changes. We investigated the direct effect of phosphorylation on the intrinsic conformational preferences of amino acids as a potential trigger of larger structural events. We conducted a comparative study of force fields on terminally capped amino acids (dipeptides) as the simplest model for phosphorylation. Our bias-exchange metadynamics simulations revealed that all model dipeptides sampled a great heterogeneity of ensembles affected by introduction of mono- and dianionic phosphate groups. However, the detected changes in populations of backbone conformers and side-chain rotamers did not reveal a strong discriminatory shift in preferences, as could be anticipated for the bulky, charged phosphate group. Furthermore, the AMBER and CHARMM force fields provided inconsistent populations of individual conformers as well as net structural trends upon phosphorylation. Detailed analysis of ensembles revealed competition between hydration and formation of internal hydrogen bonds involving amide hydrogens and the phosphate group. The observed difference in hydration free energy and potential for hydrogen bonding in individual force fields could be attributed to the different partial atomic charges used in each force field and, hence, the different parametrization strategies. Nevertheless, conformational propensities and net structural changes upon phosphorylation are difficult to extract from experimental measurements, and existing experimental data provide limited guidance for force field assessment and further development.

■

INTRODUCTION

At the atomistic level, the effect of phosphorylation is rather complex. The modified residues possess novel physicochemical properties and new potential to form electrostatic interactions such as very strong salt bridges with arginine. 13 A phosphoserine−lysine salt bridge inside a helix is the strongest helix-stabilizing interaction identified to date.14 Phosphoserine is also a strong helix former if located at the N-terminal end of the helix; otherwise, its general effect inside a helix is destabilizing.15 Phosphorylation can trigger disorder−order transitions (e.g., in the ETS-1 domain)16 and decrease conformational entropy,17 and even induce folding of the translational inhibitor 4E-BF2.18 The study of IDPs and their PTMs is a field in which experimental techniques and molecular modeling must complement each other to address research questions.19 A typical protocol for structural interpretation of sparse experimental data involves generation of a plausible structural ensemble that reproduces experimental measures such as chemical shifts, scalar couplings, paramagnetic relaxation enhancement (PRE), nuclear spin relaxation, residual dipolar couplings (RDC), small-angle X-ray scattering (SAXS), and Forster resonance energy transfer (FRET).20−22

Post-translational modifications (PTMs) of proteins enrich the repertoire of 20 amino acids. The functions of both structured and disordered proteins are extensively modulated by PTMs such as phosphorylation, acetylation, methylation, and Nlinked glycosylation.1 Serine/threonine phosphorylation is among the most common PTMs.2 However, other residuessuch as tyrosine, histidine, arginine, and lysineare also frequently modified by O- and N-phosphorylation. An estimated one-third of all proteins in eukaryotic cells, including human cells, undergo phosphorylation, which is mediated by approximately 500 putative kinases.3 Roughly 150 phosphatases are responsible for dephosphorylation.4 A large portion of phosphorylation sites are located in disordered regions,5 and intrinsically disordered regions of proteins (IDPs) are preferred substrates of protein kinases.6 The reversible character of phosphorylation makes this PTM perfectly suited for regulation of protein activity. The allosteric response and conformational changes that occur in the target protein upon covalent binding of a phosphate group can trigger signal transduction pathways, amplification, and integration during the cell cycle, and reactions to external effectors.7 Moreover, specific phosphorylation controls protein turnover,8 localization (e.g., nuclear transport),9 chromatin remodeling,10 and the assembly of protein complexes.11,12 © 2018 American Chemical Society

Received: July 12, 2018 Published: November 23, 2018 665

DOI: 10.1021/acs.jctc.8b00715 J. Chem. Theory Comput. 2019, 15, 665−679

Article

Journal of Chemical Theory and Computation Molecular dynamics simulations with atomistic force fields offer a straightforward approach to generate a structural ensemble under experimental constraints. However, the ability of current force fields to reproduce the properties of IDPs is limited and depends heavily on parametrization.23−27 Additionally, amino acids with PTMs are rarely represented in computational studies of proteins. In the past, parameters for PTM residues were often developed independently of the parental force fields, mostly relying on the transferability of parameters from model compounds. Today, parameters for phosphorylated serine, threonine, and tyrosine are available in all major biomolecular force fields. A recent edition of the AMBER ff9x family (available in AmberTools17)28 contains partial atomic charges and bonded parameters described by Homeyer et al.29 Novel LennardJones parameters of phosphate groups have been adopted for different anionic states to better describe the hydration of bioorganic phosphates. 30 Analogous parameters for the AMBER03 branch are available from Force field_PTM31 with charge distribution calculated according to the approach of Duan et al.32 with optimized torsion parameters. However, optimized Lennard-Jones parameters have not been adopted. In CHARMM (protein force field in versions 27, 36, and 36m),33−35 parameters for phosphorylated amino acids were transferred from the nucleic acid part of the force fields,36 with the exception of phosphotyrosine, which was parametrized separately.37 Gromos, another commonly used biomolecular force field, also supports phosphorylated amino acids. The “43a1p” parameter set adopts force field parameters for mono- and dianionic forms of phosphoserine, phosphothreonine, and phosphotyrosine based on the work of Åqvist et al.38 Alternative parameters for phosphorylated residues can be obtained from the Vienna PTM server,39 which supports versions 45a3,40 54a7,41 and 54a842 of Gromos force fields. Parameters for phosphorylated amino acid residues were transferred directly from analogous groups of nucleotides. The latter two versions (54a7 and 54a8) involve recently optimized parameters for phosphate and other charged PTM residues.43 A recent computational study of phosphorylation and glycosylation (namely, O-GlcNAcylation) of serine and threonine suggested that these PTMs have opposite effects on structural trends.44 The researchers demonstrated these changes in structural preferences by surveying the Protein Data Bank (PDB) and performing quantum chemistry calculations and molecular dynamics of model systems, including amino acid dipeptides, helices, and fragments of the IDP tau. In our study, we scrutinized the conformational preferences of serine, threonine, and tyrosine upon phosphorylation and subsequent (de)protonation of the phosphate group. The latter stems from the fact that the pKa of the phosphate group in the dianionic (-PO32−) and monoanionic (-PO3H−) states lies in the pH range of approximately 5.5−6.5.45 Therefore, both variants may coexist under biologically relevant conditions. We followed an approach established in our previous study that extensively investigated the effect of a force field on the propensities of amino acids to adopt different conformations.46 Analogously, in this work we used capped amino acids (dipeptides) as the simplest molecular models for IDPs. Enhanced-sampling molecular dynamics revealed a vast heterogeneity of conformational ensembles. We analyzed the net effects of phosphorylation on backbone conformers and

side-chain rotamers, as well as their dependence on the particular force fields involved in our study. Additionally, we discuss the reasons underlying and the relevance of the inconsistencies we identified, as well as outlooks for future force field development.

■

METHODS Systems and Force Fields. The model dipeptides were built with tleap utility from the Ambertools package (version 17)28 as peptides of sequence Ace-Xxx-Nme, where Ace stands for the N-terminus-capping acetyl group, Nme for the Cterminus-capping N-methylamide, and Xxx for the amino acid of interest (serine, threonine, tyrosine, or their phosphorylated variants). The following force fields were used in our study: (1) AMBER ff99SB47 extended by parameters for phosphorylated residues developed by Homeyer et al.29 and optimized Lennard-Jones parameters30 (as contributed in the phosaa10 extension in Ambertools; throughout this work, we refer to this force field as “a99”). (2) AMBER ff0332 with parameters for PTM amino acids from Force field_PTM31 (referred to as “a03”; however, this force field lacks parametrization for a singleprotonated phosphate group (monoanion), and therefore, only the fully deprotonated (dianionic) forms of phosphorylated residues were modeled by a03). (3) CHARMM36m35 provided by the MacKerell group as a port for Gromacs (version jul2017; parameters for phosphorylated amino acids were taken directly from the force field, which is referred to in this work as “c36”). Preparation Details. The topologies of peptides for AMBER force fields were generated with tleap and converted to Gromacs topologies using the Parmed tool.48 The CHARMM port was processed directly by Gromacs. All peptides were solvated in cubic boxes (with approximately 3.6 nm edges) by approximately 1,500 water molecules. The recommended water models were used for each force field: TIP3P for a99 and a03 and TIP3Ps for c36. The peptide charges were neutralized by counterions (sodium and chloride ions) in counts providing a 150 mM salt concentration. Simulation Details: Molecular Dynamics. All simulations were performed in Gromacs version 5.1.149,50 enhanced by PLUMED plugin version 2.2.51 All systems were optimized and equilibrated by standard techniques, namely, the steepest descent optimization algorithm followed by short dynamics (1 ns) with restrained peptide atom positions. The production runs used 2 fs time steps allowed by constraining bonds involving hydrogen atoms according to the LINCS algorithm. Periodic boundary conditions were applied, and electrostatic interactions were treated by the particle mesh Ewald (PME) method. Lennard-Jones interactions were truncated behind 10 Å for AMBER force fields. A force-shifting function was applied in the 10−12 Å range for c36, as recommended. A temperature of 298.15 K was maintained by the v-rescale thermostat, and a pressure of 1 bar was controlled by the Parrinello−Rahman barostat. Metadynamics. The peptides were simulated by biasexchange metadynamics52 with five replicas, each 100 ns long. Exchanges between replicas were attempted every 1,000 molecular dynamics steps. The first four replicas were biased, and the last were kept neutral. The biased collective coordinates involved φ, ψ, χ1, and χ2 torsions, each acting 666


Article

Journal of Chemical Theory and Computation

structures with resolution lower than 2.0 Å. This subset consisted of 56,849 structures, which were not filtered in advance for any kind of homology. Each PDB file was parsed, and continuous fragments of seven residues with an amino acid of interest (coded as SER, THR, TYR, SEP, TPO, and PTR) at the central position were extracted. If there were multiple fragments with the same sequence, only one representative with the lowest B-values and complete central residue side chain was counted in the final statistics. This procedure effectively suppressed bias due to multiple occurrences of a single protein chain in single or multiple PDB files, as well as occurrences of sequentially identical regions in homologous proteins. However, it might have introduced a bias toward well-resolved structures. Statistical Tests and Estimation of p-Values. Nonparametric bootstrap analysis58 was performed to test a null hypothesis assuming an equality of means of two populations. A joint distribution of samples was created from all measures (replicas and simulations) under assumption of the null hypothesis. Two groups of N samples were then chosen randomly with replacement from the joint distribution. This procedure was repeated 1,000,000 times, and statistics on differences in means of both randomly sampled groups were collected. The actual difference of means was compared with the synthetic distribution, and the corresponding p-value was evaluated as the probability of the actual (and more extreme) results under assumption of the null hypothesis. The size of a single group was typically 10 (2 × 5 replicas). However, the number of drawn samples (N) used in the test was 2, as the five replicas in a single bias-exchange simulation cannot be considered truly independent measures. Visualizations. Visualizations were created using matplotlib59 and Pymol.60

in an individual replica. The atoms used to calculate the torsion are listed in Table 1. One-dimensional bias potential Table 1. List of Atoms Defining χ1 and χ2 Torsions Used throughout This Study amino acid

χ1

serine phosphoserine threonine phosphothreonine tyrosine phosphotyrosine

N−Cα−Cβ−Oγ

Cα−Cβ−Oγ−P

χ2

N−Cα−Cβ−Cγ

Cα−Cβ−Cγ−Cδ1

was collected in all biased replicas according to the welltempered metadynamics algorithm (BIASFACTOR = 11) in the form of gaussians with an initial height of 0.5 kJ/mol, width (σ parameter) of 0.2 rad, and frequency of 1 ps. All metadynamics simulations were doubled to check convergence of results from two independent simulations and assess their statistical significance. The final analyses were performed on aggregate data. Calculation of Hydration Free Energies. Hydration free energies (HFEs) were calculated using Adaptive Poisson− Boltzmann Solver (APBS) software.53 Nonlinear Poisson− Boltzmann (PB) equations for electrostatic potential were solved on a equidistant grid of 65 × 65 × 65 points using spacing of 0.25 Å, multiple Debye−Huckel sphere boundary conditions, a solvent dielectric constant of 78.4, a solute dielectric constant of 2, a solvent probe radius of 1.4 Å, and an ionic strength of 150 mM. The same settings were used for all molecules. Non-electrostatic contributions to the HFE were neglected. Partial charges and radii of atoms were transferred from each force field. Calculations of HFE were conducted for a representative subset of conformers derived from neutral replicas (800 structures per bias-exchange metadynamics simulation). Analysis of Simulations. The biased trajectories were reweighted by the algorithm described by Bonomi et al.54 to reconstruct the canonical ensemble for further analysis. The first 20 ns were always discarded; only the last 80 ns were analyzed. Tools from the Gromacs package were used for analysis, mainly for calculation of molecular dipoles and interatomic distances. The torsion-based conformational analysis employed backbone conformers defined in our previous study.46 The Ramachandran-like plot is occupied by regions pp (−120° < φ < 0° and 60° < ψ < 180° or −180° < ψ ← 120°), αR (−120° < φ < 0° and −120° < ψ < 60°), αL (0° < φ < 150° and −180° < ψ < 180°), αr (−180° < φ ← 120° or 150° < φ < 180° and −120° < ψ < 60°), and ex (the remaining area). The scalar dipolar couplings were calculated by the Karplus equation parametrized as 3J(HN,HA) = 7.97 cos(2(φ−60°)) − 1.26 cos(φ−60°) + 0.63, where φ refers to the backbone φ torsion.55 Hydrogen bonds were identified according to the combined distance and angle criteria of Baker and Hubbard.56 Namely, the distance between a hydrogen atom and an acceptor atom must be less than 2.5 Å and the angle formed by donor− hydrogen−acceptor atoms greater than 90° for a hydrogen bond to be recognized. Analysis of PDB Structures. We analyzed a subset of the PDB57 (accessed December 2017) that was restricted to X-ray

■

RESULTS The conformational preferences of capped serine, threonine, and tyrosine and the two phosphorylated variants (mono- and dianionic) of each were investigated by metadynamics to assess the effect of phosphorylation on the behavior of the protein constituents. We assessed the convergence of our simulations by the magnitude of instantaneous fluctuations in bias potential over the course of bias-exchange metadynamics (Figures S1−3 in the Supporting Information). After 20 ns of simulation, the fluctuations reduced dramatically, indicating exploration of the whole conformational space of the studied dipeptides. From this point, production data were collected. Another convergence check involved comparison of the populations of major backbone conformers and side-chain rotamers from independent simulations (Figures S4−S6 in the Supporting Information). No major improvement in convergence would be achieved by extending the simulations for several tens of nanoseconds. Effects of Phosphorylation on the Conformational Preferences of Amino Acids. We found that phosphorylation had a significant effect on the net backbone conformation preferences. Figure 1 illustrates the observed changes in the three most populated areas in the Ramachandran plot. Populations modeled by the three different force fields had different frequencies and showed inconsistent shifts upon phosphorylation. Introduction of mono- and dianionic phosphate groups mostly led to trends in the same directions; nevertheless, the magnitudes of the changes showed 667


Article


All χ1 torsions had sharp maxima at −60°, 60°, and 180°, corresponding to the gauche+, gauche−, and trans rotamers, respectively. The location of maxima near the optimal values (−60°, 60°, and 180°) reflects the presence of little or no strain. The only exception was phosphorylated threonine in the a03 force field, in which both gauche maxima shifted noticeably. The effect of phosphorylation appears to be most pronounced for threonine, for which the trans conformer population increases significantly upon phosphorylation. On the other hand, the population of tyrosine rotamers appears nearly unchanged. There is also a clear effect of the force field. For example, c36 prefers a higher population of trans conformer for serine than the AMBER force fields (Figure 2A−C), but penalizes the gauche− conformation for tyrosine (Figure 2G−I). The χ2 torsion populations showed diverse behavior (Figure 3). Tyrosine (Figure 3G−I), the most consistent case, undergoes minimal changes upon phosphorylation and also has the smallest dependence on the force field. Unlike the χ1 torsion, χ2 in (phospho)tyrosine shows 2-fold symmetry with maxima around −90° and 90° due to the symmetry of the adjacent aromatic ring. On the other hand, based on simple stereochemical assumptions, the corresponding phosphoester or hydroxyl oxygen in (phospho)serine and (phospho)threonine could provide qualitatively similar rotational patterns as the χ1 torsion. However, this premise was confirmed only for phosphoserine (Figure 3A−C) and unphosphorylated threonine (Figure 3D−F), which populate recognizable gauche and trans conformers. In all force fields, phosphothreonine yields ambiguous and multimodal distributions in which the canonical rotamers cannot be easily assigned. Because the dependence of side-chain conformations (rotamers) on backbone torsions is well-documented for proteins, we analyzed the coupling between these two properties (Figure S10 in the Supporting Information). For each χ1 rotamer, the corresponding joint histogram of backbone φ and ψ torsions is shown as a separate panel in Figure S10. As expected, the populated regions overlap the sterically allowed areas of the Ramachandran plot. There is good distinction between helical, polyproline-helix-II-like, and extended regions, justifying our previous conformational analysis. Separation of individual χ1 rotamers introduced an additional level of resolution. This analysis confirmed our previous findings: the trans rotamers of (phospho)serine and threonine are depleted in AMBER force fields, and the gauche− rotamer of (phospho)tyrosine is depleted in c36. Moreover, Figure S10 illustrates that different rotamers might favor distinct backbone conformers and vice versa (compare panels A1−A3 and D1−D3 in Figure S10). This finding is even more pronounced when the individual panels are normalized independently (see Figure S11 in the Supporting Information). To analyze the complex form of some χ2 distributions in Figure 3, we prepared joint, two-dimensional χ1/χ2 histograms (Figure S12 in the Supporting Information). The shift of peaks from ideal positions was more pronounced for the phosphorylated side chains, in which the torsional angles of gauche− and trans χ2 conformers approach each other depending on the values of χ1 torsion (see Figure S12H,L−P). These joint χ1/χ2 histograms also illustrate the strong effect of the force field on preferred side-chain conformation. C36 favors the trans/ gauche− rotamer of the (phospho)serine side chain, which is negligibly populated in both AMBER force fields. Similarly, regions on the χ1/χ2 map that are prominent in a99 and a03

Figure 1. Effects of phosphorylation on backbone conformer populations. The propensities of extended, polyproline-helix-II-like (ppII-like), and helical backbone conformers (regions defined in Ramachandran-like plots A−C) are shown for the three force fields. The blue, green, and olive bars represent the unphosphorylated forms; the dotted and dashed bars correspond to the phosphorylated forms (monoanionic and dianionic, respectively). The green and orange arrows demonstrate a shift in preferences upon phosphorylation to the monoanionic and dianionic forms, respectively. The symbols “*” and “N” indicate statistically significant (p-value < 0.01) and nonsignificant changes in preferences, respectively. The monoanionic form is not available for a03; therefore, it is not presented.

inconsistent behavior in the different force fields. In several cases, the monoanionic phosphate group induced stronger effects than the fully deprotonated phosphate anion and vice versa (see Figure 1D−L). Tyrosine exhibited the smallest effect of phosphorylation with low statistical significance in the CHARMM force field; however, the same system modeled by the a03 force field manifested the largest changes detected for backbone conformational preferences (Figure 1J−L). Bootstrap tests were used to estimate the statistical significance of individual population shifts in Figure 1. Most of the changes were statistically significant with associated pvalue < 10−6. In other cases, the corresponding p-values are listed in Table S7 in the Supporting Information. Introduction of the bulky phosphate group is expected to change the conformational preferences of the side chains, measured as populations of individual rotamers. For all of our model dipeptides, the conformation of the side chain is described mainly by two adjacent torsions, χ1 and χ2. Distributions of torsion values obtained from simulations are plotted as histograms in Figures 2 and 3. 668


Article


Figure 2. Distribution of χ1 torsions. Histograms show the distribution of the χ1 torsion calculated from simulations with a99 (A, D, G), a03 (B, E, H), and c36 (C, F, I) for unphosphorylated and phosphorylated variants of serine (A, B, C), threonine (D, E, F), and tyrosine (G, H, I). The unphosphorylated, monoanionic,, and dianionic forms are plotted in blue, green, and orange, respectively. The monoanionic form is unavailable for a03 and therefore is not presented in panels B, E, and H. Each panel is normalized arbitrarily. The uncertainty of the plots is comparable to the width of the lines and is not shown for the sake of clarity. All differences between profiles were statistically significant (p-value < 0.01) except those for phosphorylated tyrosine in a99 (panel G). The corresponding p-values are listed in Table S8 in the Supporting Information.

Figure 3. Distribution of χ2 torsions. Histograms show the distribution of χ2 torsions calculated from simulations with a99 (panels A, D, G), a03 (B, E, H), and c36 (C, F, I) for unphosphorylated and phosphorylated variants of serine (A, B, C), threonine (D, E, F), and tyrosine (G, H, I). The unphosphorylated, monoanionic, and dianionic forms are plotted in blue, green, and orange, respectively. The monoanionic form is unavailable for a03, and therefore is not presented in panels B, E, and H. Each panel is normalized arbitrarily. The uncertainty of the plots is comparable to the width of the lines and is not shown for the sake of clarity. All profiles for serine and threonine in individual force fields are significantly different with p-values < 10−6. On the other hand, the differences upon phosphorylation of tyrosine are not statistically significant (p-value > 0.01). The corresponding p-values are listed in Table S9 in the Supporting Information.

The atoms defining χ1 and χ2 torsions do not involve any atoms of the phosphate group attached on the opposite side of the aromatic ring. We achieved the most detailed view of the propensity of amino acids to adopt different conformations by involving all four torsions (φ, ψ, χ1, and χ2) in analysis. This approach

are often underrepresented in c36. A similar conclusion can be drawn for threonine. While phosphothreonine populates a limited number of basins, minor discrepancies between force fields can be identified. The χ1/χ2 plot did not provide any additional insights for tyrosine and phosphotyrosine (data not shown) due to the spatial separation of the phosphate group. 669


Article


Figure 4. Change in propensities of amino acids to populate individual conformers upon phosphorylation. Correlations between populations of unphosphorylated and phosphorylated variants of individual conformers are shown for serine, threonine, and tyrosine (in rows) simulated by a99, a03, and c36 (in columns). The Pearson correlation coefficient (r) between unmodified and phosphorylated populations is listed in each plot. Monoanionic and dianionic forms are plotted in green and orange, respectively. Dotted lines indicate a perfect match. More details for the most frequently occurring conformers (circled) are provided in Table 2. An alternative graphical version of Table 2 picturing the conformers is included in the Supporting Information (Table S13).

tions complement our previous findings that the populations of conformers in the two AMBER force fields generally correlate better with each other (average Pearson’s correlation coefficient of 0.75) than with the c36 force field (average correlation 0.3). Serine, threonine, and their phosphorylated variants have the potential to form intramolecular hydrogen bonds between their side-chain and backbone moieties, which could have a significant effect on the stability of particular conformers. We focused primarily on hydrogen bonds introduced by the phosphate group, as shown schematically in Figure 5A. Hydrogen bonds were detected separately for amide hydrogens at the acetylated or N-methylated termini (Hac and Ham, respectively). Similarly, the terminal oxygens of the phosphate group (emphasized in Figure 5A) and the γ oxygen participating in the phosphoester bond were treated independently. An overall view of the interactions between these moieties is presented in Figure S15 in the Supporting Information. The terminal oxygens and both amide hydrogens (Figure S15A,B) had preferred distances and orientations complying well with hydrogen-bonding criteria: distance between hydrogen and oxygen less than 2.5 Å and angle between amide nitrogen, hydrogen, and oxygen higher than 90°. In fact, the second condition was fulfilled automatically in all cases if the distance criterion was accepted. On the other hand, phosphoester oxygens did not provide distributions

resulted in 45 recognizable conformers for (phospho)serine and (phospho)threonine and 15 conformers for (phospho)tyrosine (because of symmetry in χ2). The populations of these conformers before and after phosphorylation helped reveal the character of the ensembles sampled in our simulations, as shown in the correlation plots in Figure 4. These plots also show that several conformers dominate over other structures in the individual ensembles. Although the most favored conformer reached frequencies up to 50%, dozens of other conformers were populated in non-negligible amounts. Upon phosphorylation, a large shift in the population of particular conformers was observed for serine and threonine in all force fields. For tyrosine, a99 and c36 provided ensembles with more subtle shifts than those seen in the a03 force field. The three most frequently occurring conformers for each force field are listed in Table 2. Very often, the leading conformers kept their positions upon phosphorylation, with the notable exception of the (ex,t,g-) conformer massively populated by phosphothreonine in both AMBER force fields. Additionally, the list of prominent conformers in the AMBER force fields largely overlaps, but much less similarity was detected between c36 and AMBER (see Table 2). To investigate these differences between force fields, we performed an analogous analysis of the frequencies of individual conformers for all residues in a99, a03, and c36 (Figure S14 in the Supporting Information). These observa670


Article

Journal of Chemical Theory and Computation Table 2. Top Three Conformers Populated in Individual Force Fields amino acid Ser

Ser-PO3H− Ser-PO32−

Thr

Thr-PO3H− Thr-PO32−

Tyr

Tyr-PO3H− Tyr-PO32−

force field

rank no. 1 conformer



a99 a03 c36 a99 c36 a99 a03 c36 a99 a03 c36 a99 c36 a99 a03 c36 a99 a03 c36 a99 c36 a99 a03 c36

(pp, g+, t) [13%] (αR, g+, t) [20%] (pp, t, g-) [14%] (αR, g+, t) [25%] (pp, t, g-) [26%] (αR, g+, t) [37%] (αR, g+, t) [21%] (pp, t, g-) [24%] (αR, g-, t) [21%] (αR, g-, t) [15%] (pp, t, g-) [13%] (ex, t, g-) [31%] (pp, g+, t) [19%] (ex, t, g-) [26%] (ex, t, g-) [46%] (pp, t, g-) [34%] (pp,t)[25%] (αR, t) [29%] (pp, t) [29%] (pp, t) [39%] (pp, g+) [33%] (pp, t) [29%] (pp, t) [33%] (pp, g+) [37%]

(αR, g+, t) [10%] (αR, g-, t) [15%] (ex, t, g-) [9%] (αR, g-, t) [18%] (pp, g+, g+) [21%] (αR, g-, t) [14%] (αR, g-, t) [20%] (pp, g+, g+) [15%] (pp, g-, t) [10%] (pp, g+, t) [13%] (pp, g-, g-) [12%] (αR, g-, t) [18%] (pp, g-, t) [19%] (αR, g-, t) [22%] (αR, g-, g-)[15%] (pp, g+, g-) [16%] (ex, t) [17%] (pp, t) [17%] (pp, g+) [29%] (ex, t) [14%] (pp, t) [26%] (ex, t) [19%] (pp, g+) [18%] (pp, t) [23%]

(αR, g-, t) [7%] (αR, g+, g-) [9%] (pp, g+, g-) [8%] (pp, g+, t) [2%] (ex, t, g-) [9%] (pp, g+, t) [10%] (pp, g+, t) [18%] (pp, g+, t) [11%] (αr, g-, t) [8%] (pp,g-,t) [10%] (ex, t, g-) [12%] (pp, t, g-) [10%] (pp, t, g-) [12%] (pp, t, g-) [8%] (αR, g+, g-) [14%] (pp, g+, t) [14%] (pp, g+) [10%] (αR, g+) [12%] (αR, g+) [11%] (pp, g+) [8%] (αR, g+) [10%] (pp, g+) [8%] (ex, t) [11%] (αR, g+) [11%]

The conformers are listed in the following formats: backbone state, χ1, χ2 for serine and threonine; backbone state, χ1 for tyrosine. The number in square brackets indicates the percentage of molecules adopting the conformation. Uncertainty is below ±2%. a

We inspected the individual phosphoserine and phosphothreonine conformers sampled in our simulations for their ability to form intramolecular hydrogen bonds. The most prominent conformers with Hac participating in hydrogen bonding were as follows: (αR, g+, g-), (αR, g-, g+), (pp, g+, g), and (pp, g-, g+), according to the shorthanded notation “(backbone state, χ1, χ2).” Analogously, hydrogen bonds with Ham were allowed mostly for conformers (ex, t, g-) and (pp, t, g-). A detailed list of prominent conformers with an internal hydrogen bond is provided in Table S18 in Supporting Information. Comparison of hydrogen-bonding capabilities (Table S18), populations of conformers (Table 2), and observed hydrogen bonds (Figures 5 and S16) is not straightforward. For phosphoserine analyzed in AMBER force fields, the top three conformers do not show hydrogen-bonding tendencies, in agreement with Figure 5B,D. On the other hand, we identified a higher number of Ham-bonding conformers than Hac-bonding conformers at the top positions with the greatest impact on the ensemble. This is in apparent contrast with the findings shown in Figure 5E−G, which suggest a comparable frequency of Hacbonding. We conclude that hydrogen-bond-compatible conformers do not form hydrogen bonds all the time, the frequency of internal (Ham or Hac) hydrogen-bonding is small in general (at most ∼50% for phosphothreonine in a03 and typically half of this value or less; see Figure S16), and there is a non-negligible contribution from less frequent conformers that are not listed in Table 2. As internal hydrogen bonding does not seem to be the major determinant of the structural ensemble, we investigated the hydration free energy (HFE) of representative samples. Distribution of HFE for individual dipeptides and force fields is shown in Figure S19 in the Supporting Information. All

preferring optimal hydrogen-bond geometries, although they were partially compatible with the criteria used. This indicates the presence of strong geometric constraints that prohibit formation of optimal hydrogen bonds with strong energetic contributions. Therefore, hydrogen bonds involving the phosphoester oxygen were neglected in further analysis. This behavior was observed for all individual phosphopeptides and force fields (see Figures S16 and S17 in the Supporting Information). A compact version of the complete analysis for terminal phosphate oxygens is presented in Figure 5B−G. In both AMBER force fields (Figure 5B,D), we observed little or no formation of hydrogen bonds in the phosphoserine molecule. Despite the low overall presence of hydrogen bonds, participation of Hac was more frequent than Ham. In comparison, in the c36 force field, phosphoserine showed a much greater propensity to form hydrogen bonds (Figure 5F). In particular, we observed increased participation of Hac in the dianionic form of phosphoserine, but other hydrogen bonds involving Ham also formed significantly more frequently in than in the AMBER force fields. In contrast to phosphoserine, phosphothreonine populated more conformers fulfilling the criteria for hydrogen bonding of both amide hydrogens in a99 (Figure 5C). The amide hydrogens of phosphothreonine showed the strongest tendencies to interact with the phosphate group in a03 (for the dianionic form, Figure 5E). In c36, the dianionic form formed both types of hydrogen bonds (those with Hac and those with Ham) more frequently than the monoanionic form (Figure 5G). On the other hand, the most preferred binding motif in a99 involved Ham and the monoanionic molecule (Figure 5C). Phosphotyrosine did not form any hydrogen bonds between the phosphate group and amide backbone in any of our simulations (data not shown). 671


Article


more similarity with each other than with (phospho)tyrosine. Ff03 systematically produced the lowest (more favorable) hydration energies for serine and threonine and highest for tyrosine. C36 provided the lowest HFE for dianionic species, but this trend was not reproduced for monoanionic molecules. An analysis presented in Figure 6 focuses on phosphoserine and phosphothreonine, relating HFE of the three most frequent conformers (listed in Table 2) and hydrogen-bonded conformers. Figure 6 provides important insight into the impact of internal hydrogen bonding on HFE. In all cases, the Hac- and Ham-bonded structures were penalized and populated in the less favorable tail of the total histogram (most obvious in plots B, E, F, G, H, and I of Figure 6). Furthermore, the three most prominent conformers do not contribute preferentially to the energetically favorable tail of HFE distribution; rather they mostly follow the main peak of the distribution (except ThrPO3H− in the c36 force field, Figure 6G). Differences in HFE in individual force fields likely originate from different charge distributions modeled by partial atomic charges. To test this hypothesis, we investigated the net dipole moments of the dipeptides. The net dipole moment is the second most important term in expansion of molecular electrostatic potential, and unlike a net charge (the most prominent term), it is conformation-dependent. Distributions of net dipole moments are shown in Figure S20 in the Supporting Information. The plots in Figure S20 reveal a strong force field effect on molecular dipoles, similarly to those for HFE. Furthermore, both dipole and HFE distributions mostly share similar but inverted features (the larger the dipole, the lower the HFE). Correlations of both properties for all systems are presented in Figure S21 in the Supporting Information. In general, the measured correlations were weaker for neutral species (starting at Pearson’s r = −0.45 for threonine in a99, Figure S21I) and stronger for charged molecules (up to r = −0.91 for phosphothreonine, Figure S21N). In general, the partial atomic charges of individual dipeptides in the three force fields correlate very well (Pearson’s correlation coefficient 0.92 and higher; see Figure S22 in the Supporting Information), but large deviations in the charges of particular atoms can be observed. The most prominent cases are CZ in Tyr-PO32− and CB in Thr-PO32−, which differ more than 0.7e in charge between a99 and c36. For the sake of completeness, all partial charges are listed in Tables S23−S31 in the Supporting Information. Comparison with Experimental Observations. NMR analysis of capped phosphoserine and phosphothreonine at various pHs suggested that these molecules have different internal hydrogen-bonding potencies.61 According to these data, the dianionic phosphate group of phosphothreonine can bind both backbone hydrogens, whereas only the amide hydrogen of the acetylated terminus (Hac) is accessible in phosphoserine. Such strict discrimination was not observed in our simulations, although some of our observations agree partially with these experimental conclusions (see Figure 5). We compiled conformation-sensitive scalar dipolar couplings 3 J(HN,HA) from experimental work led by Avbelj et al.62 and Kim et al.45 The coupling constant values predicted from simulations were calculated according to the Karplus equation and compared with the experimental data (Figure 7). There was reasonable numerical agreement between the calculated coupling constants from the a99 force field and the measured values, except for the dianionic forms (Figure 7A). However,

Figure 5. Analysis of intramolecular hydrogen bonding. A schematic drawing of the phosphoserine dipeptide with hydrogen-bond donors and acceptors circled is shown in panel A. The amide hydrogens on the acetyl end (Hac) and N-methyl (Ham) capping group are shown explicitly. Internal hydrogen bonding is indicated by histograms of distances between terminal phosphate oxygens and Hac (blue and green curves for monoanionic and dianionic forms, respectively) and Ham (orange and red curves for monoanionic and dianionic forms, respectively) in panels B−G. Data for monoanionic forms are unavailable for a03; therefore, the corresponding distance distributions are not plotted in panels D and E. The dotted lines represent the 2.5 Å cutoff distance for hydrogen bonding. The y-axis units in panels B−G are arbitrary but the same for all panels. The insets in panels B and D magnify the sparsely populated regions. The uncertainties are not shown for the sake of clarity.

distributions had the form of bell-shaped curves with one prominent peak, and they were not strongly biased toward favorable (more negative) values. The HFE of neutral, singlecharged, and double-charged molecules differed in magnitude, reflecting the highest contribution from solvation of an electric charge. Nevertheless, the distributions of HFE corresponding to individual force fields revealed significant shifts. We observed several trends and inconsistencies in the behavior of the force fields. In particular, the shifts of distributions for (phospho)serine and (phospho)threonine dipeptides share 672


Article


Figure 6. Hydration free energy (HFE) of phosphoserine and phosphothreonine in mono- and dianionic forms (in rows) in the a99, a03, and c36 force fields (in columns). The black line delineates the distribution of the entire ensemble. The dotted lines show the contributions of the three most frequent conformers in a cumulative (non-overlapping) fashion. The filled area highlights the fraction of conformers with internal bonding of amide hydrogens (Hac and Ham). The frequency on the y-axis is in arbitrary units for each plot. Vacancies are due to the unavailability of monoanionic forms in a03.

Distribution of backbone torsions φ and ψ varies significantly between the libraries, but rotameric torsions show more consistent features. Distribution of χ1 for phosphoserine follows the distribution for serine, both slightly preferring the gauche− rotamer over the others. On the other hand, phosphothreonine preferentially adopts the gauche+ rotamer, a noticeable shift in its conformational preferences upon phosphorylation. Preferences for these particular rotamers were not observed in our simulations (compare with Figure 2); however, all force fields predicted dramatic changes upon phosphorylation. Tyrosine did not show any major changes upon phosphorylation, and the χ1 distributions mostly resemble those produced by the c36 force field. Our PDB analysis indicated that phosphoserine strongly prefers the trans conformer for χ2, as predicted by AMBER force fields (see Figure 3). Analogously, phosphothreonine shows a distinct peak around 120° (Figure 8H), in accordance with force field predictions.

the trends predicted by a99 for nonphosphorylated and monoanionic and dianionic phosphorylated forms of serine and threonine were in opposition to the experimental findings. On the other hand, the experimentally observed trends were reflected by the a03 and c36 force fields (Figure 7B,C, respectively). However, the predicted values provided by c36 were systematically underestimated. Based on the limited amount of experimental evidence and weak agreement with experimental data in general, no ranking of force fields was attempted. No force field in our study substantially outperformed the others when compared with experimentally observed hydrogen-bond selectivity and measured 3J(HN,HA) dipolar couplings. Comparison with PDB Structures. To provide insight into the behavior of phosphorylated residues within proteins, we compiled conformational preferences obtained from structures deposited in the PDB. A direct comparison between the PDB and simulations of model peptides is complicated by the limited number of well-resolved phosphorylated residues in the PDB, inability to distinguish protonation of the phosphate group, bias from many sources, and the unspecified effect of other residues in the protein chain. Nevertheless, our PDB analysis highlighted some general trends. Distributions of φ, ψ, χ1, and χ2 torsions gleaned from the PDB data set are plotted in Figures 8 and S32−S36 for three structural libraries (all residues, residues outside helices and sheets, and regions without any recognizable structural motif).

■

DISCUSSION In this work, we studied the direct effect of phosphorylation on the intrinsic conformational and hydrogen-bonding preferences of the three most important residues that undergo in vivo phosphorylation: serine, threonine, and tyrosine. To avoid force-field bias, we employed three different force fields for biomolecular simulations: a99, a03, and c36. There have been continual efforts to increase the reliability of force fields for 673


Article


fitting strategy described by Duan et al.,32 which is used both in the original a03 and Force field_PTM parameters. Unlike in the other force fields (a99 and c36), the partial atomic charges on backbone atoms were not constrained, but rather optimized during a two-stage RESP procedure.66 Therefore, the partial charges on the backbone are sensitive to the different electronic structures of the side chains. We previously demonstrated that partial charges on backbone atoms contribute significantly to the backbone conformational preferences of amino acids in AMBER force fields.46 The mobility of the phosphate group in phosphoserine and phosphotyrosine is highly restricted, unlike the hydroxyl group in serine and threonine. Therefore, the resulting profiles along the χ2 torsion show a strong preference for particular values (Figures 3 and S12). However, this limitation does not directly prevent population of any other backbone conformer or χ1 rotamer. Indeed, the steric hindrance due to the phosphate group contributes to the energetics, but it does not seem to have a prohibitive effect on the majority of leading conformers in the simulated ensembles (Figure 4 and Table 2). The actual effects stabilizing particular leading conformers are difficult to determine and are likely force-field-dependent. For example, hydrogen bonding between phosphate group oxygen atoms and amide hydrogen atoms plays an unimportant role in phosphoserine modeled by a99 and a03 (Figure 5B,D), suggesting that these atoms are preferentially solvated by water molecules. On the other hand, phosphoserine behaves differently in c36, and we observed more frequent internal hydrogen bonds. In all force fields, phosphothreonine has a tendency to form more internally bonded conformers than phosphoserine. In the literature, this effect has been explained by partial descreening of the phosphate group by the extra methyl group of the threonine side chain. The stronger effective charge of the partially desolvated phosphate group consequently helps form stronger bonds with available hydrogen-bond donors.67 The competition between hydrogen bonding and solvation is illustrated in Figure 6. In distributions of HFEs, the fraction of conformations stabilized by internal hydrogen bonds always occupies the energetically unfavorable region. Even the most frequently occurring conformers were not predominant in the most favorable HFE regions. These observations indicate that neither of these opposing forces completely prevails over the other. Additionally, the fact that individual conformers were not entirely determined by the presence or absence of a particular type of internal hydrogen bond suggests that other internal terms in the force fields, such as torsion potential, play an important role. The range of HFE differs between force fields (Figure S19), drawing attention to the balance between solvent−solute and solvent−solvent interactions. Indeed, the delicate balance between hydration, internal hydrogen bonding, and torsion potential depends on the force field parameters. The calculated HFEs of individual structures correlated significantly with their molecular dipole moments (Figure S21), indicating that the partial atomic charges are of particular importance. Furthermore, hydrogen bonds in AMBER and CHARMM force fields are modeled mainly by electrostatics, and both force-field families also include electrostatic terms as part of the torsion (1−4) interaction. Thus, the fact that a99, a03, and c36 obtained partial atomic charges in different ways might explain the source of the observed inconsistencies between the force fields. Although we

Figure 7. Comparison of experimental and calculated 3J(HN,HA) dipolar couplings. The experimental 3J(HN,HA) values (x-axis) are plotted against their counterparts calculated from simulations (y-axis) in the a99 (A), a03 (B), and c96 (C) force fields. The dotted lines represent perfect correlation between the experimental and simulated results. The missing data points for a03 are due to the unavailability of parameters for monoanionic phosphorylated forms.

extraordinary protein forms such as IDPs, resulting in various optimized parameter sets.63−65 Because the parameters for phosphorylated residues predate these recent efforts to improve force-field balance, we deliberately avoided any previously unreported, untested, or unintended combinations of parameters, particularly in the AMBER force-field family. Unsurprisingly, we found significant differences in the intrinsic conformational preferences of serine, threonine, and tyrosine across the force fields. This finding was expected in light of results from our previous study that scrutinized the preferences of all basic proteinogenic amino acids in four different force fields.46 Therefore, instead of evaluating absolute values of propensities for different conformations, we focused more on general trends caused by phosphorylation and the introduction of a negatively charged bulky substituent. We originally hypothesized that such a change in the amino acid structure may be accompanied by a very strong shift in its intrinsic conformational preferences, overriding differences in the force fields. The conformational preferences of serine, threonine, and tyrosine responded differently to phosphorylation. In general, we observed the subtlest changes for tyrosine (see Figures 1−4); however, in the a03 force field, the preferences of tyrosine were influenced markedly (see Figure 1K and Figure 4H). We speculate that these changes are due to the charge 674


Article


Figure 8. Distribution of side-chain torsions in PDB structures. Panels A−I show distributions from the whole data set filtered only for identity of sequences; i.e., helices, sheets, and turns are all included in the statistics. Phosphorylated amino acids are labeled according to PDB nomenclature: SEP (phosphoserine), TPO (phosphothreonine), and PTR (phosphotyrosine). Data from the PDB do not distinguish the individual protonation states of the phosphate group. The distributions of χ2 torsion are not shown for unmodified serine and threonine because PDB does not record the position of hydrogen atoms. The distribution of χ2 for tyrosine was omitted for its high symmetry and low information content.

include all these effects, studies on longer sequences are necessary, but the number of relevant sequences grows rapidly with peptide length. However, the intrinsic conformational preferences of individual amino acid residues generate the primary energy landscape of a polypeptide, which is further modified by inter-residue interactions. Therefore, we conclude that correct description of the intrinsic conformational preferences of amino acids is important for the performance of a force field in situations in which inter-residue interactions are weak or transient, including IDPs and unfolded ensembles of structured proteins. Unfortunately, no experimental method provides unambiguous and straightforward insights into the detailed conformational state of a single amino acid residue. Although NMR chemical shifts and scalar couplings, Fourier transform infrared spectroscopy (FTIR), and vibrational circular dichroism (VCD) provide valuable pieces of information about particular structural details, they cannot match the atomic resolution obtained from molecular modeling. Agreement between predicted and measured observables is the desired output of a validation process, but assessment of a force field in cases of numerical disagreement with observed values remains challenging. Similarly, it is difficult to obtain feedback for direct improvement of the force field. Molecular modeling results also can be compared with distributions derived from experimentally resolved 3D protein structures. However, extracting the intrinsic conformational preferences of individual amino acids is a nontrivial task due to the ubiquitous inter-residue interactions. The closest approximation of a single isolated residue can be found in carefully constructed coil libraries,77,78 which have been used to

observed strong overall correlation of partial charges among the force fields (Figure S22), the differences between individual atoms were as large as 0.7e (Tables S23−S31). On the other hand, only the net electrostatic distribution has a physical meaning, and the distribution of partial charges on individual atoms is ambiguous. However, these assumptions are not valid for one to four electrostatic interactions, which are highly sensitive to particular values of atomic charges and contribute significantly to conformation energies. Phosphorylation can induce large conformational changes in both structured proteins and IDPs. Many of these changes have been successfully investigated by computational methods employing the force fields compared in this work. This includes stabilization of transient helices, disordered parts of the sodium−proton exchanger,68 and the C-terminal part of the polybasic domain of ARNO;69 formation of an “arginine claw” motif in desmoplakin;70 long-distance rearrangements in STAT5 signaling proteins;71 modulation stability of the p53TAD:TAZ2 complex;72 allosteric activation of Aurora kinase A;73 partial unfolding of the CDK inhibitor p19;74 and order−disorder transitions and aggregation of τ peptides.75,76 The success of the force fields at reproducing the conformational changes induced by phosphorylation suggests a more robust mechanism than a simple change in the intrinsic preferences of the phosphorylated residue. Such mechanisms may originate from interactions of the phosphate group in salt bridges or analogous charge-driven interactions, which may sequentially stabilize (or destabilize) transient structural elements in ensembles of IDPs. The strength of these effects may override the underlying intrinsic conformational preferences investigated in this study. To properly 675


Article

Journal of Chemical Theory and Computation calibrate force fields with promising results.79 On the other hand, well-resolved structures of PTM residues are sparse in databases, and the repertoire of proteins represented is biased toward those of research interest. This approach may therefore have limited applicability for PTM residues. Most of the recently reported improvements in force fields have relied on quantum chemical calculations. Beginning with reparametrization of backbone torsions in a9947 and the CMAP correction for CHARMM,80 researchers recognized a necessity for correction of side-chain torsions. Improved rotamer distributions were targeted by the ILDN correction81 for a99 and in the c36 force field.33 Our results reveal that various combinations of backbone and side-chain conformers are populated in solution. This highlights a need for a thorough and balanced description of the potential energy surface that correctly captures various couplings between the backbone and side-chain degrees of freedom. Such an elaborate reparametrization has been attempted in the AMBER force field family, resulting in ff15ipq82 and the “force-balanced” fb15 variant.83 Development of both state-of-the-art force fields required a large amount of ab initio reference calculations (MP2 level of theory) and extensive optimizations. However, these force fields still do not agree quantitatively on populations of sampled conformers for serine, threonine, and tyrosine (see Figure S37). This example illustrates the ambiguity in parametrization schemes, as well as our insufficient understanding of the condensed phase effects that must be included implicitly during construction of a fixed-point-charge force field. These effects are likely even more critical for polarizable anionic species such as phosphorylated amino acids. The fact that these force fields perform well on complex systems might seem surprising in light of our study. This may indicate either that small details such as intrinsic conformational preferences are irrelevant in large polypeptides or, more likely, that the well-calibrated strength of stability- and foldingdriving interactions compensate for imprecise terms. On the other hand, the inconsistently described microscopic details might be partly responsible for the suboptimal performance of force fields on IDPs or denatured proteins. More research is needed to address these questions. Studies on longer peptides might shed more light on this topic, especially if they are supported by detailed and conclusive experimental evidence.

Detailed analysis of the simulations revealed that a complicated balance of hydration, hydrogen bonding, and likely other force field terms (such as torsions) influence populations of individual structures in the sampled ensembles. Therefore, partial atomic charges are considered responsible for at least part of the observed inconsistencies. No force field showed satisfying agreement with experimental observations or with conclusions drawn from experimental data. In addition to suboptimal reproduction of experimental observations, our finding that individual force fields inconsistently portray the microscopic details of phosphorylation reveals their current limits and exposes areas for further improvements. However, no progress is feasible without detailed high-quality reference data.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.8b00715.

■

Additional simulations information; figures showing bias potential instantaneous changes, calculated population differences, Ramachandran-like plots, joint χ1/χ2 distributions, population correlations, various statistical plots, HFE and dipole moment distributions, various correlations, and backbone torsion distributions; and tables listing estimated p-values, top three conformer characteristics, subensemble compositions, and atomic charges (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Jiří Vymětal: 0000-0002-0165-8707 Jiří Vondrásě k: 0000-0002-6066-973X Present Address †

Swiss Federal Institute of Technology in Lausanne, BCH 5121, Ave. F.-A. Forel 2, CH-1015 Lausanne, Switzerland. Funding

The work was supported by the ERDF/ESF project “Chemical biology for drugging undruggable targets (ChemBioDrug)” (Grant No. CZ.02.1.01/0.0/0.0/16_019/0000729).

■

CONCLUSION We studied the effect of phosphorylation on the conformational preferences of serine, threonine, and tyrosine in the simplest model of an IDP consisting of a single capped amino acid (dipeptide). We analyzed ensembles sampled by a99, a03, and c36 force fields accompanied by parameters for phosphorylated residues, either bundled with a particular force field or developed in a compatible manner. Metadynamics simulations generated diverse structural ensembles in which we detected quantitative changes of different magnitudes in the populations of individual conformers upon phosphorylation. However, the ensembles rendered by individual force fields differed significantly, hindering analysis of the phosphorylation effects themselves. General trends extracted from our simulations revealed that the conformational preferences are at least somewhat influenced by phosphorylation. The strongest effect was observed for introduction of a phosphate group to the Cβ atom of threonine. These trends were confirmed by our PDB survey.

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Access to computing and storage facilities owned by parties and projects contributing to the National Grid Infrastructure MetaCentrum, provided under the program “Projects of Large Research, Development, and Innovations Infrastructures” (CESNET LM2015042) is greatly appreciated. Computational resources were provided by the ELIXIR-CZ project (LM2015047), part of the international ELIXIR infrastructure. Computations were performed on resources provided by UNINETT Sigma2the National Infrastructure for High Performance Computing and Data Storage in Norway.

■

ABBREVIATIONS IDP, intrinsically disordered protein; PDB, Protein Data Bank; PTM, post-translational modification; PTR, phosphotyrosine; 676


Article


(19) Burger, V. M.; Gurry, T.; Stultz, C. M. Intrinsically Disordered Proteins: Where Computation Meets Experiment. Polymers (Basel, Switz.) 2014, 6 (10), 2684−2719. (20) Jensen, M. R.; Zweckstetter, M.; Huang, J.; Blackledge, M. Exploring Free-Energy Landscapes of Intrinsically Disordered Proteins at Atomic Resolution Using NMR Spectroscopy. Chem. Rev. 2014, 114 (13), 6632−6660. (21) Kikhney, A. G.; Svergun, D. I. A Practical Guide to Small Angle X-Ray Scattering (SAXS) of Flexible and Intrinsically Disordered Proteins. FEBS Lett. 2015, 589 (19), 2570−2577. (22) Schuler, B.; Soranno, A.; Hofmann, H.; Nettels, D. SingleMolecule FRET Spectroscopy and the Polymer Physics of Unfolded and Intrinsically Disordered Proteins. Annu. Rev. Biophys. 2016, 45 (1), 207−231. (23) Palazzesi, F.; Prakash, M. K.; Bonomi, M.; Barducci, A. Accuracy of Current All-Atom Force-Fields in Modeling Protein Disordered States. J. Chem. Theory Comput. 2015, 11 (1), 2−7. (24) Rauscher, S.; Gapsys, V.; Gajda, M. J.; Zweckstetter, M.; De Groot, B. L.; Grubmüller, H. Structural Ensembles of Intrinsically Disordered Proteins Depend Strongly on Force Field: A Comparison to Experiment. J. Chem. Theory Comput. 2015, 11 (11), 5513−5524. (25) Gao, Y.; Zhang, C.; Zhang, J. Z. H.; Mei, Y. Evaluation of the Coupled Two-Dimensional Main Chain Torsional Potential in Modeling Intrinsically Disordered Proteins. J. Chem. Inf. Model. 2017, 57 (2), 267−274. (26) Henriques, J.; Cragnell, C.; Skepö, M. Molecular Dynamics Simulations of Intrinsically Disordered Proteins: Force Field Evaluation and Comparison with Experiment. J. Chem. Theory Comput. 2015, 11 (7), 3420−3431. (27) Chong, S.-H.; Chatterjee, P.; Ham, S. Computer Simulations of Intrinsically Disordered Proteins. Annu. Rev. Phys. Chem. 2017, 68 (1), 117−134. (28) Case, D.; Cerutti, D. S.; Cheatham, T.; Darden, T.; Duke, R.; Giese, T. J.; Gohlke, H.; Götz, A.; Greene, D.; Homeyer, N.; et al. Amber 2017; University of California: San Francisco, 2017. (29) Homeyer, N.; Horn, A. H. C.; Lanig, H.; Sticht, H. AMBER Force-Field Parameters for Phosphorylated Amino Acids in Different Protonation States: Phosphoserine, Phosphothreonine, Phosphotyrosine, and Phosphohistidine. J. Mol. Model. 2006, 12 (3), 281−289. (30) Steinbrecher, T.; Latzer, J.; Case, D. A. Revised AMBER Parameters for Bioorganic Phosphates. J. Chem. Theory Comput. 2012, 8 (11), 4405−4412. (31) Khoury, G. A.; Thompson, J. P.; Smadbeck, J.; Kieslich, C. A.; Floudas, C. A. Forcefield_PTM : Ab Initio Charge and AMBER Force Field Parameters for Frequently Occurring Post-Translational Modifications. J. Chem. Theory Comput. 2013, 9, 5653−5674. (32) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; et al. A PointCharge Force Field for Molecular Mechanics Simulations of Proteins Based on Condensed-Phase Quantum Mechanical Calculations. J. Comput. Chem. 2003, 24 (16), 1999−2012. (33) Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E. M.; Mittal, J.; Feig, M.; MacKerell, A. D. Optimization of the Additive CHARMM AllAtom Protein Force Field Targeting Improved Sampling of the Backbone ϕ, ψ and Side-Chain χ 1 and χ 2 Dihedral Angles. J. Chem. Theory Comput. 2012, 8 (9), 3257−3273. (34) MacKerell, a D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102 (18), 3586−3616. (35) Huang, J.; Rauscher, S.; Nawrocki, G.; Ran, T.; Feig, M.; De Groot, B. L.; Grubmüller, H.; Mackerell, A. D. CHARMM36m: An Improved Force Field for Folded and Intrinsically Disordered Proteins. Nat. Methods 2017, 14, 71−73. (36) MacKerell, A. D.; Banavali, N. K. All-Atom Empirical Force Field for Nucleic Acids: II. Application to Molecular Dynamics Simulations of DNA and RNA in Solution. J. Comput. Chem. 2000, 21 (2), 105−120.

SEP, phosphoserine; SER, serine; THR, threonine; TPO, phosphothreonine; TYR, tyrosine

■

REFERENCES

(1) Bah, A.; Forman-Kay, J. D. Modulation of Intrinsically Disordered Protein Function by Post-Translational Modifications. J. Biol. Chem. 2016, 291 (13), 6696−6705. (2) Khoury, G. A.; Baliban, R. C.; Floudas, C. A. Proteome-wide Post-Translational Modification Statistics: Frequency Analysis and Curation of the swiss-prot Database. Sci. Rep. 2011, 1, 90. (3) Manning, G. The Protein Kinase Complement of the Human Genome. Science (Washington, DC, U. S.) 2002, 298 (5600), 1912− 1934. (4) Jackson, M. D.; Denu, J. M. Molecular Reactions of Protein Phosphatases–Insights from Structure and Chemistry. Chem. Rev. 2001, 101 (8), 2313−2340. (5) Iakoucheva, L. M.; Radivojac, P.; Brown, C. J.; O’Connor, T. R.; Sikes, J. G.; Obradovic, Z.; Dunker, A. K. The Importance of Intrinsic Disorder for Protein Phosphorylation. Nucleic Acids Res. 2004, 32 (3), 1037−1049. (6) Gsponer, J.; Futschik, M. E.; Teichmann, S. A.; Babu, M. M. Tight Regulation of Unstructured Proteins: From Transcript Synthesis to Protein Degradation. Science (Washington, DC, U. S.) 2008, 322 (5906), 1365−1368. (7) Cohen, P. The Regulation of Protein Function by Multisite Phosphorylation–a 25 Year Update. Trends Biochem. Sci. 2000, 25 (12), 596−601. (8) García-Alai, M. M.; Gallo, M.; Salame, M.; Wetzler, D. E.; McBride, A. A.; Paci, M.; Cicero, D. O.; De Prat-Gay, G. Molecular Basis for Phosphorylation-Dependent, PEST-Mediated Protein Turnover. Structure 2006, 14 (2), 309−319. (9) Nardozzi, J. D.; Lott, K.; Cingolani, G. Phosphorylation Meets Nuclear Import: A Review. Cell Commun. Signaling 2010, 8 (1), 32. (10) Horn, P. J.; Carruthers, L. M.; Logie, C.; Hill, D. A.; Solomon, M. J.; Wade, P. A.; Imbalzano, A. N.; Hansen, J. C.; Peterson, C. L. Phosphorylation of Linker Histones Regulates ATP-Dependent Chromatin Remodeling Enzymes. Nat. Struct. Biol. 2002, 9 (4), 263−267. (11) Nishi, H.; Hashimoto, K.; Panchenko, A. R. Phosphorylation in Protein-Protein Binding: Effect on Stability and Function. Structure 2011, 19 (12), 1807−1815. (12) Pawson, T.; Nash, P. Assembly of Cell Regulatory Systems through Protein Interaction Domains. Science (Washington, DC, U. S.) 2003, 300 (5618), 445−452. (13) Mandell, D. J.; Chorny, I.; Groban, E. S.; Wong, S. E.; Levine, E.; Rapp, C. S.; Jacobson, M. P. Strengths of Hydrogen Bonds Involving Phosphorylated Amino Acid Side Chains. J. Am. Chem. Soc. 2007, 129 (4), 820−827. (14) Errington, N.; Doig, A. J. A Phosphoserine−Lysine Salt Bridge within an α-Helical Peptide, the Strongest α-Helix Side-Chain Interaction Measured to Date †. Biochemistry 2005, 44 (20), 7553− 7558. (15) Andrew, C. D.; Warwicker, J.; Jones, G. R.; Doig, A. J. Effect of Phosphorylation on α-Helix Stability as a Function of Position. Biochemistry 2002, 41 (6), 1897−1905. (16) Pufall, M. A.; Lee, G. M.; Nelson, M. L.; Kang, H.-S.; Velyvis, A.; Kay, L. E.; McIntosh, L. P.; Graves, B. J. Variable Control of Ets-1 DNA Binding by Multiple Phosphates in an Unstructured Region. Science 2005, 309 (5731), 142−145. (17) Xiang, S.; Gapsys, V.; Kim, H. Y.; Bessonov, S.; Hsiao, H. H.; Möhlmann, S.; Klaukien, V.; Ficner, R.; Becker, S.; Urlaub, H.; et al. Phosphorylation Drives a Dynamic Switch in Serine/Arginine-Rich Proteins. Structure 2013, 21 (12), 2162−2174. (18) Bah, A.; Vernon, R. M.; Siddiqui, Z.; Krzeminski, M.; Muhandiram, R.; Zhao, C.; Sonenberg, N.; Kay, L. E.; Forman-Kay, J. D. Folding of an Intrinsically Disordered Protein by Phosphorylation as a Regulatory Switch. Nature 2015, 519 (7541), 106−109. 677


Article

Journal of Chemical Theory and Computation (37) Feng, M. H.; Philippopoulos, M.; MacKerell, A. D.; Lim, C. Structural Characterization of the Phosphotyrosine Binding Region of a High-Affinity SH2 Domain-Phosphopeptide Complex by Molecular Dynamics Simulation and Chemical Shift Calculations. J. Am. Chem. Soc. 1996, 118 (45), 11265−11277. (38) Hansson, T.; Nordlund, P.; Aqvist, J. Energetics of Nucleophile Activation in a Protein Tyrosine Phosphatase. J. Mol. Biol. 1997, 265 (2), 118−127. (39) Margreitter, C.; Petrov, D.; Zagrovic, B. Vienna-PTM Web Server: A Toolkit for MD Simulations of Protein Post-Translational Modifications. Nucleic Acids Res. 2013, 41 (W1), W422−W426. (40) Schuler, L. D.; Daura, X.; van Gunsteren, W. F. An Improved GROMOS96 Force Field for Aliphatic Hydrocarbons in the Condensed Phase. J. Comput. Chem. 2001, 22 (11), 1205−1218. (41) Schmid, N.; Eichenberger, A. P.; Choutko, A.; Riniker, S.; Winger, M.; Mark, A. E.; Van Gunsteren, W. F. Definition and Testing of the GROMOS Force-Field Versions 54A7 and 54B7. Eur. Biophys. J. 2011, 40 (7), 843−856. (42) Reif, M. M.; Hünenberger, P. H.; Oostenbrink, C. New Interaction Parameters for Charged Amino Acid Side Chains in the GROMOS Force Field. J. Chem. Theory Comput. 2012, 8 (10), 3705− 3723. (43) Margreitter, C.; Reif, M. M.; Oostenbrink, C. Update on Phosphate and Charged Post-Translationally Modified Amino Acid Parameters in the GROMOS Force Field. J. Comput. Chem. 2017, 38 (10), 714−720. (44) Rani, L.; Mallajosyula, S. S. Phosphorylation versus OGlcNAcylation: Computational Insights into the Differential Influences of the Two Competitive Post-Translational Modifications. J. Phys. Chem. B 2017, 121 (47), 10618−10638. (45) Kim, S. Y.; Jung, Y.; Hwang, G. S.; Han, H.; Cho, M. Phosphorylation Alters Backbone Conformational Preferences of Serine and Threonine Peptides. Proteins: Struct., Funct., Genet. 2011, 79 (11), 3155−3165. (46) Vymětal, J.; Vondrásě k, J. Critical Assessment of Current Force Fields. Short Peptide Test Case. J. Chem. Theory Comput. 2013, 9 (1), 441−451. (47) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins: Struct., Funct., Genet. 2006, 65 (3), 712−725. (48) Shirts, M. R.; Klein, C.; Swails, J. M.; Yin, J.; Gilson, M. K.; Mobley, D. L.; Case, D. A.; Zhong, E. D. Lessons Learned from Comparing Molecular Dynamics Engines on the SAMPL5 Dataset. J. Comput.-Aided Mol. Des. 2017, 31 (1), 147−161. (49) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1−2, 19−25. (50) Páll, S.; Abraham, M. J.; Kutzner, C.; Hess, B.; Lindahl, E. Tackling Exascale Software Challenges in Molecular Dynamics Simulations with GROMACS. Lect. Notes Comput. Sci. (including Subser. Lect. Notes Artif. Intell. Lect. Notes Bioinformatics) 2015, 8759, 3−27. (51) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New Feathers for an Old Bird. Comput. Phys. Commun. 2014, 185 (2), 604−613. (52) Piana, S.; Laio, A. A Bias-Exchange Approach to Protein Folding. J. Phys. Chem. B 2007, 111 (17), 4553−4559. (53) Jurrus, E.; Engel, D.; Star, K.; Monson, K.; Brandi, J.; Felberg, L. E.; Brookes, D. H.; Wilson, L.; Chen, J.; Liles, K.; et al. Improvements to the APBS Biomolecular Solvation Software Suite. Protein Sci. 2018, 27 (1), 112−128. (54) Bonomi, M.; Barducci, A.; Parrinello, M. Reconstructing the Equilibrium Boltzmann Distribution from Well-Tempered Metadynamics. J. Comput. Chem. 2009, 30 (11), 1615−1621. (55) Vögeli, B.; Ying, J.; Grishaev, A.; Bax, A. Limits on Variations in Protein Backbone Dynamics from Precise Measurements of Scalar Couplings. J. Am. Chem. Soc. 2007, 129 (30), 9377−9385.

(56) Baker, E. N.; Hubbard, R. E. Hydrogen Bonding in Globular Proteins. Prog. Biophys. Mol. Biol. 1984, 44 (2), 97−179. (57) Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E. The Protein Data Bank. Nucleic Acids Res. 2000, 28 (1), 235−242. (58) Efron, B.; Tibshirani, R. J. An Introduction to the Bootstrap; Chapman & Hall/CRC Monographs on Statistics & Applied Probability, Vol. 57; Taylor & Francis: Boca Raton, FL, USA, 1994. (59) Hunter, J. D. Matplotlib: A 2D Graphics Environment. Comput. Sci. Eng. 2007, 9 (3), 90−95. (60) PyMOL Molecular Graphics System, Version ∼ 1.8; Schrödinger: Cambridge, MA, USA, 2015. (61) Lee, K. K.; Kim, E.; Joo, C.; Song, J.; Han, H.; Cho, M. SiteSelective Intramolecular Hydrogen-Bonding Interactions in Phosphorylated Serine and Threonine Dipeptides. J. Phys. Chem. B 2008, 112 (51), 16782−16787. (62) Avbelj, F.; Grdadolnik, S. G.; Grdadolnik, J.; Baldwin, R. L. Intrinsic Backbone Preferences Are Fully Present in Blocked Amino Acids. Proc. Natl. Acad. Sci. U. S. A. 2006, 103 (5), 1272−1277. (63) Piana, S.; Donchev, A. G.; Robustelli, P.; Shaw, D. E. Water Dispersion Interactions Strongly Influence Simulated Structural Properties of Disordered Protein States. J. Phys. Chem. B 2015, 119 (16), 5113−5123. (64) Piana, S.; Lindorff-Larsen, K.; Shaw, D. E. How Robust Are Protein Folding Simulations with Respect to Force Field Parameterization? Biophys. J. 2011, 100 (9), L47−L49. (65) Best, R. B.; Hummer, G. Optimized Molecular Dynamics Force Fields Applied to the Helix-Coil Transition of Polypeptides. J. Phys. Chem. B 2009, 113 (26), 9004−9015. (66) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. A WellBehaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The RESP Model. J. Phys. Chem. 1993, 97 (40), 10269−10280. (67) Wong, S. E.; Bernacki, K.; Jacobson, M. Competition between Intramolecular Hydrogen Bonds and Solvation in Phosphorylated Peptides: Simulations with Explicit and Implicit Solvent. J. Phys. Chem. B 2005, 109 (11), 5249−5258. (68) Hendus-Altenburger, R.; Lambrughi, M.; Terkelsen, T.; Pedersen, S. F.; Papaleo, E.; Lindorff-Larsen, K.; Kragelund, B. B. A Phosphorylation-Motif for Tuneable Helix Stabilisation in Intrinsically Disordered Proteins − Lessons from the Sodium Proton Exchanger 1 (NHE1). Cell. Signalling 2017, 37 (May), 40−51. (69) Srinivasaraghavan, K.; Nacro, K.; Grüber, G.; Verma, C. S. Effect of Ser392 Phosphorylation on the Structure and Dynamics of the Polybasic Domain of ADP Ribosylation Factor Nucleotide Site Opener Protein: A Molecular Simulation Study. Biochemistry 2013, 52 (41), 7339−7349. (70) McAnany, C. E.; Mura, C. Claws, Disorder, and Conformational Dynamics of the C-Terminal Region of Human Desmoplakin. J. Phys. Chem. B 2016, 120 (33), 8654−8667. (71) Langenfeld, F.; Guarracino, Y.; Arock, M.; Trouvé, A.; Tchertanov, L. How Intrinsic Molecular Dynamics Control Intramolecular Communication in Signal Transducers and Activators of Transcription Factor STAT5. PLoS One 2015, 10 (12), e0145142. (72) Ithuralde, R. E.; Turjanski, A. G. Phosphorylation Regulates the Bound Structure of an Intrinsically Disordered Protein: The P53TAZ2 Case. PLoS One 2016, 11 (1), e0144284. (73) Ruff, E. F.; Muretta, J. M.; Thompson, A. R.; Lake, E. W.; Cyphers, S.; Albanese, S. K.; Hanson, S. M.; Behr, J. M.; Thomas, D. D.; Chodera, J. D.; Levinson, N. M.; et al. A Dynamic Mechanism for Allosteric Activation of Aurora Kinase A by Activation Loop Phosphorylation. eLife 2018, 7, 1−22. (74) Löw, C.; Homeyer, N.; Weininger, U.; Sticht, H.; Balbach, J. Conformational Switch upon Phosphorylation: Human CDK Inhibitor P19 INK4d between the Native and Partially Folded State. ACS Chem. Biol. 2009, 4 (1), 53−63. (75) Lyons, A. J.; Gandhi, N. S.; Mancera, R. L. Molecular Dynamics Simulation of the Phosphorylation-Induced Conformational Changes 678


Article

Journal of Chemical Theory and Computation of a Tau Peptide Fragment. Proteins: Struct., Funct., Genet. 2014, 82 (9), 1907−1923. (76) Prokopovich, D. V.; Whittaker, J. W.; Muthee, M. M.; Ahmed, A.; Larini, L. Impact of Phosphorylation and Pseudophosphorylation on the Early Stages of Aggregation of the Microtubule-Associated Protein Tau. J. Phys. Chem. B 2017, 121 (9), 2095−2103. (77) Fitzkee, N. C.; Fleming, P. J.; Rose, G. D. The Protein Coil Library: A Structural Database of Nonhelix, Nonstrand Fragments Derived from the PDB. Proteins: Struct., Funct., Genet. 2005, 58 (4), 852−854. (78) Swindells, M. B.; Macarthur, M. W.; Thornton, J. M. Intrinsic φ,ψ Propensities of Amino Acids, Derived from the Coil Regions of Known Structures. Nat. Struct. Biol. 1995, 2 (7), 596−603. (79) Jiang, F.; Han, W.; Wu, Y. D. The Intrinsic Conformational Features of Amino Acids from a Protein Coil Library and Their Applications in Force Field Development. Phys. Chem. Chem. Phys. 2013, 15 (10), 3413−3428. (80) MacKerell, A. D.; Feig, M.; Brooks, C. L. Improved Treatment of the Protein Backbone in Empirical Force Fields. J. Am. Chem. Soc. 2004, 126 (3), 698−699. (81) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved Side-Chain Torsion Potentials for the Amber Ff99SB Protein Force Field. Proteins: Struct., Funct., Genet. 2010, 78 (8), 1950−1958. (82) Debiec, K. T.; Cerutti, D. S.; Baker, L. R.; Gronenborn, A. M.; Case, D. A.; Chong, L. T. Further along the Road Less Traveled: AMBER Ff15ipq, an Original Protein Force Field Built on a SelfConsistent Physical Model. J. Chem. Theory Comput. 2016, 12 (8), 3926−3947. (83) Wang, L. P.; McKiernan, K. A.; Gomes, J.; Beauchamp, K. A.; Head-Gordon, T.; Rice, J. E.; Swope, W. C.; Martínez, T. J.; Pande, V. S. Building a More Predictive Protein Force Field: A Systematic and Reproducible Route to AMBER-FB15. J. Phys. Chem. B 2017, 121 (16), 4023−4039.

679


AMBER and CHARMM Force Fields Inconsistently ... - ACS Publications

Recommend Documents