Optimal control design of turbo spin-echo sequences with applications to parallel-transmit systems.

PURPOSE: The design of turbo spin-echo sequences is modeled as a dynamic optimization problem which includes the case of inhomogeneous transmit radiofrequency fields. This problem is efficiently solved by optimal control techniques making it possible to design patient-specific sequences online. THEORY AND METHODS: The extended phase graph formalism is employed to model the signal evolution. The design problem is cast as an optimal control problem and an efficient numerical procedure for its solution is given. The numerical and experimental tests address standard multiecho sequences and pTx configurations.
RESULTS: Standard, analytically derived flip angle trains are recovered by the numerical optimal control approach. New sequences are designed where constraints on radiofrequency total and peak power are included. In the case of parallel transmit application, the method is able to calculate the optimal echo train for two-dimensional and three-dimensional turbo spin echo sequences in the order of 10 s with a single central processing unit (CPU) implementation. The image contrast is maintained through the whole field of view despite inhomogeneities of the radiofrequency fields.
CONCLUSION: The optimal control design sheds new light on the sequence design process and makes it possible to design sequences in an online, patient-specific fashion. Magn Reson Med 77:361-373, 2017.

INTRODUCTION

Turbo spin‐echo (TSE) sequences 1 based on the Carr–Purcell–Meiboom–Gill (CPMG) condition 2, 3 are extensively used in current clinical MRI exams. Originally intended as a rapid succession of refocusing pulses 4, the spin‐echo technique has been applied with trains of lower refocusing angles 5, considerably reducing the power deposition in the tissue and also allowing for much longer echo trains.

Previous work has addressed the design of optimal, low refocusing angle trains 6, 7, 8. In particular, Ref. 6 shows how in some cases the angles can be analytically derived for a number of multiecho sequence types such as hyperechoes 9 and transition between pseudo steady states (TRAPS) 10. Other work has focused on the optimization of the tip angle trains, usually in a heuristic fashion 11, 12, 13.

An important extension of the standard sequence design is presented in 14, 15, where the parallel transmit (pTx) system characteristics are taken into account for patient specific calculations. At high and ultra high field MRI, the Larmor frequency is such that interference and penetration effects for radiofrequency (RF) waves occur in the body. As a consequence, the RF transmit field is no longer homogeneous, leading to inconsistent contrast and signal loss over the field of view. pTx systems consist of independent and simultaneously superimposed RF fields. Many researchers have proposed methods for using the new degrees of freedom available from pTx 16, 17, 18. In the context of improving signal homogeneity, these may be divided into those which aim to improve uniformity of the RF field directly like RF shimming 19; those which aim to improve uniformity of properties of individual RF pulses using multidimensional RF pulses [e.g., “spokes” 20, “k T‐points” 21, 22, 23 or “SPINS” 24]; or finally those which seek to directly consider the resulting signals 14, 15. The latter methods have been proposed for TSE sequences since they do not require multidimensional pulses that may increase the achievable interecho spacings. Instead, signal evolution over the entire pulse sequence is considered. This approach differs from other recent work on TSE sequences 23 which must take care to design pulses whose phases are matched to obey the CPMG conditions as Massire et al. mentioned in the discussion of 23. By directly considering the resulting echo amplitudes, this is implicitly accounted for.

Previous work on signal control has used simple numerical optimization to compute optimized settings leading to long computation times. This is a weakness since these calculations must be carried out while the subject is in the scanner. In particular, use of gradient based optimization with finite difference approximations is computationally costly. Whilst later work has employed a nongradient based approach 15, the convergence of such methods remains slow since the structure of the minimization problem is not exploited.

In this work, a fundamentally different approach is taken, in which we analyze the structure of the extended phase graph (EPG) algorithm and provide insights in the modeling of the sequence design process. The EPG is seen as a discrete‐time dynamical system. The design problem is cast as the minimization of a smooth function under smooth equality and/or inequality constraints. As all functions involved in the model are differentiable, the resulting numerical optimal control problem can be solved by efficient derivative‐based methods. As a consequence, the computation time to determine the optimal sequence settings is such that online clinical implementation becomes feasible.

