Unipolar MR elastography: Theory, numerical analysis and implementation.

In MR elastography (MRE), zeroth moment balanced motion-encoding gradients (MEGs) are incorporated into MRI sequences to induce a phase shift proportional to the local displacement caused by external actuation. To maximize the signal-to-noise ratio (SNR), fractional encoding is employed, i.e., the MEG duration is reduced below the wave period. Here, gradients encode primarily the velocity of the motion-reducing encoding efficiency. Thus, in GRE-MRE, T2 * decay and motion sensitivity have to be balanced, imposing a lower limit on repetition times (TRs). We propose to use a single trapezoidal gradient, a "unipolar gradient", to directly encode spin displacement. Such gradients cannot be used in conventional sequences as they exhibit a large zeroth moment and dephase magnetization. By time-reversing a spoiled SSFP sequence, the spoiling gradient becomes an efficient unipolar MEG. The proposed "unipolar MRE" technique benefits from this approach in three ways: first, displacement encoding is split over multiple TRs increasing motion sensitivity; second, spoiler and MEG coincide, allowing a reduction in TR; third, motion sensitivity of a typical unipolar lobe is of an order of magnitude higher than a bipolar MEG of equal duration. In this work, motion encoding using unipolar MRE is analyzed using the extended phase graph (EPG) formalism with a periodic motion propagator. As an approximation, the two-transverse TR approximation for diffusion-weighted SSFP is extended to incorporate cyclic motion. A complex encoding efficiency metric is introduced to compare the displacement fields of unipolar and conventional GRE-MRE sequences in both magnitude and phase. The derived theoretical encoding equations are used to characterize the proposed sequence using an extensive parameter study. Unipolar MRE is validated against conventional GRE-MRE in a phantom study showing excellent agreement between measured displacement fields. In addition, unipolar MRE yields significantly increased octahedral shear strain-SNR relative to conventional GRE-MRE and allows for the recovery of high stiffness inclusions, where conventional GRE-MRE fails.

apparent diffusion coefficient

balanced SSFP

contrast‐enhanced‐FAST

displacement encoding with stimulated echoes

double‐echo steady‐state

discrete Fourier transform

diffusion tensor imaging

diffusion‐weighted‐SSFP

echo‐planar imaging

extended phase graph

fractional anisotropy

Fourier‐acquired steady‐state

free induction decay

fast imaging with steady‐state precession

field of view

gradient‐recalled echo

motion‐ encoding gradient

magnetic resonance elastography

magnetic resonance fingerprinting

multi‐slice

octahedral shear strain‐SNR

time‐reversed FISP

region of interest

spin echo

statistical cost, network‐flow algorithm for phase unwrapping

signal‐to‐noise ratio

steady‐state free precession

T2 contrast‐enhanced fast field echo

transistor‐transistor logic.

gyromagnetic ratio

time

wave frequency, angular wave frequency

number of wave phases, ristretto wave, ristretto delay

number of slices, number of averages

repetition time, echo time

T1/T2/T2* relaxation

total acquisition duration

gradient direction matrix

zeroth gradient moment, spatial frequency

resolution/voxel size

real‐valued displacement, displacement phase, complex‐valued displacement

complex encoding efficiency, gradient strength, encoding direction

gradient duration, gradient slope time, time‐offset between TTL and MEG

encoding fraction, fraction of wave period lost to gradient slope

Image motion phase accrual

magnetization vector

complex transverse magnetization

configuration states: transverse and longitudinal

EPG: relaxation operator

EPG: gradient shift operator

EPG: RF pulse operator

EPG: elastography motion propagators for longitudinal and transverse states

dirac delta distribution

INTRODUCTION

Magnetic resonance elastography (MRE) allows the in vivo determination of the local tissue shear stiffness by time‐resolved measurement of the 3D displacement field induced by an external actuator.1 MRE utilizes the principle of phase contrast to encode the cyclic spin motion in the image phase.2 Here, zeroth moment balanced gradients, ie, motion‐encoding gradients (MEGs), are introduced before the readout to sensitize the image phase to displacements of the magnetization occurring while the gradient is switched on. Most MRE experiments currently rely on either spin echo (SE) or gradient‐recalled echo (GRE) sequences.1, 3, 4 Here, SE‐based MRE benefits from the long echo time (TE) that allows maximization of the accumulated phase through matching of the MEG period to the wave period. Optimal motion sensitization with multiple gradient lobes and placement of MEGs before and after the echo pulse is possible.4 However, fast readout strategies including echo‐planar imaging (EPI) or spiral imaging need to be employed at the cost of T2* or frequency offset‐related spatial blurring, eddy current‐induced ghosting, or B0‐inhomogeniety‐induced geometric distortions.5 In GRE‐based MRE, a short Cartesian readout can be used due to the much shorter repetition time (TR). However, the need for a shorter MEG and reduced SNR loss due to T2* decay has to be addressed, which leads to the use of fractional encoding, where the duration of the MEG is shorter than the wave period.6, 7, 8 In addition, the loss of longitudinal magnetization due to the short TR as well as the use of RF‐spoiling to reduce spin‐history effects leads to additional reduction in SNR.9, 10 Recent developments in MRE have been driven by the shortcomings of the SE‐ and GRE‐based sequence concepts and have focused on new transducers,11 sequence timing,7, 12 encoding schemes,13 encoding approaches,14, 15, 16, 17 as well as excitation and readout schemes.15, 18, 19, 20

Instead of spoiling magnetization after the readout, refocusing gradients can be used to effectively reuse transverse magnetization in subsequent repetition periods. This balanced steady‐state free precession (bSSFP) approach has been exploited for MRE, which greatly benefits from its high SNR.6, 21, 22 However, the continuous refocusing of transverse magnetization in bSSFP is dependent on phase accrual during the TR interval, which leads to the typical banding artifacts in the presence of off‐resonance. In the context of MRE, the interference of motion phase and inhomogeneity leads to nonlinearity in the encoding of the vibratory motion, rendering bSSFP‐based MRE difficult for in vivo applications.6, 21, 22 Nonetheless, bSSFP‐based MRE has been reported for hepatic and cardiac MRE.6, 23

In addition to spoiled and balanced SSFP sequences, a generalized type of GRE sequence has recently regained interest for rapid relaxometry and diffusion‐weighted SSFP (DW‐SSFP) imaging.24, 25, 26, 27 These methods exploit additional refocusing gradients with integer multiples of the spoiling area before the readout of a conventional GRE sequence with consecutive balancing after the readout to refocus previously spoiled magnetization states while keeping the total spoiling area constant. These spoiling states, termed “configuration states”, are carried over to the next TR period, where the subsequent RF pulse then partly flips the magnetization back to the longitudinal axis (storage), leaves it unaffected (dephasing configuration), or inverts it (rephasing configuration).28 This leads to a plethora of discrete dephasing and rephasing configuration states, which can be refocused by a gradient that counteracts their accumulated spoiling.

This work proposes a new MRE sequence, which is based on the readout of the SSFP SE using a time‐reversed spoiled SSFP sequence.29 Here, the spoiling gradient acts as a highly efficient MEG, which allows considerable reduction in TRs as well as encoding fractions. In addition, 3D slab encoding leads to increased SNR compared with conventional consecutive or interleaved multi‐slice (MS) acquisitions.30 In contrast to conventional GRE‐ or SE‐based MRE sequences, which can be solely understood on the basis of phase contrast, the proposed technique depends on spin history effects and necessitates the analysis of echo pathways and their motion sensitivities. The objective of this work is 5‐fold: first, a qualitative theoretical understanding of echo formation and motion encoding in unipolar MRE is provided; second, a periodic motion propagator in the presence of gradient fields is derived for the extended phase graph (EPG) formalism, which allows for the numerical, quantitative analysis of the motion encoding in unipolar MRE using EPGs; third, a complex encoding efficiency is proposed, which takes phase differences in encoding into account and hence allows direct comparison of the complex displacement field in both magnitude and phase of conventional and unipolar MRE despite their vast difference in encoding concept and sequence timing; fourth, a simplification of the full EPG‐based simulation is attempted through adaptation of the two‐transverse TR approximation, which was originally derived for DW‐SSFP, allowing for fast encoding efficiency estimation27, 31; and fifth, a comprehensive phantom study comparing unipolar MRE and conventional GRE‐MRE is provided in order to demonstrate the feasibility of unipolar MRE and to validate theoretical findings.

