Terahertz absorption of lysozyme in solution

TL;DR

This study investigates terahertz radiation absorption in lysozyme solutions, revealing that traditional electrostatics models fail to accurately describe protein-water interface polarization, with simulations aligning better with experimental data.

Contribution

The paper introduces an analytical and simulation-based approach showing the inadequacy of Maxwell electrostatics and highlights the significance of the protein hydration shell in terahertz absorption.

Findings

Lorentz void model better fits experimental data

Simulation results agree with experimental measurements

Hydration shell acts as a distinct sub-ensemble

Abstract

Absorption of radiation by solution is described by the solution dielectric constant and can be viewed as a specific application of the dielectric theory of solutions. For ideal solutions, the dielectric boundary value problem separates the polar response into the polarization of the void in the liquid created by the solute and the response of the solute dipole. In the case of a protein as a solute, its nuclear dynamics do not project on significant fluctuations of the dipole moment in the terahertz domain of frequencies and the protein dipole can be viewed as dynamically frozen. Absorption of radiation then reflects the interfacial polarization. Here we apply an analytical theory and computer simulation to absorption of radiation by ideal solutions of lysozyme. Comparison with experiment shows that Maxwell electrostatics fails to describe the polarization of the protein-water interface and the "Lorentz void", which does not anticipate polarization of the solute void by the external field (no surface charges), better represents the data. An analytical theory for the slope of the solution absorption against the volume fraction of the solute is formulated in terms of the cavity field response function. It is calculated from molecular dynamics simulations in good agreement with experiment. The protein hydration shell emerges as a separate sub-ensemble, which collectively is not described by the standard electrostatics of dielectrics.

Equations

$$\Delta \chi (\omega )\propto {\chi }_{0s}(\omega )$$

$$M_0^{\text{M}}=-{\Omega }_{0}P\frac{3{ϵ}_s}{2{ϵ}_s+1}.$$

$$M_0^{\text{ind}}=M_0^e+⟨M_0{⟩}_E.$$

$${\alpha }_{\text{abs}}(\omega )=\frac{\omega }{c}\frac{{ϵ}^{\prime \prime }(\omega )}{\sqrt{{ϵ}^{\prime }(\omega )}}.$$

$$M_x=P(V-N_0{\Omega }_0)-N_0{\Omega }_{0}P\frac{{ϵ}_s-1}{2{ϵ}_s+1}.$$

$$\begin{array}{rl}{M}_x(\omega )& =P(\omega )(V-N_0{\Omega }_0)\\ & -N_0{\Omega }_0\alpha (\omega )P(\omega )\frac{{ϵ}_s(\omega )-1}{2{ϵ}_s(\omega )+1}.\end{array}$$

$$\frac{\Delta {\chi }_{\text{v}}(\omega )}{{\chi }_s(\omega )}=-{\eta }_0\left[1+\alpha (\omega )\frac{{ϵ}_s(\omega )-1}{2{ϵ}_s(\omega )+1}\right],$$

$$\frac{\Delta {\chi }_{\text{v}}(\omega )}{{\chi }_s(\omega )}=-{\eta }_0\frac{3{ϵ}_s(\omega )}{2{ϵ}_s(\omega )+1}.$$

$${\chi }_c^{\text{M}}(\omega )=\frac{3}{2{ϵ}_s(\omega )+1}$$

$$1+\alpha (\omega )\frac{{ϵ}_s(\omega )-1}{2{ϵ}_s(\omega )+1}=-\frac{3{ϵ}_s(\omega )}{2({ϵ}_s(\omega )-1)}\frac{{\chi }_{0s}(\omega )}{{\chi }_{00}(\omega )}.$$

$$4\pi \Delta {\chi }_{\text{v}}(\omega )=\frac{3}{2}{\eta }_0{ϵ}_s(\omega )\frac{{\chi }_{0s}(\omega )}{{\chi }_{00}(\omega )}.$$

$$S_{0a}(t)={\left[⟨\delta {\mathbf{M}}_0\cdot \delta {\mathbf{M}}_a⟩\right]}^{-1}⟨\delta {\mathbf{M}}_0(t)\cdot \delta {\mathbf{M}}_a(0)⟩,$$

$${\chi }_{0a}(\omega )=(\beta /3{\Omega }_0)⟨\delta {\mathbf{M}}_0\cdot \delta {\mathbf{M}}_a⟩\left[1+i\omega {\tilde{S}}_{0a}(\omega )\right].$$

$${\tilde{S}}_{0a}(\omega )=\sum _i\frac{A_i{\tau }_i}{1-i\omega {\tau }_i}$$

$${\chi }_{0a}(\omega )\propto \sum _i{A}_i{\left[1-i\omega {\tau }_i\right]}^{-1}.$$

$${\chi }_0^d(\omega )={\chi }_{00}(\omega )+{\chi }_{0s}(\omega )$$

$$y_0(\omega )=(4\pi /3){\chi }_{00}(\omega )+y_e,$$

$$y_e=\frac{1}{{\eta }_0^p}\frac{n_p^2-1}{n_p^2+2}.$$

