FULL PAPER Magnetic Resonance in Medicine 73:1514–1525 (2015)

A Framework for Optimization-Based Design of Motion Encoding in Magnetic Resonance Elastography Guy Nir,1* Ramin S. Sahebjavaher,1 Ralph Sinkus,2 and Septimiu E. Salcudean1 Purpose: In conventional three-dimensional magnetic resonance elastography, motion encoding gradients (MEGs) synchronized to a mechanical excitation are applied separately in each direction to encode tissue displacement generated by the corresponding waves. This requires long acquisition times that introduce errors due to patient motion and may hinder clinical deployment of magnetic resonance elastography. In this article, a framework for MEGs sequence design is proposed to reduce scanning time and increase signal-to-noise ratio. Theory and Methods: The approach is based on applying MEGs in all three directions simultaneously with varying parameters, and formulation of the problem as a linear estimation of the wave properties. Multidirectional MEGs sequences are derived by setting the problem in an experimental design framework. Such designs are implemented and evaluated on simulation and phantom data. Results: Estimation error of the displacement using the proposed MEGs designs is reduced up to a factor of two in comparison with a unidirectional design with a same number of acquisitions. Alternatively, for the same error, scanning time is reduced up to a factor of three using the multidirectional designs. Conclusion: The proposed framework generalizes acquisition of magnetic resonance elastography, and allows quantification of design performance, and optimization-based derivation of C 2014 designs. Magn Reson Med 73:1514–1525, 2015. V Wiley Periodicals, Inc. Key words: magnetic resonance imaging; elastography; linear estimation; experimental design

INTRODUCTION The characterization of mechanical properties of biological tissues is important in the detection and diagnosis of diseases such as cancer, in which affected tissues typically alter their viscoelastic properties. Manual palpation provides a nonquantitative estimation of such properties

1 Department of Electrical and Computer Engineering, University of British Columbia, Vancouver, Canada. 2 Division of Imaging Sciences and Biomedical Engineering, King’s College London, London, UK. Grant sponsor: Natural Sciences and Engineering Research Council of Canada (NSERC); Grant sponsor: Canadian Institutes of Health Research (CIHR); Grant sponsor: BC Innovation Council NRAS Program. *Correspondence to: Guy Nir, M.Sc., 2332 Main Mall, Vancouver, BC, Canada V6T 1Z4. E-mail: [email protected]

Additional Supporting Information may be found in the online version of this article. Received 3 September 2013; revised 19 March 2014; accepted 13 April 2014 DOI 10.1002/mrm.25280 Published online 6 May 2014 in Wiley Online Library (wileyonlinelibrary.com). C 2014 Wiley Periodicals, Inc. V

that is limited to the surface of the body. Magnetic resonance elastography (MRE), first introduced in (1), is a noninvasive technique for characterizing the mechanical properties of tissue based on magnetic resonance imaging (MRI). The technique have been applied to study tissue changes associated with diseases of, e.g., the liver, breast, prostate, brain, and heart (2). The most prevalent dynamic MRE technique involves the application of a harmonic mechanical excitation that generates shear waves in the tissue. The resulting displacements are estimated by using an MRI pulse sequence with motion encoding gradients (MEGs). These MEGs are bipolar gradients (typically in shape of a sinusoid or trapezoid) that are synchronized to the mechanical excitation. The phase accumulated by each isochromat is encoded in the phase component, also known as the phase image, of the MRI signal. Conventional motion encoding in MRE comprises unidirectional MEGs, in which each component of the displacement vector is encoded separately by applying the gradients in the corresponding direction. The measured phase accumulation, i.e., each voxel of the phase images, is a sinusoid with respect to a phase shift between the MEG and the mechanical wave. Thus, a sufficient number of consecutive acquisitions with varying phase shifts are required to sample the temporal characteristics of the wave. The amplitude and phase of the displacement at each voxel are then estimated through a discrete Fourier transform (DFT) of the samples. Elastograms, i.e., images that depict tissue elasticity or stiffness can be reconstructed from the displacements by solving an inverse problem using, e.g., (3–7). In (8), a generalized sampling method is proposed to achieve a sub-Nyquist sampling rate of two acquisitions per direction using a least-squares approach, thus reducing acquisition time of unidirectional sequences. In practice, however, due to noise in the measured phase images, MRE sequences typically require four or eight acquisitions in each direction to provide a displacement map and, in turn, an elastogram, with a reasonable signal-to-noise ratio (SNR). This yields a longer acquisition time that hinders usage of MRE in some cases due to time constraints and limits on image quality. Longer acquisition times also introduce artifacts in the images due to patient movement. There are approaches that have been proposed to improve SNR or, alternatively, reduce the acquisition time of MRE (2). Most approaches propose improved imaging pulse sequences, such as echo-planar imaging (3) and gradient-echo (9), as opposed to the basic spin-echo sequence. Other approaches deal directly with the design of the MEGs, e.g., fractional encoding (10,11), in which the mechanical frequency is taken as a fraction of the MEGs oscillation frequency. In Refs. (12,13), sequence

1514

Optimization-Based Design of Motion Encoding in MRE

1515