The exact derivatives of all functions playing a role in the model are presented below. For the functions which implicitly depend on the design parameters, the adjoint states method (ASM) is used. This algorithmic differentiation method makes it possible to efficiently calculate the exact derivatives of a dynamical process 25. Note that ASM was previously applied to large‐tip angle RF pulses design for pTx systems in 26.

We will give several examples of optimal control EPG‐based design. First, we will show how standard, analytically derived sequences are recovered by the algorithm. Among these cases, the derivation of the CPMG condition through numerical optimization is included. Next, we will design new sequences and finally we will apply the same framework to patient specific, pTx configurations. In the latter case, application of optimal control with ASM allows online, patient specific computations as the total computation time is in the order of 10 s on a standard desktop PC with a single CPU implementation. Finally, the optimized pTx TSE sequence will be experimentally tested on a phantom and on a volunteer's head with an eight‐channel transmit system at 7T MRI.

The new optimal control framework is a generalized and flexible approach to sequence design and can efficiently calculate sequences on a patient‐specific basis; it is useful as a support for MR scientists and as an online tool for high‐field pTx MRI exams.

This code is made available online at https://github.com/mriphysics/optimal-control-EPG.

THEORY

The Extended Phase Graph and Its Extension to pTx Configurations

We will start by giving a basic description of the EPG. More detailed information can be found in the review article of Weigel 27.

The EPG describes the evolution of the signal in terms of configuration states . For sequences with equidistant timing (fixed T R and T E), they are defined as the Fourier coefficients of the complex magnetization states : with:

In the above equations, quantifies the dephasing due to off‐resonance caused by applied gradient fields and .

The echo value is given by the two F 0 states. Note the complex conjugate relationship: .

The Bloch equation dictates the dynamics of , which can be decomposed into rotation, relaxation and dephasing effects. In the EPG formalism, these three phenomena are given respectively by: the rotation operator, , the relaxation operator , and the shift operator , which are defined as:

In the previous equations, α, , T 1, T 2, and τ represent, respectively, the tip angle amplitude, the corresponding phase, the longitudinal and transverse relaxation rates, and the echo spacing. In the case of variable refocusing angles, we write α and where the subscript n indicates the pulse number.

The EPG calculation is obtained by recursively applying the three operators to the configurations states .

Spatially Resolved EPG and pTx Configurations

In the case of a L‐channel pTx system, the effective amplitude α and phase are derived from the sum of all complex RF pulse weights multiplied by the sensitivity of the corresponding transmit coil: where denotes the complex transmit field of the ‐th channel and and represent, respectively, the real and imaginary part of the flip‐angle applied to the ‐th transmit channel at the n‐th pulse. Note that is dimensionless. For implementation reasons (see the Appendix), it is convenient to use the real/imaginary part representation for the complex RF weights of the pTx system. Needless to say, the and (x, y) representations are interchangeable. RF pulse amplitudes are given in terms of nominal flip angle: that is, the input to the system. The flip angle effectively achieved by the RF system is the nominal value multiplied by the spatial distribution of the corresponding channel.

The transmit fields are usually spatially dependent, a fact that is reflected in the spatial response of the echo train. Therefore, the EPG needs to be evaluated in each voxel of the region of interest, leading to the spatially resolved version of the EPG 14.

Sequence Design as an Optimization Problem

The EPG can be formulated as a discrete‐time dynamical process: where f is the vector of all concatenated states at the n‐th time point:

P represents the transition matrix, which depends on the sequence and tissue parameters and T 2. K is the maximum number of configuration states.

For example, represents dephasing followed by decay and refocusing. The matrices S, E and R are constructed by concatenating the refocusing, decay and rotation matrices from the previous section into block‐diagonal matrices (see the Appendix for more implementation details). For fixed and T 2, the free parameters are and . Finally, b represents the initial state (typically all zero components but the one corresponding to ).

