Shear shock waves are observed in the brain

TL;DR

This study uses advanced ultrasound imaging to directly observe shear shock waves in the brain, revealing their spontaneous formation and potential role in traumatic injuries.

Contribution

First direct observation of shear shock waves in the brain using a novel ultrasound technique, highlighting their significance in traumatic brain injury mechanisms.

Findings

Shear shock waves are spontaneously generated in the brain during impacts.

Acceleration at the shock front is amplified up to 8.5 times.

Shear shock waves may be a fundamental mechanism for traumatic injuries.

Abstract

The internal deformation of the brain is far more complex than the rigid motion of the skull. An ultrasound imaging technique that we have developed has a combination of penetration, frame-rate, and motion detection accuracy required to directly observe, for the first time, the formation and evolution of shear shock waves in the brain. Experiments at low impacts on the traumatic brain injury scale demonstrate that they are spontaneously generated and propagate within the porcine brain. Compared to the initially smooth impact, the acceleration at the shock front is amplified up to a factor of 8.5. This highly localized increase in acceleration suggests that shear shock waves are a fundamental mechanism for traumatic injuries in soft tissue.

Equations

$$\frac{\partial v}{\partial x}-\frac{\beta }{c_t^3}{v}^2\frac{\partial v}{\partial \tau }=\delta \frac{{\partial }^{2}v}{\partial {\tau }^2},$$

Full text

Traumatic brain injuries (TBI’s) are a major source of death and disability worldwide. Falls and motor vehicle related accidents are the largest contributors. Current biomechanical predictive criteria for TBI are based on measurements of skull motion such as linear and rotational acceleration Rimel et al. (1981); Rowson and Duma (2013). Although the relationship between skull motion and injury has been extensively tested Rowson and Duma (2013); Broglio et al. (2009); Crisco et al. (2011); Daniel et al. (2012); Laksari et al. (2015); Wright and Ramesh (2012), mechanisms relating the two have not been conclusively established Guskiewicz et al. (2007), due to the complexity of the deformation of the brain Bayly et al. (2005); Moss et al. (2009); Goriely et al. (2016) which behaves as a nonlinear viscoelastic medium. In situ measurements of the rapid transient motion of the whole brain during a traumatic event may establish a more accurate biomechanical description of injury.

Shear vibration experiments on small brain samples have shown that the stress-strain relationship behaves nonlinearly for amplitudes as low as $1\%$ Darvish and Crandall (2001); Donnelly (1997), which is well below the strain threshold for injury Wright and Ramesh (2012); Margulies and Thibault (1992); Sullivan et al. (2015). Therefore, even mild injuries typically occur within the nonlinear elastic regime. Furthermore, the brain’s shear wave speed ($c_{t}\sim$2 m/s) is three orders of magnitude smaller than the compressional wave speed ($c_{p}\sim$1500 m/s) therefore the deformation from an impact is almost entirely in shear mode.

Outside of ultrasound, current methods that measure brain motion cannot capture this nonlinear behavior of shear waves due to limitations in either, frame rate, penetration, or motion detection accuracy. For instance, magnetic resonance elastography (MRE) has been used to investigate shear wave propagation in the brain, typically at frame rates on the order of tens of images/second. Furthermore the temporal sampling of MRE is fundamentally limited to a few ms by the finite integration time over the spin relaxation Hamhaber et al. (2010); Johnson et al. (2013). The harmonics generated by nonlinear shear wave propagation exceed these sampling capabilities. Optical methods have larger frame rates and they can measure motion in optically transparent materials Razani et al. (2012); Song et al. (2013); Razani et al. (2014). However, they are limited to shallow penetration depths in soft tissue ($\sim 2$ mm) Park et al. (2015).

The only experimental corroboration of nonlinear shear waves and its characteristic harmonic signature was made over a decade ago in a homogeneous gelatin phantom Catheline et al. (2003). These experimental methods were derived from ultrasound elastography Tanter et al. (2002) and performed with a high frame-rate (3000 images/second) ultrasound scanner and a cross-correlation tracking algorithm that measured the inter-frame displacement. The high imaging rate was required to capture the rapid, transient, and broadband nature of the shocks. At these high frame-rates, the movement between successive acquisitions is small ($\sim 1~{}\mu$m) and highly accurate motion measurement is required. Although these gelatin phantoms were calibrated to have the same linear elastic properties as soft tissue (Young’s modulus) they did not mimic the nonlinear and attenuating properties of brain, which are key components needed to characterize the shock wave physics.