THEORY

In GRE‐based MRE, a spoiling gradient is played after the readout of the free induction decay (FID) to dephase transverse magnetization (Figure 1A,C). We propose to read out the F−1 configuration state (SSFP‐echo) instead, which is refocused by time‐reversing the spoiled GRE sequence (Figure 1B), typically referred to as contrast‐enhanced‐FAST (CE‐FAST), PSIF (Siemens; reversed FISP) or T2‐CE‐FFE (Philips).34, 35 Contrary to conventional GRE sequences, the readout in the first TR period does not acquire any signal, as the spoiling gradient fully dephases the transverse magnetization. The second RF pulse partly leaves the transverse magnetization unaffected (0° component, dephasing configuration) and partly inverts it (180° component, rephasing configuration). The latter is refocused by the spoiling gradient and forms the SSFP‐echo.32 Since the spoiling gradient is played before the readout, the sequence becomes highly motion‐sensitive, which is typically regarded as a drawback of time‐reversed spoiled SSFP.24, 25, 27, 36, 37, 38, 39, 40 Here, this increased motion sensitivity is used for highly efficient elastography acquisitions. Since the bipolar MEG is replaced by a single gradient lobe, we will refer to the sequence as unipolar MRE.

Figure 1

Depiction of echo formation in (A) spoiled (GRE) and (B) time‐reversed spoiled SSFP. While the readout in spoiled GRE sequences is dominated by the FID, the first readout in time‐reversed SSFP does not acquire any signal as the gradient spoils the signal of the FID. The second RF pulse splits the magnetization into a rephasing and a dephasing contribution. The former is refocused by the spoiling gradient and leads to the formation of the SSFP‐echo (figure adapted from Scheffler, Ref. 32). Comparison of (C) conventional GRE‐MRE with fractional bipolar motion‐encoding gradients (MEGs) and (D) the proposed unipolar MRE sequence that uses the spoiling gradient as the MEG and reads out the SSFP‐echo instead. The gray‐shaded gradients are zeroth‐moment balanced over one TR period. In (E), the global sequence timing of unipolar MRE is shown, which follows the 3D Ristretto scheme.12, 33 The innermost squares denote acquisitions according to (D), which are wave‐phase interleaved (NP = 5, ND = 1). After a set of dummy shots to reach steady‐state, successive k‐lines are acquired in both slice‐ and phase‐encoding direction (3D, slab excitation). After the full volume is acquired, the encoding direction is changed according to the Hadamard encoding scheme. A combination of TTL triggers and a correct choice of TR is used to ensure synchronization of the external wave generator and the unipolar MRE sequence

In order to acquire wave‐phase offsets (with ) to extract amplitude and phase of the displacement field, the sequence timing follows the Ristretto MRE scheme33 for 3D encoding (Figure 1E). Contrary to conventional 3D GRE‐MRE schemes,41 Ristretto interleaves the wave‐phase acquisition allowing minimization of the TR. The wave phases are acquired in an interleaved fashion in the innermost sequence loop. The repetition time is constant and given by with being a co‐prime of ( ). Each successive repetition is thus offset relative to the transistor‐transistor logic (TTL) trigger played every repetitions, where the wave phase offsets amount to Since RF pulse and net dephasing from imaging gradients is equal in each TR and the underlying motion is periodic in , a steady‐state of distinct echoes forms.42

In the proposed sequence, the unipolar gradients for spoiling and motion sensitization coincide, necessitating the switching of the spoiling direction during sequence acquisition. To correct for phase contributions from imaging gradients, which are independent of the encoding direction, a reference phase image without motion sensitization is typically acquired in GRE‐MRE.7, 8 However, since motion sensitization and spoiling is facilitated with the same gradient, an encoding scheme that encodes with equal total gradient moment for all encoding directions needs to be employed to ensure equal image contrast as well as sufficient spoiling. To this end, Hadamard encoding is incorporated, which encodes along the four diagonals of a regular cube.8, 43 The gradient direction matrix is given by where the rows correspond to the three orthogonal gradient axes, while the columns denote successive encoding directions. Since the Hadamard scheme exhibits identical gradient moments in each of the three orthogonal gradient axes, isotropic image resolution is required to fulfill the gradient spoiling condition30 (p. 353) for each spatial dimension separately. To ensure proper spoiling of the SSFP signal, the spoiling direction is changed only after the full k‐space for all wave phase offsets has been acquired and additional dummy shots are played for each encoding direction separately to reestablish the steady‐state (see Figure 1E).

Lowest order encoding efficiency of unipolar MRE: SSFP spin‐echo contribution

In phase contrast‐based MRE, the phase accrual in each voxel due to the periodic motion is proportional to the time‐varying displacement field Here, is the complex displacement amplitude, the imaginary part, and the angular wave frequency. The proportionality factor is denoted by the encoding efficiency ϵ of the gradient waveform, which is given by3 where is the gyromagnetic ratio and denotes the gradient waveform. While the phase of the complex encoding efficiency can often be neglected in conventional MRE, we specifically account for it to allow for absolute displacement field comparisons of conventional and unipolar MRE.

The sampling of the equidistant phase offsets along the m‐th direction can be stated as where accounts for the phase accumulated by imaging gradients and the notation denotes the mth gradient direction. Thus the projection of the estimated displacement including imaging gradient bias along the mth encoding direction is obtained by inverse discrete Fourier transformation which can then be solved for the unbiased using the pseudo‐inverse of the encoding matrix.8

Bipolar MEG

The complex encoding efficiency for a bipolar MEG (Figure 2A), with duration , gradient strength , and slope time played at relative to the TTL‐trigger is given by

Figure 2

Timing graphs for (A) conventional GRE‐MRE, (B) unipolar MRE, as well as the equivalent phase‐contrast sequences (C)‐(E) of unipolar MRE, where the first unipolar MEG is flipped to account for the complex conjugation of the echo signal through the second RF pulse. The phase‐contrast model is equivalent to the lowest order echo pathway contributing to the unipolar MRE signal, which we term the “SSFP‐SE” here. For comparability of the acquired displacement fields, absolute timing of the MEGs needs to be taken into account, where T0 is the time offset between MEG and transistor‐transistor logic (TTL) trigger, T is the total duration of the MEG, td is the gradient slope time and TR is the repetition time. The reference time‐point to which the displacement field is calculated is set to coincide with the second RF pulse for unipolar MRE and the RF‐pulse in conventional MRE, respectively. Approximations in (D) finite slew and (E) infinitesimal gradients are included in the analysis for better understanding of the encoding principle of unipolar MRE

Infinite‐slew unipolar MRE

The encoding efficiency of a single unipolar lobe (Figure 2D) with constant gradient strength , infinite slew, and finite duration can be obtained by solving Equation 6 and can be written as where denotes the encoding fraction,6 is the gradient area of the unipolar lobe, and is the sinc‐function.

The primary contribution to the SSFP‐echo is given by the spin echo (SSFP‐SE) of two successive TR periods, which is sensitized to motion twice (Figure 2B). The part of transverse magnetization, which is refocused, is complex conjugated at the time of the second RF pulse (Figure 2B). Thus, the motion sensitization of the first unipolar gradient is effectively inverted (Figure 2C). To obtain the correct encoding phase, the reference time with respect to the TTL trigger is set to be equal to the second RF pulse, thus the initial magnetization is excited at time . The lowest order approximation to the encoding efficiency of unipolar MRE assuming infinite slew is then given by where the timing of the Ristretto scheme as specified in Equation 1 is implied. In the limit of an infinitesimal encoding gradient with finite moment (Figure 2E), the encoding efficiency does not vanish but approaches

Finite‐slew unipolar MRE

Generalization of the previous infinite‐slew model to the real situation of a trapezoidal gradient lobe with finite slope time (Figure 2C) leads to where is the fraction of the wave period that is lost to the gradient slewing. Taking the limit , the infinite‐slew encoding efficiency Equation 11 is directly recovered. By replacing and as well as temporal shifting by , the encoding efficiency of the bipolar MEG with finite slew is obtained.

