Cancellation between auto- and mutual correlation contributions of protein/water dynamics in terahertz time-domain spectra.

Terahertz time-domain spectra (THz-TDS) were investigated using the results of molecular dynamics (MD) simulations of Staphylococcal nuclease at two hydration states in the temperature range between 100 and 300 K. The temperature dependence of THz-TDS was found to differ significantly from that of the incoherent neutron scattering spectra (INSS) calculated from the same MD simulation results. We further examined contributions of the mutual and auto-correlations of the atomic fluctuations to THz-TDS and found that the negative value of the former contribution nearly canceled out the positive value of the latter, resulting in a monotonic increase of the reduced absorption cross section. Because of this cancellation, no distinct broad peak was observed in the absorption lineshape function of THz-TDS, whereas the protein boson peak was observed in INSS. The contribution of water molecules to THz-TDS was extremely large for the hydrated protein at temperatures above 200 K, in which large-amplitude motions of water were excited. The combination of THz-TDS, INSS and MD simulations has the potential to extract function-relevant protein dynamics occurring on the picosecond to nanosecond timescale. 2019 © The Biophysical Society of Japan.

The temperature dependence of terahertz time-domain spectra (THz-TDS) is shown to differ significantly from that of incoherent neutron scattering spectra (INSS) using molecular dynamics simulations of Staphylococcal nuclease at two hydration states in the temperature range between 100 and 300 K. The mutual correlation between distinct atomic fluctuations and water dynamics to THz-TDS, which are negligible in INSS, contributes to this difference. This work provides the first theoretical and computational analysis of the dynamical components that contribute to THz-TDS.

Conformational dynamics plays an important role in protein function. Collective motions in proteins slower than pico-seconds in timescale are of particular interest. Incoherent neutron scattering experiments are among the successful methods to study protein dynamics on the picosecond to nanosecond timescale. The protein boson peak is a broad peak found in the low-frequency region (16–32 cm−1, corresponding to 1–2 ps) of the incoherent neutron scattering spectra (INSS) [1-4]. Frequencies of the protein boson peak shift higher upon hydration, indicating that solvent molecules are associated with the origin of the peak [2,4]. As the temperature increases, the boson peak for a hydrated protein shifts to lower frequency and becomes buried in the quasi-elastic contribution [1]. Another temperature-dependent phenomenon with regard to protein dynamics is a dynamical transition characterized as an abrupt increase in atomic fluctuations at a temperature above ~200 K [1,5-7]. This increase has been interpreted as a result of a transition from harmonic to diffusive anharmonic motions of the protein. Such diffusive anharmonic motions are crucial to protein function [1,7]. Moreover, the transition is strongly suppressed in dry proteins, indicating that solvent molecules are implicated in the activation of the anharmonic motions [1,4].

INSS can be simulated using the results of molecular dynamics (MD) simulations [8-12]. Our previous studies have revealed that the protein boson peak is observed when protein motions distributed over the entire protein are trapped within narrow local energy minima with higher frequencies [11,12]. This trapping causes a frequency shift in the protein motions from the frequency range lower than the boson peak to higher, resulting in lowering of the spectral density in this range and enhancing the boson peak. Our simulations have also shown that most water molecules in the hydrated state (hydration ratio in weight, h=0.49) behave like bulk water above the transition temperature and act as a lubricant for protein dynamics, which excites diffusive anharmonic motion of the protein [13]. In contrast, water molecules in the dry state (h=0.09) tend to form more hydrogen bonds with the protein even at temperatures above ~200 K, which restricts the fluctuation of these water molecules to the fluctuation level of the protein. The dynamical transition in the hydrated state is significant in loop and terminal regions exposed to solvent, whereas these regions are restricted to protein-protein contacts in the dry state [13]. Consequently, restricted protein fluctuations in the dry state cause emergence of the protein boson peak even at 300 K [12].

