Sub-THz acoustic excitation of protein motion

Proteins are an essential part of complex machinery that drives biological processes in living organisms. According to the sequence–structure–function paradigm raised by Anfinsen,¹ the function of a protein is governed by the unique three-dimensional structure adopted by the protein on the basis of its amino-acid sequence.^1,2 However, in order to understand how proteins work, it is also necessary to shed light on their behavior from a dynamic point of view.^3–6 Proteins at physiological temperatures exhibit a wide range of intrinsic motions,⁷ which span over different spatial and temporal scales and are governed by the (highly multidimensional) energy landscape.^4,8 Conformational transitions of proteins from one substate to another, separated by high energy barriers, take place at the microsecond timescale (or even slower) and are associated with the interconversion between open and closed conformations, enzyme catalysis, signal transduction, and protein–protein interaction. Looking at a more local scale, proteins exhibit faster dynamics on the picosecond-to-nanosecond timescale, which is related to local atomic fluctuations (e.g., side chain rotations) and collective fluctuations of small groups of atoms (e.g., loop motions).^4,9 In addition, these local protein motions have been shown to be involved in the initiation and modulation of slower dynamic motions, such as those associated with protein conformational changes.^10–12 

Collective motions associated with protein function can (often) be adequately described by (only) a few low-frequency normal modes.^13–17 To computationally access functionally relevant collective motions of biomolecules, normal mode analysis (NMA) is widely used. NMA is based on the harmonic approximation of potential energy^18,19 and assumes that large-scale motions relevant to protein function are described by normal modes with the lowest frequencies.²⁰ However, proteins at physiological temperatures exhibit anharmonicity in their motion. To this end, the principal component analysis (PCA), also implemented as the quasi-harmonic analysis (QHA),²¹ is carried out on a molecular dynamics trajectory covering full range of motions.^22–26 Employing PCA, collective modes of internal motion are obtained as eigenvectors of the covariance matrix of atomic fluctuations.²⁷ The eigenvalues of the covariance matrix represent the variance along each of the corresponding collective modes. In PCA, only a limited number of atoms, i.e., C^α atoms, or selected groups of atoms can be used for analysis.²² 

Both computational methods (i.e., NMA and PCA) serve as an important tool for the interpretation of spectra obtained, for example, by terahertz (THz) time-domain spectroscopy,²⁸ Raman spectroscopy,^29,30 and neutron scattering.^31–35 As these experimental techniques allow access to the picosecond-to-nanosecond timescale, they are particularly suitable for detecting, identifying, and examining functionally relevant motions in biological molecules,^{30,33,36–46} as the underlying modes of collective dynamics vibrate in the THz region, which spans from 0.1 to 10 THz (or expressed in wavenumbers from 3.33 to 333 cm⁻¹).^47–49 Moreover, THz radiation can be used for activating the process without causing denaturation or aggregation of the protein, as shown for the actin solution during polymerization.⁵⁰ In addition, THz radiation with several orders of magnitude higher peak power was displayed to cause changes in the morphology of actin filaments that were attributed to the THz-light-induced acoustic waves,^51,52 whose generation in water was also confirmed.⁵³ 

THz excitations can thus be used as a tool to modify the structure of the biomolecule or to manipulate cellular functions by directly affecting protein dynamics.^50,51,54 As illustrated, these are low-frequency internal vibrations that are key to understanding and predicting the collective dynamic behavior^16,17,55 and should be targeted by external perturbations. To investigate protein dynamics under mechanical stress induced by acoustic waves, we select ubiquitin as the benchmark protein and model it using the coarse-grained (CG) Martini 3 model.⁵⁶ We conduct nonequilibrium simulations using the open-boundary molecular dynamics (OBMD) method^57–60 that enables imposition of external boundary conditions (e.g., acoustic wave) by an additional external force.

