Real-Time Multifrequency MR Elastography of the Human Brain Reveals Rapid Changes in Viscoelasticity in Response to the Valsalva Maneuver

Modulation of cerebral blood flow and vascular compliance plays an important role in the regulation of intracranial pressure (ICP) and also influences the viscoelastic properties of brain tissue. Therefore, magnetic resonance elastography (MRE), the gold standard for measuring in vivo viscoelasticity of brain tissue, is potentially sensitive to cerebral autoregulation. In this study, we developed a multifrequency MMRE technique that provides serial maps of viscoelasticity at a frame rate of nearly 6 Hz without gating, i.e., in quasi-real time (rt-MMRE). This novel method was used to monitor rapid changes in the viscoelastic properties of the brains of 17 volunteers performing the Valsalva maneuver (VM). rt-MMRE continuously sampled externally induced vibrations comprising three frequencies of 30.03, 30.91, and 31.8 Hz were over 90 s using a steady-state, spiral-readout gradient-echo sequence. Data were processed by multifrequency dual elasto-visco (MDEV) inversion to generate maps of magnitude shear modulus | G∗| (stiffness) and loss angle φ at a frame rate of 5.4 Hz. As controls, the volunteers were examined to study the effects of breath-hold following deep inspiration and breath-hold following expiration. We observed that | G∗| increased while φ decreased due to VM and, less markedly, due to breath-hold in inspiration. Group mean VM values showed an early overshoot of | G∗| 2.4 ± 1.2 s after the onset of the maneuver with peak values of 6.7 ± 4.1% above baseline, followed by a continuous increase in stiffness during VM. A second overshoot of | G∗| occurred 5.5 ± 2.0 s after the end of VM with peak values of 7.4 ± 2.8% above baseline, followed by 25-s sustained recovery until the end of image acquisition. φ was constantly reduced by approximately 2% during the entire VM without noticeable peak values. This is the first report of viscoelasticity changes in brain tissue induced by physiological maneuvers known to alter ICP and detected by clinically applicable rt-MMRE. Our results show that apnea and VM slightly alter brain properties toward a more rigid-solid behavior. Overshooting stiffening reactions seconds after onset and end of VM reveal rapid autoregulatory processes of brain tissue viscoelasticity.


INTRODUCTION
A balance of intracranial mechanical properties is of crucial importance for normal brain function (Linninger et al., 2005;Guyton and Hall, 2006;Wagshul et al., 2006;Schmid Daners et al., 2012). Shear modulus and bulk modulus of brain tissue influence cerebrovascular compliance and pulsatility as well as intracranial pressure (ICP) (Giulioni et al., 1988;Greitz et al., 1992;Wagshul et al., 2011;Parker, 2017). While shear modulus can be measured non-invasively by magnetic resonance elastography (MRE) (Hirsch et al., 2017), there is currently no method for direct ICP measurement without an intervention or without making model assumptions (Guyton and Hall, 2006). In complex multiphasic mechanical systems such as the brain, shear modulus and pressure are linked through poroelastic interactions between the fluid and solid spaces (Bilston, 2002;Tully and Ventikos, 2011;Parker, 2014). Thus, it is likely that regulation of ICP, which is one of the most important vital functions of intracranial mechanics, also affects shear viscoelasticity (Perrinez et al., 2009;McGarry et al., 2015;Lilaj et al., 2020). However, this mechanical component of cerebral autoregulation is largely unstudied due to a lack of imaging techniques that can measure cerebral shear modulus in vivo with high spatial and temporal resolution.
Faster techniques including time-harmonic ultrasound elastography (Tzschatzsch et al., 2018;Kreft et al., 2020) and real-time MRE (rt-MRE) (Schrank et al., 2020a) have been introduced recently. While cerebral ultrasound elastography is limited by acoustic windows and cannot generate detailed maps, rt-MRE has the potential to map viscoelasticity with both high spatial resolution and high frame rates. However, feasibility of rt-MRE has as yet only been demonstrated with a small field of view in the lower extremities (Schrank et al., 2020a) and has never been tested in the brain.
Therefore, we here introduce real-time multifrequency MRE (rt-MMRE) for applications in the human brain. Multifrequency extension of rt-MRE was motivated by previous work on multifrequency wavefield inversion promising higher stability and consistency of parameter maps than single-frequency direct inversion (Papazoglou et al., 2012;Hirsch et al., 2014). Moreover, rt-MRE builds on continuous stroboscopic sampling of harmonic vibrations (Schrank et al., 2020b), which can be spectrally decomposed into multifrequency vibrations without extra scan time. As such, rt-MMRE is a natural extension of rt-MRE that yields, at no extra cost, consistent viscoelasticity maps at relatively high frame rates in the order of 6 Hz depending on the repetition time (TR). Since rt-MRE does not require gating and provides multiple viscoelasticity maps per second, we consider this method as a real-time imaging technique.
Using rt-MMRE, we investigate rapid viscoelastic changes during cerebral autoregulation associated with the Valsalva maneuver (VM). The VM is a standard maneuver to voluntarily increase ICP by forceful breathing against the closed airway with abdominal muscle contraction at the same time. VM will be compared with normal breath-holds in inspiration (BH-in) and expiration (BH-ex). To address frequency dispersion and to test the overall consistency of the values measured in association with the VM, the experiment is repeated with a second set of drive frequencies.
Overall, this study has two aims: first, we introduce rt-MMRE based on three simultaneous excitation frequencies to acquire hundreds of viscoelasticity maps within less than 1 min of scan time. Second, we explore cerebral autoregulation with the unprecedentedly high spatiotemporal resolution offered by rt-MMRE.