Despite this experimental report of shear shock waves in a homogeneous gelatin, there have been no reports of shear shock waves in the human body, including the brain. This is due to the difficulty of observing nonlinear shear displacements at depth and at high frame-rates in the brain which include: 1) the more challenging imaging environment, arising from the heterogeneous acoustical properties of its biological structures, and 2) to the higher attenuation of the shear waves. To measure the shear shock wave propagation in ex vivo porcine brain presented in this letter, we address these two challenges. We have designed a focused high frame-rate ultrasound imaging sequence that has better penetration and is more accurate compared to previous plane-wave methods and it is described here for the first time. This data is then processed with an adaptive tracking algorithm that can accurately characterize the almost discontinuous velocity profile at the shock fronts (see appendix C and Pinton et al. (2014)). With these advancements, motion is detected at high frame-rates and with high accuracy for all points within an imaging field of view as deep as the brain. It is shown that shear shock waves are spontaneously generated in the brain, i.e. a smooth excitation at the brain surface develops into a shock as it propagates within the brain. We show that the strong gradients, at the shock front, dramatically increase the local acceleration, strain and strain-rate in soft tissue which suggests that this is a primary injury mechanism for a broad range of TBI’s.

The experimental design used to generate and measure shear shock waves in the brain is illustrated in Fig. 1. Porcine brains (8 week old Yorkshire cross pigs), were extracted 1 hour post-mortem to minimize tissue degradation and were embedded in gelatin. Graphite powder (10-100 $\mu$m) was added to the gelatin with a concentration of 1 % to generate acoustical contrast for the ultrasound images. The shear wave speed of the gelatin was calibrated by setting the concentration of 225 bloom gelatin to 5.0% by weight. This yielded a shear wave speed equivalent to measurements performed in porcine brains ($2.14\pm 0.06$ m/s , see appendix A]). The speed of sound was calibrated by adding 5.0% isopropyl alcohol to obtain a speed of $1490$ m/s Duck (1990). Before the gelatin cured, a 6 mm thick acrylic plate was embedded near the surface of the brain, as shown in Fig. 1. Then, the plate was connected to an electromechanical shaker and the gel was allowed to cool and cure in an ice water bath.

By matching the shear and compressional wave speeds between the brain and gelatin a high transmission coefficient was obtained at the gelatin-brain interface. This is a non-biofidelic interface since in vivo the impact energy has to be transmitted through layers of skull, meninges, cerebrospinal fluid, etc. The effect of the skull’s shape is also ignored, although simulations indicated that the spherical nature of the skull facilitates the formation of shear shock waves deep in the brain by virtue of its natural focal geometry Giammarinaro et al. (2014). However, the complex boundary conditions at the skull-brain interface would generate wave-fields that are currently not solvable by either theoretical or numerical methods in the nonlinear regime. Instead, narrow band planar waves provide conditions for which theoretical results are known. Typically impacts that cause injury would not generate a pure polarized shear wave. Furthermore they have a short duration and a broad power spectrum.Here we generated a 75 Hz wave train, which is within the range of frequencies measured during head impacts (10-300 Hz) Laksari et al. (2015); Greenwald et al. (2008). It has 5 cycles and a -80 dB Chebyshev envelope. This spectrally narrow signal has good separability in the frequency domain which is necessary for the observation of the harmonic development and its comparison to theoretical predictions of shear shock wave development. Note that the physics of shock wave steepening and formation remains applicable to other loading conditions, such as impulse-loading. This pulse was sent to a power amplifier to generated amplitudes up to 450 m/s2 as measured by a linear accelerometer attached to the plate. This setup generates linearly polarized shear plane waves with a full control of its frequency content. The shear wave induces a movement in the $z$-axis and it is propagated through the $x$-axis.

Shear wave propagation through the brain and gelatin was monitored with a research ultrasound scanner (Verasonics). A 5 MHz commercial diagnostic ultrasound imaging probe (ATL Philips L7-4) was positioned above the phantom, as shown in Fig. 1, using a high precision (10 $\mu$m) robotic arm to monitor different sections of the brain.