To inspect the protein response to induced mechanical stress, the CG protein (shown by the gray colored beads in Fig. 1) is placed at the center of the simulation box and surrounded by CG water (represented by the blue colored beads in Fig. 1). Employing OBMD, on one side of the simulation domain (i.e., the buffer region indicated by the orange rectangle in Fig. 1), sub-THz acoustic waves are generated, while on the other side (i.e., the buffer region depicted by the yellow rectangle in Fig. 1), only the constant normal load is applied.

Next, we present details on the imposition of external boundary conditions (i.e., the constant normal load and a mechanical pressure wave) using OBMD and on the computation of collective modes of the protein employing PCA.

To simulate acoustic waves, we resort to OBMD,^57–60 which has already been applied in particle-based ultrasound simulations.^61,62 OBMD opens the molecular system, which is able to exchange mass and momentum with its surroundings.⁶³ It allows the imposition of external boundary conditions without altering Newton’s equations of motion in the bulk. The simulation box is opened in one direction, i.e., in the x-direction, while in the remaining directions, i.e., y- and z-directions, periodic boundary conditions are implemented (see Fig. 1). The external boundary condition of the constant normal load p is imposed on the region of interest (ROI) from both sides of the simulation domain (in Fig. 1 depicted by the black arrows), whereas the acoustic wave is introduced from only one side (in Fig. 1 shown by the orange rectangle) as an additional oscillatory pressure contribution to the normal load p, expressed as p + Δp sin(2πνt), where Δp, ν, and t are the pressure amplitude, frequency, and time, respectively.

The dissipative particle dynamics (DPD) thermostat is applied, as its equations conserve linear momentum and correctly reproduce the hydrodynamic behavior.^64–66 

Performing out-of-equilibrium OBMD simulations of the propagating acoustic wave through the protein system in water allows us to obtain an insight into variations of the protein’s internal dynamics.

Employing OBMD allows us to impose acoustic waves with frequencies matching those of the protein’s collective modes. The simulation box of dimensions 12.3 × 6.2 × 6.2 nm³ is divided into two buffer regions, surrounding the ROI. The protein is placed at the center of the ROI and enclosed by water molecules (see Fig. 1). The CG Martini 3 model of ubiquitin is obtained after treating the initial atomistic structure (PDB entry 1UBQ⁷¹) with Martinize2, available at https://github.com/cmbi/dssp/releases/tag/2.3.0 (for illustration, see Fig. 2). The information on the classification of the secondary structure of the protein backbone from the structure is provided from the DSSP database available at https://github.com/cmbi/dssp/releases/tag/2.3.0. The protein’s secondary and tertiary structures are maintained by additional harmonic bonds added between the backbone particles. The corresponding elastic bond constant is set to 550 kJ mol⁻¹ nm⁻², where the lower and upper elastic bond cutoffs are fixed at 0.5 and 0.9 nm, respectively. The geometry of the protein molecule is constrained using RATTLE.⁷² The equations of motion are integrated using the velocity Verlet algorithm⁷³ with an integration time step δt = 0.02 ps at a temperature of 300 K. All simulations are performed using the Martini v3.0.0 force field⁵⁶ and ESPResSo++ simulation package.⁷⁴ Non-bonded interactions are described using the Lennard-Jones 12-6 potential energy function with a cutoff value of 1.1 nm, while the long-range interactions of charged groups are treated by the Coulomb energy function with a relative permeability ɛ_r = 15 and a cutoff distance of 1.1 nm. To sample a wide range of protein’s conformational space in and out-of-equilibrium, several parallel simulations of length 500 ns are conducted, for which only the initial position of the protein is varied. The dynamic behavior of the biomolecule is analyzed over the entire trajectory and averaged over all trajectories sampled in equilibrium simulations or at a given frequency of the imposed acoustic wave. In addition, to prevent the protein from diffusing through the open ends of the simulation box or over the edges where the periodic boundary conditions are implemented, an additional harmonic constraint (k = 3.5 kg s⁻²) is added to hold the protein CoM close to the center of the ROI. As we are only interested in internal vibrational dynamics of ubiquitin, implementation of the harmonic constraint is not expected to affect the results.