The above encoding efficiency Equation 13 is only valid under the assumption that each received signal is solely given by the SSFP‐SE contribution, which is excited and refocused by the current and previous RF pulses. However, since no recovery period is included in time‐reversed spoiled SSFP sequences compared with conventional SE sequences and the sequence TR is usually much shorter than T2, spin history effects need to be taken into account. These will affect the encoding efficiency by constructive and destructive interference.

EPG formalism for unipolar MRE

In order to describe dephased magnetization states, EPG formalism is used. A comprehensive introduction to EPGs is given in the review article by Weigel,28 which we follow here in terms of notation and conventions; and which should be consulted for explicit mathematical expressions of the operators. EPGs describe spoiled magnetization as a superposition of complex configuration states and . There, the central assumption is uniformity of space so that spatial dependencies in the magnetization distribution are solely arising due to gradient fields. Each configuration state corresponds to a Fourier component of the transverse and longitudinal magnetization distribution. The formal definitions are RF pulses, gradient fields, motion and relaxation are described in terms of operators acting on these configuration states. RF pulses distribute magnetization between transverse and longitudinal configuration states independent of the spatial frequency , where denotes the flip angle and the pulse phase. Relaxation leads to attenuation of transverse and longitudinal states and the recovery of longitudinal magnetization in the configuration state only. Gradients shift the population of transverse configuration states by the respective ‐value of the gradient. Finally, motion and diffusion can be incorporated via motion propagators , which are typically k‐value dependent. Bulk motion leads to phase accrual (Fourier shift theorem) and random motion to attenuation (phase dispersion). A pulse sequence is described as a concatenation of these operators.

In Figure 3, EPGs for unipolar MRE in the presence and absence of periodic motion are shown. Since the spoiling area is constant for all imaging shots, the y‐axis counts the spatial frequencies in units of the spoiling area in the direction of the spoiling gradient. The unipolar spoiling gradient, which is played before the readout in the time‐reversed spoiled SSFP sequence, leads to dephasing of configuration states with and to respective rephasing of configurations with . Due to the presence of periodic motion, the configuration states accumulate a phase depending on their k‐value. This is depicted by the bending of the configuration states, which is proportional to the relative phase accrual with respect to the previous RF pulse. This should not be confused by an actual change in spatial frequency of the respective states.

Figure 3

(A) Pictorial description of the extended phase graph in steady‐state for the case of three encoded wave phases (NP = 3) and a delay of one (ND = 1). Solid lines denote transverse magnetization states, whereas broken lines depict longitudinal magnetization. The overlaid thin solid and dashed lines as well as the graph in (B) denote the EPG in the absence of motion. The unipolar MRE sequence is shown in abbreviated form in (D). The change in configuration state phase in the presence of periodic motion is depicted by the bending of the EPG. Here, the change in phase is depicted relative to the time of the RF pulse. Motion sensitivity of the configuration states increases linearly with|k|resulting in nonlinear encoding of displacements if configuration states are mixed. The bending of the EPG is for visualization purposes only and should not be confused with a change in k. In (C), diffusion attenuation is shown, which attenuates high spatial frequency states. This is also depicted in the EPG by the lighter color of the graph. Diffusion and T2‐decay limit the mixing of different configuration states reducing overall sensitivity to relaxation and sequence parameters. (E) shows the shortened EPG propagator ordering, which can be used to simulate the steady‐state magnetization

In addition to periodic motion, diffusion attenuates configuration states with high spatial frequency.28 This is illustrated 2‐fold: first, by reduced line thickness and opacity in the EPG with increasing and, second, by the diffusion attenuation plot to the right. At the bottom of the graph, the unipolar MRE sequence is denoted together with its simplified operator notation.

In the present work, two EPG‐based analyses are proposed: (i) the EPG representation of unipolar MRE can be used for full numerical encoding efficiency calculations by driving a forward calculation according to the operator sequence denoted in Figure 3 to steady‐state; and (ii) a reverse analysis of the echo signal contributions can be conducted, describing the lowest order perturbation theory solution of the EPG in terms of T2‐decay. This approach is commonly known as the two‐transverse TR‐approximation, which has been proposed before for the understanding of diffusion contrast in DW‐SSFP sequences, and which is extended here to incorporate periodic motion.31, 44

Motion propagator for unipolar MRE

The motion propagator for coherent bulk motion subject to an arbitrary k‐space trajectory is given by28 with the usual definition of Here, is the initial wave vector of the configuration state. Assuming periodic motion (Equation 5), the general motion propagator for a single constant gradient lobe with infinite slew and finite duration is given by and where is again the area of the gradient lobe. The propagator for longitudinal magnetization as well as transverse magnetization in the absence of a gradient field can be found by setting to zero and is thus given by only. Since a unipolar lobe causes a change in configuration state population, the full unipolar propagator is given by a concatenation of motion sensitization and subsequent shifting of the configuration state population according to the gradient operator . The motion propagator for the general case of a trapezoidal lobe with finite slew is stated in Appendix A (see the supporting information).

In addition to motion sensitization, the dephased magnetization states are subject to diffusion attenuation.28 The 1D diffusion attenuation operator is given by Thus, for transverse magnetization in the case of an infinite‐slew gradient lobe, the operator becomes with Again, is the propagator in the absence of the gradient as well as the corresponding propagator for longitudinal magnetization states. The diffusion attenuation can be incorporated in the general motion propagator for unipolar MRE by applying the replacement and , respectively.

Diffusion attenuation is an important factor in unipolar MRE, as it leads to attenuation of high spatial frequency states since the exponential scales proportional to . This ultimately reduces the sensitivity of the sequence to T1 and T2 relaxation, flip angle, as well as leading to the reduction of the nonlinearity of motion encoding.

Extending the two‐transverse TR approximation in the presence of periodic motion

The two‐transverse TR approximation is an approximate solution to the problem of diffusion weighting in SSFP‐based sequences.31, 44 It determines the signal based on the assumption that only echoes with no more than two transverse TR‐periods contribute to the signal amplitude of the SSFP‐echo. This corresponds to the lowest order Taylor series expansion of the analytical solution in .

In Figure 4, the EPG is drawn in reverse, ie, starting from an echo and tracing back the signal contributions in lowest order perturbation theory in . This leads to the graphical depiction of the three distinct echo pathways of the two‐transverse TR approximation. Solid paths are transverse magnetization states, whereas dotted lines denote longitudinal magnetization. Each path is preceded by infinitely many longitudinal magnetization‐only paths to account for the equilibrium magnetization. Thus, every pathway can be characterized by the number m of longitudinal storage periods. Transverse magnetization is subject to motion sensitization according to the transverse motion propagators and . corresponds to the motion sensitization of the refocusing pathway before readout, while determines motion sensitization after excitation. Longitudinal storage in elevated configuration states (here, only contributes) are modulated according to , respectively. The definition of the motion propagators follows from careful evaluation of Equation 18 under consideration of the number and timing of the transverse and longitudinal echo pathways. The respective terms are given in Appendix B (see the supporting information).

Figure 4

Pictorial description of the two‐transverse TR approximation for NP = 5 and ND = 1 tracing back the origin of the echo. Equilibrium magnetization undergoes an infinite series of longitudinal storage periods before the first flip into the transverse plane. In the stimulated echo case, m > 0 longitudinal storage periods are assumed before flipping the magnetization back into the transverse plane, where it refocuses to form an echo. For m = 0, no longitudinal storage is found and the magnetization is directly flipped into its conjugate state forming the spin‐echo of the spoiled SSFP sequence. The unipolar MRE sequence is depicted in the last row, where the arrows denote RF pulses, trapezoids represent the spoiling gradients and circles denote echo time‐points

Under consideration of the periodicity of the motion propagators the derivation of the two‐transverse TR approximation can be carried out by following the derivation of McNab and Miller,44 which is laid out explicitly in Appendix C (see the supporting information). In conclusion, the echo signal is given by where the SE and STE contributions are found to be with and The final step, determination of the motion phase accumulated in steady‐state, is performed by Fourier transformation of the argument of the complex steady‐state signal (Equation 24) according to Here, the integration is normalized such that the Fourier transformation reproduces amplitude and phase of a sinusoidal oscillation according to To the knowledge of the authors, no closed‐form solution of Equation 29 exists. Hence, the final solution is obtained numerically using MATLAB.

METHODS

EPG simulations

