Molecular Dynamics Simulations of Shockwave Affected STMV Virus to Measure the Frequencies of the Oscillatory Response

: Acoustic shockwaves are of interest as a possible means of the selective inactivation of viruses. It has been proposed that such inactivation may be enhanced by driving the virus particles at frequencies matching the characteristic frequency corresponding to acoustic modes of the viral structures, setting up a resonant response. Characteristic frequencies of viruses have been previously studied through opto-mechanical techniques. In contrast to optical excitation, shockwaves may be able to probe acoustic modes without the limitation of optical selection rules. This work explores molecular dynamics simulations of shockwaves interacting with a single STMV virus structure, in full atomistic detail, in order to measure the frequency of the response of the overall structure. Shockwaves of varying energy were set up in a water box containing the STMV structure by assigning water molecules at the edge of the box with an elevated velocity inward—in the direction of the virus. It was found that the structure compressed and stretched in a periodic oscillation of frequency 65 ± 6.5 GHz. This measured frequency did not show strong dependency on the energy of the shockwave perturbing the structure, suggesting the frequency is a characteristic of the structure. The measured frequency is also consistent with values predicted from elastic theory. Additionally, it was found that subjecting the virus to repeated shockwaves led to further deformation of the structure and the magnitude of the overall deformation could be altered by varying the time delay between repeated shockwave pulses.


Introduction
Laser-generated acoustic shockwaves or "transients" have been previously used to break apart viruses in order to study their structural bonding forces [1,2]. Elsewhere, selective inactivation of viruses by pulsed laser irradiation has been observed [3], with impulsive stimulated Raman scattering as the suggested mechanism. Recent work proposes acoustic shockwaves may play a role in pulsed laser inactivation of viruses with plasmonic enhancement of the laser absorption, from the addition of metal nanoparticles [4]. Other works have proposed the inactivation of viruses via resonant absorption of ultrasonic [5] or microwave energy [6,7]. Resonant absorption would involve driving the virus at its characteristic frequencies, such that the vibrational response matches the frequency of driving energy, leading to increasing amplitude of oscillations and eventual failure of the structure. For this purpose it is useful to have an estimate of the characteristic frequencies of a target virus.
Early analytical models used elastic theory to predict vibrational frequencies of viruses [8][9][10][11][12]. More recent computational models have been developed to study such frequencies using elastic network normal mode analysis [13][14][15]. Another study [16] models the subunits of the structure at the atomic level, applying constraints due to the inherent symmetries of the icosahedral structure, to predict the vibrational frequencies of a virus Acoustics 2022, 4 269 in the absence of explicit water. In contrast to these, the present work treats the entire structure of a fully solvated virus in full atomistic detail, by manipulating the surrounding water molecules to perturb the system with an acoustic shock wave, and measures the resulting oscillations in the structure.
Opto-mechanical studies involving Brillouin and Raman scattering techniques have previously been conducted to measure the optically active modes of viruses [17][18][19][20][21][22][23]. Optically active modes in these cases are those which obey selection rules for light coupling with the acoustic modes of the material [24]. In contrast to this, shock wave perturbation is a direct mechanical displacement of molecules via Coulombic and other near-field forces between the water molecules and those of the virus structure. Thus, the modes available for probing by shockwave perturbation are not restricted by the optical selection rules.
The Satellite Tobacco Mosaic Virus (STMV) has been studied extensively, as a smallsized virus~17 nm in diameter. Its structure is well-known from crystallography studiesits geometry is icosahedral as are many other viruses, and thus it is a representative case study for such viruses. The capsid is comprised of 60 identical copies of a single protein with molecular weight of 14,500. The structure has been studied via molecular dynamics (MD) simulations [25,26]. In fact, because it is a relatively large structure for such MD simulations, comprised of millions of molecules, it is often used as a benchmark to test the performance of computing systems [27,28]. Shockwaves and their effects on biological systems have been simulated in previous work [29][30][31][32][33], using different methods to set up the shockwaves in the simulation. For example Surdutovich [29][30][31] initiates the shock waves using a "hot cylinder" of atoms in a specified radius with average velocities higher than those of the surrounding media, which in turn expand and interact with neighboring atoms and thus propagate outward. These simulated shockwaves then passed through DNA structures, in order to study whether such waves would damage DNA. Since molecular dynamics simulations treat covalent bonds as harmonic oscillators, and do not take into account such bonds breaking, Surdutovich considered any bonds exhibiting vibrational energies above a certain threshold as being broken during the simulation. Other methods for setting up shockwaves such as that of [32,33] also involve selecting a region of water molecules surrounding the biological system of interest, and giving these water molecules a velocity in a particular direction. This is also the approach used in the present work, described in more detail below.

