High-Resolution Mechanical Imaging of Glioblastoma by Multifrequency Magnetic Resonance Elastography

Objective To generate high-resolution maps of the viscoelastic properties of human brain parenchyma for presurgical quantitative assessment in glioblastoma (GB). Methods Twenty-two GB patients underwent routine presurgical work-up supplemented by additional multifrequency magnetic resonance elastography. Two three-dimensional viscoelastic parameter maps, magnitude |G*|, and phase angle φ of the complex shear modulus were reconstructed by inversion of full wave field data in 2-mm isotropic resolution at seven harmonic drive frequencies ranging from 30 to 60 Hz. Results Mechanical brain maps confirmed that GB are composed of stiff and soft compartments, resulting in high intratumor heterogeneity. GB could be easily differentiated from healthy reference tissue by their reduced viscous behavior quantified by φ (0.37±0.08 vs. 0.58±0.07). |G*|, which in solids more relates to the material's stiffness, was significantly reduced in GB with a mean value of 1.32±0.26 kPa compared to 1.54±0.27 kPa in healthy tissue (P = 0.001). However, some GB (5 of 22) showed increased stiffness. Conclusion GB are generally less viscous and softer than healthy brain parenchyma. Unrelated to the morphology-based contrast of standard magnetic resonance imaging, elastography provides an entirely new neuroradiological marker and contrast related to the biomechanical properties of tumors.


Introduction
Despite recent advances in operative and postoperative treatment, glioblastoma (GB) still remains one of the most malicious and aggressive and malignant forms of cancer [1,2]. It accounts for .50% of all primary neuroepithelial tumors and approximately 20% of all brain tumors [3]. In developed countries, the incidence of GB is 3.5 per 100,000 population per year [4,5,6]. The term GB was first introduced in 1926 by Percival Bailey and Harvey Cushing and refers to the cellular origin from glioblasts and the histological heterogeneity of this brain tumor [7]. The classification of the World Health Organization (WHO) ranks GB as a grade IV tumor due to its histological characteristics with aggressive and infiltrative growth and overall poor prognosis [8]. Despite aggressive surgical resection, radiotherapy, and chemotherapy the prognosis for patients newly diagnosed with GB remains poor with a 2-year survival rate of only 13-26% and a mean survival time of 12-15 months [2].
Neuroradiological assessment of GB and differentiation from solitary intracranial metastases or lymphomas is challenging due to the tumor's heterogeneous composition resulting from the presence of cysts, necrosis, and hemorrhage [9,10]. As a consequence, diagnostic biopsy remains inevitable for a definitive diagnosis despite possible complications [11]. Advanced magnetic resonance imaging (MRI) methods such as diffusion tensor imaging provide structural information related to water mobility in white matter tracts but cannot reveal the consistency and mechanical constitution of biological tissue [12,13]. Targeting the mechanical properties of GB potentially provides information about the tumor's structural heterogeneity as well as its perifocal tissue infiltration which is of relevance for diagnosis, therapy planning and treatment monitoring.
Today, the viscoelastic properties of the brain can be assessed noninvasively by magnetic resonance elastography (MRE) [14]. Combining time-harmonic vibrations in the low audible range with motion-sensitive MRI, cerebral MRE [15,16,17] has proven sensitive to mechanostructural changes in the human brain associated with aging [18,19] and diseases [20,21,22,23,24,25]. Recent results in mouse models suggest that cerebral MRE is sensitive to demyelination, inflammation, and extracellular matrix alterations [26,27,28] and may thus provide new information about structural changes of cerebral tissue in brain tumors. Indeed, initial findings in different intracranial tumor entities [29] including meningeomas [30,31] indicate the feasibility of MRE for the presurgical assessment of neuronal tumor consistency.
However, previous studies were limited by low spatial resolution of the mechanical parameter maps achievable by MRE at a single harmonic drive frequency. Recent advances in fast image acquisition schemes [32,33] and wave field reconstruction methods [34,35,36,37] enabled us to acquire 3D wave fields at multiple vibration frequencies, generating cerebral MRE maps with a spatial resolution comparable to that of normal MRI [38,39].
In this study we applied multifrequency MRE (MMRE) for in vivo high resolution mechanical imaging of GB tissue including perifocal brain parenchyma as a surrogate for infiltrative tumor growth, and distant tumor edema in comparison to tumor and reference tissue. The lack of clearly delineable solid-type tissue within many tumors led us to further specify regions with homogeneous appearance in standard T2-and contrast-enhanced T1-weighted MR images.
By these regions, we address for the first time whether highresolution multifrequency MRE can measure the heterogeneity of GB tissue mechanics which is likely linked to the well known microstructural heterogeneity of GB including neovascularisation and central hemorrhage. Despite the fact that any spatial averaging over-simplifies the tumor's intrinsic heterogeneity, we will tabulate mechanical property values of GB as a starting point for quantification, diagnostic assessment and therapy planning of GB by MRE.