The EPG simulations were based on the Stanford EDU MATLAB implementation45 and run in MATLAB 2018a (MathWorks, Natick, Massachusetts). The unipolar motion propagator was implemented according to Equation 16 including diffusion attenuation and using the full trapezoidal gradient waveform. The steady‐state signal behavior was simulated by repeated application of the EPG formalism to drive an initial equilibrium magnetization configuration to steady‐state. The propagator ordering is given in abbreviated form in Figure 3 and was evaluated using the appropriate timing constraints given by the Ristretto MRE scheme (Equation 1). The analytical two‐transverse TR approximation model was numerically solved using MATLAB as well.

In order to determine the encoding efficiency in the steady‐state, the unwrapped signal phase was obtained for the last elements, which was then reordered according to Equation 2 and subsequently Fourier‐transformed. The first Fourier component—corresponding to the actuation frequency—was divided by the input wave amplitude yielding the encoding efficiency. If not stated otherwise, all calculations were performed with the following set of parameters to be compliant with the phantom experiment: T1/T2 = 2400/190 ms, diffusion coefficient D = 1.7 x 10−3 mm2/s, number of wave phases NP = 5, Ristretto delay parameter ND = 1, wave frequency f = 35 Hz, gradient duration T = 2 ms (encoding fraction q = 7%), gradient strength G = 20.9 mT/m, gradient slope time td = 0.13 ms, displacement amplitude d = 10 μm and flip angle = 20°, resulting in a sequence TR of 5.71 ms. EPG simulations were run with a fixed number of 700 EPG states to account for spin history and were iterated to steady‐state with 500 iterations of the phase loop corresponding to 500 x NP = 2500 TR periods.

Phantom comparison

Unipolar MRE as well as conventional GRE‐MRE were performed on a gel phantom with four cylindrical inclusions (CIRS Inc., Norfolk, Virginia) using a 3 T Philips Ingenia system (Philips Healthcare, Best, the Netherlands). The signal was received using a 15‐channel head coil, and 35 Hz electromagnetic actuation was employed.4 The phantom was positioned upright in the head coil with the actuator placed on top with primary actuation being in the top‐down direction. Slice selection direction was chosen anterior–posterior (see also Figure 8). In the following, scan parameters are stated in pairs of “conventional/unipolar MRE”; 20/10 repeat measurements were acquired to investigate the dependence of the shear stiffness estimates with increasing SNR. All measurements were performed using fractional encoding6 with encoding fractions of 17.5%/7%, 11.9/20.9 mT/m bipolar/unipolar MEGs (duty cycle limitation in conventional MRE) resulting in an encoding efficiency of 4.01/10.8 rad/mm, 3 mm3 isotropic resolution and 64 × 48 matrix size acquired using a 30°/20° flip angle, Cartesian readout, and third/first in‐phase water‐fat echo‐time (TE 6.9/2.3 ms). Conventional GRE‐MRE was performed as a MS acquisition of 12 interleaved slices resulting in a TR per slice of 110.7 ms, while unipolar MRE was performed using slab‐excitation and phase encoding of 36 slices with a TR of 5.71 ms. The increased FOV in slice direction was chosen to eliminate fold‐over artifacts due to production impurities with strong susceptibility mismatch in the bottom of the phantom. The conventional GRE‐MRE sequence was accelerated using the Ristretto MRE12, 33 scheme acquiring eight wave‐phase offsets in reverse order (ND = 7) and excitation and readout of the 12 slices for one wave‐phase offset within NW = 3 wave periods. Only five wave‐phase offsets were acquired with unipolar MRE to conform with sequence‐timing constraints. The spoiling gradient was chosen to result in 5 × 2 phase accumulation along each spatial direction, separately. The acquisition of all wave‐phase offsets and encoding directions for a single slice took 14.3/7.0 seconds, respectively. Synchronization to the external wave generator was achieved using a TTL trigger signal. In order to achieve a magnetization and wave steady‐state, 62 dummy shots were used for unipolar MRE, totaling to 1.8 seconds per encoding direction (conventional MRE: 10 dummies, 1.1 seconds/encoding direction).

Figure 8

(A) Comparison of the magnitude images of conventional GRE‐ and unipolar MRE at 35 and 70 Hz frequency for one phase‐offset, encoding direction and slice. Compared with conventional GRE‐MRE, pronounced intra‐voxel phase‐dispersion and ‐dependent contrast is visible. (B) Comparison of the real part of the complex shear displacement field in through‐plane direction with and without averaging. The displacement fields are qualitatively and quantitatively in good agreement. At 70 Hz, unipolar MRE shows slightly reduced wave displacements throughout the phantom compared with conventional MRE. Conventional MRE shows considerable displacement noise in the single scan compared with both unipolar MRE and averaged conventional MRE, respectively. The top center sketch depicts the phantom from a side view including the alignment of the imaging plane as well as the actuator placement

A second set of phantom scans was performed at 70 Hz actuation frequency. Unipolar MRE with five wave‐phase offsets was acquired with ND = 2 (36 slices), while conventional MRE was acquired with NW = 7 and ND = 3 (12 slices). Due to the shorter wavelength and thus better conditioned inversion at 70 Hz, 10 averages were acquired for both techniques. Remaining scan parameters were equivalent to the 35 Hz acquisition, leading to an encoding fraction of 35%/14%, an encoding efficiency of 7.5/10.2 rad/mm, and a scan duration of 13.6/7.0 seconds per slice and average, for conventional/unipolar MRE, respectively.

In addition to the acquisition of the MRE datasets, quantitative parameter mapping was performed to obtain local estimates for T1, T2, the apparent diffusion coefficient (ADC) and the fractional anisotropy (FA) of the phantom in one representative slice in order to calculate the local, complex encoding efficiency for the unipolar MRE acquisition. T1 and T2 estimates were obtained on the same 3 T Philips Ingenia system using a spoiled GRE‐based MR fingerprinting (FISP‐MRF)46 sequence with a constant TR of 16 ms, fully sampled multi‐shot spiral acquisition with 15 interleaves, a repetition delay of six seconds to ensure relaxation to equilibrium, and 1000 variable excitations after adiabatic inversion resulting in a total scan duration of 5.5 minutes. The flip angle train is given in Figure S1. Dictionary‐based matching was performed to determine T1 and T2 locally using MATLAB 2018a (MathWorks). The dictionary was created based on the Stanford EDU MATLAB EPG implementation.45 The dictionary resolution is given in Table S1.

SE‐based diffusion tensor imaging (DTI) with unipolar encoding gradients, six encoding directions and five b‐values (linearly interpolated, 0 to 1000 ms/mm2) was performed to obtain an estimate of the ADC and FA in a single slice.47, 48 The scan was performed on a 1.5 T Philips Achieva system (Philips Healthcare). The signal was received using an eight‐channel head coil. Sequence parameters were 3 x 3 mm2 in‐plane resolution, 8 mm slice thickness, 10 averages, single‐shot EPI readout, matrix size 64 x 45, FOV 196 x 136, flip angle 90°/180°, 60 ms echo time, and a TR of 1000 ms, resulting in 4 minutes 11 seconds scan duration.

For comparability, the FOV of the unipolar MRE scan was cropped to the 12 central slices to match the conventional MRE scan prior to postprocessing. The image phase was spatially unwrapped using the statistical cost, network‐flow algorithm for phase unwrapping (SNAPHU),49, 50, 51, 52 temporally aligned, and decoded using the pseudo‐inverse of the encoding matrix.8 Cumulative averaging of the decoded wave phase images was performed prior to Fourier‐transformation. The complex displacement fields were estimated from the first Fourier component and corrected for the encoding efficiency of the conventional and unipolar MRE scans. The phase of the complex displacement field of the conventional MRE scans was corrected for inter‐slice phase shifts according to the Ristretto MRE scheme.12, 33 The encoding efficiency of the unipolar MRE acquisitions was determined using the two‐transverse TR approximation assuming homogenous relaxation and diffusion (T1 = 2.4 s, T2 = 190 ms and D = 1.7 x 10−3 mm2/s). The complex displacement fields were then further processed without smoothing using a FEM‐based reconstruction algorithm53 to obtain the complex shear stiffness, where we report the averaged magnitude over the innermost six slices. In addition, the wave‐phase images were used to estimate the displacement noise in each voxel, allowing for the determination of the octahedral shear strain‐SNR (OSS‐SNR).54 OSS‐SNR values are reported for the inclusions and for the background as ROI‐averaged values.