designs with combination of MEGs were applied to phase-contrast MRI, a field related to MRE, to reduce scan time and increase sensitivity for estimating the magnitude of flow in magnetic resonance angiography. Recently, (14,15) have proposed approaches for encoding the three-dimensional displacements in MRE by applying the MEGs in all directions simultaneously. These multidirectional sequences have the potential of reducing the number of scans by three, while achieving the same SNR as in unidirectional sequences. In (14), MEGs are applied simultaneously with different frequencies, and the displacement is encoded for each direction by utilizing the filtering condition in (16). In (15), the authors propose a sample interval modulation MRE, in which MEGs are applied simultaneously with different phase shift intervals, such that the phase accumulation in each direction can be recovered from the total accumulation using DFT. In this article, we focus on optimization of the MEGs design. We parameterize general multidirectional sequences, and formulate the estimation of the threedimensional displacement as a linear system in the timedomain. This allows derivation of a least-squares estimator that is optimal under a given design and noise model. We then characterize and set the problem in an experimental design framework, by which we quantify the performance of MEGs sequences, and derive parameters that lead to optimal multidirectional designs under a given optimality criterion. We fit the conventional MEGs sequence in the proposed framework, and compare the sequences on simulation and phantom data. The remainder of this article is organized as follows. First, we present the generalized MEGs and corresponding phase accumulation, and propose a least-squares estimation of wave properties. We then present the framework for optimal design of motion encoding, and derive optimal sequences. We continue with a description of the experiments that were performed to evaluate our approach, and summarize the results. Finally, we conclude with discussion and future research directions. A summary of the notations used throughout this work is provided in a Supporting Information table.

THEORY Parameterized MEGs We denote the spatial coordinate vector as r ¼ ðx; y; zÞT and the time as t. We represent a steady-state harmonic mechanical wave through tissue by its three-dimensional  ¼ ðCx ; Cy ; Cz ÞT . Without loss of displacement vector C generality, we write Ci ðr ; tÞ ¼ =mfjji ðr Þjejwi ðr Þ ejvm t g ¼ jji ðr Þjsinðvm t þ wi ðr ÞÞ ; i ¼ x; y; z [1] The wave frequency vm is the same as the known mechanical excitation frequency. The spatially varying amplitudes and phases in each direction, jji j and wi , are unknown and have to be determined at each imaged voxel based on a number of phase images acquired in different acquisitions. Note that this representation can

accommodate locally any type of wave (e.g., traveling or standing), wavefront (e.g., plane or spherical), polarization (e.g., linear or circular), and attenuation. We parameterize general multidirectional sequences of  ¼ ðGx ; Gy ; Gz ÞT . Specifically, for each sinusoidal MEGs, G direction, we allow control of the (angular) frequency vi , the amplitude jGi j, the start phase bi , the phase shift relative to the mechanical wave ai , and the number of periods the gradients are applied Ni 2 N. We therefore have Gi ðt;vi ;Ni ;jGi j;ai ;bi Þ [2] ¼ jGi jsinðvi ðt  ai =vi Þ þ bi ÞWNi Ti ðt  ai =vi Þ; i ¼ x;y;z where WNi Ti ðtÞ is a window function with a ½0;Ni Ti  support, and Ti ¼ 2p=vi is the time period of each MEG. The effect of the parameters ai and bi on the MEGs is demonstrated on a single acquisition in a Supporting Information figure. Geometrically, this may be seen as if the gradients are rotating in a (piecewise) helical fashion as a function of time in gradient space. Parameterized Phase Accumulation In this article, we limit our discussion to MEGs with the same frequency as the mechanical frequency, i.e., vi ¼ vm , 8i ¼ x; y; z. However, the theory can be extended to fractional encoding in each direction. The phases accumulated by the precession of each moving isochromat during the acquisition sequence, as result of the MEGs, are encoded in the measured phase images of the MRI signal. As derived in Appendix, the phase accumulation at each voxel as a function of the MEGs parameters is given by  ; jGj; a  ¼ Fðr ; N  ; bÞ

X gpNi jGi j jji ðr Þjcosðwi ðr Þ þ ai  bi Þ vm i¼x;y;z [3]

where c is the gyromagnetic ratio. This equation is a generalization of the phase accumulation that was originally presented in (1). To formulate the problem as a linear one, for every voxel, we define Ai and Bi to be such that jji ðr Þj ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A2i ðr Þ þ B2i ðr Þ;

  Bi ðr Þ wi ðr Þ ¼ arctan ; Ai ðr Þ

i ¼ x; y; z [4] and, using the harmonic addition theorem, substitute in Eq. [3] to obtain  ; jGj; a  ¼ Fðr ; N  ; bÞ

X gpNi jGi j ½Ai ðr Þcosðai  bi Þ vm i¼x;y;z þ Bi ðr Þsinðai  bi Þ

[5]

Thus, for every voxel, estimation of the amplitudes jji j and phases wi in each direction, is given by estimation of the six unknowns Ai and Bi, i ¼ x; y; z. This requires a minimum of six acquisitions, i.e., six acquisitions of the phase image with different MEG parameters.

1516

Nir et al.

Optimal Estimation of Wave Properties In general, we acquire K  6 acquisitions with, possibly, ðkÞ ðkÞ ðkÞ ðkÞ different ai , bi , Ni , and jGi j in each direction 0 B B B B B B B B gp B B M¼ vm B B B B B B B B @

ð1Þ

ð1Þ Nx jGð1Þ x jcosðdx Þ

i ¼ x; y; z, and for each acquisition k ¼ 1; 2; . . . ; K. We ðkÞ ðkÞ ðkÞ denote di ¼ ai  bi as the net phase of each acquisition, and define the K  6 (tall) acquisition design matrix

ð2Þ

ð2Þ Nx jGð2Þ x jcosðdx Þ

ðKÞ Nx jGðKÞ x jcosðdx Þ

ð2Þ



ðKÞ Nx jGðKÞ x jsinðdx Þ

ð2Þ



ðKÞ Ny jGðKÞ y jcosðdy Þ



ðKÞ ðKÞ Ny jGðKÞ y jsinðdy Þ

ð2Þ



ðKÞ Nz jGðKÞ z jcosðdz Þ

ð2Þ



ðKÞ Nz jGðKÞ z jsinðdz Þ

ð1Þ

ð2Þ Nx jGð2Þ x jsinðdx Þ

ð1Þ Ny jGð1Þ y jcosðdy Þ

ð1Þ

ð2Þ Ny jGð2Þ y jcosðdy Þ

ð1Þ ð1Þ Ny jGð1Þ y jsinðdy Þ