Methods
Twenty-two patients with histologically proven GB (mean age 64.5615.1 years; 10 women) were included in this study. Each patient underwent clinical MRI and MRE prior to further diagnostic (e.g., biopsy) or therapeutic procedures. In addition to MRE, a clinical protocol including T1-, T2-, and proton-densityweighted sequences was applied before and after administration of gadolinium-based contrast agent for further evaluation of the lesion. Imaging slices for MRE were positioned in transverse orientation according to tumor site.
The study was approved by the ethics committee of the Charité -Universitä tsmedizin Berlin (EA1/261/12). All patients gave informed written consent prior to MRE.

MMRE
A custom-designed nonmagnetic driver based on piezoelectric ceramics [33] was mounted at the end of the patient table. The vibrations were transmitted by a carbon fiber rod connected to a custom-designed head cradle located inside the head coil. Eleven of the experiments were performed on a 1.5T MRI scanner (Magnetom Sonata; Siemens Erlangen, Germany), the remaining eleven experiments were carried out on a 3T MRI scanner (Trio; Siemens Erlangen, Germany) using a 4-channel (1.5T) and a 12channel (3T) head coil. The imaging sequence parameters for both systems are listed below. After acquisition of a localizer and a 3D T1-weighted sequence for anatomical reference, a single-shot spinecho echo-planar imaging sequence with trapezoidal flowcompensated motion-encoding gradients (MEG), consecutively applied along all three axes of the scanner coordinate system, was used for rapid motion field acquisition [40]. A custom-made head cradle was used to generate mechanical shear waves inside the brain by inducing a gentle nodding motion of the head. To allow the mechanical waves to propagate into the tissue, the vibration was initiated through a trigger pulse by the scanner at least 100 ms before the start of the MEG. The vibration frequencies (f) used in this experiment were 30,35,40,45,50,55, and 60 Hz. The trigger pulse was delayed in consecutive time-resolved scans by increments of 1/(8*f), yielding 8 dynamics of a wave cycle.

Data Postprocessing
Wave image postprocessing followed the strategy outlined in [42]. In brief: First, the complex MR images were smoothed using a 2D Gaussian filter with a kernel of 5 pixel edge length and sigma = 0.65 for noise reduction. Subsequent gradient-based unwrapping was performed as described by [43]. First-order in-plane derivatives along the image coordinate axes x k (k = 1,2) (x 1 is the phase-encoding direction and x 2 the read-out direction) of the spin phase Q j (j = 1,2,3) were calculated by: Factor j scales the spin phase w j to the physical displacement component u j (in meters) according to [40]. After Fourier transformation in time, equation (1) yields six complex-valued strain images u Ã j,k (v l ) at angular drive frequency v l , resulting in a total of 42 images at each slice invoked by the reconstruction algorithm. These images were further smoothed by a 2D Butterworth lowpass filter with a threshold of 100 m 21 . Low wave numbers as resulted by compression waves were considered sufficiently suppressed by the derivative operators. Other than in previous work [33,38,39], we abandoned curl components for wave inversion since interslice phase artifacts, as addressed by [44], impair the derivative operator in the x 3 direction. Instead of three curl components for each frequency, six independent strain wave images were obtained using eq.(1). All six images were used for stabilizing the wave inversion as described in the following. This strategy is further outlined in [42].
We applied multifrequency dual elasto visco (MDEV) inversion [33,38,39]. This algorithm provides two independent mechanical constants, |G*| and Q, corresponding to the magnitude and phase of the complex shear modulus G*. |G*| provides an indication of the softness or firmness of the tissue while Q provides an indication of the viscous, i.e. dissipative, tissue properties. Of note, both parameters are model-free and provide just another representation of the storage and loss modulus usually parameterized in MRE. It is well known that the loss and storage modulus of brain tissue display dispersion within the examined frequency range [18]. In MDEV inversion-based MRE, we sacrifice the information provided by frequency-resolved complex shear moduli for generating spatially highly resolved maps of |G*| and Q [33]. As a result, |G*| and Q refer to the amplitude and phase angle of the oscillatory response to a harmonic stress, respectively. The effective harmonic frequency of |G*| and Q is given by the mean of all vibration frequencies weighted by the wave amplitudes they produced.
Accounting for complex-valued shear strain images u Ã j,k (v l ) and making the usual assumptions in MRE such as homogeneity, linear viscoelasticity, and isotropy, we obtain the following inversion equations: Given that u Ã j,k represents in-plane strain components and D denotes the 2D-Laplacian, our inversion is entirely 2D-based. By these equations we implemented the method proposed in [39] where data and data derivatives are projected onto the ones vector instead of derivative vector as done in classical least squares solutions of the wave equation [45]. The ones-vector model refers to the almost trivial regression of repeated measurements by computing the observational average (see eqs. 2.81 and 2.82 in [46]).