The time index n runs from 0 to N and represents the echo numbers, in particular: n = 0 is the initial or equilibrium state, n = 1 is the state right after the excitation and n = 2 is the first echo. For the largest coefficient, K, in Eq. (8) we have K = N. Note that Malik et al. show in 15 that for large n, the effect of the large coefficients k is small. In this work, we have chosen to truncate the series in k such that . This allows for a considerable acceleration in the computations.

The signal at echo‐time n – 1 is given by the first and the second components of f, which are the conjugate of one another. In this work we will thus consider only the second component.

The sequence design will be cast as a minimization problem of the standard form: where s, p, and q represent, respectively, the objective, the equality, and inequality constraints functions. I and J denote, respectively, the number of equality and inequality constraints. The functions s, p, and q are required to be differentiable.

For example, to design a sequence which maximizes the signal intensity over the whole echo train, we can set:

In this equation, C is a diagonal weighting matrix which selects only the signal component of f (f contains all states), that is:

For example, setting c = 0 for means that the signal from the first M echoes is not taken into account. As default, c = 1, but it can be convenient to assign larger weights to the central k‐space echoes. We can easily implement this by setting, for instance, c = 2 for these particular k‐space indices.

As another example, suppose we want to design a sequence whose signal closely follows a predefined target response, t, and the total RF power has to be minimized. In this case we set where denotes the maximum allowed deviation from the desired target. We only specify the echo component of t and the others are zero.

Upper bounds for the tip angle amplitude can be set in the following way: where denotes the maximum values allowed.

In case of pTx configurations, the objective and/or constraint functions are evaluated over a set of R spatial positions. Making use of the real/imaginary representation, the design problem is easily extended.

The building block functions for the optimal control design are schematically presented in Table 1. Clearly, the spatially resolved version of the design problem can be applied also to single channel transmit configurations.

Table 1

Building Block Functions for the Optimal Control Problem

	Standard setup	pTx configuration
Error function	12∑n=2N(f_n−t_n)^HC_n(f_n−t_n)	12∑r=1R∑n=2N(f_n,r−t_n,r)^HC_n(f_n,r−t_n,r)
Signal amplitude	−12∑n=2NfnHC_nf_n	−12∑r=1R∑n=2Nfn,rHC_nf_n,r
Total RF power	−12∑n=0N−1αn2	−12∑n=0N−1∑ℓ=1Lxn,ℓ2+yn,ℓ2
Peak RF power	α_n−α_maxfor n=0,…,N−1	xn,ℓ2+yn,ℓ2−αmax2for n=0,…,N−1 and ℓ=1,…,L

The subscript indicates the spatial position.

Calculating the Exact Derivatives. The Adjoint States Method

To efficiently solve the optimal control problem given by Eq. (9), the derivative of all functions involved with respect to the design parameters is needed. Calculating the derivatives of the functions which explicitly depend only on the parameters is rather straightforward. For example, we have:

Conversely, the derivatives of the functions which implicitly depend on the design parameters can not be directly determined. Examples thereof are given in Table 1.

In these cases, the derivative could be approximated by finite differences schemes. For example: for a small value of δ. The advantage of this solution is the simplicity of its implementation. However, each derivative requires an EPG simulation. When dealing with a large amount of parameters and voxels (for the spatially resolved EPG), the repeated evaluation of s becomes prohibitively long.

An extremely efficient alternative to solve this problem is offered by the ASM. ASM is a very powerful tool, extensively used in numerical solutions of optimal control problems. It consists of a forward/backward simulation approach to compute exact derivatives. We will show how ASM can be applied to the problem of finding by a backward recurrence strategy.

Suppose we want to calculate the derivatives of

First of all, note that the derivative of s with respect to the last parameter, , is given by: where we used the fact that . Note that the derivatives of are analytically defined.

Defining the ‐th adjoint state, :

The next step consists of calculating : where we used the fact that . Rearranging the right‐hand side terms:

Defining the ‐th adjoint states:

Continuing this process, we obtain the backward recurrence relation for the adjoint states:

Once the adjoint states are found, the derivatives can be easily calculated: and, in an analogous way, we observe that

The derivatives of the transition matrix P in Eqs. (18) and (19) are analytically obtained in the Appendix.

Piecewise Constant Sequence Parameters and Reduction of Degrees of Freedom

