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. © 2016 The Authors Magnetic Resonance in Medicine published by Wiley Periodicals, Inc. on behalf of International Society for Magnetic Resonance in Medicine


INTRODUCTION
Turbo spin-echo (TSE) sequences (1) based on the Carr-Purcell-Meiboom-Gill (CPMG) condition (2,3) are exten-sively used in current clinical MRI exams. Originally intended as a rapid succession of 180 o 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 EPGbased 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.

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 ðF þ k ; F À k ; Z k Þ T . For sequences with equidistant timing (fixed T R and T E ), they are defined as the Fourier coefficients of the complex magnetization states ðM þ ; M À ; M z Þ T : Z k e iku : In the above equations, u 2 ½0; 2p quantifies the dephasing due to off-resonance caused by applied gradient fields and i ¼ ffiffiffiffiffiffi ffi À1 p . The echo value is given by the two F 0 states. Note the complex conjugate relationship: F þ k ¼ ðF À Àk Þ Ã . The Bloch equation dictates the dynamics of ðM þ ; M À ; M z Þ T , which can be decomposed into rotation, relaxation and dephasing effects. In the EPG formalism, these three phenomena are given respectively by: the rotation operator, R, the relaxation operator E, and the shift operator S, which are defined as: S : F k 7 !F kþ1 and Z k 7 !Z k : [5] In the previous equations, a, f, T 1 , T 2 , and s 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 a n and f n where the subscript n indicates the pulse number.
The EPG calculation is obtained by recursively applying the three operators to the configurations states ðF þ k ; F À k ; Z k Þ T .

Spatially Resolved EPG and pTx Configurations
In the case of a L-channel pTx system, the effective amplitude a n and phase f n are derived from the sum of all complex RF pulse weights multiplied by the B þ 1 sensitivity of the corresponding transmit coil: B ' Ã ðx n;' þ iy n;' Þ and B ' Ã ðx n;' þ iy n;' Þ ( ) [6] where B ' denotes the complex transmit field of the '-th channel and x n;' and y n;' represent, respectively, the real and imaginary part of the flip-angle applied to the '-th transmit channel at the n-th pulse. Note that B ' 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 ða; fÞ 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 B þ 1 < 1 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 n is the 3ðK þ 1Þ Â 1 vector of all concatenated states at the n-th time point: P n represents the transition matrix, which depends on the sequence and tissue parameters a n ; f n ; t; T 1 and T 2 . K is the maximum number of configuration states.
For example, P n ¼ Rða n ; f n ÞEðt; T 1 ; T 2 ÞS 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 t; T 1 and T 2 , the free parameters are a ¼ ða 0 ; . . . ; a NÀ1 Þ and f ¼ ðf 0 ; . . . ; f NÀ1 Þ. Finally, b represents the initial state (typically all zero components but the one corresponding to Z 0 ¼ 1).
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 K ¼ minfN; 25g. 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 n , 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: Minimize sða; fÞ s:t: p i ða; fÞ ¼ 0 ði ¼ 1; . . . ; IÞ q j ða; fÞ 0 ðj ¼ 1; . . . ; J Þ 8 > > < > > : [9] 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: f H n C n f n : [10] In this equation, C n is a diagonal weighting matrix which selects only the signal component of f n (f n contains all states), that is: C n ðj; jÞ ¼ c n 6 ¼ 0 if j ¼ 2 and C n ðj; jÞ ¼ 0 otherwise: For example, setting c n ¼ 0 for n M means that the signal from the first M echoes is not taken into account. As default, c n ¼ 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 n ¼ 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 ðf n À t n Þ H C n ðf n À t n Þ " # À s 2 [12] where s 2 denotes the maximum allowed deviation from the desired target. We only specify the echo component of t n and the others are zero. Upper bounds for the tip angle amplitude can be set in the following way: where a max 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.

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: @q @x n;' ¼ 2x n;' and @q @y n;' ¼ 2y n;' [14] 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: @s @a j % sða 0 ; a 1 ; . . . ; a j þ d; a jþ1 . . .Þ À sða 0 ; a 1 ; . . . ; a j ; a jþ1 . . .Þ d for a small value of d. 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 @s @aj by a backward recurrence strategy.
Suppose we want to calculate the derivatives of ðf n À t n Þ H C n ðf n À t n Þ: First of all, note that the derivative of s with respect to the last parameter, a NÀ1 , is given by: where we used the fact that f N ¼ P NÀ1 f NÀ1 . Note that the derivatives of P NÀ1 are analytically defined.
Defining the ðN À 1Þ-th adjoint state, l NÀ1 : The next step consists of calculating @s @aNÀ2 : where we used the fact that f N ¼ P NÀ1 P NÀ2 f NÀ2 . Rearranging the right-hand side terms: Defining the ðN À 2Þ-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: @s @a n ¼ l n @P n @a n f n [18] ðf n;r À t n;r Þ H C n ðf n;r À t n;r Þ Signal amplitude Peak RF power a n À a max for n ¼ 0; . . . ; N À 1 x 2 n;' þ y 2 n;' À a 2 max for n ¼ 0; . . . ; N À 1 and ' ¼ 1; . . . ; L The subscript r ¼ 1; . . . ; R indicates the spatial position.
364 Sbrizzi et al. and, in an analogous way, we observe that @s @f n ¼ l n @P n @f n f n ; @s @x n;' ¼ l n @P n @x n;' f n and @s @y n;' ¼ l n @P n @y n;' f n : [19] The derivatives of the transition matrix P n 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 a n and f n 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: 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 a ¼ Qa: Analogously, we can write f ¼ Qb: The derivatives with respect to the new variables are easily calculated from the derivatives with respect to a by application of the chain rule: @s @a ¼ @s @a @a @a ¼ @s @a Q and @s @b ¼ @s @f @f @b ¼ @s @f Q: 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 2 Â N p þ 1 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 2 Â 8 Â N p . The time reduction factor of ASM with respect to finite difference is thus ð2 Â 8 Â N p Þ=4 ¼ 4N p . 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: f H n C n f n : [22] We choose ðT 1 ; T 2 Þ ¼ ð1000; 150Þ ms, N ¼ 60. As a starting guess, we set a 0 ¼ 90 o and a n ¼ 60 o for n ¼ 1; . . . ; N.
The starting phases, f n , obey the CPMG condition, that is f 0 ¼ 0 o (excitation pulse) and f n ¼ 90 o for n ! 1 (refocusing pulses). We expect to recover the standard turbo-spin-echo sequence given by a ¼

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: ðf n À t n Þ H C n ðf n À t n Þ À s 2 0 8 > > > > > < > > > > > : [23] where t n ð2Þ ¼ I c and r a small positive number which controls the accuracy of the obtained signal intensity. We solve the problem for three values of I c , namely: I c 2 f0:3; 0:6; 0:9g and we set s ¼ 10 À3 . We omit relaxation by setting the decay operator E ¼ I. As a starting guess, we set a 0 ¼ 90 o ; a n ¼ 18 o for n ¼ 1; . . . ; N À 1 and f n ¼ 0.

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 60 o and we are interested only in the signal after the fifth echo. This can be easily implemented as: f H n C n f n s:t: a n 60 o n ¼ 1; . . . ; N À 1
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 a n and f n for n ! 1 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 B þ 1 maps for the eight channels are derived from low flip-angle gradientecho images, as described in (30). In addition, an absolute B þ 1 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).
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: ðf n;r À t n;r Þ H C n ðf n;r À t n;r Þ ð signal-homogenizationÞ s:t: where P max denotes the maximum allowed total power. The target signal t n;r is given by the EPG response of the echo train given by (32): The value 60 o 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 a is repeated until the conclusion of the echo train.
To reduce the degrees of freedom, we force ðx n;' ; y n;' Þ to be constant in the three intervals 23 n 32; 33 n 42 , and 43 n 52. See also Figure 1 for the corresponding diagram. The number of design parameters is then reduced to 25 Â 8 Â 2 ¼ 400 and the transmit maps are undersampled into a 3D tetrahedral lattice (15) to yield a total of R ¼ 91 spatial locations.
We set a max ¼ 270 o and P max as twice the total power obtained by the nonoptimized sequence, a standard . This choice is empirically defined so as to remain within spe-cific absorption rate constraints dictated by the power monitoring unit of the scanner. The starting guess is a standard with 0 o and 90 o phase, respectively, for the excitation and refocusing pulses (CPMG condition). The optimal control algorithm is halted after 50 iterations.

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 twodimensional (2D) TSE sequence for a volunteer's brain scan. The B þ 1 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.
The weights on the echoes (i.e., the nonzero entry of C n ) 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 B þ 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 a standard with 90 o phase off-set between the excitation and refocusing pulses. The other sequence parameters are: ðT E ; T R Þ ¼ ð60; 3500 Þ ms; echo-spacing: 12 ms; echo-train duration: 108 ms; resolution: 0.3 mm 2 ; total scan duration: 7 min and 45 s.
We will compare the proposed method with the standard sequence run in circularly polarized mode.