Materials and Methods
The structure of the Satellite Tobacco Mosaic Virus (STMV), solvated in a cube of water of cell dimensions 22 nm, was prepared based on files available from the website of the Theoretical and Computational Biophysics group, NIH Center for Macromolecular Modeling and Bioinformatics, at the Beckman Institute, University of Illinois at Urbana-Champaign. The Protein Structure File (PSF) contains the structure of the virus capsid and RNA, along with water molecules of the solution. The Protein Data Bank (PDB) file contains information about the coordinates of the atoms from the PSF file, for the start of the simulation. Using Nanoscale Molecular Dynamics (NAMD) simulation software [34], the system was then allowed to equilibrate under a constant pressure of 1 atm, and a constant temperature of 298 K for 15 ps. The output from this step was used as the input for the subsequent step, in which a shockwave was introduced in the following way. Along the edge of the water box, perpendicular to the Z-axis, atoms falling within a certain distance of the edge of the box (arbitrarily chosen to be~3 nm) had their velocities altered. The X and Y velocities were left unchanged, but the Z velocity was changed to be a fixed value for these specific atoms. This method for producing simulated shockwaves is similar to the one employed in [33]. The forces on all atoms in the system are then computed via the force field resulting from interaction with neighboring atoms within the cutoff distance-in this case 12 Å. For a given time-step, the velocity and new position of each atom are then numerically computed from the velocity and position of the prior time-step and the overall force on the atom from the current time-step. This process is repeated for each time-step, of duration 1 fs in this case, and a trajectory for all atoms is computed. Since the selected water molecules were given an elevated velocity inward, towards the rest of the water box containing the virus, they move closer to other water molecules, resulting in forces acting on these neighboring atoms, so that they too gain an elevated velocity. In this manner, in subsequent time-steps the shockwave propagates through water molecules until it eventually encounters the molecules of the virus structure.
Perhaps a more realistic model of the initial velocities which would occur in an actual physical shockwave would instead have a distribution of velocities centered around this chosen value, but this method of using a uniform value was simply chosen for convenience, since it made identifying the affected atoms simple for debugging purposes and to easily vary the parameter of energy input. To justify the use of this simplistic approach to creating shockwaves in the water box, it should be noted that the water molecules which have been subjected to this artificial instantaneous Z-direction velocity, will quickly (by the next processing step) encounter other water molecules outside the affected region. As they interact with these unaltered water molecules, there will be an exchange of momentum through the force field, and the random distribution of the velocities of those unaffected water molecules will once again smooth out and randomize the velocities of those molecules which underwent the artificial velocity change. Further processing steps will lead to more interactions with randomly moving water molecules and the strong pulse of velocities in the Z-direction will disperse into other directions, but the overall momentum will be biased in the Z-direction and the pulse will propagate through the rest of the system. Thus, provided the molecules of the biological system (the virus in this case) are far enough away from the edge of the box, by the time this pulse reaches the molecules of interest, it should resemble a planar acoustic wave transient, such as that which could be produced by an ultrashort laser pulse. Indeed, measuring the distribution of Z-direction velocities at various times shows that after a few time-steps, the velocities of the water molecules fall in a skewed Gaussian distribution. In addition to the aforementioned convenience of being able to quickly pick out the atoms which were subject to the artificial velocity change (since they are all assigned the same Z-direction value for the initial step), there is the added benefit of being able to inject a very high initial Z-directed energy pulse, while avoiding the velocity limits of the NAMD simulation system. If more energy is still needed, this can be achieved by simply choosing a larger slice of atoms along the boundary. Another approach would be to instead simply reverse the Z-direction velocities of these same atoms, since they should already have a realistic distribution of velocities after the equilibrium step; this approach is known as a "momentum mirror." [35] It has the advantage of not altering the overall kinetic energy of the system-however, this would not be representative of a real laser-generated shock wave as in this case kinetic energy is added to the system via absorption of the laser energy. In another, similar approach, a slowly moving boundary could be applied, and atoms reaching this boundary would have their velocities reversed-like a "hard edge." This is the "piston" approach of [32], and it has the advantage of avoiding clashes of atoms which may arise when the sudden change in velocities causes atoms to come too close together. Such clashes can cause the simulation to fail, since the calculated forces become so large that they are beyond the ability for the computer to store and compute. However, if a sufficiently small time-step is chosen, clashes can be avoided since the atoms have time to react to the forces generated when they approach one another, thus altering their direction and avoiding unmanageably large force values. For this reason, for some simulations performed in this work-those involving the highest energies-a very short timestep of 1 fs was used for parts of the simulation.