In practice, it is useful to reduce the degrees of freedom of the sequence design by constraining α and to be constant over a certain time interval. The consequence is that the design variables are reduced in number, improving convergence performance of the algorithm and making it possible for the scanner interface to deal with fewer input values.

In particular, suppose that we set: where is a set of indices indicating the j‐th range. For example, setting means that

The design parameters for the flip angle amplitudes are then reduced from 14 to 8. In matrix form: which can be more compactly be written as

Analogously, we can write The derivatives with respect to the new variables are easily calculated from the derivatives with respect to α by application of the chain rule:

The same method can be applied when the flip angles are expressed in terms of real and imaginary components.

Computational Complexity of the ASM for EPG

The main computational burden in the optimal control method lays in: (a) the forward recursion scheme for the configuration states; (b) the backward recursion scheme for the adjoints states given by Eq. (17); and (c) the two loops to calculate the products in Eqs. (18) and (19). Note that the last two loops can be combined into a single one, which is twice as costly in computational terms. Each of these four loops requires a similar amount of floating point operations as an EPG simulation. It follows that the computational burden for each iteration approximates four EPG calculations.

As shown in 25, the ASM becomes very attractive when the number of parameters is large. To give an example, in the case of N = 50 excitations train, with both amplitude and phase as variable gives 100 derivatives to be calculated for each iteration step of the minimization algorithm. The finite difference approximating method would require EPG simulations, where N p denotes the possibly reduced number of design parameters (see also the example from the previous paragraph, where N p = 4). If all degrees of freedom are maintained, that is N = N p, then the EPG simulation has to be repeated 101 times. The computational cost for ASM still approximates four EPG simulations.

For a pTx configuration, the signal response is calculated for several voxels in the region of interest, increasing the simulation time of the EPG. For example, in 14 approximately 75 spatial points were used, increasing the number of EPG calculations by a factor of 75. The scale of the problem increases dramatically and the minimization problem needs to be solved online. In the pTx case, the design variables are the real and imaginary component of the pulse for each n and each transmit channel. For an eight‐channel pTx system, the total number of parameters is . The time reduction factor of ASM with respect to finite difference is thus . The application of ASM is fundamental to compute the optimal sequence parameters online.

Details about the practical implementation are given in the Appendix.

METHODS

In this section, we test the optimal control design approach. We run the minimization problem for some sequence designs whose solution is analytically known and compare the outcomes. Afterwards, we show how to easily extend the requirements on the sequence by adding constraints. Finally, we address the pTx, spatially resolved EPG configurations.

In all tests, the minimization problem is solved by the built‐in Matlab implementation of the interior‐point method (function fmincon) with user‐defined gradients of the objective and constraint functions (see also the Appendix). The Hessian is approximated by the Broyden–Fletcher–Goldfarb–Shanno method 28. The spatially resolved EPG simulator is implemented in C++ by making use of the linear algebra library Armadillo 29 on a Linux PC with Intel Xeon CPU, 3.50 GHz. The Matlab and C++ functions are interfaced by means of file streams. The code does not make use of parallelization thus only one CPU core is employed.

Test 1. Maximizing the Signal on a Standard TSE

In the first test, we wish to design a sequence which maximizes the signal, that is:

We choose ms, N = 60. As a starting guess, we set and for . The starting phases, , obey the CPMG condition, that is (excitation pulse) and for (refocusing pulses). We expect to recover the standard turbo‐spin‐echo sequence given by .

Test 2. Constant Signal Intensity with Minimum RF Power

For the second test, we wish to calculate the minimum total power tip angle train for a sequence which returns constant signal intensity given by I c. The analytical solution, omitting relaxation, is given by 6:

The design problem can be cast as: where and σ a small positive number which controls the accuracy of the obtained signal intensity.

We solve the problem for three values of I c, namely: and we set . We omit relaxation by setting the decay operator . As a starting guess, we set for and .

Test 3. Exploiting Flexibility. Peak RF‐Constrained TSE and Recovery of the CPMG Condition