ð2Þ ð2Þ Ny jGð2Þ y jsinðdy Þ

ð1Þ Nx jGð1Þ x jsinðdx Þ

ð1Þ

ð2Þ Nz jGð2Þ z jcosðdz Þ

ð1Þ

ð2Þ Nz jGð2Þ z jsinðdz Þ

ð1Þ Nz jGð1Þ z jcosðdz Þ ð1Þ Nz jGð1Þ z jsinðdz Þ

By Eq. [5], for each voxel r we have a (typically overdetermined) linear system of equations  ¼ MA þ w  F

[7]

 ¼ ðFð1Þ ; . . . ; FðKÞ ÞT is a vector that comprises where F the measured phase accumulations in each acquisition, A ¼ ðAx ; Bx ; Ay ; By ; Az ; Bz ÞT is the unknown coefficient  is a vector that we need to solve the system for, and w random equation error that depends on the type of noise in the system. We can rewrite Eq. [7] as the estimation problem ^ ¼ arg min jjF   MAjj ; given M; F  A

[8]

A

 , which can be interpreted as the The distribution of w noise introduced into the measured phase images during each acquisition, determines the type of norm in Eq. [8], ^ . We model w  as a white and the choice of estimator A Gaussian noise (see Discussion section) with a variance of s2 , and assume that M is full column rank. Under these assumptions, the Gauss–Markov theorem (17) states that the ordinary least-squares estimator is optimal (best linear unbiased estimator) in the sense of minimizing the mean squared error (MSE) in the estimation ^  Ajj2 . The least-squares estimator is defined by EjjA ^ ¼ ðMT MÞ1 MT F  ¼MF  A y

ðKÞ



[9]

y

where M is the Moore–Penrose pseudoinverse of the matrix M. Optimal Design of Motion Encoding The proposed MEGs parametrization provides several degrees of freedom that allow us to modify the acquisition design matrix in Eq. [6]. The parameters should be diverse among acquisitions to yield linearly independent rows, such that M is full column rank. Clearly, the more acquisitions we acquire (larger K), the lower is the estimation error that the least-squares estimator obtains.

ðKÞ

ðKÞ

ðKÞ

ðKÞ

1T C C C C C C C C C C C C C C C C C C A

[6]

However, this would increase acquisition time, which is undesirable as explained above. Another option to reduce error is to increase the MEGs amplitude jGi j, the number of periods that they are applied Ni, or to decrease their frequency vi . However, the allowable range for the parameters above is rather small due to hardware limitations, acquisition time, and the size of the scanned target (such that enough wavelengths are imaged). Therefore, we focus on finding the parameters ai and bi that optimize the acquisition design matrix in the sense of minimizing the estimation error by the least-squares estimator in Eq. [9]. The problem above can be considered as a special case of the experimental design problem in statistics (18,19). A popular criterion for an optimal design is the so called D-optimality (20,21), in which the determinant of MT M, also known as the information matrix of the design, is maximized. Such D-optimum design minimizes the (generalized) variance of the best linear estimator (22). Thus, maximizing the determinant of the information matrix is equivalent to minimizing the estimation error. We can, therefore, formulate the D-optimal design problem in MRE as ðkÞ

ðkÞ

ðai ; bi Þ ¼ arg max MT M ; ðkÞ

i ¼ x; y; z ; k ¼ 1; 2; . . . ; K

di

given ðkÞ ðkÞ ðkÞ ðkÞ vm ; Ni ; jGi j; and possible conditions on ai ; bi [10] As we typically do not favor any particular direction or acquisition, unless stated otherwise, we choose the ðkÞ ðkÞ same Ni ¼ N and jGi j ¼ jGj, 8i ¼ x; y; z, k ¼ 1; . . . ; K. Although designs that yield the same determinant are theoretically equivalent in their estimation error, in practice, issues such as echo time (TE) and pulse repetition time (TR) increase, finite slew rate, and different flow compensation characteristics affect SNR. Thus, in some cases, to minimize the TE of the sequence in each acquisition, we set the ðkÞ same ai ¼ aðkÞ , 8i ¼ x; y; z. In other cases, to prevent disðkÞ continuities (jumps) in the MEGs, we set bi ¼ f0; pg.

Optimization-Based Design of Motion Encoding in MRE

1517

FIG. 1. RME design. Comprises a minimum number of unidirectional acquisitions.

Looking at Eq. [6], we have a constant multiplied by a ðkÞ ðkÞ matrix with elements cosðdi Þ and sinðdi Þ. In general, finding an L  L matrix that maximizes the (absolute) determinant among all L  L matrices with elements with an absolute value less or equal to one, is the wellknown Hadamard’s maximum determinant problem (23). It turns out that such (nonunique) matrix HL , called a Hadamard matrix, is known to exist on the real field for L ¼ 1; 2, L  0ðmod4Þ, and its elements take the extreme values 61. For cases of L in which a Hadamard matrix does not exist, there may still be a D-optimum design matrix (24). However, it is typically of more value to use the closest H-matrix (25), which is a non-square matrix K  L with K > L that comprises elements of 61, and satisfies HTKL HKL ¼ KI L , with I L being the L  L identity matrix. Such matrices are D-optimum with respect to a K  L design with detðHTKL HKL Þ ¼ K L .

that a minimal error with respect to sampling of the phase images is attained with evenly distributed acquisitions (in each direction) over the ½0; pÞ range for two acquisitions per direction. Thus, we limit our choice to ðkÞ ai ¼ f0; p2 g. Finally, we choose zero start phases ðkÞ bi ¼ 0. We, therefore, define a unidirectional design with a minimum of six acquisitions as jGð1;...;6Þ j ¼ ðjGj; jGj; 0; 0; 0; 0Þ; jGð1;...;6Þ j ¼ ð0; 0; jGj; jGj; 0; 0Þ ; x y j ¼ ð0; 0; 0; 0; jGj; jGjÞ jGð1;...;6Þ z ð1;...;6Þ

