Quantum Mechanics of Proteins in Explicit Water: The Role of Plasmon-Like Solute-Solvent Interactions

TL;DR

This study employs a quantum-mechanical approach to reveal the significant role of many-body dispersion interactions and plasmon-like solute-solvent effects in protein stability and dynamics within water, challenging traditional pairwise models.

Contribution

It introduces a density-functional tight-binding method with many-body dispersion formalism to explicitly account for collective dispersion effects in proteins and water, highlighting their importance.

Findings

Many-body dispersion effects decrease native state stability without water.

Aqueous environment stabilizes native conformations via collective dispersion interactions.

Long-range electron correlation persists through water, influencing protein behavior.

Abstract

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 dynamics of proteins in aqueous solution. Typically, these interactions are only treated phenomenologically via pairwise potential terms in classical force fields. Here, we use an explicit quantum-mechanical approach 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 dispersion interactions for protein energetics and the protein-water interaction. In contrast to commonly employed pairwise approaches, many-body effects significantly decrease the relative stability of the native state in the absence of water. In an aqueous environment, the collective character of the protein-water van der Waals interaction counteracts this effect and stabilizes native conformations and transition states. This stabilization arises due to a high degree of delocalization and collectivity of protein-water dispersion interactions, suggesting a remarkable persistence of long-range electron correlation through aqueous environments. Our findings are exemplified on prototypical showcases of proteins forming $\beta$-sheets, hairpins, and helices, emphasizing the crucial role of plasmon-like solute-solvent interactions in biomolecular systems.

Equations

$${\alpha }_{\mathrm{eff}}^{(I)}=\frac{{\alpha }_{\mathrm{free}}^{(I)}}{Z_I}\cdot \sum _{i\in I}{\mathbf{P}}_{ii},$$

$${\tilde{\alpha }}_{\mathrm{eff}}^{(I)}\left(i\omega \right)=\sum _J{\mathbf{B}}_{IJ},\mathrm{with}\mathbf{B}={\left[{\mathbf{A}}^{-1}(i\omega )+{\mathbf{T}}_{\mathrm{gg}}\right]}^{-1},$$

$${\alpha }_{\mathrm{eff}}^{(I)}(i\omega )=\frac{{\alpha }_{\mathrm{eff}}^{(I)}}{1+{\left(\frac{\omega }{{\omega }_I}\right)}^2}\mathrm{and}{\omega }_I=\frac{4}{3}\frac{{\alpha }_{\mathrm{eff}}^{(I)}}{{\alpha }_{\mathrm{free}}^{(I)}}\frac{C_{6,\mathrm{free}}^{(II)}}{{\alpha }_{\mathrm{free}}^{(I)}},$$

$${\mathcal{H}}_{\mathrm{MBD}}={\mathcal{T}}_{\zeta }+\frac{1}{2}{\zeta }^T\mathcal{V}\zeta ,\mathrm{with}{\mathcal{V}}_{IJ}^{(k,l)}={\omega }_I{\omega }_J\left[{\delta }_{IJ}+\left(1-{\delta }_{IJ}\right)\sqrt{{\tilde{\alpha }}_I{\tilde{\alpha }}_J}{\mathcal{D}}_{IJ}^{(k,l)}\right].$$

Full text

1. 1 Computational Details and van der Waals Models

2. 1.1 Density-Functional Tight-Binding calculations

3. 1.2 van der Waals Dispersion Models

4. 2 Summary of Gas-phase Energetics for Fip35 Hpin1 WW-domain

5. 3 van der Waals Energetics in Detail

6. 3.1 Fip35 Hpin1 WW-domain

7. 3.2 Chignolin variant “cln025”

8. 3.3 Villin headpiece 35 Nle/Nle mutant

9. 4 Correlation and Rescaling of van der Waals Solvation Energies

10. 5 Total Electronic Energy of Solvation

1 Computational Details and van der Waals Models

1.1 Density-Functional Tight-Binding calculations

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 Fip35 Hpin1 WW-domain and the Nle/Nle mutant of villin HP35 we employed the second-order (self-consistent charges) DFTB method with recent mio-1-1 parameters [41]. For the calculations on the Chignolin variant cln025 third-order DFTB [68] (DFTB3) with recent 3ob parameters [69–71] has been used. Effective atomic polarizabilities within the DFTB3 framework, as further used in our study, have been tested to comply with the results from the second-order approach. In all calculations we employed a self-consistency criterion of 10-5 elementary charges and $\Gamma$-point sampling.

1.2 van der Waals Dispersion Models

In our study on many-body dispersion effects on protein energetics and the collectivity of van der Waals (vdW) interactions in solvated biosystems, we have employed a combined approach of DFTB and the Many-Body Dispersion (MBD) formalism [13,46]. The results from MBD have been compared with common pairwise approaches to vdW dispersion interactions, as they are commonly employed in state-of-the-art simulation techniques. Pairwise models are represented by the electronic structure-based vdW(TS) [42] and Grimme’s D2 [43] and D3 [44,45] method. The central starting point of MBD and vdW(TS) is the definition of effective atomic dipole polarizabilities, $\alpha_{\mathrm{eff}}^{(I)}$, according to the chemical environment. As presented in Reference [39], these can be derived from DFTB as,

[TABLE]

where $\mathbf{P}$ is the Mulliken population matrix as obtained from DFTB. In vdW(TS), the obtained effective polarizabilities then define effective C${}_{\mbox{6}}$ interaction coefficients, which then enter pairwise additive C{}_{\mbox{6}}$$\cdot R_{IJ}^{-6} potentials. The same functional form is used in Grimme’s D2 and D3 schemes. Here, one relies on 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 involving D2 or D3 have been performed using the DFTD3 module [74]. For D3 calculations, Becke-Johnson damping as proposed in Reference [45] has been used.