The numerical optimal control approach has the advantage of being very flexible. For instance, assume that we wish to design a maximum signal sequence with refocusing angles smaller than and we are interested only in the signal after the fifth echo. This can be easily implemented as:

We set ms and the echo‐train length is N = 60. The starting candidate is for and .

This test is divided in two parts. In the first part, we require the refocusing pulses to be constant (see the paragraph on piecewise constant sequence parameters in the Theory section). This means that α and for will be represented by a unique numerical value. The excitation pulse is independent from the refocusing pulses. The aim of this numerical experiment is to investigate whether the CPMG condition is recovered by the optimal control approach.

In the second part of this test, we leave complete freedom to each refocusing pulse to vary in time.

Test 4: Eight‐Channel pTx System at 7 Tesla. 3D ROI Optimization on a Phantom

In the next two tests we investigate the sequence design for pTx systems. We consider an eight‐channels transmit system at 7.0 Tesla. A transceive headcoil (Nova Medical, Wilmington, DE) and a 7T scanner (Philips, Best, NL) are used. In the first scanner test, a spherical phantom of 10 cm diameter is used. It contains the following compounds: acetate, ethanol, phosphoric acid, and arquad solution. The approximated T 1 and T 2 values are 500 ms and 400 ms, respectively. The relative maps for the eight channels are derived from low flip‐angle gradient‐echo images, as described in 30. In addition, an absolute mapping measurement is performed with the DREAM method 31 and is used to recover quantitative information about the flip angle; the resulting maps are plotted in Figure 5a. We design optimal amplitude and phase settings for a three‐dimensional‐turbo‐spin‐echo scan (3D‐TSE, also known as SPACE, or VISTA), TSE factor = 113. The sequence should be able to compensate for the RF inhomogeneities in such a way that the signal follows the dynamics of a predefined ideal echo train over the whole region of interest (ROI). By the word “ideal” we indicate a pulse train which is designed under the assumption of a perfectly homogeneous transmit field. In particular, we look at pseudosteady‐state sequences as derived in 32.

Figure 5

Test 4. (a) Amplitude‐transmit maps in the phantom for the eight‐channel system at 7T. Central transverse, coronal and sagittal slices are shown. b: Optimized amplitude and phases for the 3D TSE 113‐echo train (first 52 pulses are shown). The phase values correspond to the off‐set with respect to the standard (circularly polarized) excitation mode. c: Relative standard deviation of the signal intensity over the ROI for the two echo trains.

Constraints on the RF total and peak power are also required. The resulting problem in terms of the real, x, and imaginary, y, parts of the pulse settings is: where denotes the maximum allowed total power. The target signal is given by the EPG response of the echo train given by 32:

The value is repeated 48 times to yield a total of 51 refocusing pulses (plus the excitation pulse). As the magnetization reaches the steady states during the first 51 pulses, there is no need to design also the remaining part of the sequence. The last value of α is repeated until the conclusion of the echo train.

To reduce the degrees of freedom, we force to be constant in the three intervals , and . See also Figure 1 for the corresponding diagram. The number of design parameters is then reduced to and the transmit maps are undersampled into a 3D tetrahedral lattice 15 to yield a total of R = 91 spatial locations.

Figure 1

Diagram indicating the distribution of the optimization settings for test 4. The first 22 pulses are driven individually. The remaining 30 pulses are subdivided in three groups of equal length and each pulse group is represented by a single parameter. Note that the amplitude and phase for each of the eight channels are independently driven, resulting into optimization parameters.

We set and as twice the total power obtained by the nonoptimized sequence, . This choice is empirically defined so as to remain within specific absorption rate constraints dictated by the power monitoring unit of the scanner. The starting guess is with and phase, respectively, for the excitation and refocusing pulses (CPMG condition). The optimal control algorithm is halted after 50 iterations.

The other sequence parameters are: ms; echo‐spacing: 2.9 ms; echo‐train duration: 374 ms; resolution: 1 mm3; total scan duration: 10 min and 11 s.

Test 5: In Vivo Scan with Eight‐Channel pTx System at 7 Tesla. 2D ROI Optimization