ai

ð1;...;6Þ

bi

¼ 0;

i ¼ x; y; z

i ¼ x; y; z

[11]

This sequence design, referred to as reduced motion encoding (RME) MRE in (8), is illustrated in Figure 1, and by Eq. [6] it yields a simple diagonal design matrix

Unidirectional Designs Here we fit a unidirectional MEGs sequence designs into the proposed framework. This is done to compare and evaluate the multidirectional designs that are presented next. This special case is also important as it shows that our framework allows estimation of the displacement map using a unidirectional design with a minimum of two acquisitions for each direction, as suggested in (8). In order to apply only one MEG direction in each acquisition, we toggle the MEGs amplitudes as ðkÞ jGi j ¼ f0; jGi jg. We follow the result of (8), which states

¼ ð0; p=2; 0; p=2; 0; p=2Þ ;

Mrme ¼

gpNjGj I 6; vm

 1 with detðMTrme Mrme Þ 2 ¼



gpNjGj vm

6

[12] For unidirectional designs that comprise more than two acquisitions per direction, (8) shows that ai should be evenly distributed over the ½0; 2pÞ range in each direction. Similar to Eq. [11], we can construct such sequences and compute their design matrices. We consider a conventional design of 24 acquisitions (eight in each direction) for our evaluation and denote its design as

1518

Nir et al.

FIG. 2. MD-RME design. Comprises a minimum number of multidirectional acquisitions, without discontinuities in the MEGs or increase in TE.

Mcon . We note that by explicitly writing the proposed least-squares solution in Eq. [9] for Mcon , it can be shown that the estimated wave properties, jji j and wi , using our time-domain approach coincide with the results achieved by a standard DFT processing of the phase images, in which jji j and wi are estimated by the amplitude and phase of the fundamental frequency. Multidirectional RME Design Here, we are interested in generating a multidirectional MEGs sequence that takes advantage of the parameterized approach to reduce estimation error. We note that using different phase shifts ai 6¼ aj , in two directions or more, increases the TE of the acquisition so as to cover the entire MEGs sequence, which decreases SNR due to a weaker signal. The TE is measured based on the union of the time the MEGs are applied in a single acquisition, rather than from t ¼ 0, therefore TE is still minimal even if ai 6¼ 0, as long as it is the same in all directions. In addition, using a start phase bi 6¼ f0; pg produces discontinuities in the MEGs waveform. In practice, gradient amplifiers impose a finite slew rate, and limitations on specific absorption rate restrict excessive dB/dt. This introduces errors to the phase accumulation calculation that lead to estimation errors in the displacement. The effect of a finite slew rate on the sequence can be observed in a Supporting Information figure. Thus, in this example we add to the optimal design problem in Eq. [10] the constraints bi ¼ f0; pg and ai ¼ a, 8i ¼ x; y; z. We propose the following design with sign-flips

ð1;...;6Þ

ai

¼ ð0; p=2; 0; p=2; 0; p=2Þ ; i ¼ x; y; z bð1;...;6Þ ¼ ðp; p; 0; 0; 0; 0Þ x bð1;...;6Þ ¼ ð0; 0; p; p; 0; 0Þ y bð1;...;6Þ ¼ ð0; 0; 0; 0; p; pÞ z

[13]

As it comprises the minimum number of acquisitions, we refer to it as a multidirectional RME (MD-RME) design. The sequence is illustrated in Figure 2, and yields a Toeplitz design matrix 0

1

0

1

0

1

0

1

0

1

0

0

1

C 1 C C C 1 0 1 0 1 0 C C Mmdrme C; 0 1 0 1 0 1 C C C 1 0 1 0 1 0 C A 0 1 0 1 0 1    12 gpNjGj 6 T with detðMmdrme Mmdrme Þ ¼ 16 vm [14] B B B B gpNjGj B B ¼ B vm B B B B @

Thus, for the same minimum number of acquisitions, the determinant value of the estimation covariance is reduced by a factor of 16 in comparison with the unidirectional RME sequence. We note that, since we are interested in the estimation error of each couple of coefficients, Ai ; Bi , we should look at the elements in the diagonal of the inverse information matrix ðMTmdrme Mmdrme Þ1 that correspond to the variances of

Optimization-Based Design of Motion Encoding in MRE

1519

FIG. 3. Optimal design with varying start phases. Notice the discontinuities in the MEGs. Using phase shifts, we achieve discontinuity of only one direction per acquisition, without increase in TE.

each coefficient. We find that these elements are half of their value in the unidirectional RME case, thus we expect the standard deviation of the estimation error to pffiffiffi be reduced by 2 using the MD-RME sequence. Optimal Design with Varying Start Phases Here, rather than proposing a sequence and evaluating its performance, we take the inverse approach where we find a design that produces an optimal acquisition matrix. In our case L ¼ 6  2ðmod4Þ, a Hadamard matrix does not exist, and the closest H-matrix is for K ¼ 8. Any such (nonunique) matrix can be chosen, and, using the technique in (24), we constructed 0

H8;6

1

B B1 B B B1 B ¼B B1 B B B1 @ 1

1

1

1

1 1

1

1

1

1 1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1T

C 1 C C C 1 1 C C C 1 1 C C C 1 1 C A 1 1

8i ¼ x; y; z. It can be proved that no combination with ðkÞ only bi ¼ f0; pg can produce H8;6 . Thus, we are forced to have discontinuities in the MEGs in this case, and to ðkÞ minimize the slew rate effect we use bi ¼ f0; pg when possible. We propose the following sequence ð1;...;8Þ

ai

