Quantum Mechanics of Proteins in Explicit Water: The Role of Plasmon-Like Solute-Solvent Interactions
QQuantum Mechanics of Proteins in Explicit Water: The Role ofPlasmon-Like Solute-Solvent Interactions
Martin St¨ohr and Alexandre Tkatchenko ∗ Physics and Materials Science Research Unit, University of Luxembourg,L-1511 Luxembourg, Luxembourg ∗ To whom correspondence should be addressed; E-mail: [email protected]
Quantum-mechanical van der Waals dispersion interactions play an essential role for both intra -protein and protein-water interactions – the two main driving forces for the structure and dynamicsof proteins in aqueous solution. Typically, these interactions are only treated phenomenologicallyvia pairwise potential terms in classical force fields. Here, we use an explicit quantum-mechanicalapproach based on density-functional tight-binding with the many-body dispersion formalism,which allows us to demonstrate the unexpected relevance of the many-body character of disper-sion interactions for protein energetics and the protein-water interaction. In contrast to commonlyemployed pairwise approaches, many-body effects significantly decrease the relative stability ofthe native state in the absence of water. In an aqueous environment, the collective character of theprotein-water van der Waals interaction counteracts this effect and stabilizes native conformationsand transition states. This stabilization arises due to a high degree of delocalization and collec-tivity of protein-water dispersion interactions, suggesting a remarkable persistence of long-rangeelectron correlation through aqueous environments. Our findings are exemplified on prototypi-cal showcases of proteins forming β -sheets, hairpins, and helices, emphasizing the crucial role ofplasmon-like solute-solvent interactions in biomolecular systems. Introduction
Water is an essential basis of life. It provides the envi-ronment in which the biomolecular machinery can existand function. By screening and stabilizing static elec-tronic multipoles, water significantly alters the structure,stability, and dynamics of biomolecules [1–3]. The fa-vorable exposure of moieties with static electronic mul-tipoles to water and the corresponding burying of non-polar residues into a hydrophobic core, is also an im-portant, two-fold driving force for protein folding: First,the (short-range) interaction with the aqueous environ-ment is optimized and, second, the disruption of the dy-namic hydrogen bond network of the surrounding wa-ter by hydrophobic residues is minimized [3–5]. Whilethe importance of this “hydrophobic effect” and the piv- otal role of water for biomolecular systems is underno dispute [4–7], the underlying fundamental physicsof solvated (bio)molecular systems is still not fully ex-plored and understood. In particular, here we focus onthe quantum-mechanical nature of solute-solvent inter-actions. It has already been shown that polarization ef-fects and the many-body character of bonded interac-tions and hydrogen bond networks play an importantrole for solvated systems [3, 8–11], but also long-rangevan der Waals (vdW) dispersion interactions form an es-sential component for water and for both intra -proteinand protein-water interactions. This vdW component,however, has not been investigated in full detail nor on afundamental level up to now. In the present study, weaddress exactly this intricate issue and find that these1 a r X i v : . [ phy s i c s . c h e m - ph ] O c t uantum-mechanical interactions can account for up to30 % of the total solvation energy. Together with theiressential role for intra - and inter-protein interactions,this calls for a more complete microscopic understand-ing of vdW dispersion forces under physiological condi-tions, which is imperative to shed light on the physics ofproteins in aqueous solvation.Non-covalent vdW dispersion interactions arise frominstantaneous, correlated fluctuations of the electrondensity and, hence, are inherently quantum-mechanicaland many-body in nature. However, current solva-tion models or molecular mechanics force fields includethem only in a phenomenological manner via pairwise-additive potentials. Recent studies demonstrate that,for a variety of systems, including polypeptides in gasphase [11, 12], such approximations can fundamentallyfail and that one has to explicitly account for the intrin-sic many-body character of vdW interactions [12–21]– even for the properties of pristine water [22] and itssurfaces [23]. Moreover, pairwise or phenomenologi-cal models only provide a very approximate descriptionof the energetic aspects of vdW forces, but no descrip-tion of the underlying quantum mechanics. Such in-sights, however, can be vital to comprehend and con-ceptually understand vdW interactions as recently illus-trated for hybrid and nanostructured systems [24, 25] or π – π stacked molecules [18]. It is important to note thatthe dielectric permittivity of water has a value around2.3 [26] at the frequencies of the electronic fluctua-tions, which are responsible for dispersion interactions, i.e. at several hundred terahertz (THz). Therefore, incontrast to static electronic multipoles, vdW interac-tions are not strongly screened by aqueous environmentsand, thus, can give rise to long-range interactions alsoin solvated systems. As such, long-range correlationforces may play an important role for the long-rangeordering often observed in biological systems or formthe quantum-mechanical basis for the emergence of co-herent molecular vibrations [16, 27]. Such collectivenuclear behavior has been proposed to play an impor-tant role in long-distance recognition among biologicalmacromolecules [28–31]. Within the conventional viewof solvated proteins, however, the basis for long-rangerecognition under physiological conditions is still con- troversially discussed. Recent studies also suggest con-nections between collective electronic fluctuations – thebasis of vdW dispersion interactions – and enzymaticaction on DNA [32, 33] or pharmaceutical activity [34].In this work, we aim at a more complete microscopicdescription of solvated proteins and report an explicitlyquantum-mechanical study, which accounts for many-body dispersion interactions as well as for electrostat-ics, polarization, and hydrogen bonding. Molecular me-chanics approaches are often found to reliably capturethe latter three, (semi-)classical interactions. This, how-ever, usually holds true only for certain conditions (typ-ically designed for the liquid phase at room temperatureand ambient pressure), which shows that the physicaldescription is incomplete. As especially the quantum-mechanical, non-local dispersion interactions are oftendescribed insufficiently and inconclusively by conven-tional approaches [12, 35], we focus on a comprehensivedescription of vdW interactions and collective electronicbehavior to highlight the role of water for the vdW ener-getics of prototypical fast-folding proteins. In this con-text, previous studies have pointed out that traditionalmolecular mechanics potentials and water models likelyprovide an unbalanced description of vdW interactionsfor proteins in water, which typically results in anover-compaction of unfolded states [36–38]. Typically,this unbalanced description is approached by adaptingthe pairwise vdW interaction coefficients for the intra -protein, water-water, or protein-water interaction. Inthis work, we seek to understand the potential funda-mental basis for the failure of the traditional models ina bottom-up approach. The calculations have been car-ried out using a combined approach [39] of the Density-Functional Tight-Binding (DFTB) method [40, 41] andaccurate ab initio dispersion models, which allows fora robust, albeit approximate, quantum-mechanical treat-ment on an atomistic level. In fact, an ab initio de-scription of vdW interactions is the only way to studythe role of the solvent, as force field methods are typ-ically strongly limited in their transferability betweengas and liquid phase due to their high degree of pa-rameterization. In our study, we focus on a compari-son of a pairwise-additive description of vdW interac-tions and an accurate, quantum-mechanical many-body2reatment. The pairwise vdW models are representedby the vdW(TS) [42] as well as Grimme’s D2 [43] andD3 [44, 45] approaches. Such an approximate, pair-wise formalism represents the basis for the standard phe-nomenological description of long-range correlation inbiomolecular simulations via Lennard-Jones potentials.For comparison, we study the vdW interaction withinthe Many-Body Dispersion (MBD) formalism [13, 46],which accounts for the many-body character of vdWdispersion interactions to infinite order in perturbationtheory within an interatomic framework and has beenproven to provide quantitative improvements and a bet-ter qualitative understanding compared to the pairwise-additive approximation in numerous studies [12–19].Yet, its computational efficiency together with modernimplementations and the ever-growing availability ofcomputational resources, allows for treatment of sys-tems consisting of several thousands of atoms as in thecase of proteins in explicit solvent. The MBD formalismalso represents a model for collective electronic fluctua-tions [18, 24, 25] – a molecular analogue to the plasmonpseudo-particle in metallic systems – which we will useto further characterize the protein-water interaction. Results
We exemplify our findings in detail for the Fip35 Hpin1WW domain (Fip35-WW) [47] and further showcasetheir general validity at hand for the de novo
Chignolinvariant “cln025” [48] and the fast-folding Nle/Nle dou-ble mutant of the villin headpiece (HP35-NleNle) [49].The folding trajectories of the three proteins have beenobtained in atomistic detail and explicit solvent in previ-ous molecular dynamics simulations by Shaw et al. [50],Lindorff-Larsen et al. [51], and Ensign et al. [52, 53], re-spectively.
Intra -protein interactions and many-body dis-persion effects
We start out by investigating the Fip35-WW trajectoryby artificially removing the surrounding solvent from amolecular dynamics trajectory in explicit water, i.e. wefirst focus on intra -protein interactions, where disper-sion forces represent one of the main sources of interac- tion within the protein core. Accordingly, we observe anincreased magnitude of the vdW energy while this coreis being formed and particularly during the hydropho-bic collapse in all applied dispersion models (see Fig. 1and Supplementary Material). Notably, in comparisonto the results obtained within the pairwise approaches,many-body dispersion effects consistently decrease therelative stability of the native state for the isolated pro-tein by 6 kcal/mol on average, cf . Fig. 1(bottom). Theoutliers of this general behavior observed around 15 and26 μ s correspond to transient, partly folded intermedi-ates. The relative destabilization by beyond-pairwisecontributions can be explained by an overestimation ofthe intra -core vdW interactions (“overcorrelation”) inthe pairwise approximation. By reducing the interactionto pairwise-additive potentials, a two-body formulationassumes ideal correlation between all pairs of atoms andwith that, neglects the complex geometrical arrangementwithin the protein core. Such geometrical constraintslimit the emergence of correlated fluctuating dipole pat-terns and thus lower the interaction energy as already ob-served for a wide variety of systems [14, 15, 17, 18, 54].For small peptides such effects have been found to bemostly negligible [35, 55]. Our findings show that forlarger biomolecules, however, a many-body treatmentof vdW interactions is indeed essential, which is in linewith the findings of Schubert et al. for 20-residue pep-tides [12].In the MBD formalism, we make use of a two-stepprocedure: We obtain effective, screened atomic polar-izabilities from self-consistent electrodynamic screen-ing to account for the presence of multiple fluctuat-ing dipoles in the system and then solve a many-bodyHamiltonian, which is defined in terms of these polar-izabilities, to capture many-body vdW interactions. Tostudy the effect of each step, we combined vdW(TS)with the self-consistent screening procedure. In thisvariant, which we refer to as vdW(TS)@SCS, screenedinteraction coefficients enter the pairwise-additive po-tentials instead of the hybridized chemical analogueused in vdW(TS). In this way, we account for the ef-fects on atomic polarizabilities due to the local field ofthe surrounding dipoles, but do not include long-rangemany-body interactions. Within vdW(TS)@SCS we al-3 ig. 1: top: van der Waals energy of Fip35 Hpin1 WW-domain in solvated geometry without solvent. bottom: Beyond-pairwise contributions, as given by difference between the many-body formalism MBD and the pairwise treatments vdW(TS)and vdW(TS)@SCS, i.e. vdW(TS) with self-consistent screening. For a more comprehensive version, see SupplementaryMaterials. ready capture some part of the destabilization of nativestates amounting to 3 kcal/mol ( cf . Fig. 1). Thus, half ofthe overstabilization in vdW(TS) is from neglecting thepresence of multiple dipoles in the system and half fromthe many-body character of dispersion interactions. Thisalso implies that for a proper description of gas phaseproteins, one has to account for the screening of polariz-abilities and the many-body nature of vdW interactions. van der Waals solvation energy Fig. 2(top) shows the vdW contribution to the solvationenergy obtained with MBD and the vdW(TS) model, asdefined by E (sol) = E vdW [ps] − E vdW [p] − E vdW [s] , (1)with ps referring to Fip35-WW in solvation, p to Fip35-WW in gas phase, and s to the pristine solvent. As an artifact of the above-mentioned over-correlation withinthe pairwise approach, we find a consistent overesti-mation of the dispersion contribution in vdW(TS). Interms of the relative solvation energy, however, pairwiseand many-body treatment show the same general trend,which qualitatively follows the inverse root-mean squaredeviation (RMSD) from the native state with a step co-inciding with the collapse of the protein into the native,more globular shape. This finding can be explained bythe removal of hydrophobic residues from the protein-water interface and thus decreasing their interaction withthe solvent. The average dispersion contribution to thesolvation energy drops by 29 kcal/mol (15 %) at the hy-drophobic collapse. The step-like behavior of the vdWsolvation energy along the trajectory is even more dis-tinct than for the intra -protein vdW interaction in thegas phase and almost resembles a two-state model offolded and unfolded states. As such, the vdW solva-4 ig. 2: Relative vdW solvation energy, E (sol)rel , during the folding process of the Fip35 Hpin1 WW domain. top : backboneroot-mean-square deviation from final conformation illustrating the hydrophobic collapse around 35 μ s (gray). The vdWcontribution to the relative solvation energy is shown for the pairwise vdW(TS) model (red) and MBD (blue). bottom :Difference in the relative stabilization by the solvent between MBD and the pairwise vdW(TS) and D3 referenced to theunfolded state. tion energy prominently captures the protein’s collapseand, thus, represents a valid descriptor for the foldingprocess. This feature has been found for all dispersionmodels considered here: the MBD formalism and thepairwise approaches vdW(TS), D2, and D3.Comparing many-body and pairwise treatment of dis-persion interactions, Fip35-WW does no longer feature aconsistent change in the relative stability of native versusnon-native conformations once embedded in an aque-ous environment. Thus, beyond-pairwise effects in theprotein-water vdW interaction stabilize folded confor-mations. Correspondingly, we see a clear increase in therelative vdW solvation energy for native and native-likestates, when comparing the pairwise models to MBD(5 kcal/mol for vdW(TS), 7 kcal/mol for D3). Thisshift is due to the lack of a systematic many-body (de-)stabilization in the total vdW energy of solvated Fip35-WW and the water box during the whole folding tra-jectory, combined with an inversion of the behavior ob- served for the isolated protein shown in Fig. 1 (see Eq. 1and Supplementary Materials). This implies that theprotein-water interaction compensates for the destabi-lization of native states via many-body dispersion ef-fects, observed in vacuo . In summary, besides screeningpermanent electronic multipoles, water also provides thenecessary environment to stabilize native conformationsvia beyond-pairwise vdW interactions, which counter-acts the destabilizing effect that such many-body termshave on the intra -protein interaction. The plasmon-like character of protein-watervdW interactions
As has been shown previously, MBD also pro-vides a model for the intrinsic electronic fluctua-tions [18, 24, 25]: In analogy to nuclear quantumeffects, the electronic fluctuations, which ultimatelygive rise to vdW interactions, can be understood asthe zero-point oscillations around the average electron5ensity. The MBD formalism gives access to anorthonormal decomposition of this zero-point fluctua-tion, which can be interpreted as “eigenmodes” of theelectron density. A detailed analysis of these electroniceigenmodes reveals, that the number of very localized,high-frequency oscillations, formerly mainly locatedon the solute, significantly decreases upon coupling tothe surrounding water. This implies a delocalization ofelectronic fluctuations and an increase of the collectivityof electronic behavior. This plasmon-like character andthe delocalization over protein and solvent form thefundamental reason for the stabilization of native statesvia many-body dispersion effects in the protein-waterinteraction. The role of the surrounding solvent can be seen as providing weakly structured polarizablematter, which attenuates the destabilizing many-bodyeffects observed in vacuo for native and partially foldedstates. To gain further insight into the characteristics ofprotein-water vdW interactions, we additionally obtainthe contribution of individual electronic fluctuationsto the solvation energy: The main contributions to thevdW solvation energy are due to highly collective,plasmon-like electronic oscillations around 450 THz( ˆ=
670 nm), such as the one depicted in Fig. 3a).Their contribution to the total protein-water interactionis 14–16 % ( ≈
41 kcal/mol) and exceeds 20 % ( ≈ Fig. 3: a)
Illustration of low-frequency plasmon-like fluctuations in solvated Fip35 Hpin1 WW domain, which show the largestcontribution to the protein-water interaction. The arrows depict the direction of simultaneous electron density deformations(eigenmode of the electron density). b) Contributions to the vdW solvation energy within the pairwise vdW(TS) approach andthe many-body dispersion formalism (MBD) as radial distribution functions.
6s demonstrated in Fig. 3a), such fluctuations com-monly feature large charge displacements along the po-larizable protein backbone coupled to electronic fluctua-tions throughout the solvent. Notably, these plasmon-like fluctuations reach from the protein backbone in-side the hydrophobic core far into the aqueous environ-ment. Comparing the radial distribution of the contribu-tions to the vdW solvation energy between the pairwisevdW(TS) and MBD models, as shown in Fig. 3b), re-veals a striking difference in the interaction range withinthe two treatments: In the pairwise model, the contri-bution of solvent atoms to the vdW solvation energysubsides beyond 6 ˚A, i.e. roughly twice the sum of thevdW radii of carbon and oxygen. Accounting for many-body dispersion, on the other hand, shows that electroniccorrelation between the protein and solvent atoms up to25 ˚A from the protein-water interface is still relevant forthe protein-water interaction. This reflects the weaknessof the screening of dispersion forces by the solvent and isin evident contrast to the often assumed locality of vdWinteractions in solvated systems. While such a range isunprecedented in the context of solvated systems, simi-lar and larger interaction ranges have already been foundfor molecular crystals [19] or nanostructures [14]. Froma different point of view, Fig. 3b) represents a radial analysis of the change in the distribution and frequencyof electronic fluctuations introduced by embedding theprotein in water. It thus demonstrates that, while theatomistic structure and the local dynamics of water typi-cally remain largely unperturbed beyond a few solvationlayers [3, 56], the instantaneous electronic structure canindeed indicate the presence of a protein over large dis-tances. Such long-range collective electronic behaviorand the ensuing non-local nuclear dynamics could beexperimentally measured thanks to recent advances inultrafast THz spectroscopy [57–59], as we will furtherdiscuss below.
Effect of secondary structure
Fip35-WW is a showcase example for the formation of β -sheets. To test the general validity of our hypothe-ses, we carried out the same analysis for the modifiedvillin headpiece, HP-35 NleNle, (formation of α -helicalentities) and the cln025 variant of the de novo proteinChignolin (plain β -hairpin formation). Our analysis con-firmed our early findings for Fip35-WW. The vdW con-tribution to the solvation energy reflects the trend of theinverse RMSD from the native structure with a drop of15 % at the hydrophobic collapse. Again, the protein- Fig. 4:
Characteristics of protein-water dispersion interactions: Independent of the secondary structure, the van der Waalssolvation energy captures the hydrophobic collapse in form of a 20–30 kcal/mol jump (“ ∆ E sol at collapse”) and many-bodyprotein-water van der Waals interactions consistently stabilize native states in solvation. Low frequency, collective electronicfluctuations contribute significantly to the relative solvation energy ( ∆ E (low ω )sol ) in all cases. ∆ E (low ω )sol ) for both, the hairpin-forming cln025and the helix-forming HP35-NleNle.Fig. 4 summarizes the above-mentioned features andhighlights the general validity of the present findings inthe biomolecular context. Independent of the secondarystructure to be formed, the vdW solvation energy cap-tures the hydrophobic collapse in the form of a sizablejump in stabilization (“ ∆ E sol at collapse”) for all con-sidered proteins. Also, the consistent increase of the rel-ative stability of native states due to the many-body char-acter of protein-water vdW interactions and the signifi-cance of low-frequency, plasmon-like electronic fluctu-ations (as characterized by their contribution to relativesolvation energies, ∆ E (low ω )sol ) turned out to be indepen-dent of the final secondary structure. It is worthwhileto note that the magnitudes of the different quantitiesare not necessarily representative for secondary struc-ture elements as also the system size varies from 6,000to 14,000 atoms. The systematic investigation of therelation between secondary structure elements and themagnitude of the present observations and characteristicfeatures of fluctuation patterns, as shown in Fig. 3a), isbeyond the scope of the present publication and subjectto ongoing studies. Discussion and Outlook
In conclusion, we have shown that many-body disper-sion effects lead to a significant relative destabilization( ≈ intra -protein interactions in general.In aqueous solvation, the vdW contribution to thesolute-solvent interaction of (small) proteins capturestheir hydrophobic collapse and thus represents a viabledescriptor for the folding process. The collapse is ac-companied by a jump of about 15 % (20–30 kcal/mol) inthe vdW solvation energy. The total electronic energy ofsolvation, for comparison, does not provide such clearinsight (see Supplementary Materials) – only the freeenergy of solvation does. The beyond-pairwise contribu-tions to the protein-water vdW interaction favor foldedstates and, thus, the many-body aspect of solvation leadsto a significant stabilization of native conformations.This stabilization is the result of a distinct many-bodycharacter of the protein-water dispersion interaction inthe form of delocalization and a high degree of collectiv-ity of electronic fluctuations across protein and solvent.This plasmon-like character of vdW interactions in sol-vated systems can also be pivotal for other quantities, asdemonstrated for the protein-water interaction range inFig. 3b). Our study shows, that these findings can begeneralized for helix-, β -sheet-, or hairpin-forming pro-teins and are, thus, independent of secondary structuremotifs. So, an accurate description of solvated proteinsin general requires capturing the subtle balance betweenbeyond-pairwise effects on the intra -protein vdW inter-action (destabilizing native states) and the highly col-lective, plasmon-like character of protein-water interac-tions (stabilizing native states). With increasing systemsize and complexity, finding this balance without explicitaccount for the quantum-mechanical many-body natureof vdW interactions is an intricate task and failure to doso can contribute to the fundamental origin of the previ-ously reported [37, 38] unbalanced description of vdWforces by pairwise molecular mechanics potentials andwater models.8n a broader perspective, our findings imply that an ef-fective pairwise-additive treatment of vdW interactionsand derivative molecular mechanics potentials can pro-vide accurate energetics for a particular application, butonly explicitly quantum-mechanical models, such as theone used here, allow to attain an unbiased microscopicunderstanding of biomolecular interactions in a moregeneral context. Based on Fig. 3b), for instance, we ex-pect that for several or larger solutes the approximationof pairwise additivity can fail on a fundamental level.The persistence and collectivity of electronic fluctua-tions through the solvent can mediate long-range cor-relation forces between individual solutes or entities ofsolvated macromolecules and an effective pairwise de-scription might not be able to reproduce the subtle bal-ance between long-range correlation on all these scales.As such, plasmon-like interactions can also substantiallyaffect molecular assembly and the formation of tertiarystructures. In addition, complex long-range fluctuationsof the electron density are less sensitive to the instan-taneous solvent structure and thus thermal fluctuations,which makes them an ideal contender for biomolecu-lar recognition. In this form of recognition, the solventprovides electron density that serves as a mediator forlong-range interaction while the actual atomistic struc-ture and the nuclear dynamics of the solvent do not nec-essarily have to be altered in the process, which has beenconcluded from a number of experiments [3, 56]. It isworthwhile to mention, however, that most of these ex-periments probed rather local interactions and dynam-ics. Recent THz-spectroscopy experiments, for instance,show that the presence of a solute can have a consid-erable effect on the long-timescale dynamics and long-range polarization of water [57–59]. The here observedlong-range persistence of electron correlation throughaqueous environments will manifest in the slow nucleardynamics of the system. Similar behavior has alreadybeen observed in crystalline molecular systems, wheremany-body dispersion effects particularly affect low-frequency (“slow”) phonon modes [16, 27]. Long-rangeelectronic correlation between (solvated) biomoleculescan also form the quantum-mechanical basis for cor-related collective nuclear motion within the respectivepartners. Such concerted motion is essential for co- ordinated enzymatic action [32, 33] or the emergenceof coherent molecular vibrations, a promising expla-nation for long-range recognition through electrody-namic interaction of the resulting oscillating molecu-lar dipoles [30, 31]. Solvent-mediated plasmon-like in-teractions can also give a quantum-mechanical founda-tion for the recent proposal by Melkikh and Meijer [62]of the existence of long-range interactions guiding pro-tein folding, assembly, and organization in cells. Alto-gether, our findings in fact apply in a broader context ofbiomolecular interactions – not just in the case of proteinfolding as exemplified here.Obviously, the stability and functionality ofbiomolecules is ultimately determined by their freeenergy. Hence, this work represents a first step towardsa more fundamental understanding of the physicsof proteins in water, but to accurately address theimplications of plasmon-like features within biologicalsystems, we need to extend our study to free energy atfinite temperature. It is already known that the many-body character of vdW dispersion interactions cansignificantly soften low-frequency vibrational modesin organic matter, which has a noticeable effect on thesystem entropy [16, 27]. In the case of polymorphs ofthe aspirin crystal, for instance, many-body dispersioneffects introduce a relative entropic stabilization of0.6 kcal/mol (per molecule in the unit cell) at roomtemperature [16]. This can be seen as an estimate forthe effect on a protein’s relative entropy per residue. Asthe dynamics and functionality of a biomolecule can bestrongly related to its eigenmodes [63, 64], the impacton vibrational modes also hints at an unrevealed role ofplasmon-like electronic fluctuations for the functionalityand coordination in the biochemical machinery. Besidesthe proposed roles in long-distance recognition, enzy-matic action [32, 33], and pharmaceutical activity [34],this further strengthens the relevance of collectiveelectronic fluctuations in biomolecular systems. OurDFTB+MBD framework provides a robust formalismfor the investigation of such intricate questions as itallows for a fully quantum-mechanical treatment oflarge-scale systems in atomistic detail.9 ethods For each snapshot along the folding trajectories, weobtain effective atomic polarizabilities, α A , as usedwithin MBD and vdW(TS), in a seamless and non-empirical formalism via net atomic populations as de-scribed in Reference [39]. Electronic structure calcu-lations have been carried out in the Density-FunctionalTight-Binding framework with self-consistent charges(SCC-DFTB) [40, 41] using a locally modified ver-sion of DFTB+ [65, 66]. In vdW(TS), the effec-tive isotropic polarizabilities define effective C inter-action coefficients for pairwise-additive interaction po-tentials [39, 42]. In MBD then, α A is subject to short-range electrodynamic screening and then enters the cou-pled fluctuating dipole model [13]. So, electron densityfluctuations in a system consisting of N atoms are mod-eled by a set of N three-dimensional, isotropic quantumharmonic oscillators (QHOs) in dipole coupling definedby the Hamiltonian, H MBD = T ζ + 12 ζ T V ζ , (2) V ( i,j ) AB = ω A ω B h δ AB + p ˜ α A ˜ α B D ( i,j ) AB i . (3)In equation (2) T is the kinetic energy, V the potentialenergy matrix, and ζ the direct sum of mass-weighteddisplacements of the individual QHOs. V AB is de-fined by the characteristic excitation frequencies, ω , thescreened atomic polarizabilities, ˜ α , and the long-rangedipole coupling tensor D AB [13, 15]. In summary,the main approximations within MBD are the coarse-graining of the system’s response to an atom-centeredframework and the dipole approximation for the cou-pling between electronic fluctuations. Unitary transfor-mation of the Hamiltonian (2) to a new set of collectivevariables, ξ = C ζ , such that C † V C = diag (cid:8) ˜ ω i (cid:9) , (4)transforms equation (2) into an uncoupled set of N one-dimensional QHOs with frequencies ˜ ω i and dis-placements ξ i . It therefore provides a model for in-trinsic collective charge density fluctuations – a molec-ular analogue to the plasmon pseudo-particle in metal-lic systems. The dispersion energy has been shown to equal the zero-point interaction energy of this set ofQHOs [46, 67], given by E MBD = 12 N X i =1 ˜ ω i − N X A =1 3 X i =1 ω A . (5)The contribution of an individual collective electronicfluctuation, ξ i , to the total vdW solvation energy, as de-fined in equation (1), is obtained as the i -th element ofthe vector, εεε int = 12 U † [ C ps ] {U [ C ps ] ˜ ω ˜ ω ˜ ω ps − ˜ ω ˜ ω ˜ ω sub − U [ ] ( ωωω ps − ωωω sub ) } . (6)The transformation matrix, U [ Y ] , is given by theelement-wise absolute square of ( C p ⊕ C s ) † Y . C cor-responds to the transformation matrix used in equa-tion (4) and ˜ ω ˜ ω ˜ ω to the vector of characteristic frequen-cies. The subscript sub always refers to the direct sumof the corresponding quantity for the isolated protein andthe pristine solvent, e.g. ωωω sub = ωωω p ⊕ ωωω s . Note thatthe above transformation preserves energy and thus thesum of all elements in εεε int equals the total solvation en-ergy. For further information on the procedure, see Ref-erence [18]. Using the above definitions, we may alsoobtain the spatial distribution of plasmon-like interac-tions relevant for the vdW solvation energy. The radialMBD interaction range, G int [MBD] , shown in Fig. 3b)is calculated via, G int [MBD] = X i ε (int) i X A δ R,R A X j ∈ A k C ( i,j )ps k (7)Here, ε (int) i is the contribution of mode ξ i to E sol , i.e. the i -th element of εεε int as defined by equation (6), and δ cor-responds to a generalized Kronecker delta on R . C ( i,j )ps denotes the elements of the transformation matrix C ps ,which define the x, y, and z components of the subvec-tor of ξ i that resides on atom A . Note that integrating G int [MBD] or G int [vdW(TS)] yields the correspond-ing interaction (solvation) energy. Further details and in-depth analysis of the physics of plasmon-like electronicfluctuations in solvated (bio)molecular systems will beprovided in a future publication. For additional theoret-ical background and computational details, see Supple-mentary Materials.10 cknowledgements The authors thank D. E. Shaw Research for provid-ing the trajectories of cln025 and Fip35-WW used inthis work and Jan Hermann for insightful discussions.M.S. acknowledges financial support from the FondsNational de la Recherche, Luxembourg (AFR PhD GrantCNDTEC). A.T. was supported by the European Re-search Council (ERC-CoG Grant BeStMo). The resultspresented in this publication have been obtained usingthe HPC facilities of the University of Luxembourg.
References [1] F. Xu, T. A. Cross,
Proceedings of the NationalAcademy of Sciences of the United States of America , 9057 (1999).[2] A. Patriksson, C. M. Adams, F. Kjeldsen, R. A. Zubarev,D. van der Spoel, The Journal of Physical Chemistry B , 13147 (2007).[3] M.-C. Bellissent-Funel, et al. , Chemical Reviews ,7673 (2016).[4] W. Kauzmann,
Advances in Protein Chemistry (Elsevier,1959), vol. , pp. 1–63.[5] D. Chandler, Nature , 640 (2005).[6] B. Honig, A.-S. Yang,
Protein Stability (AcademicPress, 1995), vol. of Advances in Protein Chemistry ,pp. 27–58.[7] A. Fersht,
Structure and Mechanism in Protein Science (World Scientific, 2017).[8] I. Shvab, R. J. Sadus,
The Journal of Chemical Physics , 194505 (2013).[9] G. A. Cisneros, et al. , Chemical Reviews , 7501(2016).[10] J. J. Dannenberg,
Peptide Solvation and H-Bonds (Aca-demic Press, 2005), vol. of Advances in ProteinChemistry , pp. 227–273.[11] M. Rossi, W. Fang, A. Michaelides,
The Journal ofPhysical Chemistry Letters , 4233 (2015).[12] F. Schubert, et al. , Physical Chemistry ChemicalPhysics , 7373 (2015).[13] A. Tkatchenko, R. A. DiStasio, R. Car, M. Scheffler, Physical Review Letters , 236402 (2012).[14] V. V. Gobre, A. Tkatchenko,
Nature Communications ,2341 (2013). [15] A. Ambrosetti, D. Alf`e, R. A. DiStasio, A. Tkatchenko, The Journal of Physical Chemistry Letters , 849(2014).[16] A. M. Reilly, A. Tkatchenko, Physical Review Letters , 55701 (2014).[17] A. M. Reilly, A. Tkatchenko,
Chemical Science , 3289(2015).[18] J. Hermann, D. Alf`e, A. Tkatchenko, Nature Communi-cations , 14052 (2017).[19] J. Hoja, A. M. Reilly, A. Tkatchenko, Wiley Interdisci-plinary Reviews: Computational Molecular Science ,e1294 (2017).[20] J. F. Dobson, T. Gould, G. Vignale, Physical Review X , 1 (2014).[21] L. Kronik, A. Tkatchenko, Accounts of Chemical Re-search , 3208 (2014).[22] A. P. Jones, F. S. Cipcigan, V. P. Sokhan, J. Crain, G. J.Martyna, Physical Review Letters , 227801 (2013).[23] V. P. Sokhan, A. P. Jones, F. S. Cipcigan, J. Crain, G. J.Martyna,
Proceedings of the National Academy of Sci-ences of the United States of America , 6341 (2015).[24] W. Gao, Y. Chen, Q. Jiang,
Physical Review Letters ,246101 (2016).[25] A. Ambrosetti, N. Ferri, R. A. DiStasio, A. Tkatchenko,
Science , 1171 (2016).[26] H. Yada, M. Nagai, K. Tanaka,
Chemical Physics Letters , 279 (2009).[27] G. Folpini, et al. , Physical Review Letters , 097404(2017).[28] H. Fr¨ohlich,
International Journal of Quantum Chem-istry , 641 (1968).[29] G. Acbas, K. A. Niessen, E. H. Snell, A. G. Markelz, Nature Communications , 3076 (2014).[30] J. Preto, M. Pettini, J. A. Tuszynski, Physical Review E , 52710 (2015).[31] I. Nardecchia, et al. , Physical Review X , 031061(2018).[32] P. Kurian, G. Dunston, J. Lindesay, Journal of Theoreti-cal Biology , 102 (2016).[33] P. Kurian, A. Capolupo, T. Craddock, G. Vitiello,
Physics Letters A , 33 (2018).[34] T. J. Craddock, et al. , Scientific Reports , 9877 (2017).[35] C. Baldauf, M. Rossi, Journal of Physics: CondensedMatter , 493002 (2015).
36] M. Knott, R. B. Best,
PLoS Computational Biology ,e1002605 (2012).[37] R. B. Best, W. Zheng, and J. Mittal, Journal of ChemicalTheory and Computation TheJournal of Physical Chemistry B , 5113 (2015).[39] M. St¨ohr, G. S. Michelitsch, J. C. Tully, K. Reuter, R. J.Maurer,
The Journal of Chemical Physics , 151101(2016).[40] D. Porezag, T. Frauenheim, T. K¨ohler, G. Seifert,R. Kaschner,
Physical Review B , 12947 (1995).[41] M. Elstner, et al. , Physical Review B , 7260 (1998).[42] A. Tkatchenko, M. Scheffler, Physical Review Letters , 73005 (2009).[43] S. Grimme,
Journal of Computational Chemistry ,1787 (2006).[44] S. Grimme, J. Antony, S. Ehrlich, H. Krieg, The Journalof Chemical Physics , 154104 (2010).[45] S. Grimme, S. Ehrlich, L. Goerigk,
Journal of Compu-tational Chemistry , 1456 (2011).[46] A. Tkatchenko, A. Ambrosetti, R. A. DiStasio, TheJournal of Chemical Physics , 074106 (2013).[47] F. Liu, et al. , Proceedings of the National Academyof Sciences of the United States of America , 2369(2008).[48] S. Honda, et al. , Journal of the American Chemical So-ciety , 15327 (2008).[49] J. Kubelka, T. K. Chiu, D. R. Davies, W. A. Eaton,J. Hofrichter,
Journal of Molecular Biology , 546(2006).[50] D. E. Shaw, et al. , Science , 341 (2010).[51] K. Lindorff-Larsen, S. Piana, R. O. Dror, D. E. Shaw,
Science , 517 (2011).[52] D. L. Ensign, P. M. Kasson, V. S. Pande,
Journal ofMolecular Biology , 806 (2007).[53] D. L. Ensign, P. M. Kasson, V. S. Pande, Molecular Sim-ulation Trajectories Archive of a Villin Variant (2007).[54] S. Yang, Y. Jiang, S. Li, W. Liu,
Carbon , 513(2017).[55] M. Rossi, S. Chutia, M. Scheffler, V. Blum,
The Journalof Physical Chemistry A , 7349 (2014).[56] F. Persson, P. S¨oderhjelm, B. Halle,
The Journal ofChemical Physics , 215104 (2018). [57] S. Ebbinghaus, S. J. Kim, M. Heyden, X. Yu,U. Heugen, M. Gruebele, D. M. Leitner, M. Havenith,
Proceedings of the National Academy of Sciences of theUnited States of America , 20749 (2007).[58] K. Meister, S. Ebbinghaus, Y. Xu, J. G. Duman, A. De-Vries, M. Gruebele, D. M. Leitner, M. Havenith,
Pro-ceedings of the National Academy of Sciences of theUnited States of America , 1617 (2013).[59] Y. Xu, M. Havenith,
The Journal of Chemical Physics , 170901 (2015).[60] A. Tkatchenko, M. Rossi, V. Blum, J. Ireta, M. Scheffler,
Physical Review Letters , 118102 (2011).[61] M. Rossi, M. Scheffler, V. Blum,
The Journal of Physi-cal Chemistry B , 5574 (2013).[62] A. V. Melkikh, D. K. Meijer,
Progress in Biophysics andMolecular Biology , 57 (2018).[63] R. M. Levy, D. Perahia, M. Karplus,
Proceedings of theNational Academy of Sciences of the United States ofAmerica , 1346 (1982).[64] W. Zheng, H. Wen, Current Opinion in Structural Biol-ogy , 24 (2017).[65] B. Aradi, B. Hourahine, T. Frauenheim, The Journal ofPhysical Chemistry A , 5678 (2007).[66] The development version of the DFTB+ code, whichhas been used in the current work, is available athttps://github.com/aradi/dftbplus.git (mbd branch).[67] J. F. Dobson, J. Wang, B. P. Dinte, K. McLennan, H. M.Le,
International Journal of Quantum Chemistry ,579 (2005).[68] M. Gaus, Q. Cui, M. Elstner,
Journal of Chemical The-ory and Computation , 931 (2011).[69] M. Gaus, A. Goez, M. Elstner, Journal of Chemical The-ory and Computation , 338 (2013).[70] M. Gaus, X. Lu, M. Elstner, Q. Cui, Journal of ChemicalTheory and Computation , 1518 (2014).[71] M. Kubillus, T. Kubaˇr, M. Gaus, J. ˇRez´aˇc, M. Elstner, Journal of Chemical Theory and Computation , 332(2015).[72] S. J. Clark, et al. , Zeitschrift f¨ur Kristallographie - Crys-talline Materials (2005).[73] S. Bahn, K. Jacobsen,
Computing in Science & Engi-neering , 56 (2002).[74] B. Aradi, GitHub - aradi/dftd3-lib: Library version of S.Grimme’s DFTD3 code. (2016). upplementary Material Quantum Mechanics of Proteins in Explicit Water:The Role of Plasmon-Like Solute-Solvent Interactions
Martin St ¨ohr and Alexandre Tkatchenko ∗ Physics and Materials Science Research Unit, University of Luxembourg,L-1511 Luxembourg, Luxembourg ∗ To whom correspondence should be addressed; E-mail: [email protected]
Contents a r X i v : . [ phy s i c s . c h e m - ph ] O c t Computational Details and van der Waals Models
Electronic structure calculations have been carried out on the Density-Functional Tight-Binding (DFTB)level of theory using a locally modified MPI-version of DFTB+ [65,66]. For the study of the Fip35Hpin1 WW-domain and the Nle/Nle mutant of villin HP35 we employed the second-order (self-consistentcharges) DFTB method with recent mio-1-1 parameters [41]. For the calculations on the Chignolin vari-ant cln025 third-order DFTB [68] (DFTB3) with recent 3ob parameters [69–71] has been used. Ef-fective atomic polarizabilities within the DFTB3 framework, as further used in our study, have beentested to comply with the results from the second-order approach. In all calculations we employed aself-consistency criterion of 10 -5 elementary charges and Γ -point sampling. In our study on many-body dispersion effects on protein energetics and the collectivity of van derWaals (vdW) interactions in solvated biosystems, we have employed a combined approach of DFTBand the Many-Body Dispersion (MBD) formalism [13,46]. The results from MBD have been comparedwith common pairwise approaches to vdW dispersion interactions, as they are commonly employed instate-of-the-art simulation techniques. Pairwise models are represented by the electronic structure-basedvdW(TS) [42] and Grimme’s D2 [43] and D3 [44,45] method. The central starting point of MBD andvdW(TS) is the definition of effective atomic dipole polarizabilities, α ( I )eff , according to the chemical en-vironment. As presented in Reference [39], these can be derived from DFTB as, α ( I )eff = α ( I )free Z I · X i ∈ I P ii , (SI-1)where P is the Mulliken population matrix as obtained from DFTB. In vdW(TS), the obtained effec-tive polarizabilities then define effective C interaction coefficients, which then enter pairwise additiveC · R − IJ potentials. The same functional form is used in Grimme’s D2 and D3 schemes. Here, one relieson fixed (D2) or geometry-dependent (D3) interaction coefficients.vdW(TS) calculations have been carried out using a standalone calculator based on ’ semp disp corr.F90 ’,originally part of CASTEP [72], within the Atomic Simulation Environment [73]. Calculations involv-ing D2 or D3 have been performed using the DFTD3 module [74]. For D3 calculations, Becke-Johnsondamping as proposed in Reference [45] has been used.In MBD, α ( I )eff is first subject to self-consistent, electrodynamic screening (SCS) to account for thepresence of the surrounding fluctuating atomic dipoles and obtain effective, screened atomic polarizabil-ities ˜ α ( I )eff . This is achieved by inverting the Dyson-like equation for the dynamic polarizabilities, ˜ α ( I )eff ( iω ) = X J B IJ , with B = (cid:2) A − ( iω ) + T gg (cid:3) − , (SI-2)where A ( iω ) is the diagonal matrix of the dynamic polarizabilities α ( I )eff ( iω ) at imaginary frequency iω and T gg is the effective short-range dipole potential constructed from the Coulomb interaction of twooverlapping Gaussian charge densities. In the MBD formalism we model the interaction of intrinsic2lectronic fluctuations in form of a set of dipole-coupled Quantum Harmonic Oscillators (QHOs). Inaccordance with this, the dynamic polarizability and effective excitation frequency, ω I , are defined as α ( I )eff ( iω ) = α ( I )eff (cid:16) ωω I (cid:17) and ω I = 43 α ( I )eff α ( I )free C ( II )6 , free α ( I )free , (SI-3)respectively. The effective, screened polarizabilities then define the MBD Hamiltonian, H MBD = T ζ + 12 ζ T V ζ , with V ( k,l ) IJ = ω I ω J h δ IJ + (1 − δ IJ ) p ˜ α I ˜ α J D ( k,l ) IJ i . (SI-4)In eq. [SI-4], T is the kinetic energy, V the potential energy matrix, and ζ the direct sum of mass-weighted displacements of the individual QHOs. V IJ is defined by the characteristic excitation frequen-cies, screened polarizabilities, and a damped dipole coupling tensor (Fermi-damping) D IJ [13,46]. AllMBD and vdW(TS)@SCS calculations have been performed using a self-written implementation.3 Summary of Gas-phase Energetics for Fip35 Hpin1 WW-domain
As can be seen from Fig. S1, all considered vdW models show comparable behavior along the foldingtrajectory. The main differences between the approaches can be found for the relative stability of nativestates with respect to unfolded conformations. Hereby, the pairwise models overestimate the stability asseen from Fig. S1(bottom). The average over-stabilization is 3 kcal/mol in Grimme’s D3 and vdW(TS)with screened C interaction coefficients, “vdW(TS)@SCS”. For D2 and conventional vdW(TS) we findan overestimation of the relative stability by 4 and 6 kcal/mol, respectively.Figure S1: van der Waals dispersion energy of Fip35 Hpin1 WW-domain in gas-phase as obtained withall considered van der Waals models (top). Beyond pairwise effects on gas-phase energetics (bottom).4 van der Waals Energetics in Detail Upon embedding in an aqueous environment, we do not observe systematic differences in the relativevdW energetics ( cf.
Fig. S2). Also in the case of the pristine solvent, the pairwise models do not showconsistent deviations from MBD. Interestingly, for the geometry-based D2 and D3, we additionally finda subtle relative over-stabilization of the pristine solvent corresponding to native state conformations.These add up with minor, yet non-systematic, discrepancies in the total vdW energy and ultimately leadto a overall underestimation of the relative solvation energy of native conformations by 12 kcal/mol forD2 and 7 kcal/mol for D3 (see main manuscript). The latter represents a surprising shortcoming whencompared to the electronic structure-based vdW(TS) as D3 outperforms all other pairwise models for gas-phase vdW energies ( cf.
Fig. S1). The spurious description of the pure solvent in the geometry-motivatedmodels can most likely be ascribed to an inaccurate description of vdW interactions involving hydrogenatoms located near the cavity and thus mainly affects the pristine solvent. We conclude that, for anaccurate description of edge effects in water, it is essential to include electronic-structure effects. In termsof the description of protein-water vdW interactions, ultimately, the pairwise approaches systematicallylack a relative stabilization of 5 to 12 kcal/mol due to neglect of beyond-pairwise interaction terms.Figure S2: Total relative van der Waals dispersion energy of Fip35 Hpin1 WW-domain in solvation asobtained using the many-body dispersion formalism and the pairwise approaches vdW(TS), D2, and D3(top). Beyond pairwise effects on relative van der Waals energetics of the Fip35 Hpin1 WW-domain inaqueous solvation (bottom). 5 .2 Chignolin variant “cln025”
For the “cln025” variant of the de novo protein Chignolin, we observe a similar behavior for the gas-phase energetics. The dispersion interaction energy is maximized when going to the native state andthe pairwise vdW model vdW(TS) overestimates the relative stability of the native state of cln025 in theabsence of water in comparison to MBD by 5 kcal/mol (see Fig. S4). Upon embedding in water, we stillfind a slight destabilization of native states via many-body dispersion effects (2 kcal/mol), while for thepure solvent we did not observe a statistically relevant change in the relative vdW energy along the foldingtrajectory. The discrepancy in the relative vdW solvation energy of cln025 between pairwise and many-body treatment, as depicted in Fig S3, is thus governed by the neglect of many-body dispersion effectson intra -protein vdW energetics during folding and a slight destabilization of native states of cln025in solvation via many-body effects. The increase of the relative stabilization through beyond-pairwiseprotein-water vdW interactions amounts to 2 kcal/mol.Figure S3: Relative van der Waals solvation energy of cln025 as obtained with the many-body disper-sion model and the pairwise vdW(TS) approach (top). Beyond pairwise contributions to van der Waalssolvation energy as given by difference between many-body and pairwise treatment (bottom).6igure S4: Relative van der Waals energy of the Chignolin variant cln025 in absence of solvent as ob-tained with the many-body (MBD) and pairwise (vdW(TS)) treatment of dispersion interaction (top).Beyond pairwise contributions to relative van der Waals energetics (bottom).
For HP35-NleNle, many-body dispersion effects destabilize the native state in absence of water, i.e. reduce the intra -protein vdW interaction, by 4 kcal/mol in comparison to the pairwise vdW(TS). Forsolvated HP35-NleNle and the pure solvent we did not observe a considerable consistent many-bodyeffect on the energetics. In sum, the pairwise formalism underestimates the protein-water vdW interactionof HP35-NleNle by about 3 kcal/mol when compared to a full many-body treatment.7
Correlation and Rescaling of van der Waals Solvation Energies
As can be seen from Fig. S1 and the main manuscript, the overestimation of vdW energies increases withthe absolute vdW interaction energy. Correspondingly, a simple rescaling of the pairwise approachesconsiderably improves the agreement with the many-body treatment. Fig. S5 shows the correlation be-tween such optimally rescaled vdW solvation energies and MBD. The obtained rescaling factors showthat relying on electronic-structure based C interaction coefficients as done within vdW(TS) providesthe best estimate for vdW solvation energies. This can mainly be attributed to the description of thepure solvent as the geometry-based D2 and D3 methods outperform vdW(TS) for gas-phase energetics( vide supra ). Despite the overall improvement, the deviation between the optimally rescaled pairwiseapproaches and MBD still regularly exceeds 4 kcal/mol. Furthermore, the optimal rescaling factors arehighly system- and method-dependent and can only be obtained as an a posteriori correction.Figure S5: Correlation of rescaled relative van der Waals solvation energies as obtained from pairwisemodels in comparison to the results obtained from many-body treatment. In contrast to the vdW solvation energy, the total electronic energy of solvation does not provide a clear-cut distinction between folded and unfolded states ( cf.cf.