In the final test, we address the design of a two‐dimensional (2D) TSE sequence for a volunteer's brain scan. The maps for the eight channels are measured in the same way as in the previous test 30, 31 and they are displayed in Figure 7a. Constraints on peak and total RF power are also required. We obtain the same design problem as in the previous test.

Figure 7

Test 5. (a) Amplitude and phase transmit maps in the volunteer's head for the 8 channel system at 7T. b: Optimized angles for the 2D TSE nine‐echo train. c: Relative standard deviation of the signal intensity for the two echo trains.

The weights on the echoes (i.e., the nonzero entry of C) are all equal to one except for the central k‐space echo, where we set it equal to 2.

The target signal intensity is given by the ideal (no inhomogeneity) EPG signal response of the echo train (TSE factor = 9):

The 2D maps are undersampled by a factor 4 in both directions and only voxels corresponding to the cerebral region are taken into account. This step is automatically performed by the Brain Extraction Tool software 33. The resulting total number of voxels for the design problem is R = 104. The (T 1, T 2) reference values are set equal to the approximate relaxation times of White matter, that is 1500 ms and 100 ms, respectively 34. As for Test 4, the optimal control algorithm is halted after 50 iterations and the starting guess is with phase off‐set between the excitation and refocusing pulses.

The other sequence parameters are: ms; echo‐spacing: 12 ms; echo‐train duration: 108 ms; resolution: 0.3 mm2; total scan duration: 7 min and 45 s.

We will compare the proposed method with the standard sequence run in circularly polarized mode.

RESULTS

The starting and optimized trains with the corresponding EPG signal responses are shown in Figure 2. As expected, the optimal control design recovers the standard spin‐echo sequence.

Figure 2

Test 1. Pulse amplitude (top) and signal intensity (bottom) for signal maximization without constraints on the RF. The phases for both pulse‐trains obey the CPMG condition. Note that the obtained train is the standard turbo spin echo.

The obtained refocusing angles for the constant signal intensity calculations are shown in Figure 3. Note that the analytical and algorithmic values coincide. We point out that the amplitudes are as shown in Figure 3 and the phases are unchanged from the CPMG phases.

Figure 3

Test 2. The refocusing angles for constant signal intensities I c obtained analytically (circles) and by optimal control design (continuous line). The two sets coincide.

Test 3. Exploiting Flexibility. Peak RF Constrained TSE and Recovery of the CPMG Condition

Setting a maximum of on the refocusing angle, we obtain the pulse trains as shown in Figure 4.

Figure 4

Test 3. Maximization of the signal when the refocusing angles are constrained by . Constant and time‐varying refocusing angles are investigated. a: Amplitude (α) and (b) phase ( ) modulation of the pulses. In part (a), the amplitude settings for the starting and the optimal constant refocusing angle setups coincide. Note the recovered CPMG condition in the optimal, constant refocusing angles setup. The excitation and refocusing angles have a phase off‐set (see arrow). c: Signal amplitude and (d) phase for each echo obtained from the three pulse trains. Note the increased signal intensity for the time‐varying with respect to the constant refocusing angle trains. The price to be paid is an additional phase modulation in the refocusing pulses, which is reflected in the phase of the signal.

For the constant refocusing angle setup, the optimized pulse train has the maximum allowed amplitude, that is, refocusing angle and a phase off‐set between the excitation and refocusing pulses (see the arrow in Fig. 4b). This is equivalent to the CPMG condition, which is thus recovered by the algorithm.

For the varying refocusing angle setup, the signal amplitude is the largest. Note in the latter case, the phase modulation of the optimized train. The consequence is an increase of the signal intensity with respect to the sequence where no phase modulation is present. The phase modulation given by the time‐varying refocusing pulses has to be accounted for prior to image reconstruction.

The computation time for the 8ch pTx 3D‐TSE sequence is about 45 s. The optimal amplitude and phase trains are shown in Figure 5b. The standard deviations of the signal intensities over the whole ROI divided by the mean value for each echo are plotted in Figure 5c. Figure 6a shows the simulated signal intensities for the central k‐space echo in four different transverse slices. The MR images obtained with this setup are shown in Figure 6b. These images are divided by the receive sensitivity profile and cropped to include only the phantom. In Figure 6c, the central vertical and horizontal image profiles are displayed. Note the improved homogeneity of the signal. Comparing Fig. 6a,b, we see that the predicted excitation patterns from the EPG simulations closely resemble the scanner images.