As a result of technical advances in generating and detecting pulsed radiation in the terahertz region, terahertz time-domain spectroscopy has been widely conducted on bimolecular systems [14-21]. Terahertz time-domain spectra (THz-TDS) studies of proteins in solution have shown that water perturbs the protein dynamics via its hydrogen bond network by both direct short-range interactions in the first hydration shell and indirect long-range interactions up to several hydration shells [17-19]. The temperature and hydration dependence of biomolecular dynamics on the picosecond timescale has also been investigated with THz-TDS of biomolecules in the powder states [14-16]. In these studies, the relation between THz-TDS and the vibrational density of states (VDOS) is discussed. The auto-correlation of the atomic fluctuations is involved in VDOS and INSS, while the auto-correlation of the dipole moment of the sample reflects THz-TDS. In this work, we show that if charge distribution is expressed by partial charges assigned on atoms, both mutual correlation between distinct atomic fluctuations and auto-correlation of atomic fluctuations affect the fluctuation of the dipole moment, indicating that protein collective motions in the sample significantly reflect THz-TDS. The contribution of mutual correlation on THz-TDS has not been examined experimentally because its estimation with only the THz-TDS experiment is difficult. However, contribution analysis can be conducted theoretically based on the results of MD simulations.

Here we show that mutual correlation between atomic fluctuations, as well as auto-correlation, significantly contribute to THz-TDS, whereas the mutual correlation does not contribute to INSS. To calculate THz-TDS, we assume that the dipole moment of the simulation system can be roughly approximated by a classical treatment using the results of all-atom MD simulations of the hydrated and dry proteins at different temperatures ranging from 100 to 300 K. In order to compare the physical characteristics of THz-TDS with those of INSS, the same MD trajectories as in our previous study [12,13] were used. We also show that water dynamics significantly affect the dipole moment of the system in this paper, whereas the contribution of water dynamics is negligible in INSS of proteins in deuterated water [12]. Although quantum treatment is shown to be important for quantitative estimation of THz-TDS in water dynamics [22], we would expect that the significance of the contributions of the mutual correlations hold true even without quantum treatment.

Methods

Classical and atomistic description of terahertz time-domain spectrum

The absorption scattering cross section function, α(ω), measured by THz-TDS or infrared spectroscopy is defined [23] as:

where ħ, c, β, n(ω) are the reduced Plank constant, the velocity of light, the thermodynamic beta and the refractive index, respectively. The absorption lineshape function, I(ω), is defined as the Fourier transform of the time-correlation function of the total dipole moment of the sample, M(t), as:

where M̃ (ω) is defined as:

In THz-TDS, the wavelength of the terahertz field pulse is ~0.1 mm, which is much larger than the length of molecular orbitals. The small change in the electron-density distribution originating from quantum effects might not significantly affect the computation of the total dipole moment. However, we expect that the overall features of THz-TDS can be well reproduced by this treatment. Therefore, we computed the dipole moment in the simulation using the partial charge defined by the AMBER ff99 force field [24] and TIP3P water model [25] as:

where q and r(t) are the partial charge and positional vector for the j-th atom, respectively. Although electronic polarizability in molecules affects the picosecond dynamics of water [26], this contribution was not considered in this paper. We assume that the significance of the contributions of the mutual and auto-correlations of the atomic fluctuations and water to THz-TDS would not change, although the contribution of electronic polarizability might quantitatively change the spectra. Technical treatment of the periodic boundary condition for the system is discussed in Supplementary Text S1.

The absorption lineshape function of THz-TDS, I(ω), can be computed from the results of molecular dynamics simulations. It is noted that the energy level of kT at 100 K corresponds to a frequency of ~70 cm−1. Therefore, quantum effects concerning the ratio ħω/kT are not dominant in protein dynamics above 100 K in the terahertz-frequency regions (~30 cm−1). To eliminate the thermal factor, the reduced absorption cross section (RACS), αR(ω), is introduced [15,27] as:

which indicates that RACS is proportional to βω2I(ω). Hereafter, we focus on the behavior of βω2I(ω) for biomolecular samples.

Relation between vibrational density of states and THz-TDS

In this subsection, we examine the relation between the vibrational density of states (VDOS) and THz-TDS. The VDOS, g(ω), for the specimen is expressed using the power spectrum of the mass weighted velocity vector, ṙ(t), as shown in [28]:

where m is the mass of the j-th atom. The power spectrum of the velocity vector is rewritten using the power spectrum of the instantaneous deviation of the position vector from its average position:

where r(t)=〈r(0)〉+Δr(t). Here, we define χ(ω) as the Fourier transform of the time-correlation function of Δr(t), i.e., the power spectrum of Δr(t):