RESULTS

In Figure 5, the encoding efficiency of unipolar MRE is shown as a function of tissue parameters T1, T2, and the diffusion coefficient within the phase‐contrast approximation (red, SSFP‐SE contribution), the two‐transverse TR approximation (yellow), and the full EPG simulation (violet). If applicable, full EPG simulations neglecting diffusion attenuation are shown in dotted violet. Sequence parameters were held constant at the values displayed by the vertical dashed lines. In Figure 5A, the encoding efficiency is shown as a function of the diffusion coefficient. The phase‐contrast model is independent of diffusion by construction. For typical diffusion coefficients in vivo (D = ~2 x 10−3 mm2/s), the encoding efficiency obtained from full EPG simulations is within 3% of the two‐transverse TR approximation and converges for strong diffusion. Neglecting diffusion, the encoding efficiency obtained from full EPG simulations deviates by ~20%.

Figure 5

Encoding efficiency comparison of unipolar MRE in the simplified phase‐contrast approximation (red), the two‐transverse TR approximation (yellow), and the full EPG simulation (violet). The dashed line denotes the default parameter assumed in all simulations, ie, f = 35 Hz, NP = 5, ND = 1 (TR = 5.71 ms), q = 7% and G = 20 mT/m, if not otherwise stated. Dependencies are plotted as a function of (A) diffusion coefficient, (B) T1 relaxation for long T2 (190 ms), (C) T2 relaxation for long T1 (2400 ms) and (D) T1 relaxation for short T2 (30 ms). In (A), the approach of the full EPG simulation towards the two‐transverse TR approximation for strong diffusion can be appreciated. For (B)‐(D), full EPG simulations for unipolar MRE neglecting diffusion (dotted violet) are shown to demonstrate deviation of the encoding efficiency due to interference of high order, motion‐sensitive configuration states. Compared with full EPG simulations, the two‐transverse TR approximation can account for more than 90% of the encoding efficiency in the presence of diffusion

In Figure 5B‐D, the encoding efficiency is shown as a function of T1 and T2 under the assumption of finite diffusion (D = 1.7 x 10−3 mm2/s). Figure 5B shows the case of long T2 relaxation (190 ms), where full EPG and the two‐transverse TR approximation show little dependency on T1 with variations below 3%. Figure 5D shows the case of short T2 relaxation (30 ms). For long T1 (>1000 ms), full EPG with and without diffusion as well as the two‐transverse TR approximation become equivalent. At very short T1 (<100 ms), maximal deviations of up to 5% are found. Contrary to the long T2 case, diffusion only plays a minor role, as shown by the near equivalence of full EPG with and without diffusion. This can also be appreciated in Figure 5C, which shows the T2 dependency for the case of long T1 (2400 ms). Below 50 ms, the two‐transverse TR approximation and the full EPG with and without diffusion are closely aligned (<1% deviation). For increasing T2, full EPG simulations, especially without diffusion attenuation, show decreasing encoding efficiency compared with the two‐transverse TR approximation. Again, the encoding efficiency of both phase‐contrast and the two‐transverse TR approximation are independent of T2 by design. The dependency of the encoding efficiency phase angle on the same set of parameters is shown in Figure S3.

In Figure 6, the encoding efficiency of unipolar MRE is compared for different sequence parameters. Figure 6A shows the encoding efficiency as a function of encoding fraction q and compares the different unipolar MRE approximations against conventional phase‐contrast MRE with bipolar MEGs (blue). The red line indicates the encoding fraction limit imposed by the sequence timing (here NP = 5, ND = 1), which limits the duration of the encoding gradient. In the small q regime (<10%), the encoding efficiency of unipolar MRE is proportional to q, whereas conventional bipolar MRE is proportional to q2. In Figure 6B, the dependency on the sequence‐timing parameters NP and ND is shown. The SSFP‐SE contribution shows its maximal encoding efficiency at NP = 2, where the phase contributions from two consecutive TRs interfere constructively. For increasing ratios, the SSFP‐SE contribution shows declining encoding efficiency. The two‐transverse TR approximation predicts approximately constant behavior, whereas the full EPG simulation shows increasing encoding efficiency with increasing NP/ND ratio above two. The behavior of the encoding efficiency is reversed if NP/ND ratios below five are considered. Figure 6C investigates the dependency of the magnitude of the complex encoding efficiency as a function of wave displacement. Both the spin echo and the two‐transverse TR approximation show constant encoding efficiency, thus displacements are encoded linearly in the signal phase. The full EPG simulation shows nonlinear behavior, especially for displacements above 15 μm. The inset depicts the maximal signal phase for the same range of displacements from 0 to 100 μm, again showing the linearity of encoding in the SSFP‐SE and two‐transverse TR approximation. The full EPG solution is nonlinear but monotonously increasing with displacement. Finally, the dependency of the encoding efficiency of unipolar MRE as a function of the flip angle is investigated in Figure 6D. With increasing flip angle, both full EPG and the two‐transverse TR approximation show an increase in encoding efficiency, while the SSFP‐SE contribution remains constant. The dependency of the encoding efficiency phase angle on the imaging parameters is shown in Figure S4.

Figure 6

Encoding efficiency comparison of unipolar MRE in the simplified phase‐contrast approximation (red), the two‐transverse TR approximation (yellow) and the full EPG simulation (violet). The dashed line denotes the default parameter assumed in all simulations, ie, T1 = 2.4 s, T2 = 190 ms, D = 1.7 x 10−3 mm2/s, f = 35 Hz and G = 20 mT/m. dependencies are plotted for (A) encoding fraction, (B) number phase offsets and delay, (C) displacement and (D) flip angle. In (a), conventional MRE with bipolar MEGs is shown in blue as well as the full EPG simulation neglecting diffusion (dotted violet). Unipolar MRE shows a linear dependency on the encoding efficiency as displacements are directly encoded, whereas bipolar MEGs encode the velocity and respectively scale proportional to q2. The refocusing of higher order configuration states leads to pronounced interference, which results in the dependency on (B) sequence timing as well as nonlinear displacement encoding (C). The inset in (C) shows the magnitude of the product of displacement and encoding efficiency. Deviations from the two‐transverse TR approximation are minor in the small displacement domain (<20 μm) showing that encoding is dominated by refocused 1st order contributions

Phantom results

In Figure 7, the results of the quantitative parameter mapping for (A) T1 and T2, and (B) ADC and FA are displayed. T1 is found to be approximately constant over the whole phantom (2425 ± 33 ms), while T2 differs significantly between background (193 ± 25 ms) and inclusions (632 ± 33 ms). The ADC is found to be approximately constant over the whole phantom ([1.68 ± 0.31] x 10−3 mm2/s). FA is close to zero, consistent with the assumption of isotropic diffusion. ROI average and standard deviation for the inclusions, background and the whole phantom are given in Table S2. In Figure 7C, the magnitude and phase of the local encoding efficiency are shown based on the determined T1 and T2 map; however, diffusion effects are neglected. Substantial encoding efficiency differences between background (|ε| = 8.77 ± 0.30 rad/mm, ∠ε = 0.192 ± 0.008 rad) and inclusions (|ε| = 5.76 ± 0.13 rad/mm, ∠ε = 0.276 ± 0.009 rad) of ~52% can be observed. Figure 7D shows the encoding efficiency taking a constant diffusion coefficient of 1.7 x 10−3 mm2/s into account. Here, the encoding efficiency of the background (|ε| = 10.605 ± 0.024 rad/mm, ∠ε = 0.006 ± 0.013 rad) and inclusions (|ε| = 10.412 ± 0.007 rad/mm, ∠ε = −0.095 ± 0.005 rad) deviates by less than 2%, allowing treatment of the encoding efficiency as homogenous over the phantom. ROI average and standard deviation of magnitude and phase of the encoding efficiency, including diffusion effects for the inclusions, background and the whole phantom, are given in Table S2.

Figure 7