We investigate the internal dynamics of the protein by subjecting it to acoustic waves with frequencies corresponding to those of its vibrational collective modes. To obtain collective modes, we perform OBMD simulations in equilibrium (where from both sides only a constant normal load is imposed on the system) and apply PCA. Given that the number of PCA modes scales with the number of beads taken into account and for the fact that we are interested in as few as possible representative modes, we perform PCA on further reduced coordinates representing CoM of six selected subunits. Accordingly, we obtain 18 modes. Figure 3 shows the lowest (left) and the highest (right) PCA modes of ubiquitin in a six-bead representation. In addition, localized motions in a few segments can be observed in Fig. 3 (right). As the low-frequency collective motions of the protein are usually related to its function,^13–17 we focus on excitation of only the lowest collective modes of ubiquitin. To this end, we excite acoustic waves with frequencies up to 0.080 THz (which covers the frequency range of the four lowest vibrational modes). Parameters p and Δp are fixed at 1 and 100 bars, respectively.

The average acoustic energy density $(w̄)$ introduced into the system is estimated using $W̄=j/c$⁠, where c and j stand for the speed of sound and energy density of the flux, respectively. The latter is expressed as j = 1/2(Δp)²/ρc. Using the expression for $w̄$ for two different pressure amplitudes, i.e., 100 and 200 bars, the value of $w̄$ is estimated to be ∼17 000 and 70 500 kg m⁻¹ s⁻², respectively. Multiplying $w̄$ by the volume of interest, the kinetic energies at Δp = 100 and 200 bars are $∼10−22$ and $∼10−21$ kg m² s⁻², respectively. The acoustic energy pumped into the system turns out to be relatively small compared to the thermal energy stored in the same volume, which is $∼10−19$ kg m² s⁻². Nevertheless, the acoustic wave propagating in one direction clearly introduces anisotropy of particle velocity distribution by enhancing the longitudinal component, as shown in Fig. 4.

Results for different frequencies of acoustic waves are shown in the top and middle panels of Figs. 5 and 6. Inspecting the $Grα(ν)$ spectrum for different beads representing protein subunits, we observe a pronounced peak at the frequency corresponding to the one used to excite the biomolecule (see the top panels of Figs. 5 and 6, where the black dotted line indicates the acoustic wave frequency to which the protein is subjected). The top panel of Fig. 5 depicts the results obtained after subjecting the protein to the acoustic wave with a frequency of 0.048 THz. As can be seen, the subunits respond distinctly at the frequency of 0.048 THz. Observed response means that exposure of the protein to acoustic wave causes the subunits to make larger deviations from their average positions. However, not all subunits respond to the same extent [for instance, compare $Grα(ν)$ of different subunits in the top panel of Fig. 5]. This observation might be explained by arguing that the protein is expected to respond strongly to the acoustic wave excitation of a certain frequency if the direction of the vibrational mode is aligned with the propagating direction of the acoustic wave. Conversely, if the collective mode is oriented perpendicular to the direction of the acoustic wave propagation, no distinct response is expected. Nevertheless, vibrational modes whose displacements are in arbitrary directions contribute to response, but their contribution depends on the deviation from the direction of the wave propagation. Similarly, the highest peak (or one of the highest peaks) in $GRg(ν)$ is observed at the excitation frequency (see the middle panels of Figs. 5 and 6). This further substantiates the sensitivity of the protein structure to the applied mechanical stress. In addition, the indicated changes in the compactness of the protein structure could also be due to the breathing-like dynamic character of the collective mode. For comparison, we add the results for the case in which the frequency of acoustic wave differs significantly from any of the mode frequencies of the protein in Fig. 7, demonstrating only a modest response in $GRg/rα(ν)$⁠, which should serve as a “good negative control.”