Conventional planar and plane-wave compounding Montaldo et al. (2009) were sensitive to off-axis clutter and did not perform optimally for ex vivo brain tracking. Sources of ultrasound image degradation, such as reverberation Pinton et al. (2011) and bright off-axis targets which generated static and dynamic artifacts, decreased the image quality and degraded motion tracking performance. This motivated the development of a new sequence which consists of 16 transmit-receive events focused at different lateral positions, each one with an F-number of 4, and a focal depth of 6 cm. Each transmit-receive event generated 1200 frames of raw acoustical data. Then, a delay-and-sum beamforming algorithm generated, from every event, 8 A-lines of focused RF data (see the highlighted area in Fig. 1). To maintain a high frame rate each transmit-receive event was synchronized to the shear wave generation. Furthermore, to increase the visualization field, the acoustical probe was mechanically displaced to a second position using the robotic arm (see appendix B). Thus, the images shown in this letter were generated from 32 synchronized transmit-receive sequences, with an effective frame rate of 6200 frames/second. The jitter error on the synchronization time was measure to be 380 ns which is equivalent to a 1/35000 of a period at 75 Hz. The experiment was performed in four different brains. For each brain, seven amplitudes were used, ranging between 122-449 m/s2. Five independent realizations were obtained for each amplitude by translating the transducer along the $y$-axis yielding a total of 35 experiments for each brain. All the experiment were performed within 15 hours post-mortem at a temperature of 1oC.

The ultrasound data was used to generate a high frame-rate B-mode movie that illustrates the brain anatomy and the movement induced by the shear deformation (Fig. 1) [25a]. We determined brain motion based on the beamformed RF data with a motion detection algorithm designed to track the large discontinuous velocities associated with shear shock waves Pinton et al. (2014). These shear excitations result in inter-frame deformations that change the relative scatterer positions and reduce the correlation between frames. The tracking algorithm uses a quality weighted three dimensional median filter to iteratively optimize the correlation values and it was calibrated to preserve the high frequency displacements necessary to characterize the sharp shock front (see appendix C for validation using simulated ultrasound imaging Pinton (2012, 2017) and Rusanov solutions of cubically nonlinear shocks).

Three snapshots of the resulting high frame-rate movie of particle velocity [25b] are shown in Fig. 2a for the case where the amplitude of the particle velocity at the plate was $v_{0}=1$m/s which produced a Mach number of $M=c_{T}/v_{0}=0.49$. A quasi-planar wave propagating from left to right can be clearly observed. Due to the high Mach number, the brain develops a nearly discontinuous velocity profile in Fig. 2b.

Describing the brain as a nonlinear elastic solid is theoretically complex due to high order tensors in the strain energy density Landau and Lifshitz (1960); Wochner et al. (2008). By assuming a linear polarization the general description can be reduced to a simpler scalar representation Zabolotskaya et al. (2004), while still retaining the nonlinearities that produce complex waves physics Pinton et al. (2010); Whitham (1974). For a plane polarized shear wave the equation of motion, for the particle velocity $v(x,t)$, can be written as:

[TABLE]

where $\tau=t-x/c_{t}$ is the retarded time respect to the physical time $t$, $x$ is the propagation distance and $\beta$ quantifies the nonlinearity. The term on the right hand side accounts for dissipation, which has a significant influence on the shock dynamics. Equation 1 has similar form to the well known Burgers’ equation Burgers (1948), except that the nonlinear term is cubic rather than quadratic. This equation is almost unstudied from the experimental point of view due to the difficulty of directly observing shear waves at depth in solids. As the wave propagates the peaks and troughs are distorted by the cubic nonlinearity which locally increases the shear wave speed as a function of amplitude. Therefore, the wave will eventually “tip” over forming shocks for both the positive and negative wave phases. The cubic nonlinearity has a recognizable impact on the spectral development that occurs with propagation because it produces only odd harmonics.

The experimental development of this nonlinearity and its very specific odd harmonic signature can be observed in the particle velocity (Fig. 2b) and power spectrum (Fig. 2c). This is direct evidence of shear shock wave development in the brain. The Rusanov scheme, along with previously measured frequency dependant attenuation, were used to fit the experimental data at $x=$6.5 mm to the equation 1. This allowed us to estimate $\beta\approx 13\pm 6$ producing an excellent agreement between the theory(black curve) and experiment Rusanov (1970); Pinton et al. (2010). In the frequency domain, consistently with theoretical predictions, energy is transferred from the fundamental frequency to the odd harmonics. Furthermore, more energy is sent to the odd harmonics as the initial wave amplitude increases (Fig. 2 d). The third harmonic amplitude as a function of the propagation distance (Fig. 2 d) indicates that the peak nonlinearity occurs at a propagation distance between 6-7 mm that is independent of initial velocity. Beyond this distance, viscosity overcomes nonlinear harmonic generation. The close match between the theoretical prediction and observation indicates that nonlinearity is a first order parameter required to accurately describe shearing mechanisms in the brain.