Morphological Tumor Assessment
MRI-based tumor morphology was classified and graded using T2-weighted (T2w) images as well as contrast-enhanced T1weighted (T1w) images. These images were used to assess tumor morphology including e.g. the presence of cysts, homogenous appearing tumor portion, and necrosis. The tumor volume was calculated using the OsiriX-Imaging software (Geneva, Switzerland) and the MiaLite plugin (SPIE medical imaging 2011, Lake Buena Vista, Florida, USA) by defining contrast enhancing tissue on T1-weighted images. Assessment of tumor morphology and Standard deviations of |G*| and Q were calculated for the tumor ROI to indicate the heterogeneity of the tumor's viscoelastic properties. Additionally, we selected a region of apparently high homogeneity within the tumor region based on morphological MRI (HAM -homogeneous appearing matter) to further compartmentalize the tumor and therewith to address the intra-tumor heterogeneity. In order to study the effect of uncertainties in tumor margins, invasion of surrounding tissue, and partial volume effects, a perifocal ROI was automatically defined by dilatation of the tumor ROI by three pixels minus the tumor region, yielding a small ring around the tumor as illustrated in Figure 1. All regions were single objects delineated in multiple slices. Similarly, tumor volume determination was based on three dimensions, i.e., area analysis was performed on consecutive sections in adjacent slices.

Statistical Analysis
The results are tabulated as arithmetic mean 6 standard deviation. The regional differences between tumor, homogenous appearing matter, perifocal region, edema, and corresponding healthy tissue in |G*| and Q were analyzed by two-tailed paired Student's t-test and Wilcoxon signed-rank test. The Gaussian distribution of the data was tested using the Lilliefors test and the Shaprio-Wilk test. Possible correlations between age and the viscoelastic properties of our regions of interest as well as between tumor volume and viscoelastic properties of the tumor were tested using linear and rank correlation. A P-value ,0.05 was considered statistically significant. All calculations were performed using the MATLAB Statistics Toolbox (MathWorks, Natick, Massachusetts, USA).

Results
The results regarding tumor morphology are summarized in   On average, healthy tissue was significantly stiffer than GB with a mean value of 1.5460.27 kPa and a range of 0.99-2.08 kPa (P = 0.001); however, in a group of 5 tumors, higher |G*| values compared to reference tissue were observed (P = 0.015).
Mean Q was 0.3760.08 for the tumor tissue and 0.5860.07 for the corresponding healthy tissue. Interestingly, this reduction in Q was seen in all tumors (P = 2.9610 210 ) regardless of their |G*| values, which suggests less dissipative (viscous) GB properties compared to healthy tissue.
Within the tumor, homogenous appearing matter showed a higher |G*| compared to full GB regions (P = 0.012) without different appearance to healthy tissue (P = 0.228), suggesting that HAM consists of less affected tissue than the remaining GB. Nevertheless, Q HAM was still lower than Q healthy (P = 0.00013) without significant difference to Q tumor (P = 0.40). This high sensitivity of Q to GB is further represented by the normalized ratios Q tumor /Q healthy which are below 1 in all patients indicating the viability of Q as diagnostic biomarker. No correlation between tumor size and |G*| or phase angle Q was seen with correlation coefficients of R = 20.391 (P = 0.072) and R = 0.101 (P = 0.655), respectively. A correlation between the viscoelastic tissue parameters (|G*|, Q) and morphological tumor staging (Table 1) was also not seen (R = 20.235, P = 0.292 and R = 20.063, P = 0.781 respectively).
|G*| in the perifocal region was not significantly different to |G*| in GB (P = 0.306), whereas Q showed a significant increase in the perifocal region (P = 0.01). The significant correlation between |G*| of tumor and perifocal region (R = 0.571, P = 0.0055) indicates the extension of the tumor's viscoelastic properties into surrounding tissue.
In 16 of 22 GB included in this study, perifocal edema was visible in T2w MRI and could be outlined for MRE parameter analysis. On average, edema tissue was significantly stiffer than GB (P = 0.004), whereas Q was not significantly altered between tumor and edema (P = 0.99). No correlation between |G*| of tumor and edema was observed (R = 0.34, P = 0.197).
Lower values for |G*| healthy at 1.5T (1.3660.21 kPa) than at 3T field strength (1.7160.22 kPa, P = 0.001) were observed, while none of the other parameters given in Table 2 was significantly different between 3T and 1.5T. Specifically, three of the five cases Table 1. Cont. with a higher |G*| tumor than |G*| healthy were measured at 1.5T, which corroborates the independence of MRE parameters from MRI field strength [47].