¼

  p 3p 3p 3p 5p 5p 5p 3p ; ; ; ; ; ; ; ; i ¼ x; y; z 4 4 4 4 4 4 4 4  p 3p ;0 ¼ 0; ; 0; p; 0; 0; bð1;...;8Þ x 2 2   p p bð1;...;8Þ ¼ 0; 0; p; ; ; p; 0; 0 y 2 2   p 3p ð1;...;8Þ bz ; p; 0 [16] ¼ 0; p; ; 0; p; 2 2

The sequence is illustrated in Figure 3, and yields the design matrix

1

[15]

Given constant gradient amplitude and number of ðkÞ periods, the matrix H8;6 suggests that di ¼ f3p=4; ðkÞ p=4; p=4; 3p=4g. The technique to determine di , for each acquisition k and direction i, is to observe the pair of elements ðk; 2i  1Þ, ðk; 2iÞ. Based on Eq. [6], these eleðkÞ ðkÞ ments are proportional to cosðdi Þ and sinðdi Þ, respecðkÞ tively. Thus, the value of di is determined by the quadrant defined by these elements [e.g., the pair ð1; 1Þ yields d ¼ p=4]. We can derive several sequences that produce H8;6 . Here, to minimize TE, we are interested in having no ðkÞ phase shift between directions, i.e., ai ¼ aðkÞ ,

h i12 gpNjGj Mopt ¼ pffiffiffi H8;6 ; with detðMTopt Mopt Þ 2 vm  gpNjGj 6 ¼ 64 [17] vm pffiffiffi Note we have a 1= 2 loss (cosine and sine of multiples of p=4), but the overall gain of 64 of the optimal design compared to the RME design is still significant. By looking at the diagonal of the inverse information matrix, we find that its elements are 0.25 of their value in the RME case, thus we expect the standard deviation of the estimation error to be reduced by two by the optimal design. Remember that we have used eight acquisitions for this design, rather than six as in the previous designs. In case we use eight acquisitions of an RME sequence (by replicating two acquisitions), we have that the optimal

1520

Nir et al.

FIG. 4. Optimal design with varying phase shifts. The maximum phase shift at each acquisition dictates the increase in TE. Using signflips, this phase shift is at most p2 .

design has a gain of 32 relative to such a design. In fact, the gain of this optimal design with only eight acquisitions equals the gain of the conventional design Mcon , which comprises 24 acquisitions.

direction is an increase of TE. Using the combinations of ðkÞ ðkÞ ai and bi above, the phase shifts between the directions in each acquisition are minimized such that the maximum phase shift is p=2, which translates into a TE increase by a quarter of the time period.

Optimal Design with Varying Phase Shifts

METHODS

In the previous optimal design, the single discontinuity in some acquisitions may hinder its usage in practice due to a limited slew rate. In addition, as pointed out in (15), different start phases may affect the SNR in in vivo environments, as MEGs with different start phase have different flow compensation characteristics. Thus, here we propose a sequence with a start phase that is limited ðkÞ ðkÞ to bi ¼ f0; pg, but with different phase shifts ai for each direction. Such multidirectional approach was proposed in (15) as sample interval modulation MRE to modulate the sampling interval such that the phase accumulation in each direction has a different frequency. As described above, to yield the same optimal design ðkÞ ðkÞ as in Eq. [17], we derive the values of ai such that di are in the corresponding quadrants determined by H8;6 . We propose the sequence

Simulation of a Plane Wave

 p p 3p 3p p p 3p ; ; ; ; ; ; ; 4 4 4 4 4 4 4 p 3p 3p p 3p p p ; ; ; ; ; ; ; ¼ að1;...;8Þ y 4 4 4 4 4 4 4 p 3p p 3p p 3p p ; ; ; ; ; ; ; að1;...;8Þ ¼ z 4 4 4 4 4 4 4 ¼ ð0; 0; 0; p; p; p; p; 0Þ bð1;...;8Þ x bð1;...;8Þ ¼ ð0; 0; p; 0; 0; 0; p; 0Þ y að1;...;8Þ ¼ x

bð1;...;8Þ ¼ ð0; p; 0; 0; 0; p; 0; 0Þ z

 3p 4 3p 4 3p 4

[18]

The sequence is illustrated in Figure 4. As discussed in (15), a disadvantage of using different phase shifts for each

We first test our approach by simulating a linearly polarized wave propagating through a three-dimensional volume of size 64  64  64 that comprises a background with a spherical inclusion. The stiffer inclusion is associated with a higher propagation speed (lower wave number), and a lower displacement amplitude than the background. We have used a shear wave frequency as in the phantom experiment (200 Hz), and chosen the propagation speeds (10 m/s), and amplitudes (10 mm) to be of the same order of magnitude as in the phantom experiment below. The goal of this simulation is to validate the least-squares wave estimation framework, and theoretical performance analysis of the proposed designs. We define a gradient sequence with the proposed parametrization of the MEGs, and compute the associated phase images for each acquisition by integrating the product of the wave and MEGs over time. Note that, in order to validate our analytical expression for the phase accumulation, we compute the integral numerically. To model the imaging noise, we add a zero mean Gaussian noise to each voxel of the phase images. We repeat the simulations for the unidirectional RME, MD-RME, and optimal designs proposed above, each with a range of noise variances. We calculate the coeffi^ using the least-squares estimator, and cient vector A derive the estimated amplitude and phase of the wave,

Optimization-Based Design of Motion Encoding in MRE

1521

FIG. 5. Simulation of a plane wave. The estimated amplitude and phase of the wave using the unidirectional RME, MD-RME, and optimal design with eight acquisitions for a standard deviation of s ¼ p=5 of the noise added to the phase accumulation images. The conventional design yields similar results as the optimal design, but with three times the number of acquisitions. Top: The amplitude in the z-axis, intensity scale is in [mm]. Bottom: The phase in the z-axis (wrapped in ½p; pÞ), intensity scale is in [rad].