Note that, in disagreement with Eq. 1, a significant but small portion of the fundamental energy is also transferred to even harmonics. This may be due to the brain structure, which is heterogeneous, and therefore violates the assumption of homogeneity. In fact it is surprising that the brain heterogeneity has so little influence on the harmonic development. We believe that this is a consequence of the shear wavelength, 28.5$\pm$0.8mm at 75 Hz, which is larger than the characteristic structures in the porcine brain such as vessels and sulci and which therefore leads to weak wave–brain structure interactions.

The particle acceleration, which is more relevant to brain injury than particle velocity Rowson and Duma (2013), was calculated with a broadband derivative in frequency Oppenheim et al. (1983). Numerical differentiation is noise sensitive, especially at the sharp shock profiles (see Fig. 2b) and due to false peak errors from of the tracking algorithm Pinton et al. (2014). A number of validation steps were therefore performed to determine the error in the derivative calculations (see appendix D). They are shown to be in a good agreement with respect to numerical simulations.

The resulting particle acceleration has a large amplitude that is finely concentrated at the shock front as it propagates into the brain (Fig. 3a and 25c]). The shear wave grows, peaks, then dissipates within the first few centimeters of propagation. As a shear wave with an initial acceleration of 449 m/s2 propagates into the brain the steepening of the shock front amplifies the velocity gradient until it reaches 4200 m/s2 at its peak, at $x=6.3$ mm of propagation (Fig. 3b). Thus the acceleration in the brain is 9.3 times greater than the acceleration at the surface of the brain. The peak acceleration into the brain consistently occurs at $x=6.3$mm for the range of initial conditions considered here (Fig. 3c). Furthermore the maximum acceleration within the brain increases rapidly as the surface excitation increases, i.e. it is strongly nonlinear (Fig. 3d). At the maximum surface excitation the average amplification factor is 8.5.

In addition to acceleration, shear strain and strain rate are often used to predict brain injury Wright and Ramesh (2012); Sullivan et al. (2015). Both are shown in Fig. 4. The shear strain has an amplitude of 0.12 and 0.28 at the surface and at 6.3mm, i.e. an amplification factor of 2.3. The nonlinear distortion of the strain with propagation, which is similar to what was observed for the velocity, generates a dramatic increase in the strain-rate at the shock front (cf. Fig. 4b and Fig. 2b). This generates an amplification factor of 6.4 with respect to the brain surface. These shear strain and strain-rate amplitudes are above the thresholds for injury estimated for TBI (0.07 of strain and 40s*-1* for strain-rate) Sullivan et al. (2015).

In a real life impact, the skull, reflective fluid-solid interfaces and non-polarized initial conditions will generate more complex and disorganized shear waves than a plane wave. In this highly complex scenario, we think that at peak locations a shock front will be locally generated, magnifying the local acceleration, shear strain and strain-rate experienced by neuronal tissue. This etiology is consistent with diffuse axonal injury, one of the most common and devastating forms of TBI, which consists of numerous small (1-15mm) lesions that are distributed deep in the brain, far from the area of impact.

In conclusion, by directly measuring brain motion with high frame-rate ultrasound techniques, we have demonstrated that the combination of nonlinearity and attenuation is sufficient to develop shear shock waves over short propagation lengths starting from a smooth, low impact, shear wave in the ex vivo porcine brain. The physical mechanism from which this shear shock wave is generated is consistent with theoretical predictions of cubic nonlinearities that have a specific and readily identifiable odd harmonic signature. The imaging techniques and shock wave tracking methods developed here are not restricted to the brain, and can be applied to soft tissue anywhere in the body that is accessible to conventional ultrasound imaging.

APPENDIX A: Mechanical characterization of the brain

Nonlinear propagation transfers energy from a signal’s fundamental frequency to higher frequencies via harmonic generation. To avoid the complex frequency dependent interaction between nonlinearity, attenuation, and velocity, measurements of the mechanical properties of the brain had to be performed in the linear regime. To characterize these linear mechanical parameters of the brain, a low amplitude $2$m/s2, $75$Hz excitation was used to generate a shear wave. The particle velocity for different propagation distances in the brain, $x=0$, $3$, $6$, $9$, $12$mm (Fig. 5) show that as the shear wave propagates in the brain, there is a decrease in amplitude but no measurable harmonic generation. This confirms that subsequent measurements of attenuation and dispersion were performed in the linear regime, and are independent of the effects of nonlinearity.