Discussion
Noninvasive characterization of GB remains a challenge in the present clinical routine. Conventional MRI provides only little information about tissue structure and intraparenchymal tissue connectivity. GB may include cystic, solid, and necrotic fractions as well as diffuse tumor infiltrations of the surrounding tissue; each fraction may alter the mechanical tissue properties measured by MRE.
This manuscript presents the first analysis of viscoelastic constants of intracranial tumors obtained with high spatial resolution MDEV inversion MMRE. To eliminate artifacts resulting from ill-posed inverse problems related to time-harmonic wave patterns and improve MRE parameter maps, we included multifrequency information in the solution of the inverse problem of time-harmonic elastography. The mechanical parameters elucidated in this study are well known in material science and provide full information on the complex shear modulus of human brain tissue. While |G*| relates to our haptic distinction between stiff and soft materials, Q represents the dispersion of the complex modulus, which is dictated by the topology of the underlying cellular network [48]. A highly elastic material such as agarose gel has a low Q value and is thus regarded as less dissipative than biological soft tissues composed of dense and irregular viscoelastic networks including energy-absorbing motile chains. This example illustrates the importance of considering both elastic and viscous terms for characterizing the mechanical properties of a material: agarose gel and biological tissue can have the same elasticity while their distinct viscous behavior may be appreciated by manual palpation.
|G*| and Q are not correlated with each other, and the two parameters convey different and independent mechanical information. In our study, this fact is illustrated by Figure 2, Q was clearly different in GB and healthy brain tissue, while |G*| was lower in only 17 of the 22 tumors. The uniform reduction in the dissipative GB properties may suggest a causal relationship between homeostatic tumor pressure and malignant growth as recently proposed [49] and may in the future be used as a neuroradiological marker of tumor malignancy similar to recent findings in liver tumors [50].
The heterogeneity of |G*| deserves further investigations in GB animal models. Since neither a correlation between |G*| and tumor size was observed (geometry bias) nor system-specific reasons may account for the higher stiffness in five of our patients (three were investigated at 1.5 T and two at 3 T) we expect that |G*| bears potentially valuable information for the characterization of GB. The large variability in morphological tumor assessment scores resulting from the fact that GB may be solid masses or contain cystic and necrotic fractions reflects the potential source of heterogeneity in |G*|. The fact that |G*| is not correlated to the morphological score underlines the novelty of information measured by MRE.
Although encouraging, our study has some limitations: since we conducted a pilot study, we investigated the feasibility of highresolution MMRE in a relatively small group of patients. Future studies should include more patients and compare the findings in different tumor entities. Furthermore, no other mechanical tests could be performed to provide reference values since MRE is unique for the in vivo assessment of tumor consistency. The subjective haptic impression of surgeons in our departments varied widely, preventing us from using their scores as a gold standard of tumor consistency. Future studies in animal models can tackle this issue by using indentation tests or other microelastography methods. Finally, we cannot draw any conclusions regarding the cause of variability and the sensitivity of |G*| to diagnostically relevant tissue changes. This information has to be gathered by MRE in animal models and in a higher number of patients including post-treatment follow-up.
In summary, using multislice MMRE in combination with MDEV inversion enabled us to characterize intracranial tumors  Table 2. The region of HAM is indicated by the dashed line. The yellow arrows indicate compartments of soft tissue properties (low |G*|) but different dissipative behavior (Q) in both tumors. |G*| was scaled from 0 to 3 kPa, Q was scaled from 0 to 2.5 rad. doi:10.1371/journal.pone.0110588.g002 by high-resolution mechanical parameter maps of the human brain. In our cohort of 22 GB patients, the mechanical tissue parameter |G*| indicated that GB are generally softer than healthy tissue, although we noted a large heterogeneity of values. A second mechanical parameter, Q, which is related to the dissipative behavior of tissue, was significantly reduced in all cases. High-resolution MRE may provide an early imaging marker sensitive to pathological changes of mechanical networks in brain tissue. Its diagnostic value, in particular concerning post-treatment follow-up, has to be verified by future studies.