Results
Using NAMD, the STMV virus structure in a water box was subjected to a single pulse of a simulated acoustic shockwave using the method described above, at various pulse energies, and the response of the system studied. Using Visual Molecular Dynamics (VMD) software [36], conformational changes were studied. The Root Mean Squared Displacement (RMSD) for the entire protein backbone of the virus capsid structure was calculated using the built-in plug-in of this software. A higher value indicates a larger overall conformational change to the entire structure throughout the trajectory. Shown in Figure 1a are the RMSD curves for simulated shockwaves at various energies. Shown in Figure 1b are RMSD curves for repeated shockwaves compared to a single pulse. For each pulse, the initial velocity of water molecules is 39.999 Å/ps-a relatively gentle pulse. By varying the time delay between pulses, the amplitude of the overall deformation of the structure is seen to vary. Shown in Figure 2 are screenshots from the VMD visualization of the virus subject to a shockwave of initial water speed 99.999 Å/ps (the maximum allowable velocity of the software-chosen to show the deformation most clearly), at time points corresponding to extrema in the RMSD curve. Figure 2a shows an illustration of a planar cross-section of a sphere with the L3 mode at maximum deformation for comparison. Figure 2b shows the virus before encountering the shockwave-at time 0 and 0 RMSD, with the spherical crosssection from Figure 2a overlaid for comparison. Figure 2c shows the virus at maximum deformation, with the cross-section of L3 mode at max deformation overlaid.

Results
Using NAMD, the STMV virus structure in a water box was subjected to a single pulse of a simulated acoustic shockwave using the method described above, at various pulse energies, and the response of the system studied. Using Visual Molecular Dynamics (VMD) software [36], conformational changes were studied. The Root Mean Squared Displacement (RMSD) for the entire protein backbone of the virus capsid structure was calculated using the built-in plug-in of this software. A higher value indicates a larger overall conformational change to the entire structure throughout the trajectory. Shown in Figure 1a are the RMSD curves for simulated shockwaves at various energies. Shown in Figure 1b are RMSD curves for repeated shockwaves compared to a single pulse. For each pulse, the initial velocity of water molecules is 39.999 Å /ps-a relatively gentle pulse. By varying the time delay between pulses, the amplitude of the overall deformation of the structure is seen to vary. Shown in Figure 2 are screenshots from the VMD visualization of the virus subject to a shockwave of initial water speed 99.999 Å /ps (the maximum allowable velocity of the software-chosen to show the deformation most clearly), at time points corresponding to extrema in the RMSD curve. Figure 2a shows an illustration of a planar cross-section of a sphere with the L3 mode at maximum deformation for comparison. Figure 2b shows the virus before encountering the shockwave-at time 0 and 0 RMSD, with the spherical cross-section from Figure 2a overlaid for comparison. Figure 2c shows the virus at maximum deformation, with the cross-section of L3 mode at max deformation overlaid.   Response of virus structure to shockwave pulses (a) RMSD curves exhibiting the oscillatory response of the virus to single shockwave pulses-each curve represents the time evolution of the deformation of the virus structure for a given initial velocity value assigned to water molecules in the water slice to initialize the shockwave. (b) RMSD curves showing response of the virus to repeated shockwaves pulses compared to a single shockwave pulse-a varying delay between pulses leads to variation in the overall amplitude of deformation of the structure. As seen in Figures 1 and 2, as well as video SV1, when the simulated shockwave front begins to hit the virus at time 200 fs, the structure begins to deform and the magnitude of this deformation increases as the wave propagates through the structure and the surrounding water medium, peaking at time 1.4 ps, when the structure is at maximum compression. As the wave reaches the end of the structure, it is stretched out, re- As seen in Figures 1 and 2, as well as video SV1, when the simulated shockwave front begins to hit the virus at time 200 fs, the structure begins to deform and the magnitude of this deformation increases as the wave propagates through the structure and the surrounding water medium, peaking at time 1.4 ps, when the structure is at maximum compression. As the wave reaches the end of the structure, it is stretched out, resulting in a secondary peak in the RMSD at time (9 ps). Finally, there is a recoiling of the overall structure as the structure is compressed from the recoil. This is reminiscent of a damped harmonic oscillator type motion, better described in the model of [8], with spherical Bessel functions. The time in between the first and second peaks ranges from 7.0 ps to 8.4 ps, for pulses with average initial velocity ranging from 99.999 Å/ps and 29.999 Å/ps. Since the primary and secondary peaks represent the times of maximal compression and stretching, respectively, the overall time for the complete oscillatory motion is double the time between first and second peak. These times correspond to 59.5 to 71.4 GHz. As seen in Table 1, there is little correlation between this measured frequency and the initial velocities of the water slice molecules used to set up the simulated shockwave.