When the dynamics of the system is harmonic, χ(ω) is proportional to kT because of the equipartition of energy. Using Eqs. (6) and (8), g(ω) is expressed as:

Next, we formulate βω2I(ω) using χ(ω). When we treat the dipole moment classically as in Eq. (4), the time-correlation function of M(t) is rewritten as:

A Fourier transform of Eq. (10) yields the absorption lineshape function from THz-TDS,

where and,

A(ω) and C(ω) originate from the auto-correlation of the atomic fluctuation and the mutual correlation (or cross-correlation) of the atomic fluctuations, respectively. It should be noted that C(ω) can adopt a negative value, though A(ω) is always positive. Although Eq. (11) is derived based on the atomic partial charge model, a similar formulation is also possible by introducing local dipoles around atoms. From Eq. (11), we can derive the following equation:

If the absolute value of C(ω) is negligibly small, the density of states can be estimated from THz-TDS, although the weighting factors for each atom j differ between density of state (m) and THz-TDS (q2). However, C(ω) is not negligible as will be shown in the Results and Discussion section. The contribution of the mutual correlation between atomic fluctuations (or βω2C(ω)) to RACS (or βω2I(ω)) has not been evaluated because it is impossible to measure it directly from the experiment. In the Results and Discussion section, we examine the contribution of βω2C(ω) to βω2I(ω) using the results of molecular dynamics simulations of crystalline proteins.

Contributions of protein dynamics, solvent dynamics and their coupling to βω2I(ω)

As described in the Introduction, solvent dynamics play an important role in protein function. In order to clarify the contribution of the solvent to THz-TDS, here we decompose the total dipole moment of the sample, M(t), into two terms, the dipole moment of the protein molecule and that of the solvent:

The time-correlation function of M(t) can be decomposed into three terms:

Eq. (15) shows that the time-correlation function of M(t) is rewritten as the sum of the contributions of protein, solvent, and coupling between protein and solvent dynamics. It should also be noted that this decomposition is possible without assuming the atomic partial charge model as shown in the first line ofEq. (14). By performing a Fourier transform of Eq. (15), βω2I(ω) is decomposed into three parts as:

where Ipro(ω), Isol(ω), and Icoupling(ω) are the Fourier transforms of the first, second, and third terms in Eq. (15), respectively. In the Results and Discussion section, we examine the contributions of protein dynamics, solvent dynamics and their coupling to βω2I(ω).

Incoherent neutron scattering spectra

Here, we summarize the formulation of INSS, S(Q,ω), for comparison with THz-TDS (see Supplementary Text S2 for details) where Q and ω represent the momentum and energy transfers between incident neutron and sample, respectively. S(Q,ω) is defined as the weighted sum of the Fourier transform of the time-correlation functions:

where 〈exp(–iQ·Δr(0))exp(iQ·Δr(t))〉 indicates the time-correlation function, and b, and Δr(t) are the incoherent atomic scattering length and the instantaneous deviation of the position vector of the j-th atom from its average position at time t, respectively. This equation shows that INSS only consists of the contribution of the atomic auto-correlations, whereas mutual correlations reflect coherent neutron scattering, which is not considered in this work. On the sphere of |Q|=Q, g(ω), S(Q,ω) and χ(ω) have the following relation:

Computational procedure

Typically, powder samples are employed to investigate the temperature and hydration dependence of THz-TDS of biomolecules [14-16]. Here, we regarded the powder state as an ensemble of micro-crystals and employed molecular dynamics (MD) simulation trajectories of a crystalline Staphylococcal nuclease (SNase) at six temperatures ranging from 100 to 300 K [12,13] to calculate THz-TDS. In order to examine the solvent effect on protein dynamics, two types of biomolecular states, minimal- and full-hydration states (hereafter, MHS and FHS, respectively), were simulated. The hydration levels for MHS and FHS are h=0.09 and h=0.49 g D2O/g protein, respectively. The simulations were described in detail in our previous paper [13]. From the results of the simulations, dynamical transitions were observed at around 200 K in both MHS and FHS, which were more significant in FHS. The characteristics of water in FHS were similar to those for bulk water above the transition temperature, which enhances protein dynamics. The protein boson peak was observed in the frequency range 16–32 cm−1 at 100 K in both states. The peak frequency in FHS shifted to higher than that in MHS. At 300 K, the boson peak was buried in the quasi-elastic contributions originating from diffusive anharmonic motions in FHS, but was still visible in MHS.