Experimental Setup
All experiments were performed in a 3T MRI scanner (Siemens MAGNETOM Prisma, Erlangen) using a 32-channel head coil. Triple-harmonic vibrations in a narrowband frequency regimen were synchronously induced by four pressurized air drivers attached to a transmission plate and placed underneath the head. The applied frequencies were: 30.03, 30.91, and 31.8 Hz (hereinafter referred to as 31-Hz regimen). The two outmost drivers were operated at the highest frequencies with alternated phases relative to each other. The two inner drivers were operated with the same frequency, again with alternated phases. This way each frequency induced mainly lateral-rotational head motion with minimized compression components. The setup is shown in Figure 1.
Vibrations and radiofrequency (RF) excitation started 5 s before data acquisition to ensure establishment of steady states of time-harmonic oscillations and magnetization before start of each experiment. The following rt-MMRE experiments were performed: i Valsalva maneuver (VM) ii deep inspiration and breath-hold (BH-in) iii expiration and breath-hold (BH-ex) The VM experiment included four consecutive phases: 30 s baseline, 5 s breath-hold in inspiration, 20 s VM and 35 s recovery (total scan time: 90 s). Prior to the experiment, subjects were trained to perform a moderate Valsalva maneuver that could be easily sustained for 20 s to prevent involuntary movement after deep breathing. This experiment was repeated with a second narrowband frequency regimen comprising 40.77, 41.67, and 42.55 Hz (hereinafter referred to as 42-Hz regimen) in order to check the overall consistency of MRE during VM and if there is a noticeable influence of frequency.
BH-in and BH-ex experiments consisted of 30-s baseline acquisition with the volunteer breathing normally, followed by a 25-s breath-hold in inspiration or expiration, and a final 35 s recovery phase (total scan time: 90 s).
A resting period of at least 30 s was observed between the experiments. Start and stop commands were given as visual signals to the volunteers. The finger pulse was continuously recorded to track changes in heart rate.
Additionally, anatomical images were acquired using a T1weighted, turbo-spin echo (TSE) sequence.