Discussion
The simulated structural response of STMV exposed to a shockwave exhibited oscillatory behavior, the periodicity of which did not show high dependency on the initial energy of the shockwave. This suggests that the frequency may be an intrinsic property of the structure. This frequency was determined to be 65 ± 6.5 GHz.
For a comparison with values predicted from elastic theory, the icosahedral virus can be approximated as a sphere, in order to employ Lamb's theory [9,37,38] for estimating vibrational frequencies. Taking values of longitudinal wave speed and Poisson's ratio to be 1800 m/s and 0.30, respectively (corresponding to a transverse wave speed of 900 m/s), consistent with values used in prior works [6,11,18,22], and with a radius of 17 nm, the L1 and L3 modes are predicted to be~60.6 GHz and 66.6 GHz, respectively, which are in good agreement with the measured value. It should not be surprising that these oddnumbered modes would be activated by a pulse energy which asymmetrically compresses the structure from one side, leaving the other side initially unperturbed. Indeed from visual inspection of the oscillatory motion using the VMD software, the shape of the compressed state appears to most resemble that of the L3 mode, as seen in Figure 2c. Thus, it seems reasonable that the predicted L3 mode frequency (66.6 GHz) should be a good approximation of the measured values (65 ± 6.5 GHz.) Note, however that the error in this measurement is quite significant (10%), so the L1 mode could be considered an equally good numeric fit. Table A1 shows the first few spheroidal modes, for comparison. Further, the real motion of the virus oscillations is more complex than that of a simple sphere, as the composition of the structure is very complex and the shockwave perturbation can activate more than one mode. Moreover, the parameters (Poisson's ratio, wave speeds of the material) used in the Lamb's theory approximation are not entirely well-established in the literature, so these should be considered estimations as well. Given all the sources of error, it is reassuring that the predicted values match the measured values so well, but this should be taken with some care. For example if a Poisson's ratio of 0.40 were instead used, this would result in a predicted value of 53.6 GHz for the frequency of the L3 mode. It is also worth noting that, although it is reassuring that elastic theory gives similar values for vibrational frequencies, the full atomistic molecular dynamics model gives additional information by taking into account atomic interactions of the biological components with one another and with the water molecules. These could be approximated by introducing a damping term to the elastic theory model, but molecular dynamics simulation gives a more comprehensive understanding of the dynamics of the system. For example, it appears that the vibrations quickly dissipate and begin to level off after the first couple of peaks in the RMSD curve, as seen in Figure 1a.
It was further observed that an enhancement of the amplitude of deformation could be achieved by subjecting the virus to repeated shock-wave pulses. By choosing a time delay between pulses, such that the deformation of the structure is at a maximum as a subsequent shockwave pulse arrives at the structure, the overall deformation increased, as seen in Figure 1b. This supports the plausibility of achieving a resonant response by driving the virus structure at the proper frequency, via acoustic waves such as shockwaves.
Future research may explore lower energy perturbations-for example, using water speeds which do not exceed the wave speed of water, in order to see if the oscillatory motions are still observable and whether the frequencies of such motions would be consistent with those measured using these super-sonic shockwave pulses. If oscillatory frequencies are ubiquitous, irrespective of excitation energy, it may suggest the possibility of driving the virus structure at such frequencies, in order to set up a resonant response.
Another interesting direction for future study would be to examine the possibility of virus inactivation, by investigating bond-breaking events which may occur during shockwave excitation. The results described here treated bonds as a simple harmonic potential, which is unrealistic as bond energies simply increase unbounded with bondlength. More interesting would be to simulate shockwave interaction with viruses using more advanced MD simulation techniques such as ab initio or reactive force field methods in order to study bond-breaking and re-formation. This would give a more realistic view as to the plausibility of viral inactivation via resonant excitation of the structure.
In spite of the relative simplicity of this model, some useful observations can be made. To wit: frequencies of oscillation due to shockwave perturbation seem to correspond to the odd-numbered spheroidal modes (L1 and L3), and the frequencies are around 65 GHz, irrespective of shockwave energy. With these observations, it may be possible to tailor a scheme whereby selective inactivation of viruses could be enhanced by driving the system at such a characteristic frequency-either via shockwaves or otherwise.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Shockwave simulations were performed with Charm22 force field [39], rigid water bonds and periodic boundary conditions. For further simulation details, the NAMD simulation instruction file is provided in the supplementary materials section.