$$\Delta \chi =\Delta {\chi }_{\text{v}}+\Delta {\chi }_0.$$

$$\frac{4\pi \Delta {\chi }_0}{{\eta }_0{y}_0}=3{ϵ}_s{\chi }_c+J_0({\eta }_0)\left[3{ϵ}_s{\chi }_c-{ϵ}_s-2\right],$$

$${\chi }_c(\omega )=\frac{{\chi }_0^d(\omega )}{{\chi }_{00}(\omega )}=1+\frac{{\chi }_{0s}(\omega )}{{\chi }_{00}(\omega )},$$

$$J_0({\eta }_0)=\frac{3}{\pi }{\int }_0^{\infty }[j_1(x/2){]}^2\left[S_0(x)-1\right]dx.$$

$$\Delta {\chi }_0=N_0\frac{\beta ⟨\delta {\mathbf{M}}_0^2⟩}{3V}\frac{3{ϵ}_s}{2{ϵ}_s+1}.$$

$$4\pi \Delta \chi (\omega )=-\frac{3}{2}{\eta }_0{ϵ}_s(\omega )\left[1-{\chi }_c(\omega )\frac{3n_p^2}{n_p^2+2}\right],$$

$$4\pi \Delta {\chi }^{\text{M}}(\omega )=-\frac{3}{2}{\eta }_0{ϵ}_s(\omega )\left[1-\frac{9n_p^2}{(2{ϵ}_s(\omega )+1)(n_p^2+2)}\right].$$

$$4\pi \Delta {\chi }^{\text{L}}(\omega )=-\frac{3}{2}{\eta }_0{ϵ}_s(\omega )\left[1-\frac{{ϵ}_s(\omega )+2}{{ϵ}_s(\omega )}\frac{n_p^2}{n_p^2+2}\right].$$

$$\Delta {\bar{\alpha }}_{\text{abs}}(\omega )={\alpha }_{\text{abs}}(\omega )/{\alpha }_{\text{abs}}^s(\omega )-1,$$

Full text

I Introduction

Inserting a solute into a molecular liquid leads to a nonlinear perturbation of its density.Chandler (1993) The density profile formed around a solute is a sensitive function of the strength of the intermolecular potential. Weak dewettingChandler (2005); Sarupria and Garde (2009) or, alternatively, collapse of hydration shellsBrovchenko and Oleinikova (2007, 2008); Friesen and Matyushov (2012) occurs for hydration when the solute-water interaction is either weaker (dewetting) or stronger (collapse) than the water-water interaction.

In contrast to the diversity of scenarios for the density profile of hydration layers, the orientational structure of the solute-solvent interface is typically viewed as a specific case of the standard dielectric boundary-value problem. The dipoles of water then orient along the lines of the solute’s electric field according to the rules of electrostatics of dielectrics.Jackson (1999) There are a number of fundamental reasons why this simplistic picture can break down,Matyushov (2014) particularly in the case of patterned interfaces of nano-scale dimension.Giovambattista et al. (2008) By the patterned interface we mean a substrate with its polarity or charge changing on the length-scale comparable to the correlation length between the dipoles in the surrounding water. The most illuminating example is a patchwork of positive and negative charges, each orienting the nearest water molecules along the field of the surface charge. If the distance between different charges is comparable with the size of the surface polarized domain, this spatial arrangement will lead to the orientational frustration of the water molecules at the domain boundaries. One can anticipate that the response of this patchwork of frustrated domains can potentially be different from the response of a homogeneous dielectric.

The picture of a patchwork of positive and negative charges distributed over a nanometer-scale surface is a good match of the protein-water interface. Proteins possess a large density of surface ionized residues, which provides the free energy required to stabilize a protein in solution. The distribution of positive and negative charges at the surface of a globular protein is nearly uniform, with the average distance between the charges $\sim 10$ Å.Barlow and Thornton (1986); Takashima (2002) This distance is comparable to the length of a chain of 3–4 water molecules. When the two waters at the ends of the chain are pinned by the protein charges, the one-two molecules in the middle must be strongly frustrated orientationally. The result is the separation of the hydration shell into polarized domains with their low-temperature behavior reminiscent of the phenomenology of relaxor ferroelectrics close to the glass transition.Martin and Matyushov (2015) The key question relevant to this observation is that of the depth of propagation of these domains into the bulk, or, in other words, of the number of water molecules involved into the hydration layer with properties distinct from the bulk. Given that $\sim(300-500)$ molecules of water can be counted only in the first hydration layer of a typical globular protein, changes in the orientational structure of even the second hydration layer will produce an ensemble of water molecules sufficiently large to be characterized as a separate mesophase ($\sim 1063$ water molecules within a 6 Å shell around lysozymeMartin and Matyushov (2015)). Such a shell will carry new properties requiring characterization.