rt-MMRE Pulse Sequence
Single-frequency rt-MRE using a 2D single-shot gradient echo MRE pulse sequence with spiral readout was recently introduced for directly mapping skeletal muscle function (Schrank et al., 2020a). For rt-MMRE we used a similar prototype-a single-shot, gradient-echo sequence with dual-density spiral readout, which samples multifrequency vibrations in a stroboscopic fashion as illustrated in Figure 2. TR was 62 ms including RF excitation with 20 • flip angle, 20 ms TE, 28 ms readout length, signal spoiling and fat suppression. For motion encoding, a single-cycle, bipolar motion-encoding gradient of 17.5-ms duration (57 Hz) and 40-mT/m amplitude was deployed within each TR according to the principle of fractional encoding (Rump et al., 2007). Images were reconstructed using the SPIRiT non-Cartesian parallel imaging technique (Lustig and Pauly, 2010). Three Cartesian motion components were encoded in an interleaved fashion within the series of consecutive TRs, yielding a sequence of 1,458 wave images. Collapsing these three components into a single viscoelasticity map resulted in a total MRE frame rate of 3 × TR = 186 ms, or approximately 5.4 Hz.
Data were acquired in a single transverse image slice with a field of view (FoV) of 192 × 192 mm 2 and 2 × 2 × 5 mm 3 voxel size. The slice was automatically positioned using the localizerbased auto-align function at the level of the basal nuclei along the largest diameter of the lateral ventricles in the sagittal plane as shown in Figure 1.

Parameter Reconstruction
The 1,458 raw, complex-valued MR images were smoothed with a Gaussian filter (σ = 0.65 px) and subsequently unwrapped using gradient unwrapping (Hirsch et al., 2017). The three vibration frequencies were decomposed by temporal Fourier transformation. Due to stroboscopic undersampling of vibrations in rt-MMRE, the frequencies appeared at aliased positions in the spectrum (see Figure 3A). The frequencies were selected FIGURE 2 | Steady-state gradient echo timing diagram with spiral readout trajectory for single-shot multifrequency real-time MRE. From top to bottom: Experimental design with timing, harmonic vibrations at three frequencies over a period of 9 × 3 TRs with stroboscopic image acquisition below, interleaved wavefield encoding over 2 × 3 TRs, simplified sequence diagram with combined RF and x-gradients to illustrate fat saturation, RF excitation, motion encoding and spiral readout by x-gradients over a single TR period.
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org by three Gaussian bandpass filters (σ = 0.1 Hz) each of which centered at the expected (aliased) frequency of the fundamental drive frequency. These filters were used for inverse Hilbert transformation to compute complex-valued wave fields (wave images) for each vibration frequency, separately yielding 4,374 (1,458 × 3 vibration frequencies) time-resolved wave images (Schrank et al., 2020b). Nine wave images of three Cartesian field components and three vibration frequencies (see Figure 3B) were fed into multifrequency dual elasto-visco inversion (Papazoglou et al., 2012), yielding 486 (4,374/3 encoding components/3 vibration frequencies) consecutive maps of stiffness (| G * |) and loss angle (ϕ) with 5.4-Hz frame rate over the entire examination time. While | G * | is a measure of stiffness, ϕ describes the ratio of elastic to viscous tissue properties indicating fluid properties as explained in Streitberger et al. (2020) | G * | and ϕ maps from the beginning and end of the series were discarded within 5s margins to minimize transient effects introduced by periodic boundary conditions of the Hilbert transform. Consequently, the final observation window was 80 s. All data processing was done in MATLAB (version 2020a). The inversion pipeline is publicly available at https://bioqic$-$apps.charite.de (Meyer et al., 2019). Main results are given in Table 1, Table 2

