Fundamental modeling of wave propagation in temporally relaxing media with applications to cardiac shear wave elastography

: Shear wave elastography (SWE) might allow non-invasive assessment of cardiac stiffness by relating shear wave propagation speed to material properties. However, after aortic valve closure, when natural shear waves occur in the septal wall, the stiffness of the muscle decreases signiﬁcantly, and the effects of such temporal variation of medium properties on shear wave propagation have not been investigated yet. The goal of this work is to fundamentally investigate these effects. To this aim, qualitative results were ﬁrst obtained experimentally using a mechanical setup, and were then combined with quantitative results from ﬁnite difference simulations. The results show that the amplitude and period of the waves increase during propagation, proportional to the relaxation of the medium, and that reﬂected waves can originate from the temporal stiffness variation. These general results, applied to literature data on cardiac stiffness throughout the heart cycle, predict as a major effect a period increase of 20% in waves propagating during a healthy diastolic phase, whereas only a 10% increase would result from the impaired relaxation of an infarcted heart. Therefore, cardiac relaxation can affect the propagation of waves used for SWE measurements and might even provide direct information on the correct relaxation of a heart.

Cardiac diseases are a major cause of death in developed countries. Early diagnoses might help prevent the development of life-threatening conditions by detecting signs of deterioration before cardiac functionality becomes compromised. Such diagnoses may be obtained by monitoring the stiffness of the cardiac muscle, which has been observed to correlate with the health condition of the heart. [1][2][3] In order to monitor the material properties of the heart, however, non-invasive techniques must be employed, as invasive measurements are highly uncomfortable and potentially harmful to patients.
Shear wave elastography (SWE) exploits wave propagation phenomena to explore the elastic properties of a material, and it has already been proven to be a viable tool in clinical applications. [4][5][6][7] Its application to cardiac settings, however, is hindered by a challenge intrinsic to the functioning of the heart: the heart cycle.
As the heart performs its pumping function, its stiffness increases and decreases cyclically to allow for the heart chambers to fill with blood and expel it. Aortic valve closure, which provides one of the sources of waves that can be employed for cardiac SWE, [8][9][10][11][12][13] takes place at the beginning of the isovolumic relaxation of the muscle; due to the muscle relaxation, the waves generated at this time could experience a change in propagation speed of %15% in just 10 ms (estimation based on the stiffness variation measured in isolated perfused rabbit hearts 14 ). While SWE performed on waves naturally occurring in the heart could be more precisely called "natural" SWE, we will refer to it simply as SWE, because the phenomena involved are similar to those of shear waves from non-natural sources (e.g., acoustic radiation force pushes).
Several models are employed in literature to describe the mechanical properties of the cardiac tissue, [15][16][17][18] however, to the best of our knowledge, all the models employed for elastography assume the mechanical properties of the medium to be constant or, at most, slowly varying 19 in the timescales of the propagating wave. This assumption may hold true for mechanically inactive organs, yet its validity is questionable in the context of SWE measurements performed during diastolic relaxation.  In fact, the existence of measurable effects of timevarying medium properties on propagating waves has already been established in the field of electromagnetism, with theoretical descriptions of media changing smoothly or instantaneously. [20][21][22][23][24][25][26][27] These studies predict that a wave that propagates at varying speed (i.e., in a medium with temporally varying dielectric or magnetic constant) is subjected to a variation in amplitude and oscillation period; additionally, reflected waves are generated at time-discontinuities of the medium, similar to what happens at the spatial interface of two media. Experimental studies 20,28 confirmed the predictions regarding amplitude and frequency changes by observing magnetic waves propagating in media subjected to an externally modulated magnetic field. To the best of our knowledge, reflected waves were not observed experimentally. It remains an open question whether these effects are also present in elastic waves, and whether they should be taken into account while performing cardiac elastography.
The goal of the present study is to observe and describe the effects that a temporally varying propagation speed has on mechanical waves in order to assess their relevance to shear wave elastography of the heart. While a live heart itself could be, in principle, used as a medium to perform these studies, it would be impractical, since its complicated geometry and material properties would make it hard to reliably isolate and identify the specific effects of temporal variations. For this reason, we have chosen to model a simplified setting in which speed variations represent the only complication to one dimensional (1D) wave propagation. We have developed an experimental setup consisting of rotating metal rods suspended by nylon wires (a wave machine) in which the tension can be controlled in real time to alter propagation speed. The rotational displacement of the rods can then travel through the setup as a 1D torsional wave with varying speed. Moreover, we developed finite difference simulations that describe these phenomena numerically. We employed the setup to obtain a first, qualitative confirmation that mechanical waves can also be affected over time by variations of medium properties. The simulation, on the other hand, allowed us to investigate quantitatively how the time-dependent effects are related to the dynamical parameters of the system, i.e., the amount of speed variation and the rate at which this happens.
We apply our findings to data on rabbit hearts 14 to predict the effects of muscle relaxation on cardiac elastography measurements, and we discuss their relevance and possible applications as a new diagnostic tool.

