Are Current Atomistic Force Fields Accurate Enough to Study Proteins in Crowded Environments?
[['Different macromolecules occupy up to 40% of the total cytoplasmic volume in typical cells , with proteins being the most abundant class of molecules [1] , [2] .', 'Importantly , such densely packed environments strongly affect various thermodynamic and kinetic properties of proteins including their structural stabilities , intermolecular binding affinities and enzymatic rates [3] , [4] .', 'Most drastically , as a consequence of high intracellular concentration , proteins can form cytotoxic aggregates that have been linked with numerous pathologies [5] , [6] .', 'Additionally , most solution-based biophysical experimental methods , such as NMR or different spectroscopies , involve samples with relatively high protein concentration .', 'For these reasons , increasing attention has recently been devoted to studying protein behavior in crowded environments .', 'For example , volume-exclusion and confinement effects in the context of crowding have been qualitatively well-understood by statistical mechanical theories and computer simulations [4] , [7] .', 'Moreover , Brownian dynamics and coarse-grained simulations have been used to provide a detailed description of multipart mixtures of biomolecules and have sometimes even matched real systems in their complexity [8]–[10] .', 'Although highly successful , however , such approaches still fall short of capturing the fully atomistic , dynamic picture of high-concentration macromolecular systems .', 'To this end , molecular dynamics ( MD ) simulations , a high-resolution computational biology tool [11] , have recently been employed to model different aspects of crowding when it comes to protein structure , dynamics and interactions as well as solvent behavior [12]–[16] .', 'Moreover , MD simulations have been applied to explore early events in the formation of protein aggregates , focusing predominantly on short peptides with high aggregation propensity [17]–[22] .', 'On the other hand , arguably due to high computational expenses , control tests to show that non-aggregating polypeptides do not aggregate are rarely performed , with only a few attempts in this direction .', 'For example , Gsponer et al . [17] and Tsai et al . [20] have shown that control mutant peptides aggregate into structures with reduced amyloid character as compared to aggregation-prone peptides .', 'However , to the best of our knowledge , not a single study so far has provided clear evidence of a known non-aggregating polypeptide as a negative control .', 'Here , we use classical MD simulations to study the behavior of the 36-residue villin headpiece mini-protein [23]–[25] at atomistic resolution with multiple copies of the protein in the simulation box .', 'In contrast to other MD simulation studies of protein-protein interactions or aggregation , our choice of the model system is primarily motivated by the fact that villin headpiece does not self-associate or aggregate at moderate protein concentrations .', 'For example , this well-studied actin-binding polypeptide remains fully soluble at protein concentration of 1–2 mM as shown by infrared spectroscopy [26] , [27] .', 'Additionally , Fourier transform infrared spectra of the peptide indicate no aggregation at 6 mM when suspended in a 50 mM sodium acetate buffer with a possibility of aggregate or dimer formation only at a higher protein concentration of ∼18 mM [25] .', 'Similarly , circular dichroism experiments suggest that there is no significant aggregation of villin headpiece at the concentration of 9 . 2 mM , although here the data was collected at −40°C in glycerol/water solution [24] .', 'Finally , NMR spectra of villin headpiece can successfully be recorded even at the concentration of 32 mM , although with changes in chemical shifts of surface residues as a consequence of increased protein-protein interactions upon crowding [28] .', 'As a control for our simulations of wild-type , folded villin headpiece , we also simulate the fully carbonylated , unfolded and aggregation-prone villin headpiece at high concentration .', 'Increased levels of protein aggregation have been repeatedly related to protein oxidative damage with highly oxidized proteins being a frequent component of potentially cytotoxic aggregates [29]–[31] .', 'In a recent study , we have shown that metal-catalyzed carbonylation , arguably one of the most important types of irreversible oxidation , drastically increases the intrinsic aggregability of villin headpiece by directly affecting its hydrophobicity , net charge and secondary structure [32] , protein properties shown to strongly influence aggregation propensity [33] .', 'Here , we study wild-type and carbonylated villin headpiece at experimentally relevant protein concentrations ( 6 mM and 9 . 2 mM ) and different sodium chloride concentrations ( 0 to 0 . 8 M ) using two different force fields ( GROMOS 45a3 [34] and GROMOS 54a7 [35] ) , SPC water model [36] and two types of electrostatics treatment ( PME - particle mesh Ewald [37] and RF - reaction-field [38] ) .', 'Moreover , in order to account for potential inaccuracies in the physical description of the system , we use four additional force fields ( AMBER94 [39] , AMBER99SB-ILDN [40] , CHARMM22-CMAP [41] and OPLS-AAL [42] ) to reanalyze the wild-type villin headpiece GROMOS trajectories and employ two of these ( AMBER99SB-ILDN and CHARMM22-CMAP ) in combination with the TIP3P water model [43] to perform additional MD simulations .', 'The principal question that we ask is how does the behavior of villin headpiece differ at high protein concentration as compared to infinite dilution and , in the process , we explore the limitations of current atomistic force fields in describing biologically and experimentally realistic protein solutions .'], ['Villin headpiece is a 3-helix bundle protein with a tightly packed hydrophobic core comprised of three phenylalanine residues ( Fig . 1A ) [23] .', 'In order to structurally characterize the simulated proteins , we use:', '1 ) atom-positional backbone root-mean-square deviation ( RMSD ) from the native experimental villin headpiece structure [23] ,', '2 ) the number of residues in α-helical conformation ( #α ) , and', '3 ) the sum of distances between the centers of mass of core phenylalanines ( ΣdPHEs ) .', 'In Figure 1B , we present aggregate averages of these three measures over all simulations performed using GROMOS 45a3 and 54a7 force fields for wild-type and carbonylated villin headpiece , with the average values for all individual simulation protocols given in Table S1 .', 'If not stated otherwise , these GROMOS simulations are used for analysis throughout this study .', 'Note also that the initial configurations for simulations with multiple copies of the protein in the simulation box were chosen randomly from the equilibrium ensembles of the last 25 ns of the previously simulated 110-ns-long , infinite dilution trajectories of wild-type and carbonylated villin headpiece [32] .', 'These same ensembles were also used for analysis of protein behavior at infinite dilution , thus ensuring that the two sets of simulations are fully comparable .', 'Overall , wild-type polypeptides remain folded at infinite dilution as previously shown [32] with <RMSD>\u200a=\u200a2 . 7±1 . 2 Å , <#α>\u200a=\u200a22 . 3±2 . 3 and <ΣdPHEs>\u200a=\u200a20 . 3±2 . 4 Å ( Fig . 1B , wild-typesingle ) .', 'On the other hand , carbonylated molecules at infinite dilution populate the unfolded state [32] , with an increase in the average RMSD ( <RMSD>\u200a=\u200a8 . 1±0 . 8 Å ) and ΣdPHEs ( <ΣdPHEs>\u200a=\u200a33 . 3±5 . 8 Å ) and a significant decrease in the α-helical content ( <#α>\u200a=\u200a7 . 0±3 . 9 ) ( Fig . 1B , carbonylatedsingle ) .', 'Importantly , both wild-type and carbonylated villin headpiece , simulated with multiple copies in the simulation box , closely resemble the structural properties of the equivalent polypeptides at infinite dilution for both GROMOS force fields used ( Fig . 1B , wild-typemultiple and carbonylatedmultiple ) .', 'The only appreciable difference concerns the α-helical content in carbonylation simulations where the average value at infinite dilution exceeds that for multiple-copy simulations by approximately 40% ( 12 . 1±5 . 1 vs . 7 . 0±3 . 9 ) .', 'Additionally , we observe major unfolding ( RMSD>5 Å ) for 4 out of 116 simulated wild-type villin headpiece molecules in multiple-copy simulations , which is consistent with the fact that villin headpiece is a marginally stable protein .', 'However , it is noteworthy that two out of these four instances of unfolding occur in simulations obtained using GROMOS54a7 force field in combination with the RF method .', 'Finally , the overall structural stability of wild-type villin headpiece in multiple-copy GROMOS simulations is very similar to that observed in AMBER99SB-ILDN and CHARMM22-CMAP simulations ( Table S1 ) .', 'Altogether , 10 different simulation setups have been examined ( Table S1 ) using GROMOS 45a3 and 54a7 force fields in combination with different types of electrostatics treatment at various ionic strengths .', 'As predicted , carbonylated villin headpiece molecules start to associate after only a few nanoseconds of free diffusion , leading to an exponential decrease in the number of free monomers in solution and a concomitant increase in the number of intermolecular atomic contacts at all conditions examined ( Fig . 2 , see Methods for more details ) .', 'Configurations in which the number of free carbonylated monomers reaches 0 are observed in 23 out of 26 simulated trajectories with aggregates comprised of at least 75% of all monomers observed in every single trajectory ( Fig . 2 ) .', 'Surprisingly , in contradiction with experimental findings [25] , [26] , the fraction of free monomers in wild-type villin headpiece simulations decreases over time in the same manner as for the carbonylated villin headpiece and the same is true for the increase in the number of inter-protein atomic contacts for both GROMOS force fields used ( Fig . 2 ) .', 'In fact , complete disappearance of free monomers and , at the same time , appearance of aggregates comprised of all simulated monomers , is observed in 22 out of 26 GROMOS trajectories of the wild-type villin headpiece , with aggregation being largely independent of the specific force field or simulation setup used .', 'To further characterize the aggregation process of the studied systems , we have analyzed the kinetics of protein-protein association and dissociation in simulations obtained using the two GROMOS force fields .', 'The average waiting time for dissociation of individual wild-type molecules from aggregates over all simulated conditions is 45 ns , which is significantly longer than the average waiting time required for their association of 17 ns .', 'Note that the observed association time is in an excellent agreement with the inverse of the diffusion-limited association rate ( ron−1\u200a= ( kon[villin] ) −1\u200a=\u200a16 ns ) as estimated following the formalism of Smoluchowski [44] from the concentration ( [villin]\u200a=\u200a6 mM ) , the size and the diffusion coefficient of wild-type villin headpiece monomers directly calculated from simulated trajectories .', 'In this framework , the rate constant is simply given by kon\u200a=\u200a4πRD , where R is the sum of radii of interacting molecules ( here taken as twice the radius of gyration of villin headpiece ) and D is the relative diffusion coefficient .', 'Similarly , the average waiting times of protein dissociation and association for the carbonylated system are 51 ns and 19 ns , respectively .', 'Importantly , approximately 65% of wild-type and 70% of carbonylated complexes that have formed never dissociate during our simulations , including all those with life-times longer than 55 ns ( Fig . S1 ) .', 'This fact clearly indicates that the actual life-times of villin headpiece in the aggregated state may be much longer than the average values obtained from simulations , which are limited by sampling .', 'On the other hand , longer simulations would not significantly influence the estimated association times , since only less than 5% of proteins in both systems remain free in solution throughout simulations .', 'This in turn suggests that even much lower protein concentration would likely not prevent wild-type villin headpiece from aggregating , but only increase the search time needed for proteins to find each other by free diffusion .', 'For both wild-type and carbonylated villin headpiece simulations performed using GROMOS 45a3 and 54a7 force fields , starting from free monomers , peptide dimers begin to form first , followed by the formation of trimers and tetramers , leading to the maximum number of tetrameric aggregates after approximately 50 ns of simulated time ( Fig . 3A ) .', 'Overall , tetrameric aggregates are the most abundant and free monomers the least abundant species for both wild-type and carbonylated villin headpiece simulations , representing on average approximately 45% and 15% of the total protein content at 50 ns in both cases , respectively ( Fig . 3B ) .', 'In order to identify specific residues involved in self-association , we have calculated the fraction of intermolecular atomic contacts formed by each residue ( side chains only ) over all wild-type and carbonylated GROMOS simulations pooled together separately , including a normalization by solvent-accessibility to account for surface exposure , i . e . , probability to interact .', 'Surprisingly , while association of hydrophobic residues is believed to be a key element in protein aggregation , our analysis reveals that interactions in both wild-type and carbonylated aggregates are dominated by either backbone or glutamine and asparagine residues ( Fig . 4 ) .', 'In addition to these largely hydrophilic moieties , aggregates are also characterized by contacts involving hydrophobic ring-containing phenylalanines and tryptophans ( Fig . 4 ) .', 'Finally , we have also analyzed the average contact maps for both wild-type and carbonylated systems ( Fig . S2 ) , which demonstrate that the N-terminal residues form fewer intermolecular contacts than the rest of the molecule in both studied systems .', 'Similarly , the region between residues 15 and 25 of the wild-type contact map is depleted in contacts for the two GROMOS force fields , whereas the carbonylated map displays a somewhat more even distribution of contacts .', 'We should emphasize , however , that no pronounced system-specific pattern among the residues exhibiting a high number of intermolecular contacts has been observed .', 'To address this , we have explored the role of enthalpic contributions in the aggregation of the studied systems , including solvent-solvent , protein-solvent and protein-protein interactions , by calculating the average difference in potential energy between the fully aggregated ( tetrameric ) and the non-aggregated ( monomeric ) conformers ( Fig . 5 ) .', 'Moreover , we have monitored the potential energy as a function of simulated time ( Fig . S3 ) .', 'Expectedly , in GROMOS simulations solvent-solvent and protein-protein interactions provide favorable contribution , while protein-solvent interactions provide unfavorable contribution to the total potential energy of aggregation , as seen in the case of simulations with the most extensive sampling ( 200 ns ) and the PME electrostatics treatment ( Fig . 5A and S3 ) .', 'However , the total potential energy of the same systems decreases with simulated time ( Fig . S3 ) and is significantly lower in the aggregated state ( by approximately 150 kJ/mol and 330 kJ/mol for the wild-type and the carbonylated system simulated by GROMOS45a3 force field , respectively , and approximately 100 kJ/mol and 270 kJ/mol for the wild-type and the carbonylated system simulated by GROMOS54a7 force field , respectively , Fig . 5A ) , suggesting that self-association may be an enthalpically driven process .', 'Further analysis of the simulations obtained using RF electrostatics treatment ( the GROMOS 45a3 parameter set only ) shows that the aggregated species are also favored at every salt concentration examined ( Fig . 5B ) .', 'In order to study other force fields in addition to GROMOS 45a3 and 54a7 , we have re-evaluated the potential energy of the simulated GROMOS trajectories by employing 4 additional widely used MD force fields ( AMBER94 , AMBER99SB-ILDN , CHARMM22-CMAP and OPLS-AAL ) and two types of electrostatics treatments ( RF and PME ) on energy-minimized snapshot configurations , where the force field used for energy minimization was also used to re-evaluate the potential energy .', 'For consistency , the same procedure was also repeated for the GROMOS 45a3 and 54a7 force fields as well .', 'Note that we have only re-evaluated potential energies for wild-type systems since parameters for carbonylated residues are only available for the GROMOS force fields .', 'Moreover , for all force fields other than the two GROMOS force fields , we have used the TIP3P water model [43] in the above procedure .', 'Remarkably , regardless of the type of electrostatics treatment , all of the evaluated force fields show the same trend in favoring the tetrameric aggregated over the monomeric non-aggregated state of wild-type villin headpiece ( Fig . 6A and S4 ) with the extent of such bias ranging from approximately −50 kJ/mol to −240 kJ/mol ( PME calculations ) for GROMOS54a7 and OPLS-AAL force fields , respectively .', 'Importantly , in order to decrease the uncertainty in estimation , we have averaged the re-evaluated differences in potential energy between the tetrameric and monomeric configurations for each force field over all trajectories no matter which simulation condition they came from ( Fig . 6A and S4A ) .', 'Under the assumption that the entropic contribution to the free energy of association is comparable for all different force fields , this result suggests that utilizing any of these force fields would most likely lead to aggregation of wild-type villin headpiece in MD simulations , and therefore be at odds with experimental observations .', 'Moreover , while different force fields exhibit different dependence on salt when it comes to aggregation , they all favor the aggregated state at most conditions examined ( Fig . 6B ) .', 'Interestingly , the difference in potential energy of tetrameric and monomeric configurations reaches a maximum at either 0 . 2 M NaCl ( CHARMM22-CMAP ) or 0 . 1 M NaCl ( all other force fields ) ( Fig . 6B ) .', 'The analysis of this effect , which is possibly related to the different chaotropic behavior of NaCl at different concentrations , is nevertheless beyond the scope of our present study .', 'Finally , we should emphasize that our energy re-evaluation procedure is only approximate: the calculated differences in the total potential energy come with error bars of approximately 200 kJ/mol for each simulation and 50 kJ/mol for the average difference ( Figs . 6 and S4 ) , as estimated using the root-mean-square errors from the potential energies derived directly from simulations , calculated for the same GROMOS force fields as used in the simulations .', 'In order to further explore the above findings , we have additionally simulated the wild-type villin headpiece using AMBER99SB-ILDN [40] and CHARMM22-CMAP [41] force fields , focusing on simulation conditions ( [protein]\u200a=\u200a6 mM and [NaCl]\u200a=\u200a0 . 05 M ) , which are most similar to the experiment described in reference [25] .', 'Note that carbonylated versions of the peptide were not modeled with these force fields since parameters for carbonylated amino acids were not currently available .', 'Similarly to trajectories obtained using GROMOS force fields , for both AMBER99SB-ILDN and CHARMM22-CMAP force fields , the number of intermolecular atomic contacts increases and the number of free monomers decreases with concomitant decrease in the potential energy of the system over time , as shown by the averages over three 200 ns-long independent simulations ( Fig . 7A and Fig . S5 ) .', 'Indeed , for the AMBER99SB-ILDN force field , the long-time fraction of free monomers drops down to 0 , while for CHARMM22-CMAP it stabilizes at a low 8 . 3% .', 'Importantly , the difference in total potential energy between the aggregated tetrameric and the non-aggregated monomeric state for both force fields is approximately −100 kJ/mol ( Fig . 7B ) , which is similar to the value observed in simulations generated using GROMOS54a7 force field and PME electrostatics treatment ( Fig . 5A ) .', 'On the other hand , these numbers are by almost 100 kJ/mol less negative than the values obtained by re-evaluating the energies of the snapshots obtained in GROMOS simulations under similar conditions ( Fig . 6A and 7B ) .', 'This further suggests that , although qualitatively likely correct , our energy re-evaluation procedure , especially in the case of individual salt conditions ( Fig . 6B ) , is quantitatively only approximate as already discussed above .'], ['By probing various simulation conditions , we have observed that simulated wild-type and carbonylated villin headpiece molecules aggregate in a similar fashion , despite the fact that they display markedly different aggregation propensities [32] .', 'This finding is even more remarkable given the fact that we have focused on conditions experimentally known to render the wild-type villin headpiece soluble .', 'Villin headpiece self-association in our simulations is accompanied by a decrease in protein-protein and solvent-solvent , and an increase in protein-solvent short-range potential energy , resulting in a significant net decrease of total potential energy upon aggregation .', 'Strikingly , by re-evaluating the potential energy of the aggregated and non-aggregated conformers from GROMOS simulations using 6 widely-used atomistic force fields , we have shown that all of them favor the aggregated state ( Fig . 6 and S4 ) .', 'Here , one should emphasize that this result was obtained by averaging over conformers originally obtained by differently defined Hamiltonians and using a limited number of energy minimized configurations , significantly increasing the error in the estimated potential energy ( see Methods for more details ) .', 'These shortcomings notwithstanding , the general trends for all force fields consistently point in the direction of favoring aggregation .', 'Moreover , this result was further confirmed by direct simulations of villin headpiece at high concentration using AMBER99SB-ILDN and CHARMM22-CMAP force fields ( Fig . 7 and S5 ) , in which spurious aggregation was also observed , albeit for a single simulation setup only .', 'While more extensive simulation efforts are undeniably needed to complete the analysis of different atomistic force fields , we believe that our present results do reveal significant deficiencies in how they capture the behavior of proteins in crowded solutions and invite caution when it comes to using them for such purposes .', 'Even though it is experimentally known that villin headpiece can begin to aggregate at protein concentration of 6 mM after a period of a few days [25] , it is highly unlikely that this is related to the nanosecond-time-scale aggregation observed here .', 'In order to address the possibility that aggregation of the wild-type system was still induced by the “borderline” villin headpiece concentrations used , we have analyzed the kinetics of protein-protein association and dissociation .', 'This analysis has revealed that a large number of protein-protein complexes never dissociate in the course of simulated trajectories , and that the average protein dissociation time is markedly longer than the average association time ( Fig . S1 ) , suggesting that aggregation in MD simulations would occur even at lower concentrations than applied here .', 'Having said this , it is important to emphasize that in our analysis we do not directly compare simulated and experimental observables reporting on villin headpiece aggregation , but rather rely on indirect interpretation of different experimental observables .', 'It is well known that such comparisons at the level of interpretation can be fraught with difficulties and are never fully unambiguous [45] , [46] .', 'For example , it is possible that the CD spectra of tetrameric aggregates are similar to those of the monomeric form and that the experimental interpretation of such spectra , which was used to suggest that villin headpiece does not aggregate at the concentration of 9 . 2 mM , was simply inaccurate [24] .', 'On the other hand , both 2D IR and NMR are extremely sensitive to local changes in the environment of individual polypeptide groups and have successfully been used to detect dimerization , multimerization and aggregation of different proteins [47]–[52] .', 'This further suggests that villin headpiece is indeed monomeric at concentrations studied herein , as determined in the experiments performed by using these methods [25]–[28] .', 'In conjunction with our analysis of residence times in different states , this in turn bolsters the claim that the aggregation observed in our simulations is indeed unphysical , but the above caveats should still be borne in mind .', 'Our analysis shows that interactions in villin headpiece aggregates are mainly dominated by polar glutamines and asparagines as well as hydrophobic ring-containing phenylalanines and tryptophans ( Fig . 4 ) .', 'Although the two former hydrophilic amide-containing amino acids are generally considered to be aggregation reducing [33] , [53] , they have been repeatedly linked to protein aggregation and deposition disorders , most notably in poly-Q diseases [5] , [20] , [54] .', 'In addition to this , the high occurrence of peptide backbone atoms among villin headpiece intermolecular contacts ( Fig . 4 ) supports the hypothesis that poly-peptide chains have a general tendency to aggregate due to intrinsic aggregability of the protein backbone [55] .', 'This is further corroborated by recent findings that poly-glycine and poly-alanine chains aggregate readily [56] .', 'Finally , Andrews and Elcock have recently used MD simulations in combination with several atomistic force fields and water models to analyze the behavior of high-concentration aqueous solutions of glycine , valine , phenylalanine and asparagine [16] .', 'In agreement with the available experimental data , their analysis shows that asparagine and phenylalanine exhibit a considerably greater density increment , a measure of interaction propensity , when compared to glycine and valine , which is also well matched by our present data .', 'Taken together , our results suggest that typical classical MD force fields possibly bias protein aggregation by overestimating protein-protein and solvent-solvent as opposed to protein-solvent interactions .', 'Parenthetically , the primary strategy in parameterizing the GROMOS 54a7 force field , which exhibits one of the lowest preference towards the aggregated state when compared to other evaluated force fields , was to reproduce experimental hydration free energies of amino-acid side-chain analogs , a property closely related to solubility in water , while the other force fields examined here significantly underestimate this property [57]–[59] .', 'Overall , the imbalance between protein-protein , protein-solvent and solvent-solvent components of total potential energy may be partly a consequence of the widely-used force field validation approaches , which frequently aim at reproducing secondary and tertiary structures of well-characterized proteins .', 'Simply put , strengthening of protein-protein and weakening of protein-solvent interactions leads to stabilization of protein structure and this can have direct repercussions on the behavior of proteins in crowded environments , as shown here .', 'This , if true , further suggests that realistic polypeptides may display more dynamics and unstructuredness than generally observed in MD-simulation studies .', 'Interestingly , a recent study exploring the limitations of MD simulations by employing CHARMM22* force field [60] and a state-of-the-art designated supercomputer to perform a 200-µs-long simulation of an intrinsically unfolded protein revealed that the modeled protein appears to be more compact and collapsed than observed by NMR [61] , further supporting this speculation .', 'Importantly , these potential flaws of current force fields may have strong implications when it comes to the accuracy of MD models in describing protein dynamics and interactions in biologically relevant crowded environment .', 'Recent validation studies of force fields show that they improve over time , but are still not able to reproduce all relevant experimental data [62] , [63] .', 'Such synergistic efforts between experiment and MD simulations should lead to improvements in computational models of biomolecular systems in the context of experimentally or biologically realistic conditions .', 'We hope that the results presented herein will provide a new source of motivation in this direction .'], ['Classical MD simulations were used to study the behavior of villin headpiece with multiple copies of the molecule in the simulation box .', 'We examined both wild-type villin headpiece ( sequence: MLSDEDFKAVFGMTRSAFANLPLWKQQNLKKEKGLF ) and its carbonylated form .', 'The seven letters in bold mark the most important carbonylable amino acids ( K , R and P ) in villin headpiece , which were all modified in the carbonylated form of the molecule .', 'Upon carbonylation , lysine is converted into aminoadipic-semialdehyde , while arginine and proline are converted into glutamic-semialdehyde , for which force field parameters were taken from refs [32] , [64] , [65] .', 'Ten simulation protocols were applied ( Table S1 ) , varying in protein ( 6 mM and 9 . 2 mM ) and salt ( 0 M , 0 . 05 M , 0 . 1 M , 0 . 2 M , 0 . 4 M and 0 . 8 M ) concentration , force field ( GROMOS 45a3 and 54a7 parameter sets ) , electrostatics treatment ( RF - reaction-field [38] and PME - particle mesh Ewald [37] ) , the number of protein molecules in a simulated box ( 4 and 8 ) , simulation time ( 50 ns , 100 ns and 200 ns ) , and the number of independent simulations ( 1 and 3 ) .', 'Additionally , the wild-type villin headpiece was simulated at 6 mM protein concentration and 0 . 05 M salt concentration using AMBER99SB-ILDN [40] and CHARMM22-CMAP [41] force fields , PME electrostatic treatment and three 200 ns-long independent replicas in each case ( Table S1 ) .', 'All MD simulations were carried out using GROMACS biomolecular simulation package [66] , keeping the system at 300 K and 1 bar using a Berendsen thermostat ( τT\u200a=\u200a0 . 05 ps ) and barostat ( τp\u200a=\u200a1 ps and compressibility\u200a=\u200a4 . 5×10−5 bar−1 ) [67] .', 'A cutoff of 1 . 4 nm for both Lennard-Jones and Coulomb potentials was used with the dielectric constant of 65 for RF , and the Fourier spacing of 0 . 1 nm for PME calculations .', 'Starting from either wild-type or carbonylated free monomers ( 4 and 8 for 6 mM and 9 . 2 mM protein concentration , respectively ) maximally separated and randomly oriented in a simulation box , polypeptides were allowed to diffuse freely and interact with each other and the solvent .', 'The initial configurations were prepared by solvating villin headpiece molecules in a cubic simulation box , with the size defined by the protein concentration , i . e . , sides of approximately 10 . 4 nm ( 6 mM ) and 11 . 3 nm ( 9 . 2 mM ) , for a total number of atoms exceeding 100 , 000 in all cases ( see Table S1 for more details ) .', 'Villin headpiece monomers were randomly selected from the ensemble of the last 25 ns of five 110-ns-long independent simulated trajectories of both the peptide in the wild-type and carbonylated form from ref [32] , only taking into account structures with atom-positional backbone root-mean-square deviation ( RMSD ) from the native villin headpiece structure [23] in the range of the average ensemble plus or minus one standard deviation .', 'After filling the simulation box with SPC water molecules [36] for GROMOS 45a3 and 54a7 simulations or TIP3P water molecules [43] for AMBER99SB-ILDN and CHARMM22-CMAP simulations , sodium chloride was added at a given concentration by replacing the equivalent number of randomly chosen water molecules , all the while ensuring charge neutrality with an excess of 2 Cl− and 4 Na+ ions per villin headpiec
['Introduction', 'Results', 'Discussion', 'Methods']
['The high concentration of macromolecules in the crowded cellular interior influences different thermodynamic and kinetic properties of proteins , including their structural stabilities , intermolecular binding affinities and enzymatic rates .', 'Moreover , various structural biology methods , such as NMR or different spectroscopies , typically involve samples with relatively high protein concentration .', 'Due to large sampling requirements , however , the accuracy of classical molecular dynamics ( MD ) simulations in capturing protein behavior at high concentration still remains largely untested .', 'Here , we use explicit-solvent MD simulations and a total of 6 . 4 µs of simulated time to study wild-type ( folded ) and oxidatively damaged ( unfolded ) forms of villin headpiece at 6 mM and 9 . 2 mM protein concentration .', 'We first perform an exhaustive set of simulations with multiple protein molecules in the simulation box using GROMOS 45a3 and 54a7 force fields together with different types of electrostatics treatment and solution ionic strengths .', 'Surprisingly , the two villin headpiece variants exhibit similar aggregation behavior , despite the fact that their estimated aggregation propensities markedly differ .', 'Importantly , regardless of the simulation protocol applied , wild-type villin headpiece consistently aggregates even under conditions at which it is experimentally known to be soluble .', 'We demonstrate that aggregation is accompanied by a large decrease in the total potential energy , with not only hydrophobic , but also polar residues and backbone contributing substantially .', 'The same effect is directly observed for two other major atomistic force fields ( AMBER99SB-ILDN and CHARMM22-CMAP ) as well as indirectly shown for additional two ( AMBER94 , OPLS-AAL ) , and is possibly due to a general overestimation of the potential energy of protein-protein interactions at the expense of water-water and water-protein interactions .', 'Overall , our results suggest that current MD force fields may distort the picture of protein behavior in biologically relevant crowded environments .']
['Protein behavior is strongly affected by highly crowded and interaction-rich environments , i . e . , typical conditions in both biologically relevant systems , such as the cellular interior , and solution-based structural experiments , including NMR and different spectroscopies .', 'On the other hand , primarily because of limited computational power , molecular dynamics ( MD ) simulations , a premier high-resolution method for analyzing structure , dynamics and interactions of proteins , have been predominantly used to study individual proteins at infinite dilution .', 'To fill this gap , we use MD simulations to study the behavior of wild-type ( aggregation-resistant ) and oxidatively damaged ( aggregation-prone ) forms of villin headpiece at high concentration , and reveal unexpected limitations and inaccuracies of modern-day MD force fields when it comes to modeling proteins at physiologically or experimentally relevant concentrations .']
['macromolecular complex analysis', 'biochemistry', 'proteins', 'protein structure', 'biology and life sciences', 'computational biology', 'biophysics', 'molecular biology', 'macromolecular structure analysis', 'biophysical simulations']