Parameter Analysis and Statistical Tests
For every time frame, | G * | and ϕ were quantified by averaging values over the same region of interest (ROI). ROIs were manually drawn based on anatomical T1-weighted images, as shown in Figure 1A. Furthermore, these ROIs were refined by empirical thresholds of 10 (time-averaged MRE signal magnitude) and of 950 Pa (time-averaged | G * | map) to remove ventricles and larger sulci similar to Shahryari et al. (2020) (see Figure 4).
The same ROI was also used to determine magnetization signal-to-noise ratio (SNR) and wave displacement SNR (WSNR) for every time frame. WSNR was derived using the blind noise estimation method proposed by Donoho et al. (1995) as outlined and previously applied to MRE data in Bertalan et al. (2019b) and Schrank et al. (2020a) This noise estimation method is suited for wave image analysis since the spatial frequencies of MRE waves and noise are well separated in the wavelet domain (Selesnick et al., 2005;Barnhill et al., 2017).
To test if multifrequency inversion yields more stable values than single-frequency inversion we determined the coefficient of variation (CV) during the baseline phase prior to VM, BH-in and BH-ex for both | G * | and ϕ in all volunteers. The same raw data was used, but for the single-frequency inversion only one frequency from the temporal Fourier spectrum was selected.
We further analyzed difference | G * | and ϕ values relative to mean baseline values given as | G * | = | G * | (t) − | G * | (baseline) (correspondingly for ϕ) in order to quantify individual parameter changes. In addition, peak viscoelasticity values and their temporal delays relative to the onset and end of VM were identified and tabulated for each volunteer.
Finally, group statistics was applied to the absolute values of | G * | and ϕ, after temporal averaging over the following experimental phases for each participant: Of note, these time intervals were given by the aforementioned study design (30-25-35 s. for baseline-breathold/VM-recovery) minus 2.5 s transition phases at the beginnings and ends of these phases including an additional late-response phase. The transition phases were discarded from our analysis in order to minimize transients resulted by the frequency bandpass filter. Also, 5 s BH (25-30 s) and 10 s of post-VM (60-70 s) were considered as transition phases and henceforth not included in our group statistical analysis. All phases are demarcated in Figure 5.
To test possible deformations of lateral ventricles due to VM as reported previously (Ertl-Wagner et al., 2001), we applied automatic segmentation of cerebral spinal fluid (CSF) to the temporal averaged MRE magnitude images of the different experimental phases using SPM12 (Penny et al., 2011; see Figure 4). CSF probability maps were thresholded at 0.5 to generate logical CSF-associated voxel masks. A linear mixedeffects model with varying intercept was employed. CSF volume was used as dependent variable and the individual phases as independent variables. Participants were assigned as random effect, and P-values were calculated using Tukey's post hoc test with Bonferroni correction for multiple comparisons. To test for significant changes in | G * | and ϕ between phases (1)-(4), a linear mixed-effects model with varying intercept was employed. | G * | and ϕ were used as dependent variables and the individual phases as independent variables. Participants were assigned as random effect, and P-values were calculated using Tukey's post hoc test with Bonferroni correction for multiple comparisons. This test does not account for inter-individual slope variations of | G * | and ϕ but analyzes the significance of temporal changes of these parameters. SNR and viscoelastic parameters were correlated using a linear mixed model with | G * | and ϕ as dependent variables and SNR or WSNR as fixed effects with subjects as random factor. All statistical analysis was done in R (version 3.6.2). Unless otherwise stated, errors are given as standard deviation (SD). Correlations between viscoelastic baseline values as well as individual peak responses and participant characteristics (see Table 1) were analyzed using Pearson's correlation coefficient. P-values below 0.05 were considered statistically significant.