II. SETUP
A. The wave machine We have built a modified wave machine 29 in which the speed of the wave can be controlled during propagation by means of tension variations. As shown in Fig. 1, the setup consists of two wooden frames, placed at 3.63 m from each other, that support three wires on which 32 aluminum rods are suspended. Each rod is 60 cm long and has a 1 cm Â 1 cm square cross-section. The central wire is made of steel and runs through a hole in the midpoint of the long side of the rod, providing a pivot around which the rods can rotate freely; bolts are fixed to the central wire before and after each rod, to prevent them from translational movements. The other two wires, symmetrically placed at both sides of the steel wire, are made of nylon and can freely slide through their holes in the rods, so that stretching of these wires does not cause translation of the rods. The nylon wires provide the restoration force that opposes rotational displacements from the mutual angular position of the rods. During experiments, the amplitude of the applied perturbation had a Gaussian-like shape, as this was the easiest to produce by hand. When a rotational perturbation is applied to one of the rods, it propagates along the setup through the nylon wires, effectively creating a discretized 1D torsional wave. This system can be seen as a discrete approximation of the continuous case in which the distance between two consecutive rods approaches zero. For our case the propagation of the torsional wave can approximately be described by the wave equation where r ¼ 0:060 6 0:001 m is the distance between steel and nylon wires, F is the variable tension in each wire, L ¼ 3:410 6 0:005 m is the distance between the two extremal rods, I ¼ 0:0049 6 0:0004 kg Á m 2 is the moment of inertia of a rod, N ¼ 32 is the total number of rods and h is their angular displacement. This yields a torsional wave speed

B. Tension control
The tension in the wires could be manually controlled, during experiments, by means of a mechanical lever connected to the nylon wires: pulling the lever would result in stretching of the wires, with a consequent increase in tension. Relaxation of a muscle can be mimicked by prestretching the wires, then releasing the lever over a transition time s during wave propagation; with a framerate of 60 frames per second, it was not possible to determine s from the video recordings.
To verify the reliability of Eq. (2), a first experiment was run with a force gauge (Force Gauge TMT-5020, OCS.tec GmbH & Co. KG, Neuching, Germany) connected to one end of the nylon wire, so as to compare the values of c derived from the equation and those measured directly. During this experimental validation, the tension read by the force gauge was F ¼ 13:75 6 0:2 N, corresponding to a calculated speed of c ¼ 1:50 6 0:14 m/s, a value comparable, within the experimental tolerance, with the speed measured directly from the video recordings, c ¼ 1:65 6 0:02 m/s. While the value of speed determined from the tension measurement suffered from a relatively high experimental uncertainty, the uncertainty of the direct speed measurement depends essentially on the framerate of the recording, allowing for more precise measurements. In our experiments, therefore, the propagation speed was measured by acquiring and analysing video-recordings of the wave. Although not necessary, the tension of the pre-stretched nylon wires could be estimated from the speed measurements via Eq. (2): with a speed c 1 ¼ 3:00 6 0:02 m/s, the tension was estimated to be 57.0 6 0.1 N. With the lever in its rest position, on the other hand, the wave traveled at a measured speed of c 2 ¼ 1:40 6 0:02 m/s, corresponding to a tension of F ¼ 12:4 6 0:1 N.