Absorption of radiation in the THz domain of frequencies has recently appeared as a novel technique to study hydration shells of proteins.Zhang and Durbin (2006); Xu, Plaxco, and Allen (2006); Heugen et al. (2006); Ebbinghaus et al. (2007); Born et al. (2009); Niehues et al. (2011); Vinh, Allen, and Plaxco (2011); Acbas et al. (2014); Yamamoto et al. (2016); Moix, Parker, and Echchgadda (2017) The absorption of light as the means of probing the solvation structure has a significant advantage compared to broad-band dielectric spectroscopy since the sensitivity of the response is improved by a factor roughly equal to the dielectric constant for the same magnitude of the liquid polarization or by a factor of dielectric constant squared for the same magnitude of the external electric field.Matyushov (2010) The difference is due to the transverse geometry of light absorption compared to the longitudinal geometry of the dielectric experiment.

Further advantage of the THz domain spectroscopy over absorption and dielectric spectroscopy at lower frequencies is that the response of the protein dipole (relaxation time $\sim 9$ ns for lysozymeMartin and Matyushov (2014)) is dynamically frozen at the THz frequency. While there are a number of vibrational modes in the THz region belonging to the vibrational density of states of a protein,Yu, Park, and Leitner (2003); Lerbret et al. (2012) they do not project on the protein dipole to produce sufficient THz absorption.Matyushov (2012) The protein is therefore mostly transparent in the THz domain and main absorption comes from the water component of the solution. The challenge is how to separate the hydration shell from the dominant bulk absorption.

The analysis of the experimental data in terms of simplistic mixing model,Zhang and Durbin (2006); Heugen et al. (2006); Ebbinghaus et al. (2007); Luong et al. (2011); Sushko, Dubrovka, and Donnan (2015) disregarding non-additivity of the interfacial polarization (see discussion after Eq. (5) below), suggested the possibility of an extended hydration shell around a protein, with absorption distinct from the bulk.Ebbinghaus et al. (2007); Luong et al. (2011) It was later shown that the change in the absorption of solutions relative to bulk water can be related to the cross-correlation between the protein dipole and the dipoles of water molecules in the hydration shell.Matyushov (2012) Specifically, the alteration of the dipolar susceptibility of the solution relative to the bulk is directly proportional to the solute-solvent dipolar correlation function $\chi_{0s}(\omega)$

[TABLE]

(see Eq. (11) below for a complete form). When the cross-correlation $\chi_{0s}(\omega)$ is calculated for hydration shells of varying thickness, its length of saturation to the bulk value turns out to be in the range of 20–40 ÅMatyushov (2012); Heyden, Tobias, and Matyushov (2012) from the protein surface. Therefore, deviations from the ideal behavior with increasing solute concentrationEbbinghaus et al. (2007) are likely related to changes of $\chi_{0s}$ from its infinite-dilution value.

As expected from the general arguments presented above, simulations have shown that orientational correlations of water dipoles in the hydration layer are significantly different from those in a homogeneous dielectric. The major question, still posed to the field, is whether these differences can be reliably extracted from absorption of THz radiation with its wavelength far exceeding the length-scale of any structural correlations in the solution. Since long-wavelength radiation does not provide spatial resolution, theory is required to interpret the data.

The first question to ask is what would the standard electrostatics predict for the absorption coefficient of solutions. This problem can in fact be addressed by Maxwell electrostatics of dielectrics, by which we mean the application of the Maxwell boundary conditions to the solution of the Poisson equation.Jackson (1999) The problem at hand is illustrated in Fig. 1. The polarization of the interface of a spherical void approximating the hydrated protein integrates to an effective dipole assigned to the void, which is anti-parallel to the external field

[TABLE]

Here, $P$ is the polarization of the bulk caused by the uniform field of the electromagnetic wave, $\Omega_{0}$ is the volume of the spherical solute and $\epsilon_{s}$ is the dielectric constant of water; the superscript “M” specifies the solution of the Maxwell boundary value problem. In contrast to this prediction, simulationsMartin, Friesen, and Matyushov (2011); Matyushov (2012) and absorption measurements for hydrated lysozymeNovelli et al. (2017) and for some amino acidsNiehues et al. (2011) have suggested that the polarization of the interface can significantly deviate from the Maxwell scenario. We indeed show below that the so-called Lorentz void,Martin and Matyushov (2008) instead of the Maxwell void (Eq. (2)), provides a better representation of the collective interfacial polarization. However, this conclusion cannot be reached without taking into account the induced dipole of the protein itself.

The negative projection of the interface Maxwell dipole is counterbalanced by a positive induced solute dipole $M_{0}^{\text{ind}}$. The latter is composed of the instantaneous (on the THz time-scale) electronic induced dipole $M_{0}^{e}$ and the average permanent dipole $\langle M_{0}\rangle_{E}$ established along the field within an ensemble of proteins in solution

[TABLE]