Figure 6

Test 4. Phantom 3D experiment. a: Simulated signal intensity for the 3D phantom experiment. Images correspond to four transverse slices separated by 20 mm. b: MR images for the 4 transverse slices. The images are scaled by the receive sensitivity profile. c: Profiles of the MR images taken along the central vertical and horizontal lines.

Test 5: Eight‐Channel pTx System at 7 Tesla. 2D ROI Optimization, In Vivo Scanner Experiment

The computation time for the eight‐channel pTx 2D‐TSE sequence is about 8 s. The optimal amplitude and phase trains are shown in Figure 7b.

The relative standard deviations of the signal intensities over all brain voxels for each echo are plotted in Figure 7c. Note the improved homogenization level. The simulated signal amplitudes are shown in Figure 8. Echo number 5 corresponds to the central k‐space line acquisition (indicated by the box). The excitation fidelity is improved by the optimal control echo train.

Figure 8

Test 5. Simulated echo amplitudes in the volunteer's head for the eight‐channel system at 7T. The echo corresponding to the central k‐space line is indicated by the red box. Top: standard quadrature excitation. Bottom: signal obtained with the optimized pTx settings. Note the improved homogeneity in the temporal regions (see arrows).

The in vivo images are shown in Figure 9. The contrast is more homogeneous over the whole brain and the signal is recovered also in the temporal regions (see arrows).

Figure 9

Test 5. In vivo MR experiment. Left: Image obtained by the standard 2D TSE sequence. The temporal lobes suffer of inhomogeneous contrast and signal loss (see arrows). Right: Image obtained with the Optimal Control pTx train. As predicted from the simulations, the contrast is maintained in the temporal lobes.

DISCUSSION

We have presented a generalized and efficient numerical approach to the design of optimal pulse amplitudes and phases in fast spin‐echo sequences. The design problem is modeled as optimal control of differentiable functions. The resulting optimization procedure can thus be solved by efficient derivative‐based algorithms such as, for instance, the interior‐points method. Convergence to the known analytical solution is obtained for the simple cases where an analytic solution is available. Furthermore, the CPMG condition is shown to be recovered in test 3. This fact illustrates the robustness and validity of the presented approach. Additionally, we have shown how the flexibility of the algorithm can be exploited to design new sequences, revealing new insights in the field of sequence optimization. In particular, the sequence derived in test 3, where the refocusing angles are constrained by , shows how a specific phase modulation can considerably increase the signal intensity. We believe that many more examples can be found where the interplay between varying amplitudes and phases lead to better sequences.

In the case of pTx systems at high fields MRI, the online optimization of the pulse train becomes necessary since the sensitivity profile of the RF transmit array is exam‐dependent. The numerical optimal control approach can work in a clinical environment thanks to the application of the ASM. This powerful algorithmic differentiation technique allows for time reduction factors in the order of 100 with respect to brute‐force finite differencing. The resulting computation time for a 3D TSE sequence and eight‐channel pTx system is in the order of 10 s, with a standard PC, single thread. The design problem can be easily solved by making use of multithreading implementation. Since the spatially resolved EPG and adjoint states simulations are voxel‐dependent, the calculations can be carried out simultaneously for different groups of voxels assigning each group to a separate thread. The acceleration factor in the computing time should be close to the number of threads used. As a result, the actual computing time for the most demanding sequences could be further reduced to few seconds. Furthermore, the derivatives obtained by ASM are exact, which means that the method is more robust to experimental imperfections and numerical errors.