The MD simulations at each temperature generated a 10-ns trajectory. Each 10-ns trajectory was divided into five 2-ns trajectories, and physical quantities such as I(ω) were calculated as the average of the results from the five 2-ns trajectories. THz-TDS presented here were smoothed by convolution with a resolution of a 0.323 cm−1, which corresponded to 0.1 ns, while INSS were smoothed by convoluting a 1.61 cm−1 (~200 μeV) resolution function, which corresponded to the LAM40 spectrometer [7].

Results and Discussion

Figure 1 shows the temperature dependence of the absorption lineshape function, I(ω), and INSS, S(Q,ω) at Q=2 Å−1, around the THz frequency (33 cm−1) in MHS and FHS (See Supplementary Figure S1 for higher resolution spectra). The temperature dependence of I(ω) was significantly different from that of S(Q,ω). No distinct broad peak was observed in I(ω) at any temperature in either FHS and MHS, while the protein boson peak was observed in INSS except for FHS at 260 and 300 K. The sum of χ(ω) weighted by the atomic scattering length was previously shown to have a peak in the frequency range of 16–32 cm−1 except for FHS at 260 and 300 K [12]. Considering this fact and Eq. (11), we expected that the contribution of C(ω) concealed other components in I(ω) having a peak in the frequency range 16–32 cm−1, which we will show later. As seen in Figure 2, the temperature-normalized spectra, βI(ω), in the frequency range ω>40 cm−1 in MHS were in good agreement at all temperatures, while those in FHS differed particularly at high temperatures. Water molecules in MHS are mainly hydrogen-bonded with the protein, resulting in restricted water movement as compared to the fluctuation of the protein [13]. The loop and terminal regions, which are exposed to solvent in FHS, are restricted in MHS by protein-protein contacts [13]. Because of these restrictions, dynamics in MHS can be considered as almost harmonic. In contrast, most water molecules in FHS behave like bulk water and act as a lubricant for protein dynamics at higher temperatures [12,13], which causes the difference in the temperature dependence of I(ω) between MHS and FHS.

Figure 1

Temperature dependence of the absorption lineshape function, I(ω), for (a) MHS and (b) FHS, and that of the incoherent neutron scattering spectra (INSS), S(Q,ω), for (c) MHS and (d) FHS at Q=2 Å−1. Error bars show standard deviations.

Figure 2

Temperature dependence of the absorption lineshape function normalized with temperature, βI(ω), for (a) MHS and (b) FHS.

Figure 3 shows βω2I(ω) and g(ω) at 100 and 300 K for MHS and FHS. In FHS, βω2I(ω) and g(ω) at 300 K were notably greater than those at 100 K. As shown in our previous paper [12], g(ω) for protein only in FHS at 300 K was similar to that at 100 K, indicating that the large difference in total g(ω) between 100 K and 300 K originated from g(ω) for water. The origin of the large temperature difference in βω2I(ω) is examined later. On the other hand, in MHS, βω2I(ω) and g(ω) at 300 K were only slightly greater than those at 100 K, due to restricted protein fluctuations in MHS [12,13]. In all the hydration and temperature conditions, βω2I(ω) increased monotonically in the frequency range ω<80 cm−1, whereas g(ω) became flat in the frequency range ω>40 cm−1. From the comparison between Eqs. (9) and (13), the contribution of C(ω) was concluded to be the origin of the difference.

Figure 3

(a) βω2I(ω) at 100 and 300 K for MHS (magenta) and FHS (cyan). (b) Density of states, g(ω), at 100 and 300 K for MHS (magenta) and FHS (cyan). Error bars show standard deviations.