RESULTS
Variation in baseline | G * | and ϕ was smaller when using multifrequency inversion (CV = 0.74%, 0.51%) than single frequencies (CV = 0.99%, 0.77%, P < 0.001). Figure 4 shows representative time-averaged MRE magnitude images, automatically segmented CSF masks, as well as | G * | and ϕ maps acquired during the four phases of the experiment. A slight increase in | G * | was visible in the late VM response, whereas no response of ϕ was apparent in individual maps. Group statistics revealed no significant change of CSF-associated voxels between the different states of the maneuver. A descriptive statistic for the individual phases of the VM experiment in the 31-Hz regimen and for each participant is given in Supplementary Table 2.
The BH-ex experiment showed no clear peak, neither in | G * | nor ϕ. | G * | increased continuously with onset of BHex and reached a maximum 2.5 ± 1.5 s after the end of BH-ex (26 ± 23 Pa, P < 0.001).
Absolute | G * | and ϕ Changes Figure 6 shows boxplots with median effects for different states of the maneuver for | G * | and ϕ. The significance levels, indicated by asterisks, were determined from a linear mixed model analysis with varying intercept and participants as random effect. For the VM, different individual effect sizes were observed; however, all subjects showed an increase in | G * | and a decrease in ϕ due to the maneuver. Averaged | G * | values changed between all phases of the experiment (range: 1,370-1,446 Pa) with significance levels indicated in the figure. Averaged ϕ values changed both from BSL to ESM and again from ESM to LRM (range: 0.784-0.800 rad).
By contrast, | G * | only changed at the start and end of the maneuver in BH-in (range: 1,338-1,372 Pa), whereas ϕ changed between BSL and ESM as well as between ESM and LRM (range: 0.783-0.792 rad). In the BH-ex experiment, | G * | changed between LRM and REC (range: 1,348-1,371 Pa) while ϕ changed between ESM and LRM (range: 0.786-0.791 rad).
Results of the second VM experiment performed using the 42-Hz regimen are presented in Supplementary Material. No significant differences in viscoelastic responses between the 31-Hz and 42-Hz regimen were observed (P = 0.24).
Descriptive statistics of | G * | in Pa and ϕ in rad for the individual phases of the VM experiment and for each participant are summarized in Table 2A. Correspondingly, statistical results  -ex). Asterisks at the top demarcate significant changes in | G*| and ϕ which were determined from a linear mixed-effects model with varying intercept. | G*| and ϕ assigned dependent variables and the individual phases as independent variables. Participants were assigned as random effect, and P-values were calculated using Tukey's post hoc test with Bonferroni correction for multiple comparisons. (*P < 0.05, ***P < 0.001).
for the breath-hold and VM experiments performed with the 42-Hz regimen are presented in Supplementary Tables 1a-c. Participant characteristics did not correlate with | G * | or ϕ .

SNR Analysis
Time-averaged SNR and WSNR values did not change significantly across volunteers and over time (P = 0.43). Mean SNR was 29 ± 2 dB across all volunteers with minor and insignificant variations of ± 0.5 dB over the course of the experiment. Mean WSNR was 36 ± 2 dB with minor and insignificant variations of ± 1 dB over the course of the experiment. Significant correlation between group mean | G * | and ϕ was observed (31-Hz regimen: R = −0.4, P < 0.001, 42-Hz regimen: R = −0.5, P < 0.001).

DISCUSSION
This paper presents a novel rt-MMRE technique for the in vivo measurement of rapid and non-periodic changes in brain viscoelasticity in humans. MRE exploiting stroboscopic sampling of multifrequency harmonic vibrations revealed the viscoelastic response of brain tissue to the Valsalva maneuver. Overall, the extension of rt-MRE to rt-MMRE by simultaneous excitation of multifrequency oscillations has increased the consistency of our measurements without adding scanning time. Probably for this reason, all subjects consistently showed an increase in | G * | and a decrease in ϕ with VM, resulting in high statistical significance. This basic finding is remarkable, since the VM is known to induce variability by subjective pressure generation. To further discuss our results we start by briefly reviewing the basic effects of VM on cerebral perfusion and ICP.