The permanent dipole of course depends on the frequency of the signal and is expected to be dynamically frozen in the THz domain. The average over an ensemble of randomly oriented protein permanent dipoles then produce a negligible net value of $\langle M_{0}\rangle_{E}$. Nevertheless, dynamical freezing depends on the radiation frequency and should not be blindly assumed. The induced dipole cannot be generally neglected, as we indeed show below. Given the mutual compensation between the anticipated negative interface dipole and the solute induced dipoles, a quantitative theory is required to incorporate both components. Here, we construct such a theoretical description focused on experimental applications and based on previous theoretical advances.Matyushov (2010, 2012) The frequency dependence of the permanent protein dipole induced by the external field is supplied by molecular dynamics (MD) simulations. These results confirm that the protein permanent dipole is dynamically frozen in the THz frequency window, allowing one to neglect $\langle M_{0}\rangle_{E}$ in Eq. (3). The main compensation therefore happens between the positive induced electronic dipole and the negative interface dipole. We compare the results of our analysis to recent THz absorption data reported for lysozyme solutions.Novelli et al. (2017)

II Formalism

The absorption coefficient for electromagnetic radiation at frequency $\omega$ is fully defined in terms of the real and imaginary parts of the dielectric constant of the solution $\epsilon(\omega)=\epsilon^{\prime}(\omega)+i\epsilon^{\prime\prime}(\omega)$

[TABLE]

The problem is therefore reduced to calculating $\epsilon(\omega)$ in terms of the complex-valued dielectric constant of the solvent $\epsilon_{s}(\omega)$ (water in our case) and the properties of the solute (protein). The general formalism for $\epsilon(\omega)$ of an ideal solution was developed in Refs. Matyushov, 2010, 2012. Here we summarize the main steps of this theory relevant to our simulations and comparison to experiment.

We first consider a fictitious solute which does not produce any dipolar response and thus can be approximated as a void in a polar liquid. This approximation is relevant to many THz experiments when the solute dipole is too slow and is dynamically frozen. Computer simulations presented below specifically address the question of whether this approximation is applicable in the case of hydrated lysozyme. Independently of the outcome, this limit is the simplest conceptual point of departure for our discussion.

II.1 Voids in water

A void removes polar material from its volume thus reducing the overall dipole moment of the liquid from its homogeneous magnitude $M^{\text{liq}}_{x}=PV$ to the value $M_{x}^{\text{liq}}-N_{0}\Omega_{0}P$. Here, we have adopted the geometry of an absorption experiment sketched in Fig. 1 in which light is propagated with the wave-vector along the $z$-axis of the laboratory frame and its electric field polarized along the $x$-axis. Further, the solution contains $N_{0}$ solutes, each having the volume $\Omega_{0}$, $\eta_{0}=(N_{0}\Omega_{0})/V$ is the volume fraction.

Removing polarized medium from the void’s volume is not sufficient to calculate the polarization of the solution. Any discontinuous interface inside a dielectric is polarized by the field, and there will be an additional interfacial polarization, which requires solving the dielectric boundary value problem.Jackson (1999) In the case of a spherical void, the transverse geometry of the absorption experiment leads to the following relationMatyushov (2010)

[TABLE]

Here, the last summand accounts for an additional polarization of the interface not accounted for by the removal of the polarized liquid.

The first term in Eq. (5) projects on the reduction of absorption by removing water from the volume of the protein. This is the water term in the two-component analysis often used in the literature.Zhang and Durbin (2006); Heugen et al. (2006); Ebbinghaus et al. (2007); Luong et al. (2011); Sushko, Dubrovka, and Donnan (2015) In this approach, the absorption coefficient of the solution is represented as the volume-fraction weighted sum of the protein (p) and water (w) components: $\alpha_{\text{abs}}=(1-\eta_{0})\alpha_{w}+\eta_{0}\alpha_{p}$. This approximation can produce negative protein absorption when applied to analyzing experimental data.Vinh, Allen, and Plaxco (2011) This seemingly astonishing result is in fact trivial when compared to Eq. (5): the additional negative contribution, which is assigned to absorption by the protein in the two-component analysis, is caused by polarization of the void in addition to the removal of water (the second summand in Eq. (5)).

The two-component analysis often leads to seemingly reasonable results, but this comes not from the soundness of this approximation, but instead from the combination of the breakdown of the Maxwell dielectric picture for the heterogeneous protein-water interface and a potential compensating effect of the dipole induced at the solute ($M_{0}^{\text{ind}}$ in Eq. (3) and in Fig. 1, also see below).Matyushov (2012) Deviations between the two-component model and observations are often assigned to the third component, the hydration water. Given that such deviations, according to Eq. (5), are predicted even within the standard dielectric picture, any conclusions based on such an analysis are unfounded. Numerically, the additional polarization of the interface [second summand in Eq. (5)] is a significant effect, amounting to about a half, at $\epsilon_{s}\gg 1$, of the polarization loss from the direct expulsion of water from the solute volume.

The appearance of the dielectric constant of bulk water $\epsilon_{s}$ in Eq. (5) is a consequence of using the standard dielectric boundary conditions at the discontinuous surface separating the void from the polar liquid. Hydration water surrounding proteins is strongly altered compared to the bulk and it is a priori not clear if the standard dielectric boundary conditions apply. In fact, we have shown beforeMatyushov (2012) and demonstrate more conclusively here that polar response of the protein-water interface is qualitatively different from that of the standard dielectric interface. In order to account for these differences, we introduced the parameter $\alpha(\omega)$, which allows us to re-write Eq. (5) in the following form