C. Data acquisition
A digital single lens reflex camera (Nikon D5300, Nikon Corporation, Tokyo, Japan) with framerate of 60 frames per second, facing the cross section of the rods, was used to record the propagating wave from one side, with a field of view of approximately 2 m that allowed the imaging of 18 rods in the center of the setup. In order to increase the contrast between the setup and the background, the tips of the rods were painted with an orange phosphorescent paint that reacts to ultraviolet (UV) light. We performed the experiments in a dark room illuminated only by 5 UV 60-W lamps, so that the bright orange glow of the extremities of the rods would be easily distinguishable from the background. This allowed us to isolate and track their motion in post-processing, using the software ImageJ (National Institute of Health, Bethesda, Maryland, U.S.) to isolate the motion of each individual bar, and then importing all data in Matlab (version r2016b, MathWorks, Natick, MA). Frames of the recorded propagating wave are shown in Fig. 2, while Fig. 3 shows a single videoframe and its numerical reconstruction.

III. RESULTS
During the experiments, the nylon wires were first prestretched to the maximum tension of 57 N. A single unipolar wave pulse was then manually generated by perturbing the first metal rod, and the tension was subsequently dropped to 12.4 N by releasing the lever during propagation. In order to avoid the superposition of boundary reflections with the waves we wanted to study, the release of the lever was timed so that the waveform would always be in the center of the system when the tension dropped. The tension drop happens in a fraction of a period. Figure 4 shows the propagation measured during an experiment: the axes represent the time elapsed and the spatial coordinates of each rod expressed by rod number; due to the limited field of view of the camera, the extremal rods were not imaged and therefore do not appear in the plot. The color scheme represents the amplitude of rotation of each bar, so that the bright yellow "bands" essentially correspond to the positive forward traveling wave. At time t % 1:5 s the effects of the tension drop can be seen. A reflected wave appears as the broad band that moves backwards from right to left, with its negative amplitude represented in dark blue. In addition, broadening of the We note a remarkable similitude with the behaviour of waves crossing the spatial discontinuity between two media, the main difference being that, in the case of a temporal discontinuity, it is wavelength that is conserved, while the period changes. In order to better understand the parallelism, let us look at a wave equation in which the wave speed depends on the spatial coordinate Next, we consider the current situation, in which the wave speed depends on time, giving the equation By defining a new parameter sðtÞ ¼ 1=cðtÞ called slowness, we can rewrite the equation as This resembles Eq. (3), but with the roles of z and t being interchanged. From the wave behaviour at a spatial discontinuity, it is known that the spatial period of the wave (i.e., the wavelength) varies proportionally to c, whereas its temporal period remains constant. By considering Eq. (5) and performing the same reasoning as above, that is inverting the roles of space and time, we can expect that at a temporal discontinuity the time period will vary proportionally to s, whereas the spatial period will remain constant.
We then proceeded to implement in Matlab a 1D first order explicit finite difference scheme to solve Eq. (1) numerically and simulate the behavior of the setup in different circumstances. In order to test the viability of the simulation to investigate these phenomena, we compared its results with experimental data: the parameters of the wave function, the space, and the time discretizations were all chosen to match those of the experimental setup, i.e., spatial steps of 11 cm (the distance between rods) and time steps of 16 ms (the time between video-frames). The wave was simulated as a Gaussian pulse with an e À1 width of 0.3 s. A comparison between this function and the excursion of the third rod (the first one to be imaged in our measurements) is shown in Fig. 5. Figure 6 shows a comparison between simulation and experiment, by plotting the amplitude over time of two rods. The rods were chosen so that one would oscillate once before the tension lever was released, while the other would be crossed by the wave only after the sudden tension drop, which was simulated to happen in 16 ms, i.e., one time-step. The results of the simulation (solid line in Fig. 6) are in qualitative agreement with the experiment in terms of amplitudes and period of the incident, transmitted and reflected wave around the transition phase. This validates our numerical approach. Having thus established its reliability, we continued our study by performing simulations only, in order to perform a systematic, quantitative study.
First, we ran simulations to determine the relation between the deceleration of the wave and the formation of transmitted and reflected waves. We computed the behavior of a one-cycle, sinusoidal wave as could be generated on a 1D string. The spatial discretization was refined to a spacing of 0.01 m to increase spatial sampling of the waves, and the time discretization was shortened to steps of 0.1 ms, ensuring stability of the numerical scheme (Courant number C num 0:028) as well as correct sampling of all phenomena. Figure 7 shows how the amplitude and period of the transmitted and the reflected waves change as a function of c 1 =c 2 , for fixed duration s ¼ 10 ms of the deceleration. We can see that both amplitude and period increase linearly with the ratio between the initial and final propagation speed (c 1 and c 2 , respectively), in agreement with the relations A T =A I ¼ ðc 2 þ c 1 Þ=2c 2 and A R =A I ¼ ðc 2 À c 1 Þ=2c 2 detailed in the Appendix, and the relation T ¼ c=k with constant wavelength k; here, A represents the amplitude of the wave, T the period, and the subscripts I, R, and T represent the initial, reflected, and transmitted waves, respectively. When referring to the period of the wave before and after the tension drop, the subscripts 1 and 2 will be used, so that the two periods will be indicated by T 1 and T 2 , respectively. Furthermore, we investigated the effects of a deceleration taking place over longer spans of time, up to about twice the period of the wave. In our numerical model, for t < t 1 the speed was 3 m/s; the speed then decreased linearly from 3 to 1 m/s, between t 1 and t 2 ¼ t 1 þ s, and was consequently kept at a constant value of 1 m/s for t > t 2 . As shown in Fig. 8, as the ratio s=T 1 between transition time and wave period increases, the amplitude of transmitted and reflected wave decreases; the phase of the reflected wave appears to be opposite of that of the incident wave.  Moreover, we can also notice that two reflected waves are actually present when s=T 1 ! 1, whereas for shorter transition times there appears to be only one, albeit distorted, reflected wave. The variation in period of the wave does not seem to be affected by s, as shown in Fig. 9. Based on the results of these simulations, we conclude that reflected waves are generated at points in time when the acceleration of the wave is discontinuous; if two such points are separated by less than T 1 , two reflected waves will be generated, but they will partially overlap with each other, appearing to be a single, distorted waveform. On the other hand, if the distance in time between the two points is greater than T 1 , both distinct waveforms will be visible.
The results of the simulations detailed above offer an overview of the general effects that the temporal variation of a medium has on a propagating wave. In order to assess the relevance of these effects for cardiac SWE, we can input in our model realistic values of stiffness variation of the cardiac muscle to predict how the relaxation would affect a propagating shear wave. Let us first consider the diastolic phase of a heart in which the muscle isovolumic relaxation phase can be modeled by an exponential decrease in stiffness with time constant a ¼ 50 ms, as was measured in Langendorff perfused rabbit hearts. 14 Let us further consider that a shear wave traveling through the muscle is imaged (e.g., by means of ultrasound scanners) for a duration of d ¼ 10 ms during this phase. Figures 7 and 9 showed that the variation in period of the wave depends only on the ratio between initial and final speed, T 2 =T 1 ¼ c 1 =c 2 . The ratio itself can be easily computed knowing a and the d, since c 1 =c 2 ¼ e d=2a . Considering a ¼ 50 ms and d ¼ 10 ms, one could expect to observe an increase in wave period of 20%. On the other hand, an infarcted heart might have a relaxation constant a ¼ 100 ms 14 ; under the same measuring conditions, one would then observe an increase in wave period of just 10%. Therefore, when reconstruction of the shear wave in space-time domain is possible, the broadening of the tracked shear wave in time could be used as an indirect measure of the relaxation that took place. The amplitude variation, however, being in the order of 3% and 2% for healthy and unhealthy hearts, respectively, would be unlikely to be detected in realistic measurements.