Physiological Effects of VM on Cerebral Blood Flow, ICP and MRE
In this study, elevation of intrathoracic pressure during VM was induced by deep inspiration following and increased abdominal pressure similar to the maneuver used in Ipek-Ugay et al. (2017). With onset of VM and elevated intrathoracic pressure, arterial blood pressure (ABP) increases (Elisberg, 1963;Smith et al., 1987). Intrathoracic pressure is communicated through the vascular tree into the cranial cavity, leading to a transient increase in ICP and obstruction of venous outflow from the brain with, thus, increased venous pressure (Prabhakar et al., 2007). Reduced venous return to the heart causes ABP to decrease. Hence, cerebral perfusion pressure is reduced, leading to a reduction in cerebral blood flow (CBF). Cerebral autoregulation is a mechanism to maintain constant CBF. For this reason, cerebral autoregulation, after the decrease in CBF, immediately FIGURE 7 | Diagram of average viscoelastic response with typical mean arterial pressure reported by Pstras et al. (2016) and heart rate variations during and after the Valsalva maneuver with inspiration before the maneuver as measured in this study.
responds to reduce vascular resistance by dilating the cerebral arteries in order to facilitate blood flow and maintain stable CBF. At the same time, the heart rate is increased through the baroreflex (Eckberg, 1980;Looga, 1997), which restores normal ABP and accumulation of blood in the brain, since venous return is still diminished. Constant influx of blood with reduced outflow steadily increases ICP. With release of intrathoracic pressure, there is a significant drop of ABP (Stone et al., 1965), and ICP returns to normal. As a result, normal venous return is restored and more blood flows back into the heart, leading to a transient increase in cardiac output and overshoot in ABP. Since vascular resistance is still low, CBF overshoots as well.
The time curves of MRE parameters presented in Figure 7 suggest that stiffness (| G * |) correlates with ICP while viscosityrelated ϕ correlates with reduced venous outflow or cerebral perfusion pressure. Perfusion pressure is proportional to CBF normalized by mean vascular diameter (Hetzer et al., 2019) and, thus, decreases upon vasodilation with constant CBF. The ramp-up of | G * | during the continuing VM phase seems to reflect the increasing heart rate and steady accumulation of blood in the brain, which drives ICP. By contrast, ϕ remains low throughout the VM phase as if viscous damping in brain tissue is lower when perfusion pressure is reduced. It is an intriguing result that possible ICP changes during VM can be indirectly monitored using rt-MMRE since non-invasive ICP measurement are still an unsolved problem. These findings could help to relate pathologically increased ICP to overall brain stiffness for clinical applications.
In previous work we observed an increase in ϕ of the brain due to hypercapnia (2% increase) (Hetzer et al., 2019) and arterial pulsation (0.5% increase) (Schrank et al., 2020b). In both studies, there was an increase in CBF with a concomitant increase in perfusion pressure while, as explained above, CBF in VM is, due to cerebral autoregulation, associated with a fairly constant CBF and reduced perfusion pressure. Together, the two rt-MMRE parameters, | G * | and ϕ, provide complementary information on the concert of physical parameters involved in ICP autoregulation.
Overall, our baseline parameters of brain viscoelasticity are in good agreement with previously reported values acquired in similar frequency ranges (Hetzer et al., 2018(Hetzer et al., , 2019Schrank et al., 2020b). We observed no significant differences in the responses of rt-MMRE parameters to VM between the 31 and 42 Hz regimens. This consistency of multifrequency data further validates the technique of rt-MMRE. Furthermore, this observation indicates that the poroelastic response of brain tissue  is similar at 30 and 40 Hz (McGarry et al., 2015). Additional validation of rt-MMRE was obtained by reference experiments performed during breath-holds but without sustained VM. BH-in induced a similar increase in stiffness and decrease in ϕ as observed during VM. Thus, from an MRE perspective, deep inspiration followed by breath-holding induces effects similar to a light VM. Otherwise, no such changes were observed in BH-ex, rendering this maneuver neutral with regard to ICP. Nevertheless, even BH-ex had some small effect on MRE parameters, which, notably, were not correlated to changes in either SNR or WSNR. Also, analysis of CSF volume and ventricle size did not reveal any significant correlation with VM. Previous work by us and others showed that total brain volume increases due to VM by approximately 3% while ventricle volume shrinks by 20% (Ertl-Wagner et al., 2001;Mousavi et al., 2014). In contrast to these studies, our subjects were instructed to perform a moderate Valsalva maneuver to minimize variations in thoracic pressure, muscle strain, and head position.
In our previous work we used ultrasound time-harmonic elastography in a temporal bone window to acquire VM-induced rapid changes in shear wave speed in the temporal lobes of healthy volunteers (Tzschatzsch et al., 2018). Effect sizes in that region were higher (10.8 ± 2.5%) than revealed by MRE in the full brain tissue slice. It should be noted that the regions covered by our current study do not correspond to the medial temporal gyrus addressed by transtemporal timeharmonic elastography, which makes a direct comparison of effect sizes between the two studies difficult. To analyze the spatial representation of viscoelasticity changes we performed automatic image segmentation using MNI-based registration as well as voxel-wise correlation analysis based on a boxcar function. No significant patterns of viscoelasticity changes could be detected. A more detailed analysis of the spatiotemporal representation of brain viscoelasticity in response to the VM is warranted.
Our study has limitations. The nature of stroboscopic sampling of vibrations by steady-state single-shot acquisitions limited our technique to 2D wave field sampling including three encoding components. This intrinsic limitation of rt-MMRE can currently not be overcome by a multishot variant because VM is a non-periodic event and cannot be repeated with enough temporal reliability. Consequently, our multifrequency inversion technique was entirely 2D, which may have led to variability due to different slice positioning and oblique intersection of 3D shear wavefields rendering our values as effective viscoelasticity parameters. Nevertheless, our conclusions are drawn from group values in two different frequency regimens. The fact that these values changed with statistical significance in VM, while neither SNR nor anatomy changed, emphasizes the robustness of the observed MRE effects. Furthermore, our data could be used for suppression of bulk waves based on the in-plane curl component. However, this curl-analysis did not provide more consistent values than our standard MDEV inversion with respect to confidence intervals and statistical power. Finally, 2D brain MRE has a long tradition in disease detection (Wuerfel et al., 2010;Streitberger et al., 2012Streitberger et al., , 2017Streitberger et al., , 2020Lipp et al., 2013Lipp et al., , 2018Fehlner et al., 2016;Gerischer et al., 2018) as well as in the study of brain physiology (Sack et al., 2009(Sack et al., , 2011Schrank et al., 2020b;Herthum et al., 2021). It remains to be determined whether single-frequency 3D MRE can provide similarly consistent clinical and physiological brain data. Instead, as shown herein, unintentional breath-holds may affect 3D MRE due to long scan times. Generally, current MRE techniques cannot account for poroelasticity, heterogeneity, hyperelasticity, anisotropy and temporal variations of brain tissue at the same time. Therefore, to date all values measured by brain MRE should be considered as effective parameters.
In summary, we studied the viscoelastic response of the human brain to breathing and the Valsalva maneuver using a novel real-time multifrequency MRE technique. Significant increases in brain stiffness and decreases in ϕ due to VM were observed with use of two different frequency regimens. Control experiments showed that breath-holds after inhalation induce a response similar to VM but with a smaller effect size. By contrast, breath-holds after exhalation had the smallest effects on cerebral MRE parameters. The time courses we report here provide a reference for the VM response in healthy subjects and might be of value for studying dysfunctional autoregulation as associated with various neurological diseases. rt-MMRE is a fast technique which can provide consistent imaging markers of brain viscoelasticity within a fraction of a minute.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics committee of Charité -Universitätsmedizin Berlin in accordance with the Ethical Principles for Medical Research Involving Human Subjects of the World Medical Association Declaration of Helsinki. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
HH carried out all experiments and contributed to all parts of the manuscript. MS assisted in interpreting the results with regard to physiological changes and statistical analysis. FS, HT, and SH helped to carry out the data processing and verified the results. CW and JP carried out the MRI sequence implementation and image reconstruction. SG and HN contributed to the experimental setup. JB helped supervise the project and constructed the actuation system. IS designed and directed the project and aided in interpreting the results. All authors provided critical feedback and helped shape the research and manuscript.

FUNDING
Funding from the German Research Foundation (GRK 2260 BIOQIC, SFB1340 Matrix in Vision, Neuro-MRE Sa901/17-2) and from the European Union's Horizon 2020 Program (ID 668039, EU FORCE -Imaging the Force of Cancer) is gratefully acknowledged.