In MBD, $\alpha_{\mathrm{eff}}^{(I)}$ is first subject to self-consistent, electrodynamic screening (SCS) to account for the presence of the surrounding fluctuating atomic dipoles and obtain effective, screened atomic polarizabilities $\tilde{\alpha}_{\mathrm{eff}}^{(I)}$. This is achieved by inverting the Dyson-like equation for the dynamic polarizabilities,

[TABLE]

where $\mathbf{A}(i\omega)$ is the diagonal matrix of the dynamic polarizabilities $\alpha_{\mathrm{eff}}^{(I)}(i\omega)$ at imaginary frequency $i\omega$ and $\mathbf{T}_{\mathrm{gg}}$ is the effective short-range dipole potential constructed from the Coulomb interaction of two overlapping Gaussian charge densities. In the MBD formalism we model the interaction of intrinsic electronic fluctuations in form of a set of dipole-coupled Quantum Harmonic Oscillators (QHOs). In accordance with this, the dynamic polarizability and effective excitation frequency, $\omega_{I}$, are defined as

[TABLE]

respectively. The effective, screened polarizabilities then define the MBD Hamiltonian,

[TABLE]

In eq. [SI-4], $\mathcal{T}$ is the kinetic energy, $\mathcal{V}$ the potential energy matrix, and $\zeta$ the direct sum of mass-weighted displacements of the individual QHOs. \mathcal{V}$${}_{\mbox{IJ}} is defined by the characteristic excitation frequencies, screened polarizabilities, and a damped dipole coupling tensor (Fermi-damping) \mathcal{D}$${}_{\mbox{IJ}} [13,46]. All MBD and vdW(TS)@SCS calculations have been performed using a self-written implementation.

2 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 folding trajectory. The main differences between the approaches can be found for the relative stability of native states with respect to unfolded conformations. Hereby, the pairwise models overestimate the stability as seen from Fig. S1(bottom). The average over-stabilization is 3 kcal/mol in Grimme’s D3 and vdW(TS) with screened C${}_{\mbox{6}}$ interaction coefficients, “vdW(TS)@SCS”. For D2 and conventional vdW(TS) we find an overestimation of the relative stability by 4 and 6 kcal/mol, respectively.

3 van der Waals Energetics in Detail

3.1 Fip35 Hpin1 WW-domain

Upon embedding in an aqueous environment, we do not observe systematic differences in the relative vdW energetics (cf. Fig. S2). Also in the case of the pristine solvent, the pairwise models do not show consistent deviations from MBD. Interestingly, for the geometry-based D2 and D3, we additionally find a 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 lead to a overall underestimation of the relative solvation energy of native conformations by 12 kcal/mol for D2 and 7 kcal/mol for D3 (see main manuscript). The latter represents a surprising shortcoming when compared 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-motivated models can most likely be ascribed to an inaccurate description of vdW interactions involving hydrogen atoms located near the cavity and thus mainly affects the pristine solvent. We conclude that, for an accurate description of edge effects in water, it is essential to include electronic-structure effects. In terms of the description of protein-water vdW interactions, ultimately, the pairwise approaches systematically lack a relative stabilization of 5 to 12 kcal/mol due to neglect of beyond-pairwise interaction terms.

3.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 and the pairwise vdW model vdW(TS) overestimates the relative stability of the native state of cln025 in the absence of water in comparison to MBD by 5 kcal/mol (see Fig. S4). Upon embedding in water, we still find a slight destabilization of native states via many-body dispersion effects (2 kcal/mol), while for the pure solvent we did not observe a statistically relevant change in the relative vdW energy along the folding trajectory. 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 effects on intra-protein vdW energetics during folding and a slight destabilization of native states of cln025 in solvation via many-body effects. The increase of the relative stabilization through beyond-pairwise protein-water vdW interactions amounts to 2 kcal/mol.

3.3 Villin headpiece 35 Nle/Nle mutant

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). For solvated HP35-NleNle and the pure solvent we did not observe a considerable consistent many-body effect on the energetics. In sum, the pairwise formalism underestimates the protein-water vdW interaction of HP35-NleNle by about 3 kcal/mol when compared to a full many-body treatment.

4 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 with the absolute vdW interaction energy. Correspondingly, a simple rescaling of the pairwise approaches considerably improves the agreement with the many-body treatment. Fig. S5 shows the correlation between such optimally rescaled vdW solvation energies and MBD. The obtained rescaling factors show that relying on electronic-structure based C${}_{\mbox{6}}$ interaction coefficients as done within vdW(TS) provides the best estimate for vdW solvation energies. This can mainly be attributed to the description of the pure 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 pairwise approaches and MBD still regularly exceeds 4 kcal/mol. Furthermore, the optimal rescaling factors are highly system- and method-dependent and can only be obtained as an a posteriori correction.

5 Total Electronic Energy of Solvation

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. Fig. S6 and S7). This is not surprising as the, in the end decisive, free energy of a solvated molecule has a large entropic component and it is known that the hydrophobic effect mainly arises from entropic contributions [5]. In the case of cln025, we do observe a slight shift in the total electronic energy of solvation when comparing folded and unfolded states, see Fig. S7. Considering the absolute spread of the solvation energy, however, this shift is less clear as for its mere vdW component (Fig. S3) and coincides with the spread of total solvation energies of the unfolded conformations.