To measure the linear mechanical properties of the brain, the amplitude ($A$) and the phase ($\phi$), as a function of distance of propagation ($x$), were extracted from the Fourier transform of the experimentally measured particle velocity for the fundamental frequency. The attenuation, $\alpha$, was obtained by fitting the expression $A(x)=e^{-\alpha x}$ to the measured amplitude. The average characteristic attenuation length was found to be $L_{A}=1/\alpha=25\pm 2$mm. The shear speed, $c_{t}$, was determined by fitting the experimentally measured phase to a linear phase shift, given by $\phi(x)=x\omega/c_{t}+\phi_{0}$, where $\omega$ is the angular frequency and $\phi_{0}$ is an arbitrary reference phase. The average shear speed was found to be $2.14\pm 0.06$m/s. Here the error represents the standard deviation across the four brain samples.

APPENDIX B: Imaging sequences and validation

To validate the proposed shock tracking ultrasound imaging sequence, referred to as flash focus, two experiments were performed in the brain with a driving amplitude of $449$m/s2. The first imaging sequence consisted of a flash focus sequence at one transducer position. The second sequence consisted of the gold-standard conventional plane-wave compounding scheme Montaldo et al. (2009), each plane wave had a different angle of propagation in the $x-z$ plane. To be comparable to the flash focus acquisition time we used $15$ angles and one transmit-receive event per angle. These angles were chosen to be equally spaced between $\pm 18^{\circ}$ with respect to the $z$-axis. Each angle was synchronized with the shear wave emission. After the multiple angles were acquired, the data was beamformed to generate $1200$ frames of RF data at $6200$ frames per second.

Fig 6 shows the B-mode image corresponding to the conventional plane-wave compounding (a) and the flash focus sequences (b). The interframe displacement was calculated from the beamformed RF data with the same adaptive tracking algorithm for both cases and the resulting particle velocity is shown in Fig. 7 a, for plane wave compounding, and b, for flash focus. It is qualitatively apparent from this figure that the flash focus sequence provides a more uniform estimate of motion that is consistent with a planar shear wave. The time plot, shown in Fig. 8 a, is smoother for the flash focus sequence, and unlike the plane-wave compounding sequence, it does not exhibit the discontinuities that are characteristic of false-peak motion tracking errors Pinton et al. (2006). The frequency spectrum for the flash focus sequence (Fig. 8b), is consistent with theoretical predictions of odd harmonic generation, whereas the plane-wave compounding angle exhibits a second harmonic artifact which is much larger than any of the odd harmonics in its spectrum.

For a more quantitative comparison of the improvement generated by the flash focus sequence, the average spectrum was computed for two regions, first in the surrounding gelatin (Fig. 9 a), second in the brain (Fig. 9 b). Fig. 9 a indicates that the two sequences generate similar results to the gelatin. The only significant improvement is a reduction of 14dB of the noise floor when the flash focus sequence was used. In contrast, in the brain region plane wave compounding underestimates the fundamental component by 7.2dB, it over estimates the second harmonic by 10.3dB and it underestimates the third harmonic by 4.4dB. Furthermore, the flash focus sequences reduces the noise floor by 14.1dB allowing visualization up to the ninth harmonic. Therefore, the flash focus sequence is a significant improvement that allows the observation and characterization of nonlinear shear waves propagated in the brain tissue.

APPENDIX C: Validation of the tracking algorithm

The shock front of a shear wave is sharp and consequently has wide support in the spectral domain. This is challenging for tracking algorithms because the numerous high frequency harmonics are at a significantly lower amplitude than the fundamental but they define the shape and rise time of the shock front Pinton et al. (2014). Our experimental setup does not afford opportunities for direct validation of the adaptive tracking algorithm since there is no a priori knowledge of brain motion. Therefore the adaptive tracking algorithm was validated with two simulation tools, one for nonlinear shear propagation Pinton et al. (2010) and another for ultrasound imaging of particle displacements Pinton (2017).

First the particle displacement was determined by solving the nonlinear elastic wave equation (Eq. 1) with a Rusanov scheme Rusanov (1970); Pinton et al. (2010). This method is able to solve nonlinear hyperbolic equations with third order accuracy and low computational time. The linear and nonlinear coefficients in this equation were assigned the values that were measured in the porcine brains, $c_{t}=2.14$m/s, $L_{A}=25$mm, and $\beta=13\pm 6$. The Chebyshev enveloped $75$Hz signal that was used in the experiments was also used as an initial condition to the numerical solution. An initial amplitude of $0.8$m/s was chosen since it represents the high range of what is currently achievable by our experimental setup.