[TABLE]

Here, we have indicated the dependence of all components on the radiation frequency $\omega$. It has to be stressed that, similarly to the dielectric constant $\epsilon_{s}(\omega)$, the parameter $\alpha(\omega)=\alpha^{\prime}(\omega)+i\alpha^{\prime\prime}(\omega)$ is a complex-valued function, as will be clear from the following discussion.

From Eq. (6), one can directly calculate the change in the susceptibility of the solution $4\pi\chi(\omega)=\epsilon(\omega)-1$ relative to the susceptibility of the bulk liquid $4\pi\chi_{s}(\omega)=\epsilon_{s}(\omega)-1$

[TABLE]

where we used the subscript “v” to stress that this susceptibility change is due to the ideal solution of voids. In the case of the Maxwell dielectric, $\alpha(\omega)=1$ and the solution susceptibility decreases linearly with the volume fraction $\eta_{0}$

[TABLE]

The slope of the linear scaling with $\eta_{0}$ is equal to the the product of $\epsilon_{s}$ with the cavity field susceptibility of dielectric theoriesBöttcher (1973); Fröhlich (1958)

[TABLE]

The susceptibility $\chi_{c}=E_{c}/E_{0}$ is the ratio of the field inside the void (cavity) $E_{c}$ to the uniform external $E_{0}$ (see below), the superscript “M” indicates the Maxwell solution for this property.Jackson (1999); Martin and Matyushov (2008)

Equation (7) does not accomplish much unless a connection between $\alpha(\omega)$ and parameters accessible by either simulations and/or experiment can be established. It was suggestedMatyushov (2012) that the effective dipole established at the surface of a spherical void polarized by the uniform external field can be alternatively calculated from the response to a dipole placed inside the void. The connection is provided by the relation

[TABLE]

Equation (7) is therefore re-written in the following form

[TABLE]

Equation (11) provides a significant advantage for the theory-experiment connection since it directly yields the THz slope ($\alpha_{\text{abs}}$ vs $\eta_{0}$), instead of absorption at a given concentration as calculated from simulations of the total dipole of the mixture.Heyden et al. (2008); Moix, Parker, and Echchgadda (2017) The required input is two susceptibility functions, $\chi_{00}(\omega)$ and $\chi_{0s}(\omega)$. They are accessible by computer simulations in the limit of a single solute in the simulation cell, thus bypassing the difficulties of simulating mixtures. The susceptibilities $\chi_{00}(\omega)$ and $\chi_{0s}(\omega)$ arise from, correspondingly, the “self”, $\propto\langle\delta\mathbf{M}_{0}(t)\delta\mathbf{M}_{0}(0)\rangle$, and “cross”, $\propto\langle\delta\mathbf{M}_{s}(t)\cdot\delta\mathbf{M}_{0}(0)\rangle$, time correlation functions of the fluctuating dipole moments of the solute, $\mathbf{M}_{0}$, and solvent, $\mathbf{M}_{s}$. Being linear response functions, $\chi_{00}(\omega)$ and $\chi_{0s}(\omega)$ represent the permanent dipole moment induced by the external electric field at the solute (subscript “00”) and the dipole moment caused in the hydration layer by the solute dipole (subscript “0s”).

The frequency-dependent susceptibilities are Laplace-Fourier $\omega$-transforms of the corresponding time correlation functions appearing in the dynamic linear response theory.Kubo (1966) One can define the normalized time correlation functions

[TABLE]

where $\delta\mathbf{M}_{a}(t)=\mathbf{M}_{a}(t)-\langle\mathbf{M}_{a}\rangle$ and $a=0,s$. These functions are Laplace-Fourier transformed to yield the functions $\tilde{S}_{0a}(\omega)$, which then enter the corresponding linear response functions

[TABLE]

The fitting procedure using a multi-exponential decay with amplitudes $A_{i}$ ($\sum_{i}A_{i}=1$) and relaxation times $\tau_{i}$ leads to

[TABLE]

and

[TABLE]

The cross correlations $\propto\langle\delta\mathbf{M}_{s}(t)\cdot\delta\mathbf{M}_{0}(0)\rangle$ are typically poorly converged and are a challenge to calculate from simulations. In anticipation of the final results used to calculate the absorption coefficients, we note that the response function

[TABLE]

enters the final expressions. This function is based on the time correlation function $\propto\langle\delta\mathbf{M}_{0}(t)\cdot\delta\mathbf{M}(0)\rangle$ between the solute dipole and the entire dipole moment of the sample $\mathbf{M}$, which is typically a faster converging property.