j^j i j, w ^ i in each direction and for each voxel. Slices from the estimated volumes are shown in Figure 5. Since we know the ground truth of the wave phases and amplitudes, we may compute the estimation error as the MSE between the estimations and the ground truth, i.e., e2j

1 ¼ jVj

Z sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 X  j^j i ðr Þj  jji ðr Þj dr ; V

i¼x;y;z

Z sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 X  1 e2w ¼ w ^ i ðr Þ  wi ðr Þ dr jVj V i¼x;y;z

[19]

where V and jVj are the volume and its size, respectively. Phantom Experiment We have implemented the RME and the MD-RME sequences on a 3 Tesla system (Achieva 3.0T, Philips, The Netherlands) using a spin-echo echo-planar imaging pulse sequence with a TE/TR of 30ms/1600 ms. A 200 Hz vibration was applied to a rectangular shaped elasticity phantom (Model 049, CIRS, VA) using a custom shielded electromagnetic transducer synchronized to the scanner. As seen in Figure 6, the phantom contains spherical inclusions of different radii and stiffness values relative to the background material that mimic different tissue types. The entire phantom was scanned in a field of view of 128 128  20 grid, with an isotropic 2  2  2 mm3 resolution. The scanning time was 364 s per acquisition. To evaluate our approach, we have oversampled the scan of each acquisition in both designs with the phase shifts of the MEGs ranging evenly over ½0; 2pÞ relative to

the mechanical excitation. A 80  60  20 region of interest on the acquired axial phase accumulation images was selected, compensated and phase unwrapped. We estimated the displacements using the proposed approach with a corresponding design. We have also generated an elastogram using a local frequency estimation technique (3). We note that these offline processing steps are performed in the same manner as in standard processing of MRE, with the exception of the displacement estimation that involves multiplications of matrices and inversion of the information matrix instead of the standard DFT computation (total processing times are similar). We first estimated the displacements using the unidirectional and multidirectional designs with a minimum of six acquisitions each (RME and MD-RME). We then repeated for the 12 and 24 acquisition versions of the designs (where we refer to the unidirectional design with 24 acquisitions as the conventional MRE sequence), and evaluated the estimation errors in the amplitude, phase, and elasticity. As we do not have the ground truth of the displacements, we have used the best estimate obtained using our least-squares estimation approach over the entire set of acquisitions (i.e., a total of 8  6 acquisitions that include both unidirectional and multidirectional designs). The MSE of each design is computed using Eq. [19]. Further evaluation of the approach is performed by comparing the elastograms to the reported elasticity values of the phantom.1 We segmented the inclusions, and 1 Values were averaged based on recent MRE and ultrasound elastography scans proposed in (26,27). Possible causes for differences to reported values by the manufacturer are wear and tear, dehydration and/or temperature difference of the phantom.

1522

Nir et al.

FIG. 6. Phantom data. a: The elasticity phantom that was used for evaluation, b: A cross-section with the factory reported elasticity values of tissue types, c, d: A cross-section of the best estimated amplitude and phase of the wave in the z-axis, as obtained by the least-squares approach over the entire set of acquisitions, and used for evaluation.

FIG. 7. Simulation results. The MSE of the estimated amplitude and phase as a function of the standard deviation r of the noise added to the phase accumulation images. The lines and bars represent the mean and half of the standard deviation of the error over all voxels, respectively.

Optimization-Based Design of Motion Encoding in MRE

1523

Table 1 Phantom Results Unidirectional designs Amplitude error [mm] Phase error [rad] Elasticity error [kPa]

Multidirectional designs

Six Acq. (RME)

12 Acq.

24 Acq. (conventional)

Six Acq. (MD-RME)

12 Acq.

24 Acq.

46.6633.9 1.460.9 29.6622.6

2.261.8 0.460.4 3.8163.3

2.161.7 0.360.4 3.763.1

11.568.3 0.960.8 14.9613.1

1.861.5 0.260.3 2.662.2

1.461.1 0.260.3 1.761.5

Estimation errors in comparison with the best estimate obtained by the entire set of acquisitions.

calculated the MSE between the reconstructed elastograms and the elasticity of each type of tissue over the corresponding region. RESULTS Simulation of a Plane Wave The error plots as a function of the standard deviation r of the noise are shown in Figure 7. We notice linear and nonlinear relationships between the standard deviation and the MSE of the amplitude and phase, respectively. As expected, the results using the conventional sequence with 24 unidirectional acquisitions coincide with those of the optimal sequence with eight acquisitions. By dividing the amplitude error plots, we find that the estimation errors of the RME design are reduced by 1:51 60:04 and 2:0360:06 using the MD-RME and optimal designs, respectively. Similar slope ratios were computed for the phase error plots in the linear region that corresponds to small values of standard deviations. pffiffiffi These results fit the theoretical error reductions of 2 and 2 in the estimation of the coefficients, to which, by Eq. [4], the amplitude and phase are proportional linearly and nonlinearly.

Phantom Experiment The quantitative results are summarized in Table 1. The ratio between the amplitude estimation errors of the unidirectional and the multidirectional with sign-flips designs for the 6, 12, and 24 acquisitions experiments were 4.05, 1.19, and 1.47, respectively. Thus, the error reductions p by ffiffiffi the multidirectional designs are close to the theoretical 2 value that was calculated from the information matrices, and in the six acquisitions case significantly surpasses it. The elastograms that were reconstructed based on the unidirectional and multidirectional designs for 6, 12, and 24 acquisitions are illustrated in Figure 8. The corresponding errors in comparison with the factory reported elasticity are presented in Figure 9. We notice that, as expected, the more acquisitions we have, the closer are the elasticity values to their factory values. Also, on average, the values generated from the multidirectional design are closer to their reported values. DISCUSSION AND CONCLUSIONS In this article, we have presented a framework for motion encoding design in MRE. The main contributions of this work are as follows:

FIG. 8. Phantom elasticity. A slice of the elastograms that were reconstructed from the estimated amplitude and phase of the wave with a different number of acquisitions. Top: Using the unidirectional design. Bottom: Using the multidirectional design. Intensity scale is in [kPa] and same as in Figure 6b.

1524

Nir et al.

FIG. 9. Phantom elasticity results. The errors between the elastograms of the unidirectional and multidirectional designs in comparison with the true elasticity on each tissue type.

i. formulation of displacement estimation using multidirectional MEGs as a linear estimation problem in the time-domain, which allows derivation of an optimal estimator (under a given design and noise model), ii. characterization of the MEGs design as an experimental design problem, which allows quantification of design performance and derivation of optimal sequences (under the optimality criterion), iii. derivation and implementation of a multidirectional sequence with a minimum number of acquisitions that does not require an increase in TE, nor discontinuities in the MEGs. The approach was tested and evaluated using simulation and experiments on a phantom. The simulation confirms that optimal sequences have the potential of reducing the number of acquisitions required by unidirectional MRE up to a factor of three while acquiring images with the same SNR. Alternatively, for the same number of acquisitions, the estimation error of displacement amplitude can be reduced up to a factor of two. The produced elastograms are also more accurate, since reconstruction is positively correlated with the accuracy of the displacement input. A simplified analogy of the problem is the YatesHotelling weighing problem, in which it was shown that using the same number of weighings, individual weights of objects can be estimated better by weighing them in combinations, rather than separately (28). Here, each acquisition is considered a “weighing,” the measured phase accumulation are the “scale readings,” and using the MEGs parameters we assign each direction a “pan on the scale” or exclude it. With scanner limitations on the slew rate of the gradients, an optimal design with varying start phases may cause errors in phase accumulation by encoding higher frequency components. Alternative implementations of optimal designs that are not affected by the slew rate, but require longer TE, can be achieved with varying phase shifts as proposed in (15). Indeed, our framework shows that the sequence in (15) yields an optimal design under the noise model and optimality criterion we assumed, and their DFT approach coincides with the proposed least-squares solution.

Thus, in practice, we have a trade-off of SNR gain using an optimal design versus SNR loss by lower gradient amplitudes or longer TE. Currently, the MD-RME design and its extension to 12 and 24 acquisitions with sign-flips provide the best of both worlds and are straightforward to implement. Specifically, the 12 acquisition sign-flips sequence decreases estimation errors by pffiffiffi 2 compared to a conventional sequence of 24 acquisitions, with half of the acquisition time. We are also looking for other approaches to achieve optimal designs within our framework by allowing additional parametrization of the acquisition design, e.g., including varying fractional encoding ratios in each direction as in (14). We limited our experiments to use spin-echo echo-planar imaging sequences. However, other pulse sequences in the literature, e.g., (9), can be used to further improve SNR and/or shorten acquisition time. Trapezoidal waveform can also be applied rather than sinusoidal, to maximize encoding efficiency under limitation of the slew rate. Note that applying multidirectional MEGs introduces concomitant-field terms that should be corrected, e.g., using gradient symmetrization (29). In the sequences we use, the MEGs are mirrored after the echo pulse and cancel such effects. We note that modeling the noise in the phase images as white Gaussian noise is justified by the fact that the least-squares solution for the conventional sequence using our approach coincides with the standard DFT processing. Thus, such noise has been implicitly assumed so far in MRE imaging. However, in practice, it might not be the case due to, e.g., change in temperature of the coils between measurements. An advantage of our framework is that we may use a different noise model with its proper optimal estimator [e.g., weighted leastsquares (21)]. In addition, optimality criteria other than D-optimality can be considered for evaluating design performance. Our approach may allow to adapt the acquisition matrix accordingly by, e.g., using nonuniform sampling of the acquisitions. The flexibility in shaping the matrix may also be important in incorporating compressed sensing in MRE. In conclusion, the formulation of motion encoding as both linear estimation and optimal design problems allows using tools from fields such as optimization, statistics, and information theory. A careful design of

Optimization-Based Design of Motion Encoding in MRE

parameterized multidirectional MEGs yields displacement maps with a better SNR and/or shorter acquisition times. Along with the ongoing improvement in transducers and reconstruction algorithms, accurate elastograms that characterize the mechanical properties of in vivo tissue can be derived. In turn, these may pave the way for a broader clinical deployment of MRE. ACKNOWLEDGMENT Support from Professor Salcudean’s C. A. Laszlo Chair is gratefully acknowledged. APPENDIX Derivation of the Parameterized Phase Accumulation Here we derive expressions of the parameterized phase accumulation. Based on the Larmor precession, we calculate the phase accumulated by each isochromat at r  under that moves according to the wave displacement C  the MEGs G (excluding the constant phase accumulated  0 ). We have that by the static magnetic field B  ; jGj; a  ¼g Fðr ; N  ; bÞ

Z

1

 v  ; jGj; a   r ðtÞ dt Gðt; ;N  ; bÞ

Z1 1

   r ; tÞ dt  v  ; jGj; a   r þ Cð Gðt; ;N  ; bÞ ¼g 1 X Fi ðr ; vi ; Ni ; jGi j; ai ; bi Þ ¼ i¼x;y;z

  r are cancelled Using Eqs. [1] and [2], the elements in G (integration of a constant multiplied by sinusoid over its time period). For i ¼ x; y; z, we have the individual phase accumulations Fi ðr ; vi ; Ni ; jGi j; ai ; bi Þ ¼ g Z

Z

ai =vi þNi Ti

jji ðr Þjsinðvm t ai =vi

þ wi ðr ÞÞjGi jsinðvi t  ai þ bi Þ dt ai =vi þNi Ti

g ¼ jGi jjji ðr Þj ½cosððvm  vi Þt þ wi ðr Þ þ ai  bi Þ 2 ai =vi  cosððvm þ vi Þt þ wi ðr Þ  ai þ bi Þ dt In case vi ¼ vm , the second cosine is cancelled in the integration and we have Fi ðr ; Ni ; jGi j; ai ; bi Þ ¼