Quantitative parameter maps for (A) T1 and T2, (B) ADC and FA. T1 and ADC are approximately constant over the phantom and are 2425 ± 33 ms and (1.68 ± 0.31) x 10−3 mm2/s, respectively. FA is 0.042 ± 0.022 over the whole phantom, which is consistent with the assumption of isotropic diffusion. T2 of inclusions is elevated (632 ± 33 ms) compared with the background (193 ± 25 ms). In (C), magnitude and phase of the encoding efficiency is shown as a map neglecting diffusion but taking T1/T2 from (A) into account. In (D), the encoding efficiency incorporating constant diffusion (1.7 x 10−3 mm2/s) is shown. Compared with (C), differences in encoding efficiency including diffusion effects are below 2%

In Figure 8, (A) signal magnitude and (B) the real‐part of the complex displacement field for one representative slice and through‐plane direction are shown for both unipolar and conventional MRE at 35 and 70 Hz actuation frequency. The top row corresponds to the magnitude images obtained from a single phase‐offset, slice and encoding direction without averaging. Pronounced intra‐voxel phase dispersion‐based signal drop can be observed in the unipolar MRE scan due to the high encoding efficiency compared with conventional GRE‐MRE. In addition, a reduction in signal magnitude can be observed at the boundaries of the phantom due to transmit field inhomogeneity and the strong flip angle dependency of the time‐reversed SSFP sequence. The central row corresponds to single scans, whereas the bottom row denotes averaged displacement fields for increased SNR. Both conventional and unipolar MRE show qualitatively equivalent wave patterns. The SNR of the conventional MS MRE scan is considerably lower than unipolar MRE, as can be directly appreciated from the graininess of the displacement field of the single scan, which is strongly reduced in the case of 20‐fold averaging. The averaged magnitude of the ratio of unipolar to conventional displacement fields is 94.4% at 35 Hz, showing good quantitative agreement. At 70 Hz actuation, displacement fields reconstructed from unipolar MRE are slightly lower, with an average ratio of 84.0%.

In Figure 9, the magnitude of the complex shear stiffness obtained through a FEM‐based inversion algorithm is shown at 35 and 70 Hz actuation frequency for conventional and unipolar MRE in the case of single and averaged scans. Generally, background as well as the soft inclusions (left inclusions) are well recovered, independent of the technique and the number of averages acquired. Conventional MRE, however, fails to obtain reasonable stiffness values at 35 Hz for the two stiffest inclusions (right inclusions), while unipolar MRE already predicts their high stiffness with only a single scan. In the 20‐fold averaged conventional MRE scan, the increased stiffness of the two right inclusions is well observed. At 70 Hz, both techniques are in good agreement, with only slight underestimation of the shear stiffness of the stiffest inclusion (bottom right) in the case of the single average, conventional MRE scan.

Figure 9

Comparison of the magnitude of the shear stiffness for conventional GRE‐MRE and the proposed unipolar MRE scheme for 35 and 70 Hz actuation frequency with and without averaging. The reported maps are average stiffness values over the six innermost slices. While conventional MRE fails to predict the stiffness of the stiffest inclusions (right inclusions) in the single scan, the increased SNR in unipolar MRE allows recovery of these inclusions without additional averaging

In Figure 10, the ROI‐averaged magnitudes of the complex shear stiffness are displayed, against the OSS‐SNR determined over the same ROI. Conventional MRE is denoted by open circles, while unipolar MRE is shown as squares. Values for single and averaged scans are given in Table S3. For 35 Hz, the OSS‐SNR of unipolar MRE without averaging is higher than the OSS‐SNR of the 20‐fold averaged conventional scan for all regions. At 70 Hz, the single average unipolar MRE acquisition is comparable with an 8‐ to 10‐fold averaged conventional MRE scan. Stiffness of inclusions 1 and 2 as well as the background are well recovered independent of the number of averages and technique. Conventional MRE qualitatively fails to predict the stiffness of inclusions 3 and 4, showing higher stiffness for the softer, third inclusion. Through averaging, the stiffness of inclusions 3 and 4 determined by conventional MRE increases significantly by 54% (inclusion 3) and 179% (inclusion 4) at 35 Hz and by 3.8% and 10.3% at 70 Hz, respectively. In conventional MRE and at 35 Hz, the stiffness of inclusion 3 reaches a plateau after 11 averages, while the stiffness of inclusion 4 increases even after 20 averages. With unipolar MRE, the change in stiffness is below 0.3% for the first three inclusions as well as the background. Change in stiffness for inclusion 4 is 6.2% at 35 Hz and 1.9% at 70 Hz, and reaches a plateau after five averages in both cases.

Figure 10

Comparison of the region of interest (ROI)‐averaged shear stiffness (bars denote the standard deviation over the ROI) obtained with conventional GRE‐MRE (solid line, circles) and the proposed unipolar MRE (dashed line, squares) at 35 and 70 Hz actuation frequency. The x‐axis denotes the average octahedral shear strain‐SNR (OSS‐SNR) for each ROI. Sample averaging was performed to increase the SNR in both unipolar and conventional MRE to demonstrate the influence of SNR on the reconstructed shear stiffness. Here, the left‐most data‐points in each series correspond to the single‐scan case, where each successive point increases the number of sample averages up to NAvg = 20. While the stiffness of soft inclusions (numbers 1 and 2) and the background are found to be independent of the OSS‐SNR, stiff inclusions (numbers 3 and 4) depend on the number of sample averages for conventional MRE at 35 Hz. Due to the increased OSS‐SNR of unipolar MRE, even without averaging, unipolar MRE can directly capture the stiffness of these inclusions

DISCUSSION

In this work, we have proposed a novel MR elastography sequence based on time‐reversed spoiled SSFP, which utilizes the spoiling gradients as highly efficient unipolar MEGs. We have motivated the fundamental encoding principle of unipolar MRE through analogies to phase‐contrast MRI and have described an extension of the two‐transverse TR approximation developed for DW‐SSFP to account for lowest order spin‐history effects.

The theoretical predictions of the proposed approximations have been compared with full EPG simulations in a comprehensive parameter study. It was found that longer T1 and T2 values as well as lower diffusion generally leads to a reduction in encoding efficiency. We attribute this finding to the destructive interference of an increasing number of refocused echo pathways. Contrary to conventional MRE, unipolar MRE encodes displacements nonlinearly into the signal phase, as echo pathways with different motion sensitivities are mixed. However, we have identified a linear encoding domain, which is approximately given by . Since the contributions of the different echo pathways are dependent on relaxation and sequence parameters, it is safe to assume that the displacement range that is mapped linearly to the signal phase will be dependent on these parameters as well. This can also be appreciated in Figure S2, which is the short T1/T2 equivalent of Figure 6 (T1 = 300 ms, T2 = 30 ms). There, the linear encoding domain is approximately twice as large and given by . Since encoding efficiency was found to be monotonically decreasing, unipolar MRE shows heightened sensitivity to low displacement amplitudes.

In the case of conventional GRE‐MRE, the encoding efficiency decreases as the duration of the encoding gradient is reduced, even if a finite moment for each of the two lobes is assumed. By contrast, unipolar MRE shows a finite encoding efficiency in the case of the infinitesimal model given by Equation 12. This demonstrates the uniqueness of unipolar MRE, which, unlike conventional GRE‐MRE, splits encoding over two (and more) TR periods, thus effectively storing spin positions in the phase of the MR signal during the first TR and subtracting the updated position with the second TR. Analogous to SE sequences with single gradient lobes47 or displacement encoding with stimulated echoes (DENSE),55 it is precisely the mixing time between the two encoding steps that leads to the high motion sensitivity of the unipolar MRE sequence, and time‐reversed SSFP sequences in general. This relationship also manifests itself in the comparison of encoding efficiencies as a function of encoding fractions in Figure 6A, where unipolar MRE shows a linear increase in encoding efficiency as a function of the MEG duration, whereas conventional bipolar MEGs follow a quadratic behavior in the small q domain.

Note that the quadratic behavior is reminiscent of the encoding of displacement velocities, while the linear relationship arises from the direct encoding of displacement.3 Since the gradient strength G is kept constant in the plot, the moment Δk vanishes for unipolar MRE with q → 0 and thus also its encoding efficiency.