The flexibility and speed of the proposed approach is exploited to homogenize the spatial response of the spin‐echo over the whole field of view in an eight‐channel pTx system at 7T. The results show a reduced standard deviation of the signal intensity for all echoes. The obtained amplitude and phase trains are limited by peak and total RF power but they still produce noticeable improvement in the in vivo image contrast as illustrated by the scanner experiments. The design problem could explicitly be formulated to maintain a uniform signal difference between tissue types. During the developing process, we did that by optimizing the response over three different tissue types. The results obtained with this approach were analogous to the ones obtained by the single tissue approach, but the computation time increased by a factor of 3. Since the contrast uniformity is a by‐product of the single‐tissue design approach, we decided to adopt this strategy. Furthermore, the amplitude sweep of each channel is rather mild, in contrast to the pulse trains obtained in 15. This should have a positive effect regarding the robustness to maps inaccuracies, as already anticipated in 6. As the scanner experiments were successful, we believe that the optimal control pulse trains are robust to mapping errors. A thorough stability analysis of the optimal control method could be the subject of further research.

We would like to point out that the procedure outlined in this article might find a local minimum of the design problem. Given the improvement in the objective, and thus in the signal homogeneity, we have shown that even a local minimum solution has attractive properties.

The approach we followed for sequence design differs from the majority of existing methods such as spokes 20, multidimensional RF pulse design 16, 17, 18, 3D spiral nonselective RF pulses 24 and k T points 21, 22, 23 which consider the flip angle created by each RF pulse separately. TSE sequences with variable flip angles have been approached by designing multidimensional refocusing pulses that maintain the CPMG condition, either scaling them in amplitude 22 or adaptively redesigning for specific target flip angles 23. Our work follows a different approach, using simple hard pulses for refocusing, with optimization focusing on the overall state of the magnetization. This approach could be seen as signal control by dynamic shimming. By avoiding multidimensional pulses, we simplify the resulting pulse sequences and shorten the achievable echo spacing. The results shown in this article indicate that major improvement in the sequence response can be achieved by acting only on the amplitude and phase settings. The approach could however be complementary to approaches employing multidimensional pulses, with the added benefit of relaxing the constraint from the CPMG condition, since pulses are directly designed to result in stable echo amplitudes.

We have shown how to apply optimal control techniques to the EPG formalism. Optimal control of the Bloch equation is a rather developed field and with this work we hope to inspire future synergies between the two domains. One example thereof could be the combination of Bloch‐equation/EPG approaches to jointly design the excitation pulse (by Bloch equation methods) and the resulting dynamic shimming settings (by EPG methods). This combination can be implemented as a unified numerical optimal control problem. We think that the hybrid design approach Bloch equation/EPG can pave the way to new developments in the sequence design field.

In terms of safety, the sequences we have designed are constrained by peak and total RF power. The same formalism can be used to include constraints on global specific absorption rate and local specific absorption rate 35, 36. We expect that the computation time will increase by adding a large amount of local specific absorption rate constraints. Implementation on parallel computing architectures will make online application possible as is the case for multidimensional RF pulse design 37.

In this work, we have focused on the design of tip angles and corresponding phases. This is not a restriction for the algorithmic approach, since other parameters, for instance the echo spacing, τ, could be subject to optimization too. The framework includes the possibility to treat also τ as a variable and the derivatives with respect to τ can be calculated by ASM. This could be exploited to increase the contrast between different tissue species. An investigation of this aspect goes beyond the scope of this article and we leave it for future work.

The recently developed quantitative MR reconstruction technique based on the EPG signal model 38 has been implemented with a derivative‐free based method 39. The insights provided in our work will make it possible to employ more efficient derivative‐based algorithms for the same reconstruction paradigm. In particular, ASM can be applied for the calculation of the derivatives with respect to parameters as , T 1, and T 2 in the same way as we have obtained the derivatives with respect to the pulses amplitude and phase. In general, we believe that ASM can pave the way for a broader application of EPG‐based reconstruction techniques in MRI.

CONCLUSION

The design of multi spin‐echo sequences is cast as a solution of a minimization problem. The functions involved are differentiable and, by application of the ASM, the exact derivatives can be quickly calculated. This allows one to solve the design problem with standard and efficient derivative‐based methods. The flexibility and efficiency of the optimal control framework can be exploited to design new sequences and to optimize the MR exam in a patient‐specific, online fashion.