Since $\chi_{00}$ is always positive, Eq. (11) clearly shows that the sign of slope of $\Delta\chi_{\text{v}}$ vs $\eta_{0}$ is determined by the sign of the cross, solute-solvent dipolar correlations. The solute dipole and the dipoles of solvent molecules are expected to be anti-correlated, thus resulting in a usually observedZhang and Durbin (2006); Xu, Plaxco, and Allen (2006); Novelli et al. (2017) negative THz slope. The reason for such negative cross correlations is illustrated in Fig. 2, which shows the dipoles of the solvent arranged around a solute dipole. While axial solvent dipoles tend to orient in parallel, the equatorial dipoles tend to orient anti-parallel. Since there are always more equatorial dipoles that axial dipoles, the overall cross-correlation is negative. This picture is, however, based on the assumption that interfacial dipoles can freely change their orientations along the field lines of the solute dipole. Hydration shells of proteins are densely packedGerstein and Chothia (1996); Svergun et al. (1998) and orientationally frustrated,Friesen and Matyushov (2012); Martin and Matyushov (2015) thus potentially affecting the ability of hydration waters to polarize along the field of the dipole. While most solutes still show a negative slope of $\Delta\chi_{\text{v}}$ vs $\eta_{0}$, solutions of some amino acids (glycine and serineNiehues et al. (2011)) produce a positive THz slope. The latter result might be, at least partially, related to the compensating positive dipole moment induced at the solute (Eq. (3)), which we consider below.

In the analysis below, we use the input for the frequency-dependent dielectric constant of water from experimentVinh et al. (2015) and model the electronic induced dipole of the protein by its refractive index. The force fields typically used in MD simulations miss polarizability of at least one of these two components. It is not a priori clear if the force fields typically used in condensed phase simulations are suitable for modeling absorption. To illustrate this point, Fig. 3 compares the real and imaginary parts of $\epsilon_{s}(\nu)$ ($\omega=2\pi\nu$) for TIP3P water (dashed lines) with experimentVinh et al. (2015) (solid lines). The differences between experiment and TIP3P water are significant. In order to establish whether these differences affect absorption of radiation, one needs a formalism operating in terms of solvent parameters accessible by both experiment and simulations. Equation (11) presents this opportunity, and we discuss this analysis below. It turns out that the use of the dielectric constant of TIP3P water significantly deteriorates the agreement between calculations and experiment within the analytical (Lorentz void) framework, but does not strongly affect the results directly obtained from MD simulations.

II.2 Solute dipole

The solute dipole can include both the permanent and electronic components. As is typically done in mean-field theories of polarizable dielectrics,Wertheim (1979); Stell, Patey, and Høye (1981) it is convenient to combine them in one polar density parameterMatyushov (2012)

[TABLE]

where $y_{e}=(4\pi/3)(\alpha_{0}/\Omega_{0})$ and $\alpha_{0}$ is the electronic polarizability of the solute. Clausius-Mossotti equation can be used to estimate $y_{e}$

[TABLE]

Here, $n_{p}$ is the refractive index of the protein and $\eta_{0}^{p}$ is the packing fraction of proteins in the powder used to measure $n_{p}$. We will assume $\eta_{0}^{p}\simeq 1$ in our estimates below, which is also the assumption made by Onsager for molecular polarizability.Onsager (1936) Corrections can be introduced if this parameter is known independently.

The frequency-dependent dipolar response of the solute adds to the solution susceptibility, which is now defined as

[TABLE]

Here, the first summand, due to the void, is given by Eq. (11). The second terms, the change in the susceptibility due to dipolar response of the solutes, is given by the following equationMatyushov (2012)

[TABLE]

where here and in Eq. (19) the dependence on frequency has been dropped for brevity. The frequency-dependent response function $\chi_{c}(\omega)$ is

[TABLE]

where $\chi_{0}^{d}(\omega)$ given by Eq. (16) is the total susceptibility of the solute dipole in response to the uniform polarizing external field.Stern and Feller (2003); Friesen and Matyushov (2011) The physical meaning of $\chi_{0}^{d}/\chi_{00}$ is that it is equalMatyushov (2012) to the ratio, $E_{c}/E_{0}$, of the cavity field $E_{c}$ experienced by the solute dipole from the polarized solvent to the field of the external charges $E_{0}$ (the result of Maxwell electrostatics is given by Eq. (9)). Although a direct measurement of the cavity field is hardly possible, it enters a number of spectroscopic observables.Toptygin (2003) Correspondingly, the response function $\chi_{c}(\omega)$ represents the ratio of the cavity field to the external field at a given frequency.

The first summand in Eq. (20) is the infinite-dilution contribution of the solute dipole to the solution susceptibility. The second summand, which includes the function $J_{0}(\eta_{0})$, represents the non-ideal effect of polarizing the void by the permanent dipoles of other solutes in the solution. Therefore, $J_{0}(0)=0$, and it is given for an arbitrary $\eta_{0}$ by an equation involving the density structure factor of the solutes, $S_{0}(k)=N_{0}^{-1}\sum_{i,j}^{N_{0}}\exp[i\mathbf{k}\cdot\mathbf{r}_{ij}]$, where $\mathbf{r}_{ij}$ is the vector connecting two solutes in the solution and $\mathbf{k}$ is the wavevector. The density structure factor determines mutual correlations of the positions of the solutes and $J_{0}(\eta_{0})$ is given asMatyushov (2010)

[TABLE]