A direct consequence of encoding displacements in unipolar MRE (as opposed to encoding displacement velocities with conventional MRE) is the dependency of the encoding efficiency on sequence timing rather than encoding fraction (assuming is kept constant). The SSFP‐SE contribution to unipolar MRE is shown to decrease with increasing ratio, as displacements between successive encoding/decoding steps become more and more similar. The full EPG simulations, on the other hand, show increasing encoding efficiency with increasing ratio. We speculate that the interference of echo pathways is reduced as successive TR periods accumulate more and more similar motion phases. The parameter study also revealed an increase in encoding efficiency with increasing flip angle. Since an increase in flip angle leads to a reduction of the population of high‐spatial frequency states, this finding is in general accordance with the dependency on T1, T2 and diffusion. Furthermore, it supports the hypothesis of destructive interference of motion encoding due to differences in motion sensitivity of echo pathways. Contrary to bSSFP‐based MRE,6, 21, 22 the encoding efficiency of unipolar MRE is independent of off‐resonance, as dephasing and refocusing configuration states experience opposite phase accrual. Thus, each echo pathway shows the same dependency on off‐resonance Δ, which is given by . Thus, in equivalence to conventional spoiled GRE‐MRE, off‐resonance leads to an image phase offset, which is independent of the wave‐phase offset and is removed once the Fourier transformation over wave phases is performed.

Since unipolar MRE uses the spoiling gradients as MEGs, the sequence TR can be considerably reduced compared with conventional GRE‐MRE. Given that time‐reversed spoiled SSFP sequences cannot be used with interleaved MS acquisitions, the sequence TR is generally much shorter than typical T2 values. A Taylor series expansion of the echo signal in is at first sight not feasible. Nonetheless, the two‐transverse TR approximation proves to be valuable in describing the encoding efficiency of the unipolar MRE acquisition, especially within the linear encoding regime. This can be appreciated from the closeness of predicted encoding efficiency of the two‐transverse TR approximation and the full EPG simulations in Figures 5 and 6, as well as the close coincidence of the displacement fields in the phantom experiments depicted in Figure 8. While differences between the full EPG and the two‐transverse TR approximation certainly persist in theory, we believe that the attenuation of configuration states is not only due to diffusion and transverse relaxation but is complemented by intra‐voxel phase dispersion (see below) additionally driving the sequence closer to the two‐transverse TR solution.

The propagator definition in Equation 16 assumes periodic bulk motion. From conventional MRE it is known that the variation of the displacement field within one voxel leads to intra‐voxel phase dispersion (IVPD) and thus signal attenuation.56, 57 In the 1D case of a propagating plane‐wave with wavelength , signal attenuation by IVPD in the conventional MRE case can be estimated to be where denotes the expected signal without IVPD and the attenuated signal.8, 57 In the case of unipolar MRE, intra‐voxel variations of the displacement field can be treated using the spatially resolved EPG, which assumes each spatial position as independent EPG simulations.58 The echo signal can then be approximately expressed as the average of sub‐voxel signals, each being a superposition of refocused echo pathways. Since summation over echo pathways and spatial averaging can be exchanged, IVPD in unipolar MRE can not only be discussed as a total effect of signal attenuation, it can also be studied for each echo pathway separately. As the motion sensitivity of configuration states increases proportionally to their spatial frequency, each echo pathway accumulates a different motion‐dependent phase. From Equation 31 it becomes clear that echo pathways involving higher order configuration states will exhibit stronger attenuation due to their increased phase dispersion. Thus, similar to diffusion, the presence of IVPD is expected to reduce the contributions of higher order configuration states to the steady‐state signal and therefore reduces the influence of tissue, experiment and sequence properties on the encoding efficiency.

In the derivation and analysis of unipolar MRE, the assumption of isotropic diffusion was made. While the DTI measurements support this assumption in the phantom used, a discussion of anisotropic diffusion is necessary when unipolar MRE is employed in vivo. McNab and Miller extended the two‐transverse TR approximation for anisotropic DW‐SSFP by replacing the exponent of the diffusion propagator by its respective anisotropic form.44 Weigel et al later introduced anisotropic diffusion in the context of the EPG formalism.59 Both approaches can be equally used here. Since the encoding direction and thus diffusion sensitization is held constant in successive TRs, the diffusion coefficient that enters into the 1D diffusion propagator (Equation 20) can be simply replaced by the projection of the diffusion tensor on the encoding direction Thus, in the presence of anisotropic diffusion, encoding direction‐dependent attenuation of configuration states, and with these direction‐dependent encoding efficiencies, are expected to occur.

In the nonlinear regime (here ), full EPG simulations show a monotonous decline in encoding efficiency, which leads to systematic underestimation of large displacements (Figure 6C). Since the dependency on the displacement amplitude is found to be monotonous within the relevant displacement domain, a bidirectional relationship is likely to exist. With exact knowledge of the field as well as tissue parameters within the ROI, full EPG simulations with different displacement amplitudes can be performed to correct the displacement fields, eg, via a look‐up‐table approach. This also holds for the variation in phase angle (as observed in Figure S4).

The sequence has been validated in a phantom study against a conventional GRE‐MRE sequence at two driving frequencies, where displacement fields were found to be in excellent agreement. We attribute the remaining differences to differences in the slice profile of MS and 3D acquisitions as well as to residual differences in encoding efficiency due to local relaxation parameter variations. Especially at low driving frequency, the displacement noise is visually higher in the conventional MRE scan than in unipolar MRE, which is also confirmed by the OSS‐SNR calculations. The low SNR of the conventional MRE scan leads to underestimation of the stiffness of the two stiffest inclusions. Here, conventional MRE without averaging even fails to qualitatively predict the elevated stiffness of inclusion 4 compared with inclusion 3 (Figure 9). The stiffness of inclusions 1 and 2 as well as the phantom background were found to be in excellent agreement between unipolar and conventional MRE.

Through averaging, the SNR of the conventional MRE scan could be successfully increased to the point of close agreement of recovered stiffness values with unipolar MRE. At 35 Hz, however, even 20‐fold averaging of the conventional MRE scan did not provide sufficient SNR to allow estimation of the stiffness of the two stiffest inclusions. Minor differences in the order of 0.2 kPa persist, however, for inclusion 3 at 35 Hz, and also the two stiffest inclusions at 70 Hz. We believe these differences to be again attributable to slice profile differences, since the high stiffness of the inclusions leads to very low spatial derivatives and thus renders the inversion especially sensitive to geometrical mismatches.

In this study, conventional MRE was performed with fewer slices compared with unipolar MRE (12 vs. 36 slices) to keep scan durations within acceptable ranges. Since unipolar MRE uses slab excitation of the whole phantom and phase encoding in the slice‐encoding direction, SNR of 3D unipolar MRE is ‐fold higher than a respective 2D scan with equal TR, where denotes the number of phase‐encoding steps in the slice direction. Accordingly, the SNR of the unipolar MRE scan was √3 higher, equaling 3‐fold averaging. An OSS‐SNR efficiency adjusted for the slice acquisition technique can be defined, which relates the OSS‐SNR to the total acquisition duration and accounts for the number of averages and slices . For conventional MRE, where MS acquisition is performed, the OSS‐SNR efficiency can be defined as The use of slab excitation (3D) in unipolar MRE leads to At 35 Hz, this results in an OSS‐SNR efficiency of 120 min−1 for conventional MRE and 226 min−1 for unipolar MRE for a single slice and average (at 70 161 min−1 and 1355 min−1, respectively). The OSS‐SNR was determined over the whole phantom for a single scan. Hence, for the acquisition of a single slice and average, unipolar MRE is approximately twice as effective at 35 Hz than conventional GRE‐MRE, and more than 8‐fold at 70 Hz.

Reconstructed stiffness in unipolar MRE converges after 5‐fold averaging compared with at least 20‐fold averaging with conventional MRE. This corresponds to a per‐slice speed‐up of at least 8‐fold in the case of 35 Hz, and 4‐fold in the case of 70 Hz, when taking the acquisition speed‐up of factor 2 per slice into account. Based on the OSS‐SNR, single‐average unipolar MRE corresponds to more than 40‐fold acceleration, as it shows considerably higher OSS‐SNR at 35 Hz than the 20‐fold averaged conventional MRE scan. At 70 Hz, unipolar MRE would be up to 18‐fold faster when compared with conventional MRE.

In the phantom study, the actuation strength was specifically chosen to remain within the linear encoding regime. While in conventional MRE the encoding efficiency can be reduced by gradient scaling or lowering the encoding fraction,6 in unipolar MRE spoiling and motion sensitivity are coupled. Thus, a reduction of encoding efficiency to remain within the linear encoding regime is not feasible. Changes to the duration of the encoding gradient at constant zeroth moment only leads to minor reductions in encoding efficiency. Thus, if linear encoding is intended, displacement amplitudes must be chosen accordingly.