IV. DISCUSSION
Our results confirm that mechanical waves show a response to temporal variations of the medium comparable to the behavior of electromagnetic waves detailed in literature: the amplitude and the period of the traveling wave increase proportionally to the decrease in propagation speed, the wavelength is unaffected, and reflected waves can be generated at points of discontinuity in the acceleration.
Based on our numerical results, the predicted effects of cardiac relaxation are of particular interest: in fact, not only would a measured wave be affected by the relaxation of the muscle, it would actually be affected by the specific relaxation curve within the measurement. In other words, if two muscles relax at different rates, they may in principle be told apart by observing the variation in amplitude and period of waves traveling through them. The variations in period, in particular, could potentially be employed to directly help diagnoses of diastolic impairment (i.e., reduced muscle relaxation), which is connected to diastolic heart failure. It should be noted that these results show a relation between initial and final period (or frequency), without any assumption about the initial value itself; therefore, the results are independent of the frequencies involved, and can apply to the frequency ranges typical of naturally occurring shear waves, as well as artificially induced ones.
As promising as these results may be, however, their practical implementation still presents challenges. First of all, the amplitude increase caused by diastolic relaxation is expected to be small and could be hard to measure during in vivo experiments, due to signal attenuation and noise. Moreover, because the relaxation of the muscle happens smoothly, there are no discontinuities to give rise to clearly observable reflected waves. Therefore, the only effect of temporal variations to be visible in clinical measurements may be the shift in wave period, and measurements with high time-resolution could be necessary to distinguish the different shifts of a healthy and a diseased heart. For example, for an initial frequency of 50 Hz, the time difference between an increase in the period by 10% and 20% is 2 ms, corresponding to a minimum imaging framerate of 500 Hz, while an initial frequency of 100 Hz would require a minimum framerate of 1000 Hz to discern the different period increases. Finally, due to the anisotropic, viscoelastic, threedimensional nature of the heart, additional effects will compound to those produced by muscle relaxation, so that more complex models will be required in order to analyze data accurately.
Another matter to be noted is that the results presented in this study only show explicitly the effects of a wave that is generated before the relaxation begins; the relaxation of the medium takes place entirely during the propagation of the wave, and the effects are determined after the entire process has ended. However, when the closure of the aortic valve generates a wave, the cardiac muscle has already started relaxing; moreover, a SWE measurement based on such a wave would be over well before the cardiac relaxation process ended. Therefore, one could question whether the predictions formulated for real hearts can truly be trusted based solely on our numerical study. We argue that this is indeed the case: Figs. 7 and 9 show that the variation in period of the wave depends only on the ratio between initial and final speeds c 1 =c 2 , where c 1 and c 2 correspond the speed of the wave at the beginning and at the end of a measurement; this means that for fixed c 1 =c 2 , the variation in period will always be the same, regardless of whether the speed varies with a step function or a smooth curve. As mentioned above, the difference with our simulations, introduced by the smoothness of the diastolic relaxation, lies in the absence of reflected waves, which are generated at discontinuities. It is important to notice that, during a fixed measuring time, different relaxation curves (as could characterize healthy and diseased hearts) will result in different values of c 1 =c 2 , therefore affecting the propagating wave differently.
Finally, we consider a wave produced by mitral valve closure, which may also be used to perform SWE. The mitral valve closes at the onset of systole, the contracting phase of the heart cycle, and, due to the stiffening of the muscle, waves traveling during this phase experience an increase in propagation speed. Based on our results, we expect that the amplitude and period of an accelerating wave will decrease proportionally to the increase in speed, and reflected waves will be once again generated at discontinuities in the acceleration.
From an experimental point of view, compared with other approaches proposed in literature to investigate waves in temporally varying media, the use of our setup shows several advantages. To begin with, the setup allows us to observe the full phenomenon in its complexity, including not only the initial and final state of the wave, but its transition phase as well. Moreover, it can be easily built and modified so as to tailor the wave parameters (speed, amplitude) to specific needs, and the tension can be varied directly and in real time during an experiment; thanks to these features, more complex interactions (e.g., combining temporal variations with viscous behaviors) could be investigated with minor modifications of the setup.
There are a few limitations in our choice of setup as well: in particular, non-ideal matching of wires and rods, non-ideal fixed boundaries, and friction, are all factors that play a role in altering the behavior of the wave, degrading the match between experiment and idealized numerical simulation, as can be observed, e.g., for the reflected wave in Fig. 6 at times greater than 1.5 s. Moreover, due to the limited resolution and spatial window with which the motion of the bars can be detected by the camera, only wave amplitudes between 3 and 20 cm could be reliably detected. This meant that only a restricted parameter space could be practically investigated with the experimental setup. Nevertheless, the setup provided useful qualitative information, as well as a validation of our numerical approach; together, experiments and simulations proved to be a robust tool to investigate waves in time-dependent media.