Here, $j_{1}(x)$ is the spherical Bessel function of the first orderAbramowitz and Stegun (1972) and integration is performed over the dimensionless variable $x=k\sigma_{0}$, where $\sigma_{0}$ is the solute diameter. While $S_{0}(k)$ is experimentally available from small-angle scattering, we estimated it here based on Percus-Yevick solution for a fluid of hard spheres.Hansen and McDonald (2003) These estimates (see $J_{0}(\eta_{0})$ tabulated in Ref. Matyushov, 2010) show that $J_{0}(\eta_{0})$ can be dropped from the final equations at the low volume fraction $\eta_{0}<0.1$. When $J_{0}=0$, one can apply Maxwell electrostatics ($\alpha=1$) as a consistency check of Eq. (20). One then gets Eq. (9) for $\chi_{0}^{d}/\chi_{00}$ and for the static ($\omega=0$) response

[TABLE]

This is the standard linear-response resultBöttcher (1973) for an ideal solution of impurities experiencing the cavity field $E_{c}=3\epsilon_{s}E/(2\epsilon_{s}+1)$ from the surrounding dielectric (cf. to Eq. (9), $E$ is the Maxwell field in the bulk).

III Simulation results

All-atom MD simulations were performed using the initial crystallographic structure of lysozyme (PDB entry 1AIK) solvated with TIP3P water. The simulation cell consisted of a total of 87050 atoms, with 28361 water molecules. The total charge of the lysozyme protein was $-7$ e. The constant pressure temperature equilibration simulations were done using the Langevin temperature-pressure control with the damping coefficient of 5 ps*-1*, a piston pressure of 1 atm, a piston decay time of 50 fs, a piston oscillation period of 200 fs, and at temperature of 300 K. A cutoff radius of 12  ̵̊Å and full electrostatics using the particle mesh Ewald technique at every simulations step were employed. NAMD 2.10Phillips et al. (2005) with the CHARMM27 force field was used to produce the MD trajectories. An initial optimization of the simulation cell was performed by conjugate gradient minimization for 2000 steps, followed by a 5 ns NPT equilibration simulation. Proceeding from equilibration run, production NVT simulations were carried out for 10 ns.

The main reason for a relatively short length of the simulation trajectory is that the simulations (both the equilibration NPT and the production NVT) were performed with a relatively short timestep of 0.5 fs and flexible hydrogen bonds of the protein. The trajectory was saved every 10 fs. This focus on the short-time dynamics has allowed us to calculate the time correlation functions with a 10 fs time resolution required for modeling the THz response. Exponential fits of the correlation function were done on the time window of 1 ns. Those were Laplace-Fourier transformed to obtain Eq. (14) for the correlation functions and Eq. (15) for the response functions.

IV Comparison to experiment

A significant result of our present and previousMatyushov (2012) computer simulations is that $(4\pi/3)\chi_{00}$ can be neglected in Eq. (17) relative to $y_{e}$ in the THz domain. The nuclear motions of the protein do not produce sufficient fluctuations of its dipole and protein absorption is negligible. Only the electronic dipole induced at the protein makes a non-negligible contribution. Note that, in contrast with the previous simulations,Matyushov (2012) the current setup allowed vibrations of protons in the protein, along with a significantly smaller integration step and more frequent savings. Nevertheless, $S_{00}(t)$ is nearly flat at $t\simeq 0.1-1$ ps (Fig. 4) and, correspondingly, $\chi_{00}(\omega)$ is very small for $\nu=\omega/(2\pi)$ in the THz domain of frequencies.

With this simplification and after dropping the term proportional to $J_{0}$ from Eq. (20), one arrives at a simple relation

[TABLE]

where the Clausius-Mossotti relation between the protein refractive index $n_{p}$ and $y_{e}$ was applied [$\eta_{0}^{p}=1$ in Eq. (18)]. In the limit of the Maxwell electrostatics (superscript “M”) one obtains

[TABLE]

At $n_{p}=1$, this relation converts to Eq. (8) for the susceptibility of a Maxwell void.

One can additionally consider the limit in which no surface charges, originating from introducing a dividing surface into the dielectric,Jackson (1999) are created at the interface of a void (no “$+$” and “$-$” at the void’s surface in Fig. 1). This limit corresponds to $\alpha(\omega)=0$ in Eq. (7) and the assumption that the only result of creating the void is the expulsion of the polarized water from the void’s volume. Such a void is known in the theory of dielectrics as the virtual, or Lorentz, cavity.Böttcher (1973) We will therefore dub this limit as the “Lorentz scenario”.Martin and Matyushov (2008) The ratio $\chi_{0}^{d}/\chi_{00}$ becomes equal to the Lorentz cavity field,Matyushov (2012) $E_{c}/E_{0}=(\epsilon_{s}+2)/(3\epsilon_{s})$, and one obtains from Eq. (24) (superscript “L” is for the Lorentz scenario)

[TABLE]