The computed $R$ values, defined in Eq. (5), are shown as the function of the collective mode frequency in the bottom panels of Figs. 5 and 6. In general, by comparing the $R$ values at the frequencies used to perturb the protein (depicted by the purple crosses) with the $R$ values of the protein in equilibrium (shown by the orange crosses), we find that the former increase. This suggests that the acoustic wave perturbation causes dynamic motion of the collective mode to become more breathing-like. As can also be observed, frequencies at which the breathing-like character of collective modes is present are correlated with frequencies at which peaks in $GRg(ν)$ occur. Therefore, part of the change in the compactness of the protein structure could be attributed to the breathing-like dynamic motion. Bead representations of vibrational modes with the frequency corresponding to the frequency used for excitation are shown in the insets of the bottom panels of Figs. 5 and 6. As can be seen, displacement vectors in modes matching the excitation frequencies of 0.048 and 0.062 THz have a pronounced radial character, which is in agreement with higher $R$ values of collective modes depicted in Figs. 5 and 6 (see the purple crosses next to the black dotted lines).

In experimental setups applying electromagnetic THz waves, the most common parameters varied are the frequency of radiation, intensity, and duration of exposure to the terahertz radiation.^50,78 For the implemented acoustic waves, we can vary all these parameters: the frequency (i.e., input parameter ν), the amplitude of the imposed wave (i.e., input parameter Δp), and the time of exposure (i.e., simulation length). By varying the frequency, we show that the protein is highly susceptible to induced mechanical stress, exhibiting larger deviations from average positions and increased breathing-like dynamic motion of modes. As we anticipate that the higher amplitude (i.e., Δp) would cause even larger displacements, we additionally perform OBMD simulations in the similar manner as before, except that we impose only acoustic waves with a frequency of 0.062 THz and set Δp to 200 bars.

Repeating the analysis, we obtain the $Grα(ν)$ spectra and compare them with those determined at the amplitude of 100 bars. The results are shown in Fig. 8. Taking these into consideration, even more pronounced peaks are observed at frequencies matching the excitation ones (indicated using the black dotted lines). This suggests that the beads describing protein subunits increasingly deviate from their average positions when subjected to mechanical pressure waves of higher amplitudes. Note that not all objects respond in the same manner due to the deviation of displacement vectors from the direction in which the imposed acoustic wave is propagated. Nevertheless, at higher amplitudes of the pressure wave, the protein is even more susceptible to mechanical stress, as supported by a higher peak in $GRg(ν)$ (see Fig. 9).

In the present study, we show that the internal motion of a protein can be manipulated by mechanical pressure waves of different frequencies. To this end, a harmonically restrained coarse-grained protein model is immersed in water and subjected to acoustic waves whose frequencies are chosen around the frequencies of the protein’s quasi-harmonic collective modes. Such a protein model is biased by its elastic character, which affects the resonant behavior of the protein compared to a real system dominated by highly anharmonic vibrations at low frequencies. Nevertheless, it allows straightforward application of the PCA to determine the collective motions of the protein. By analyzing the autocorrelation functions of the generalized velocities derived from displacements of beads that describe the protein subunits and from the radius of gyration, we show that the protein is susceptible to mechanical stress induced by acoustic waves. The response of the protein is more pronounced at frequencies corresponding to the acoustic wave. The assumption that the pressure wave with a larger amplitude would cause even larger deviations and enhanced changes in the protein compactness is confirmed and demonstrated by our simulations. Moreover, the dynamic response of the protein also depends on the correspondence between the direction of the corresponding collective mode and the propagation direction of the acoustic wave at a given frequency.

The authors acknowledge the financial support under Grant Nos. P1-0002 and J1-3027 from the Slovenian Research and Innovation Agency. The financial support through ERC Advanced Grant (Grant No. 885155) from the European Research Council is also gratefully acknowledged.

The authors have no conflicts to disclose.

Petra Papež: Formal analysis (lead); Methodology (equal); Validation (lead); Visualization (lead); Writing – original draft (lead). Franci Merzel: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Validation (equal); Writing – original draft (supporting); Writing – review & editing (equal). Matej Praprotnik: Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Methodology (equal); Supervision (equal); Validation (equal); Writing – original draft (supporting); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.