Then the output of this simulation was used to determine the position of the acoustical scatterers in a second simulation of ultrasound propagation. This second simulation tool consists of a finite difference time domain solution of the full acoustic wave equation that we have previously developed and used to model small particle displacements and to generate highly realistic ultrasound images Pinton (2012, 2017). Here it was used to model the field generated by the ultrasound imaging system with the same transducer geometry and transmit-receive imaging sequence. The acoustical maps generated for the simulation had the gross acoustical properties of soft tissue, i.e. an average density of $1000$kg/m3, a sound speed of $1540$m/s, and an attenuation coefficient of $0.3$dB/MHz/cm. To obtain accurate speckle statistics subwavelength scatterers with a size of $39\mu m$, a concentration $>20$ scatters per resolution cell, and a random speed of sound between $1540$m/s and $1502$m/s were distributed randomly in the medium. A separate acoustical map was generated for each of the $1200$ frames and in each map the scatterers were displaced according to the time dependent displacement previously generated by the numerical solution of Eq. 1.

The macroscopic view of two scatters maps, one displaced respect to the other, is shown in Fig. 10 a. In the mesoscopic scale (b) the reference scatter map is shown in red and the displaced scatter map is shown in blue. Note that the displacement changes as a function of the lateral position. Subresolution displacement is represented by shifting the relative impedance of the two element that represent a single scatter (Fig. 10c). This is known as the impedance flow method Pinton (2017).

In this scatter the displacement is less than a pixel, for that reason the displacement is imposed using an impedance flow method. A transmit-receive acoustical simulation was performed for each map to generate the raw RF data in silico. The simulated RF data was matched to experimental data by downsampling and then adding white noise to achieve the $20dB$ SNR observed in the ex vivo brain experiment. Then the RF data was beamformed using the same aperture and delay-and-sum methods that were previously used experimentally. This yielded the $1200$ frames that were finally processed with the adaptive tracking code to estimate particle motion Pinton et al. (2014).

The imposed shear wave calculated from Eq. 1, is shown as the red curve on the left plot in Fig. 11. It is almost indistinguishable from the ultrasound-based motion estimate shown as the black curve (the two are superimposed) and the root-mean-square difference between these two curves is $2.1$%. The corresponding power spectrum (Fig. 11b) is also superimposed which confirms that the ultrasound imaging sequence and adaptive tracking algorithm preserve the spectral content of the particle motion and can accurately capture the higher harmonics that are necessary for the characterization of the shock front.

The final validation step for the tracking algorithm and its ability to represent the harmonics was performed with experiments driven at a low amplitude that was assumed to be linear ($0.03m/s$). In this experiment, shown in Fig. 5, the second and third harmonics represent less that $2.0\%$ of the amplitude of the fundamental component. Thus, the waveform does not distort nonlinearly with propagation. This not only confirms the initial assumption of linearity but also eliminates the possibility that the nonlinear behaviour observed experimentally was a tracking-algorithm based artifact.

APPENDIX D: Differentiation Method

Numerical differentiation, on which our estimates of acceleration from velocity rely, can be noise sensitive especially since the derivatives are required to be accurate at the steep shock front. The acceleration, $a(t)$, was computed from the experimental data with a standard derivative in the Fourier domain Oppenheim et al. (1983).

To determine the expected error in acceleration the derivative of the known velocity determined by the simulation of Eq. 1 was compared to the derivative of the ultrasound-based velocity. Since the imposed velocity was calculated numerically it does not include experimental noise and it was sampled at high frequency (10 times the experimental sampling). Thus, its time derivative is robust and does not depend on the method used for differentiation. On the other hand, the experimental velocity was obtained by applying our tracking algorithm to the simulated RF data as described in the previous section. Therefore, the simulated tracked velocity has similar noise characteristics and the same sampling frequency ($6200$Hz) as the experimental data, where the actual velocity is unknown a priori. These acceleration calculations (Fig. 12) show that, even though this method increases high frequency noise proportionally to $\omega$, it preserves the acceleration peaks at the shock front and underestimates them by approximately $15\%$. This indicates that estimates of the shock front magnification factors reported in the main article are conservatively low.