The motion propagator in Equation 18 cannot only be used for the simulation of unipolar MRE, but for sequences with periodic motion in general. Since the propagator describes motion sensitivity of configuration states subject to a generalized trapezoidal gradient, concatenation of the propagator allows performing EPG simulations with bipolar or flow‐compensated gradients. Hence, it can be employed to study spin‐history effects in conventional 3D GRE‐MRE with short TRs or the influence of RF spoiling and its efficiency in canceling higher order echo pathways.

Previous work by Scheffler et al has described an encoding concept similar to unipolar MRE to induce an “oscillating steady‐state” in a spoiled SSFP sequence with readout of the F‐2, F‐1, F0 and F1 configuration states.60 Contrary to our objective to encode and extract the full 3D shear displacement field, Scheffler et al's concept aimed at studying the effect of alternating RF‐phases on the signal amplitude of SSFP‐based sequences. While unipolar MRE can be seen to directly build on their encoding concept, it differs from the work by Scheffler et al in the following important aspects: (i) unipolar MRE is designed to readout the F−1 state only in order to reduce the TR and accelerate MRE acquisition to the fullest extent possible. A dual‐echo approach as employed in Scheffler et al's work for aforementioned purposes could possibly be employed to (a) reduce relaxation time dependencies with combined displacement estimation using the high motion sensitivity of the F−1 state and the relaxation parameter insensitive F0 state, and (b) to combine the technique with double‐echo‐steady‐state (DESS) relaxometry.26 (ii) The work by Scheffler et al is based on the description of a steady‐state that oscillates between two TRs. Our steady‐state is composed of NP distinct signals, which is closer to the description of pseudo‐steady states by Assländer et al and Amthor et al.61, 62 (iii) Spoiling and readout gradients are uncoupled to allow for the sensitization of the sequence to the different motion components of the shear wave field. (iv) Our formulation of steady states is based on the EPG formalism. The two‐transverse TR approximation describes the signal to the second order in and thus makes assumptions on the contributing echo pathways. On the other hand, Scheffler et al developed a signal model for small perturbations of the accumulated phase during a TR period without making assumptions on relaxation properties. Transferring their approach to unipolar MRE, it becomes a model to describe motion encoding for small displacements in the case of NP = 2. An extension of their work to more than two states is beyond the scope of this work; however, it would provide a valuable addition to the understanding of encoding in unipolar MRE in the absence of diffusion or IVPD.

In the present work, Hadamard encoding has been used to allow for the extraction of the 3D volumetric displacement field. It should be noted that unipolar MRE is not restricted to Hadamard encoding alone—all encoding schemes lacking a separate reference scan without motion sensitization can be used. Clinically, MR elastography is often performed on a small set of discontinuous 2D slices with only through‐plane motion sensitization.63, 64 In unipolar MRE, slab excitation and phase encoding in the slice direction can be replaced by single slice excitation. Interleaved MS acquisitions cannot be paired with unipolar MRE as it is described here, as spoiling gradients in consecutive TR periods will accumulate and change the overall spoiling area as well as the motion sensitization scheme.

While the present work has focused on the theoretical foundations and phantom validation of unipolar MRE, future work is required to assess the utility of the method in vivo. The high intended motion sensitivity of unipolar MRE requires careful consideration and provisions depending on the target anatomy, in particular if physiological and patient motion is present. Accordingly, unipolar MRE is expected to require breath‐holding in the case of thoracic and abdominal applications. To address cardiac‐related motion, artifact reduction techniques such as averaging or gating may need to be incorporated.65 In brain and spine, the pulsatile CSF flow is expected to be encoded by unipolar MRE, potentially interfering with displacement encoding.36, 42, 66 However, since SNR is much higher in unipolar MRE than in a comparable conventional GRE‐MRE sequence, a combination of compressed sensing and self‐gating techniques is envisaged.13, 65, 67, 68 Overall, the influence of physiological and patient motion on the steady‐state signal and resulting changes in encoding efficiency remains to be understood. In addition, in vivo parametric mapping of T1, T2 and ADC would be required for the tissue of interest to estimate their impact on the local encoding efficiency. While the phantom results and theoretical considerations suggest that relaxation effects are not critical, these results need to be confirmed in vivo.

Given the considerations above, a first application of unipolar MRE could focus on brain, where mechanical attenuation by the skull results in low displacement amplitudes necessitating high encoding efficiency and SNR.

In the future, unipolar MRE might prove valuable for multi‐frequency MRE, when different actuation frequencies are acquired consecutively. This is especially the case for low driver frequencies as the SNR is significantly increased with unipolar MRE compared with conventional MRE with fractional encoding. At high frequencies, however, unipolar MRE loses its advantage over conventional MRE as encoding fractions typically become much larger and acceleration is limited by the combined duration of RF pulse and readout, not by the MEG duration.

CONCLUSION

We have proposed a novel MR elastography sequence termed “unipolar MRE”, which utilizes the spoiling gradients of a time‐reversed SSFP sequence to acquire the full 3D shear displacement field. Compared with conventional GRE‐MRE, unipolar MRE offers significantly reduced TRs with comparable or even increased displacement encoding efficiency. Accordingly, 4‐ to 8‐fold faster acquisitions at equal or higher OSS‐SNR can be achieved.

FUNDING INFORMATION

This project received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 668039.

Table S1: Dictionary resolution for the fingerprinting reconstruction. T1 resolution was decreased between 100 and 1500 ms to reduce the overall dictionary size. The dictionary was created from all combinations fulfilling T1  > T2, resulting in 35280 fingerprints.

Table S2: Mean and standard deviation of quantitative parameter maps determined for each ROIs. T1 is approximately constant over the whole phantom (2425 ± 33 ms), whereas T2 is significantly higher for the inclusions (632 ± 33 ms) compared to the background (193 ± 25 ms). The ADC has been found to be slightly lower in the phantom background then the inclusions, while fractional anisotropy (FA) is close to zero in the whole phantom (isotropic diffusion). Both magnitude and phase of the encoding efficiency in the case of finite diffusion (D = 1.7 10−3 mm2/s) closely coincide throughout the phantom.

Figure S1: Flip angle train of the spoiled MRF sequence used for T1/T2 determination. TR was kept constant at 16 ms. An inversion pulse precedes the MRF acquisition with a 20 ms inversion delay. Fifteen interleaved Archimedean spirals were used, each interleaved acquired in successive MRF experiments after a pause of six seconds to fully relax the system to equilibrium.

Figure S2: Encoding efficiency comparison of unipolar MRE as shown in FIGURE 6, however, with short T1/T2 values of T1 = 300 ms, T2 = 30 ms. Here, D = 1.7·10−3 mm2/s, f = 35 Hz, and G = 20 mT/m. In a), the Unipolar EPG model with and without diffusion now coincide with the two‐transverse TR approximation. In b), the crossing point of the three models is shifted to lower NP/ND ratios compared to the long T1/T2 case. In c), the linear encoding domain can be seen to increase to about 25 μm

Figure S3: Encoding efficiency comparison of unipolar MRE as shown in Figure 5 for different tissue properties, however, plotting the phase angle instead of magnitude. In a), the change in encoding efficiency phase is shown as a function of diffusion coefficient. For diffusion coefficients found in‐vivo, the phase change is negligible. In b) and d), the dependency on T1 is plotted for long and short T2 values. A monotonous increase in phase angle can be observed, which evaluates to approximately.4 radians over the full range of T1 values of 0 to 2500 ms. c) denotes the change as a function of T2, which drops monotonically with increasing T2 by.2 radians over the simulated range (T2 = 0 to 200 ms) for the case of finite diffusion.

Figure S4: Encoding efficiency comparison of unipolar MRE as shown in Figure 6 for various imaging parameters, however, plotting the phase angle instead of magnitude. Dependencies are plotted for a) encoding fraction, b) number phase offsets and delay, c) displacement, and d) flip angle. For increasing displacements, a non‐linear change in phase angle can be observed, which can be accounted for using a look‐up table approach. In addition, the phase angle is found to depend approximately linearly on the flip angle.

Click here for additional data file.