gpNi jGi j jji ðr Þjcosðwi ðr Þ þ ai  bi Þ vm

REFERENCES 1. Muthupillai R, Lomas DJ, Rossman PJ, Greenleaf JF, Manduca A, Ehman RL. Magnetic resonance elastography by direct visualization of propagating acoustic strain waves. Science 1995;269:1854–1857. 2. Glaser KJ, Manduca A, Ehman RL. Review of MR elastography applications, recent developments. J Magn Reson Imaging 2012;36:757– 774. 3. Manduca A, Oliphant TE, Dresner MA, Mahowald JL, Kruse SA, Amromin E, Felmlee JP, Greenleaf JF, Ehman RL. Magnetic resonance elastography: non-invasive mapping of tissue elasticity. Med Image Anal 2001;5:237–254. 4. Oliphant TE, Manduca A, Ehman RL, Greenleaf JF. Complex-valued stiffness reconstruction for magnetic resonance elastography by alge-

1525 braic inversion of the differential equation. Magn Reson Med 2001; 45:299–310. 5. Sinkus R, Tanter M, Catheline S, Lorenzen J, Kuhl C, Sondermann E, Fink M. Imaging anisotropic, viscous properties of breast tissue by magnetic resonance-elastography. Magn Reson Med 2005;53:372–387. 6. Papazoglou S, Rump J, Braun J, Sack I. Shear wave group velocity inversion in MR elastography of human skeletal muscle. Magn Reson Med 2006;56:489–497. 7. Honarvar M, Sahebjavaher RS, Salcudean SE, Rohling R. Sparsity regularization in dynamic elastography. Phys Med Biol 2012;57:5909– 5927. 8. Wang H, Weaver JB, Doyley MM, Kennedy FE, Paulsen KD. Optimized motion estimation for MRE data with reduced motion encodes. Phys Med Biol 2008;53:2181–2196. 9. Garteiser P, Sahebjavaher RS, Ter Beek LC, Salcudean SE, Vilgrain V, Van Beers BE, Sinkus R. Rapid acquisition of multifrequency multislice, multidirectional MR elastography data with a fractionally encoded gradient echo sequence. NMR Biomed 2013;26:1326–1335. 10. Rump J, Klatt D, Braun J, Warmuth C, Sack I. Fractional encoding of harmonic motions in MR elastography. Magn Reson Med 2007;57: 388–395. 11. Herzka DA, Kotys MS, Sinkus R, Pettigrew RI, Gharib AM. Magnetic resonance elastography in the liver at 3 Tesla using a second harmonic approach. Magn Reson Med 2009;62:284–291. 12. Pelc NJ, Bernstein MA, Shimakawa A, Glover GH. Encoding strategies for three-direction phase-contrast MR imaging of flow. J Magn Reson Imaging 1991;1:405–413. 13. Hausmann R, Lewin JS, Laub G. Phase-contrast MR angiography with reduced acquisition time: new concepts in sequence design. J Magn Reson Imaging 1991;1:415–422. 14. Yasar TK, Klatt D, Magin RL, Royston TJ. Selective spectral displacement projection for multifrequency MRE. Phys Med Biol 2013;58: 5771–5781. 15. Klatt D, Yasar TK, Royston TJ, Magin RL. Sample interval modulation for the simultaneous acquisition of displacement vector data in magnetic resonance elastography: theory, application. Phys Med Biol 2013;58:8663–8675. 16. Sack I, Mcgowan CK, Samani A, Luginbuhl C, Oakden W, Plewes DB. Observation of nonlinear shear wave propagation using magnetic resonance elastography. Magn Reson Med 2004;52:842–850. 17. Kailath T, Sayed AH, Hassibi B. Linear estimation. New Jersey: Prentice Hall; 2000. 854 p. 18. Kiefer J, Wolfowitz J. Optimum designs in regression problems. Ann Math Stat 1959;30:271–294. 19. Pukelsheim F. Optimal design of experiments. New York: Wiley; 1993. 454 p. 20. John RC, Draper NR. D-optimality for regression designs: a review. Technometrics 1975;17:15–23. 21. Boyd S, Vandenberghe L. Convex optimization. Cambridge, UK: Cambridge University Press; 2004. 727 p. 22. Kiefer J. Optimum designs in regression problems II. Ann Math Stat 1961;32:298–325. 23. Hadamard J. R esolution d’une question relative aux d eterminants. Bull Sci Math 1893;17:240–246. 24. Galil Z, Kiefer J. D-optimum weighing designs. Ann Stat 1980;8: 1293–1306. 25. Coppel WA. Hadamard’s determinant problem. In: Coppel WA, editor. Number theory. New York: Springer; 2009. pp 223–259. 26. Baghani A, Salcudean SE, Honarvar M, Sahebjavaher RS, Rohling R, Sinkus R. Travelling wave expansion: a model fitting approach to the inverse problem of elasticity reconstruction. IEEE Trans Med Imaging 2011;30:1555–1565. 27. Baghani A, Eskandari H, Wang W, Da Costa D, Lathiff MN, Sahebjavaher RS, Salcudean SE, Rohling R. Real-time quantitative elasticity imaging of deep tissue using free-hand conventional ultrasound. In: Ayache N, Delingette H, Golland P, Mori K, editors. Medical image computing and computer-assisted intervention—MICCAI. Berlin-Heidelberg: Springer; 2012. pp 617–624. 28. Mood AM. On Hotelling’s weighing problem. Ann Math Stat 1946;17: 432–446. 29. Bernstein MA, Zhou XJ, Polzin JA, King KF, Ganin A, Pelc NJ, Glover GH. Concomitant gradient terms in phase contrast MR: analysis and correction. Magn Reson Med 1998;39:300–308.