Not surprisingly, there is no change to the solvent response when the virtual void is filled with the same material as the solvent and $\Delta\chi^{\text{L}}=0$ at $\epsilon_{s}=n_{p}^{2}$. Note that the Lorentz cavity field is used in deriving the Clausius-Mossotti equation, and that fact leads to the complete cancellation of the interface and electronically induced solute dipoles at $\epsilon_{s}=n_{p}^{2}$. This does not happen in the Maxwell case in Eq. (25) since the cavity field of dielectric theories,Böttcher (1973) instead of the Lorentz field, is used in the Maxwell electrostatics.

In application to analyzing experimental results, the complex-valued dielectric constant of bulk water was taken from recent measurementsVinh et al. (2015) covering the range of frequencies $5.9-1120$ GHz. The experimental dataNovelli et al. (2017) for the absorption coefficient of lysozyme were obtained by averaging absorption over the interval of frequencies $0.38-0.92$ THz. There is therefore a sufficient frequency overlap between the data for bulk waterVinh et al. (2015) and solutions.Novelli et al. (2017) The refractive index of the protein, $n_{p}\simeq 1.55$, is from Ref. Cervelle, Cesbron, and Berthou, 1974.

Figure 5 presents the calculations of the change in the absorption coefficient of the solution with increasing the volume fraction of the solute according to Eq. (4). Specifically, we normalize the change in absorption by the absorption of bulk water

[TABLE]

where the absorption of water $\alpha_{\text{abs}}^{s}(\omega)$ is calculated from Eq. (4) with the water dielectric constant $\epsilon_{s}(\omega)$ used in place of the solution dielectric constant $\epsilon(\omega)$. The solution dielectric constant is calculated as $\epsilon(\omega)=\epsilon_{s}(\omega)+4\pi\Delta\chi(\omega)$, with three scenarios presented by Eqs. (24)–(26) used for $\Delta\chi(\omega)$.

Both the Maxwell and Lorentz voids are sensitive to the choice of the protein refractive index. The dashed lines in Fig. 5 refer to the calculations neglecting the induced dipole of the protein ($M_{0}^{\text{ind}}=0$ in Eq. (3)), which is achieved by putting $n_{p}=1$. There is, however, no physical reason to adopt this assumption and this limit is shown here as a mere warning that a “good agreement” with experiment can follow from unphysical approximations (the Lorentz scenario is in perfect agreement with experimental data in this case, blue dashed line). It is also clear that MD simulations of $\chi_{c}(\omega)$ provide a good account of the cavity field dynamics and are consistent with experiment. We additionally show in Fig. 6 the comparison of the calculations employing the experimentalVinh et al. (2015) dielectric constant $\epsilon_{s}(\omega)$ (solid lines) and the same function obtained from MD simulations of TIP3P water (dashed lines). While there is little sensitivity to the choice of the water model in simulations, which access the dynamics of the cavity field directly from $\chi_{c}(\omega)$, the Maxwell and Lorentz equations for the cavity field are obviously sensitive to the choice of the frequency-dependent dielectric constant function. The agreement with experiment is much worse when $\epsilon_{s}(\omega)$ for TIP3P water is used Eqs. (25) and (26). This is not surprising given large discrepancy in the dielectric properties between the experimental and TIP3P water shown in Fig. 3.

V Conclusions

Absorption of radiation provides access to the orientational structure of the hydration shell in terms of the dipolar solute-solvent (cross) correlation function [Eqs. (1) and (11)]. The structure of water can be disrupted to a different extent depending on the solute. The extent of the disruption and the resulting structure of the hydration shell affect the decay of the solute-solvent correlations into the bulk.Matyushov (2012) The decay length can be in the range of $20-40$ Å. The extended hydration shell anticipated in the pastEbbinghaus et al. (2007) in fact implies an extended cross-correlation. Nuclear motions of the protein do not contribute to the signal and protein itself is effectively transparent in the THz domain. Therefore, it is the protein-water interface that determines the change of THz absorption by solutions relative to the bulk.

The analytical modelMatyushov (2010, 2012) for the absorption of radiation by solutions gives the slope of the absorption coefficient vs the solute volume fraction in terms of the cavity response function $\chi_{c}(\omega)$ [Eq. (24)], which is the ratio of the field inside the solute to the external field of radiation measured at the radiation frequency. Simulations of a single lysozyme in solution performed here are capable of producing this function and lead to a good agreement with experimentNovelli et al. (2017) (Fig. 5).

In the absence of direct simulations, two analytical limits, corresponding to the Maxwell and Lorentz voids, can be produced [Eqs. (25) and (26)]. We find, in agreement with previous calculations,Heyden, Tobias, and Matyushov (2012) that the Maxwell void gives a poor representation of the data, while the Lorentz void gives a better account of experiment. One still has to realize that the appearance of the Lorentz limit for the void polarization is not consistent with the standard electrostatics, but is an emergent consequence of the orientational and density restructuring of the protein hydration layer compared to bulk water. It is therefore the deviation from the Maxwell limit that should be considered as the specific effect of the protein hydration shell. The Lorentz scenario is an effective-medium representation of the complex structure of the protein-water interface, which requires modification of the standard dielectric boundary-value problem.

Acknowledgements.

This research was supported by the National Science Foundation (CHE-1464810) and through XSEDE resources (TG-MCB080116N).