Figure 4 shows contributions of the auto-correlation (the first term in Eq. (13)) and the mutual correlation (the second term in Eq. (13)) to total βω2I(ω) at 100 and 300 K in MHS and FHS. Similar to the ω dependence of g(ω), the auto-correlation term, βω2A(ω), becomes almost flat above 40 cm−1 in all cases. In contrast, while the ω dependence of the mutual correlation term, βω2C(ω), had a shape similar to that for the auto-correlation term, the sign was negative, almost a mirror image of the auto-correlation with respect to the abscissa. The negative value of βω2C(ω) almost canceled the positive value of βω2A(ω), which resulted in a positive value of the total function βω2I(ω) and a monotonic increase of RACS. In other words, the physical nature of RACS is significantly different from VDOS because C(ω) is not negligibly small. If the motion of the j-th atom is independent of that of the k-th atom, 〈Δr(0)·Δr(t)〉 in Eq. (10) becomes zero. Our results indicate that the collective motions of the system appear in C(ω). As theoretically shown above, THz-TDS includes information on both auto- and mutual correlations between atomic fluctuations, whereas INSS originates only from the auto-correlations.

Figure 4

Auto-correlation (broken line) and mutual correlation (dotted line) contributions to total βω2I(ω) (solid line) at (a) 100 and (b) 300 K for MHS (magenta) and FHS (cyan). This decomposition is defined in Eq. (13). Error bars show standard deviations.

Finally, we examined the contribution of the solvent to THz-TDS. As seen in Figure 1, the intensity of I(ω) for MHS differed entirely from that for FHS at all temperatures. On the other hand, INSS for MHS and FHS were in good agreement above 40 cm−1 at all temperatures [12], because the contribution of deuterated water to INSS was negligibly small in the higher frequency region. Figure 5 shows the contributions of protein dynamics, water dynamics and protein-water coupling to βω2I(ω) at 100 and 300 K for MHS and FHS. βω2Ipro(ω) in MHS and FHS was comparable at both temperatures. βω2Isol(ω) in MHS was lower than that in FHS at both temperatures, and that in FHS was comparable to βω2Ipro(ω) at 100 K and became almost double that of βω2Ipro(ω) at 300 K, which indicated that solvent dynamics in FHS dominantly contributed to βω2I(ω) at 300 K. In all cases, βω2Icoupling(ω) was small and negative. Because βω2Isol(ω) in MHS can nearly cancel out βω2Icoupling(ω) at both temperatures, βω2Ipro(ω) and βω2I(ω) look similar. We should note that the contribution of the solvent to βω2I(ω) and RACS should be extremely large for a fully hydrated protein. In order to extract information on protein dynamics from THz-TDS of the hydrated proteins, its combination with other methods such as INSS and MD simulation is essential.

Figure 5

Contributions of βω2Ipro(ω) (thin solid line), βω2Isol(ω) (broken line) and βω2Icoupling(ω) (dotted line) to βω2I(ω) (thick solid line) at (a) 100 and (b) 300 K for MHS (magenta) and FHS (cyan). This decomposition is defined in Eq. (16). Error bars show standard deviations.

Conclusion

We calculated THz-TDS of SNase in different hydration states using the results of MD simulations conducted at different temperatures. The temperature dependence of the absorption lineshape function was very different from that of INSS because THz-TDS significantly depended on both auto- and mutual correlations between atomic fluctuations. Since the decomposition of THz-TDS into auto- and mutual correlations is experimentally impossible, this work provides the first theoretical and computational analysis of motional components that contribute to THz-TDS. Since the contribution of the mutual correlation of the atomic fluctuations is not negligibly small, RACS is significantly different from VDOS, in contrast to INSS which depends only on the auto-correlation of the atomic fluctuations. Our results clearly demonstrated that the mutual correlation term has a negative intensity comparable to the auto-correlation term in absolute value. Therefore, the auto- and mutual correlations nearly cancel each other out, and the resulting positive residue can be observed as a monotonically increasing feature of RACS. Because of this cancellation, no distinct broad peak is observed in the absorption lineshape function, whereas the protein boson peak is observed in INSS except for FHS at 260 and 300 K. The behavior of water molecules in FHS is similar to that of bulk water above ~200 K [13], which enhances protein dynamics and results in the large contribution of water molecules to THz-TDS. Therefore, the contribution of water dynamics to THz-TDS at temperatures above 200 K is extremely large in FHS.

As shown above, THz-TDS can be decomposed into many contributions, such as the auto- and mutual correlations of the atomic fluctuations and solvent effects on the picosecond to nanosecond timescale, which are related to protein function. Although it is difficult to extract such information from THz-TDS only, the combination of THz-TDS, INSS and molecular simulations has the potential to allow the analysis of function-relevant protein dynamics.