V. CONCLUSIONS
We conclude that the variation over time of the stiffness of a medium, such as the beating heart, produces an inversely proportional variation in amplitude and period of a wave traveling through it. Furthermore, if the stiffness variation presented discontinuities, reflected waves would be generated. Based on the results of our numerical simulations and literature data, we predict that a healthy diastolic relaxation will affect propagating waves differently than an impaired process, producing a wave period increase twice as large as the one that might be observed in a dysfunctional organ. This difference may potentially be exploited to directly obtain information on diastolic functionality with one SWE measurement.

ACKNOWLEDGMENTS
This work is part of the STW-Dutch Heart Foundation partnership program "Earlier recognition of cardiovascular diseases" with project number 14740, which is financed (in part) by The Netherlands Organization for Scientific Research (NWO). Moreover, the authors would like to thank Henry Den Bok, from the Delft University of Technology, and Geert Springeling, from the Erasmus MC, University Medical Center Rotterdam, for technical support with the construction of the experimental setup.

APPENDIX: REFLECTION AND TRANSMISSION COEFFICIENTS
We show here how the reflection and transmission coefficients can be calculated for a 1D transverse wave travelling on a string that undergoes an instantaneous material property variation (e.g., a decrease in tension). This situation can be considered to be the temporal equivalent of a wave travelling across the interface between two different, ideally bonded, strings. Let us assume that the instantaneous variation happens at time t 0 ¼ 0, and let us consider an initial wave function w i ¼ f I ðx À c 1 tÞ for all t < 0. As we have seen, after the temporal discontinuity, in the setup two waves are present, a transmitted wave w T ¼ f T ðx À c 2 tÞ and a reflected wave w R ¼ f R ðx þ c 2 tÞ. Notice that the transmitted and reflected waves travel at the same speed c 2 , with opposite signs. We can then consider two instants in time, before (B) and after (A) the discontinuity, and define two waveforms and representing the deformation of the string in these two instants. Let us now make two considerations: (1) The deformation of the string has to evolve continuously from w B to w A at all coordinates x (a discontinuity of w between two consecutive instants requires movements at infinite speeds). Therefore, it follows that for all x: (2) The vertical speed of each individual particle in the string also needs to be continuous over time (a discontinuity of @wðx; tÞ=@t in time would require infinite acceleration). We can write this condition as @w B ðx; tÞ @t t¼0 ¼ @w A ðx;tÞ @t We have thus shown that the ratios of the amplitudes of initial, reflected and transmitted waves (A I , A R , and A T , respectively) are equal to A R =A I ¼ c 2 À c 1 2c 2 (A8) and A T =A I ¼ c 2 þ c 1 2c 2 : (A9)