Test 1. Maximizing the Signal on a Standard TSE
The starting and optimized ða; fÞ trains with the corresponding EPG signal responses are shown in Figure 2. As expected, the optimal control design recovers the standard 90 o =180 o spin-echo sequence.

Test 2. Constant Signal Intensity with Minimum RF Power
The obtained refocusing angles for the constant signal intensity calculations are shown in Figure 3. Note

Test 3. Exploiting Flexibility. Peak RF Constrained TSE and Recovery of the CPMG Condition
Setting a maximum of 60 o on the refocusing angle, we obtain the pulse trains as shown in Figure 4.
For the constant refocusing angle setup, the optimized pulse train has the maximum allowed amplitude, that is, 60 o refocusing angle and a 90 o 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.

Test 4: Eight-Channel pTx System at 7 Tesla. 3D ROI Optimization on a Phantom
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.
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.
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).

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 60 o , 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 voxeldependent, 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 spinecho 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 B þ 1 maps inaccuracies, as already anticipated in (6). As the scanner experiments were successful, we believe that the optimal control pulse trains are robust to B þ 1 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 FIG. 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. 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, s, could be subject to optimization too. The framework includes the possibility to treat also s as a variable and the derivatives with respect to s 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 B þ 1 , 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 EPGbased 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.

IMPLEMENTATION ASPECTS
The transition matrix, P, is given by the product

P ¼ ðI REÞS
where S denotes the shift operator matrix, I is the N Â N identity matrix and denotes the kronecker tensor product. R and E the rotation and decay matrices as given in Eqs. [3][4][5]. S is the shift operator matrix.