Accepted Manuscript Real-time 2D spatially selective MRI experiments: Comparative analysis of optimal control design methods Ivan I. Maximov, Mads S. Vinding, Desmond H.Y. Tse, Niels Chr. Nielsen, N. Jon Shah PII: DOI: Reference:

S1090-7807(15)00065-8 http://dx.doi.org/10.1016/j.jmr.2015.03.003 YJMRE 5619

To appear in:

Journal of Magnetic Resonance

Received Date: Revised Date:

9 September 2014 12 March 2015

Please cite this article as: I.I. Maximov, M.S. Vinding, D.H.Y. Tse, N. Chr. Nielsen, N. Jon Shah, Real-time 2D spatially selective MRI experiments: Comparative analysis of optimal control design methods, Journal of Magnetic Resonance (2015), doi: http://dx.doi.org/10.1016/j.jmr.2015.03.003

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Real-time 2D spatially selective MRI experiments: Comparative analysis of optimal control design methods Ivan I. Maximov1a,1,∗, Mads S. Vindingb,1,∗, Desmond H. Y. Tsea , Niels Chr. Nielsenb , N. Jon Shaha,c a Institute of Neuroscience and Medicine 4, Forschungszentrum J¨ ulich GmbH, 52425 J¨ ulich, Germany for Insoluble Protein Structures (inSPIN), Interdisciplinary Nanoscience Center (iNANO), and Department of Chemistry, Aarhus University, DK-8000 Aarhus, Denmark c Department of Neurology, Faculty of Medicine, RWTH Aachen University, JARA, 52074 Aachen, Germany

b Center

Abstract There is an increasing need for development of advanced radio-frequency (RF) pulse techniques in modern magnetic resonance imaging (MRI) systems driven by recent advancements in ultra-high magnetic field systems, new parallel transmit/receive coil designs, and accessible powerful computational facilities. 2D spatially selective RF pulses are an example of advanced pulses that have many applications of clinical relevance, e.g., reduced field of view imaging, and MR spectroscopy. The 2D spatially selective RF pulses are mostly generated and optimised with numerical methods that can handle vast controls and multiple constraints. With this study we aim at demonstrating that numerical, optimal control (OC) algorithms are efficient for the design of 2D spatially selective MRI experiments, when robustness towards e.g. field inhomogeneity is in focus. We have chosen three popular OC algorithms; two which are gradient-based, concurrent methods using first- and second-order derivatives, respectively; and a third that belongs to the sequential, monotonically convergent family. We used two experimental models: a water phantom, and an in vivo human head. Taking into consideration the challenging experimental setup, our analysis suggests the use of the sequential, monotonic approach and the second-order gradient-based approach as computational speed, experimental robustness, and image quality is key. All algorithms used in this work were implemented in the MATLAB environment and are freely available to the MRI community. Keywords: RF pulse design, optimal control, 2D spatially selective excitation, MRI, blOCh

1

2 3 4 5 6 7 8 9 10 11 12

1. Introduction Magnetic resonance imaging (MRI) is a very powerful, non-invasive imaging modality [1]. By exploiting different MRI contrasts, one can visualise not only the anatomical structure [2], but also the function [3], metabolism and tissue organisation of living cells and/or organs [4]. MRI has found numerous applications in medicine and the natural sciences as a brilliant diagnostic tool and real-time surgery assistant. However, some of the main problems of MRI still remain; low signal-to-noise ratio (SNR) and poor contrast due to multiple artefacts originate from either the subject (low natural polarisation, relaxation losses, physiological interference), or from imperfect hardware (magnetic field distortions, insensitive coil constructions, imprecise field gradients etc.). MRI safety is an inescapable issue to address as well; to avoid peripheral nerve stimulation and dangerous tissue heating, field gradient switching rates and specific absorption rates are limited [5–7]. Higher magnetic field strengths [8], hyperpolarisation [9], improvements in coil sensitivities and the strength of the magnetic field gradients [10] have been made to surmount the signal and contrast 1 Present

address: Experimental Physics III, TU Dortmund University, 44221 Dortmund, Germany authors Email addresses: [email protected] (Ivan I. Maximov), [email protected] (Mads S. Vinding) 1 These authors contributed equally to the work

∗ Corresponding

Preprint submitted to Journal of Magnetic Resonance

March 19, 2015

41

issues. Advanced radio-frequency (RF) pulse design has been used to overcome various problems such as field inhomogeneities [11, 12], susceptibility artefacts [13–15], hardware imperfections [16, 17], dissipative relaxation mechanisms [18–20], as well as to limit pulse power by taking hardware capabilities, and the local and global specific absorption rates [6, 7] into account. One group of interesting advanced RF pulses is 2D spatially selective RF pulses [21, 22]. They have applications of clinical relevance including for example: spinal cord imaging [23], MR spectroscopy [9], human brain diffusion [24] or functional [25] imaging. To obtain 2D spatially selective RF pulses encompassing such a large number of conditions often demand numerical methods which can be time consuming. The clinical sessions, however, are bound with a tight schedule. Hence, retrieving the 2D spatially selective RF pulses efficiently is imperative. In this study, our focus is fast optimisation of 2D spatially selective RF pulses for in vivo experiments. In this context, optimal control (OC) is a suitable candidate amongst many numerical approaches or theoretical approximations, since OC can handle multiple constraints and vast controls efficiently [21, 22]. Often OC is applied in engineering control tasks [26], but also in optical spectroscopy [27], solid- [28–31] and liquid-state [32–35] NMR, to design experiments for protein structure determination [36, 37], and in quantum information processing [38–40]. The first application of OC in RF pulse design for MRI was made more than 25 years ago [21, 22]. Decades later, OC found an increasing number of applications in MRI due to the development of parallel-transmit systems and ultra-high-field magnets [13, 41, 42]. We have previously given a detailed numerical and algorithmic account for design of 2D spatially selective RF pulses with OC efficiently [43]. Here we focus on the experimental application and concentrate on three popular OC candidates; the gradient-ascent pulse engineering (GRAPE) algorithm [34]; the Broyden-Fletcher-Goldfarb-Shanno (BFGS), quasi-Newton algorithm [44]; and the monotonic convergent (MC), sequential algorithm [47], which is very distinct from the other two in its numerical implementation. We considered their ability to adapt in experiments with two typical model systems: a water phantom, and an in vivo human head. We designed these RF pulses to excite square patterns and analysed their robustness towards B0 and B1 inhomogeneities. The paper is arranged as follows: the theory section gives a brief description of the function and subtleties of each OC algorithm seen from the MRI RF pulse design viewpoint. The experimental section explains the protocols we used and the details of the measurements. The results section describes the outcome of our measurements. We round off the paper with a discussion and conclusion of our work.

42

2. Theory

13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

T

43 44

45

The dynamics of the magnetisation vector M(r) = [Mx (r), My (r), Mz (r)] equation which without relaxation is written as     M (r, t) Mx (r, t) d  x My (r, t)  = Ω(r, t)  My (r, t)  , dt Mz (r, t) Mz (r, t) where the matrix Ω(r, t) incorporates all spatial-temporal terms:  0 γG(t) · r + ∆ω0 (r) 0 Ω(r, t) =  −γG(t) · r − ∆ω0 (r) β(r)ωy (t) −β(r)ωx (t) T

46 47 48 49 50 51 52 53 54

is described by the Bloch

 −β(r)ωy (t) β(r)ωx (t)  . 0

(1)

(2)

Here the vector r = [x, y, z] is a spin position vector. Spins at location r are irradiated by a resonant RF T field γB1 (r, t) = β(r) [ωx (t), ωy (t), 0] and are at the same time subject to a time-varying magnetic field T gradient G(t) = [Gx (t), Gy (t), Gz (t)] . The normalised factor β(r) corresponds to the spatial sensitivity of the B1 field. Ideally in the completely homogeneous case this factor is 1. The off-resonance condition is described by the ∆ω0 (r) term and should ideally be 0. In 2D spatial selective excitation, the goal is to transfer the magnetisation from its initial, equilibrium state to a state, with non-zero transverse components in a specific 2D pattern, while leaving the magnetisation outside the pattern unperturbed. Generally, we can describe this total target or desired state with a T spatial dependent magnetisation vector, Mdes (r) = [Mx,des (r), My,des (r), Mz,des (r)] . 2

55 56 57 58

The objective is to find the RF pulse waveform ωx (t) and ωy (t) which, in a finite time T , maximise the transfer efficiency from the initial state to the target state, in other words, minimise the difference between the final state at time T and the target state. In terms of OC, we maximise a performance functional J(r) in each spatial position ∫

T

J(r) = Φ [M(r, T )] +

Θ [M(r, t), ωx (t), ωy (t), t] dt,

(3)

0

73

where Φ [M(r, T )] corresponds to a final fidelity (expressed through the target and final magnetisations), and Θ[· · · ] is a regularisation term typically incorporating implicit constraints such as RF power. The OC problem Eq. (3) and its possible solutions are very well described in the literature [26, 46]. Often it is solved using Lagrange multipliers and the Euler-Lagrange equation [43, 46]. We have previously in Ref. [43] described how Eq. (3) is tackled with respect to achieving the goal, and described in detail how the numerical procedure for this is implemented. In short, formulae deduced from Eq. (3) are used to iteratively update the RF pulse waveforms from an initial guess to a newly improved waveform. The differences between the different algorithms lie in the performance of this step, the robustness during consecutive steps, and the number of steps required in total for an acceptable solution. Frequently, there is a trade-off between robustness and convergence rate [43]. One can also observe a difference between the expected theoretical performance and the numerically achieved performance. The discretisation of the continuous formulas in terms of computer interpretable functions is a crucial step that can compromise the outcome. We should therefore stress that observations within one field cannot be assumed general and are therefor not necessarily applicable in other fields. Our algorithms were developed for in vivo MRI applications and their computational speed is the main objective.

74

2.1. GRAPE

59 60 61 62 63 64 65 66 67 68 69 70 71 72

75 76

77 78 79 80 81 82

The gradient-ascent method GRAPE was initially realised for NMR by Khaneja and colleagues [34]. Trailing Ref. [43], it becomes apparent that an optimisation of Eq. (3) follows from these simple steps: 1. 2. 3. 4.

Make an intelligent initial guess for controls in the pulse waveform. Propagate the magnetisation forward in time from its known initial state using the given set of controls. Check if the final state is acceptable. If so quit, otherwise continue. Propagate the Lagrange multiplier backwards in time from the target state using the given set of controls. 5. Update the controls (see Eq. (21) in Ref. [43]) and go to step 2.

91

Qualitatively, the difference between the forward and backward propagated paths is used to evaluate the new set of controls in Step 5. Once the two paths coincide or are close to each other, the optimum is achieved. Usually, the GRAPE algorithm in our implementation finds a good solution after a few iterations. This, however, depends on how aggressively one attempts to pursue in the new directions, controlled by a step size parameter. Aggression in the beginning can bring results, but once the optimum is approached, GRAPE can overstep the goal and become unstable. Since the method is gradient-based it can only be guaranteed to find a local optimum. Multiple attempts may in some cases be necessary, but it is our experience that GRAPE almost always finds an acceptable solution for 2D spatially selective RF design. Providing GRAPE with an educated initial guess as proposed, a positive result is almost guaranteed.

92

2.2. BFGS

83 84 85 86 87 88 89 90

93 94 95 96 97

The Newton method exploits both first- and second-order derivatives, i.e. the gradients and Hessian matrix, respectively. The Hessian matrix is computationally expensive to evaluate and it is not always possible to deduce an analytic expression. Thus, the quasi-Newton schemes attempt to build the Hessian matrix approximately using iterative methods. One popular quasi-Newton method is BFGS [44]. As demonstrated in the literature [45], BFGS can keep up the pace towards an optimum in the later iterations, where GRAPE 3

98 99 100 101 102 103 104 105 106 107 108 109 110

111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126

slows down. For this promised convergence of BFGS, it is important that the gradients are exact. Otherwise approximation errors accumulate in the Hessian, eventually yielding a suboptimal convergence. It is, in fact, possible to analytically calculate the exact gradient of our system; however, it is evident that these expressions are computationally more demanding than first-order approximations. Moreover, in our situation it is more efficient to follow the proposal of Floether et al. [40], where the exact gradient computation is very efficient and calculated with machine precision, but still longer than that of the firstorder approximations. We have shown a comparison in the “Supplementary Materials”. Thus, for our BFGS implementation [43], a compromise existed between a fast, but suboptimal, and an improved, but slower variant. For the sake of speed in in vivo experiments, we chose the faster version, which in our experience yields comparable results within a tolerable time frame. Our implementation is similar to that of GRAPE, and the Hessian computation is done via the MATLAB (MathWorks, Natick, MA US) function fminunc. Hence, this function actually controls the flow, but the gradient calculation is done externally just as in our GRAPE implementation. 2.3. MC The MC algorithm is required to be monotonic convergent, i.e. improves the objective function with each iteration step until it reaches a desired precision. We used a general implementation of an MC algorithm [47] with two sets of controls, for the magnetisation and Lagrange multipliers, respectively. In the update procedure, the sets are entangled. Each set receives both a fraction of gradient-based information and another fraction of the other set. That fraction is individual to each set, and is determined by the operator. Another feature of the implemented MC algorithm is the sequential update. Instead of a full forward and full backward propagation as in the concurrent methods above, only one time step is done before the present control is updated. More thorough descriptions of this peculiar behaviour can be found in Maday and Turinici [47] and Vinding et al. [43]. The main goal of MC is a very fast convergence in early iterations; afterwards one observes the convergence rate goes down. In our experience the MC algorithm seems to be rather invariant to the initial starting point, making it more adaptable than GRAPE and BFGS. If for some reason the MC algorithm fails and primarily a loss of monotonicity is observed, e.g., moving from well-adjusted optimisation under ideal conditions to conditions with field inhomogeneities, this can often be repaired by adjusting the two fraction parameters slightly. The same may also be necessary with the gradient-based methods and their step size parameter.

138

2.4. Other methods We would like to mention that there are also other OC methods, for instance, genetic algorithms [53], unified computational approaches [54], exact derivatives in NMR computations [55], Krotov-based optimisation [56], hybrid-algorithms [43, 57] and many others (see, for example, Ref. [58] for an extensive comparison). Some of these novel methods demonstrate enormous efficiency for a given problem and may be used as an alternative to the algorithms described above. Unfortunately, a comparison of all algorithms is not possible and is outside of the scope of this manuscript. A proper adoption of one OC algorithm from one field to another requires extensive work. Unless needed, a quantum mechanical description, say of an NMR OC algorithm, may be too cumbersome to run in an MRI OC algorithm only requiring a Bloch simulation. Also, the lack of spatial information in NMR algorithms requires additional rearrangements for an algorithm to suit MRI applications. However, we would like to mention that we consider several of the aforementioned algorithms for our next projects due to the many interesting features they provide.

139

3. Experimental

127 128 129 130 131 132 133 134 135 136 137

140 141 142 143 144 145

All datasets were acquired with a 32-channel phase array head coil on a Siemens 3T Tim-Trio MAGNETOM scanner (Siemens Medical Systems, Erlangen, Germany) with a maximum gradient strength of 40 mT/m and a maximal slew rate of 200 T/m/s. A turbo spin echo (TSE) [1] was modified to include volume-selective pulses in place of the standard excitation pulse. The imaging parameters for each test case, such as repetition time (TR), echo time (TE), acceleration factor (AF), acquisition time (TA), resolution, field-of-view (FoV) etc. are specified below. 4

146

3.1. Workflow

Figure 1: Graphical representation of the real-time OC experiment workflow. 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168

169 170 171

A schematic workflow of the real-time 2D spatially selective experiments is presented in Fig. 1. After the volunteer or the sample was prepared in the scanner, the B0 field was shimmed, and then a localiser scan was acquired. Next, we ran an anatomical scan Protocol (either a conventional magnetisation prepared rapid gradient echo (MPRAGE) [1] or our implemented TSE pulse sequence). In the case where B0 and B1 maps are used, they were measured subsequently. The B0 map was in each case captured from a standard gradient echo sequence [1] with two echoes at echo times of 4.9 and 7.4 ms, respectively. For one slice image TA = 3.7 s with 5 mm3 isotropic resolution. For the B1 map, we used the actual flip-angle imaging (AFI) sequence [61, 62] with TE = 3.06 ms, TR1 = 16.7 ms, TR2 = 83.3 ms and a nominal flip angle of 60 degrees, TA = 25 s, with 5 mm3 isotropic resolution. In both cases, the maps were acquired with dimensions matching each other in order to avoid a resampling step in the simulations. The anatomical scan data were then sent to the optimisation computer with the appropriate image analysis software for localising the region of interest (ROI). We used ITK-SNAP [60] (www.itksnap.org) with easy ROI drawing and export facilities. Our export procedure was from the NifTI to the MATLAB format (using Jimmy Shens scripts of www.mathworks.de/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image). The optimisation package “blOCh” has a few options for importing a ROI, for example via an image file saved in png format. When the B0 and B1 maps were measured, they were also transferred to the optimisation computer. The results of the optimisation (RF pulse and gradient values) were saved in binary format and sent back to the scanner’s pulse sequence using a network connection. Depending on the optimisation parameters, the pulse design procedure which takes place before the actual 2D spatial-selection experiment can take around 5– 7 min. The bottleneck is the mapping procedure, hence, to optimise the remaining procedure, we recommend to use two operators: one at the scanner workstation and one at the optimisation computer. It allows one to perform the routine measurements and optimisations in parallel. 3.2. Optimisation The OC methods are implemented in a framework called “blOCh” which runs in MATLAB and is freely accessible on request under the General Public Licence. The 2D spatially selective RF pulses were designed 5

172 173 174 175 176 177 178 179 180 181 182 183 184 185

186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202

203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218

to make a flip angle of 90◦ degrees inside the ROI (see below); that is, to complete the transfer Mt=0 = [0, 0, M0 ]T to Mt=T = [M0 , 0, 0]T inside the ROI, and outside leave the magnetisation at equilibrium. We used a spiral gradient scheme during the excitation RF pulse period with the amplitude and slew rate limited design by Glover [59]. In order to stay safely within the hardware limits, we used 50% and 70% of the maximum allowed gradient and slew rate strengths, respectively. The numerical spatial grid was 128×128 spanning 280 mm by 280 mm in the sagittal plane, and the spiral had an in-plane resolution of 2 mm and a duration of 11.67 ms. The minimal discretisation time was set to 10 µs which corresponds to the scanner’s minimum gradient discretisation time; hence the total number of controls was 2 × 1167 (phases and amplitudes, or x- and y-amplitudes are individual controls). The optimisation parameters were adjusted with fast optimisation and field inhomogeneity compensation as the primary concern. One of the optimisation success criteria for each approach was similar numerical performance to the third decimal. Gaining the necessary experience to choose “good” parameters of simulation is essential. For the stationary samples, the time spent on searching the parameter space is not important. Thus through these experiments we could find parameters that would be feasible in the in vivo case. 3.3. Water phantom We used a bottled, 1900-mL water phantom from Siemens Medical Systems doped with 3.75 g of NiSO4 ×6H2 O and 5 g of NaCl per litre of distilled water. This water phantom was designed to exhibit the RF inhomogeneities similar to those measured in vivo in the human head. Protocol Ia: TR/TE = 200 ms / 17 ms, FoV = 280 mm, AF = 1, TA = 64 s, and the voxel size was 5 × 1 × 1 mm3 in one slice; Protocol Ib: TR/TE = 400 ms / 17 ms, FoV = 280 mm, AF = 1, TA = 128 s, and voxel size was 5 × 1 × 1 mm3 in 11 slices (slice spacing ∆x = 5 mm). The B0 and B1 maps are shown below in Figs. 2e,f, respectively. Protocol Ia was used to evaluate the effect of 1) ideal field maps as input to the optimisation; 2) the use of a measured B0 map in the optimisation; 3) the use of a measured B1 map in the optimisation; and 4) the use of both B0 and B1 maps in the optimisation. In all four cases, B0 -field shimming was performed on the imaging slice. Protocol Ib was used to estimate the robustness of the pulses. The RF pulses were optimised as in 1) to 4) above, for field maps acquired at the central slice location, but the data were gathered from all the slices. Hence, the slightly changing field inhomogeneities as we move away from the central slice, are used to statistically estimate the robustness of the pulses. Pulse sequences obtained from each algorithm and used in Protocols Ia,b are presented in the Supplementary Materials. 3.4. In vivo measurements A male, 33 years old, volunteer gave written informed consent prior to participation. The volunteer had no self-reported history of any psychiatric or neurological diseases and had no brain injuries. The anatomical localisation (Protocol IIa) with the TSE sequence had these parameters: TR/TE = 250 ms / 17 ms, FoV = 240 mm, TA = 77 s and the voxel size was 2.5 × 1 × 1 mm3 for one slice. Selective excitation was done using Protocol IIb: TR/TE = 250 ms / 17 ms, FoV = 240 mm, AF = 1, TA = 8 min 13 sec and the voxel size was 2.5 × 1 × 1 mm3 for one slice and 8 repetitions. The ROI, a 35 by 35 mm square was as illustrated, shown overlaid on the anatomical localisation in Fig. 5e. The B0 and B1 maps are shown in Fig. 5c,d. Protocol IIb was used to evaluate the case 1) when the slice was manually shimmed and only ideal field maps were used in the optimisation, and 2) the pulse was optimised with the B0 and B1 maps of the manually shimmed slice. In order to probe the robustness of the designed pulses accounting B0 and B1 field distortions we applied Protocol IIc: TR/TE = 350 ms / 17 ms, FoV = 240 mm, AF = 1, TA = 5 min 54 sec and the voxel size was 2.5 × 1 × 1 mm3 for 9 slices (slice spacing ∆x = 2.5 mm), and 4 repetitions. The pulse sequences computed by all algorithms and used in in vivo measurements are presented in the Supplementary Materials.

6

219

220 221 222 223 224 225 226

227 228 229 230 231 232 233

4. Results The robustness and excitation/suppression quality of designed RF pulses by the three different approaches were analysed. Our criterion for successful optimisation of the RF pulses initially was a visual inspection of the simulated excitation profile and reading of the numerical fidelity numbers. The latter can vary significantly depending on the algorithm and conditions, and we used our experience in this regard. However, we have to note the resulting numerical fidelities for all approaches differ only on the third decimal in all experimental cases. The numerical fidelity for each case is listed in the Supplementary Materials. Experimentally we have considered SNR and tissue-contrast as comparative criteria. 4.1. Water phantom The results from Protocol Ia are shown in Fig. 2. Histograms from the excitation and suppression areas for the images in Fig. 2a, b, c and d are shown in Fig. 2g, h, i and j, respectively. The histograms were calculated over the normalised images taking all cases from a to d together. Quantitative changes in the histograms of excitation quality are presented in Tab. 1. Using the smoothed and cubic spline interpolated histograms we estimated a maximum peak position (maximal signal) and a full width at half maximum (FWHM) for each case. Table 1: Quantitative histogram comparison of excitation quality in Fig. 2

g h i j

234 235 236 237 238 239 240

241 242 243 244 245 246 247 248 249 250 251 252 253 254

FWHM maximum FWHM maximum FWHM maximum FWHM maximum

GRAPE 0.388 0.35 0.201 0.60 0.190 0.65 0.202 0.64

BFGS 0.453 0.41 0.227 0.69 0.208 0.67 0.203 0.67

MC 0.378 0.38 0.212 0.62 0.202 0.65 0.188 0.64

The excitation and suppression distribution results from Protocol Ib in the form of scatter plots are shown in Figs. 3 and 4, respectively. For clarity, we have coloured those points that belong to the central imaging slice red, and all other points from neighbouring slices blue. The corresponding linear regression parameters of the points in the scatter plots, Fig. 3, are summarised in Table 2 for the central slice (x = 0) and the neighbouring slices (x = ±∆x, 2∆x, 3∆x, ...). The fitting model was ordinate = slope × abscissa + shif t and corr is the Pearson correlation coefficient assuming that all methods are linearly correlated. In all cases, the correlation significance p < 10−5 . 4.2. In vivo measurements Here our proposed real-time workflow (Fig. 1) was put to the test, in terms of thrifty labour. The volunteer had to remain still in the scanner until all three algorithms was tested for all field variations. It took 2 hours 30 minutes. Therefore, we limited the measurements to full FoV scans, Protocol IIb, shown in Fig. 5. Lacking the high resolution of “zoomed” FoV scans, we limited the histogram analysis to excitation and suppression signals only (see Fig. 5f and g, respectively). The histograms were normalised over all images similar to the water phantom case. In order to provide quantitative parametrisation of the histograms we summarised the FWHM and maximum signal values in Table 3. Again we tested each algorithm with the “uncompensated” pulses (with ideal B0 and B1 field maps in the optimisation), and “compensated” pulses (with measured B0 and B1 field maps in the optimisation). A robustness analysis was done for the “compensated” pulses in order to estimate the experimental stability of the applied approaches to the field variations (see Fig. 6 showing the results from Protocol IIc). The statistical assessments of the excitation quality are presented in Table 4. 7

Figure 2: The results with the water phantom and 2D spatially selective RF pulses obtained from the three algorithms GRAPE, BFGS, and MC. The Protocol Ia was used. Panel a) shows responses with a manually shimmed B0 field as the only means to improve field homogeneity, thus assuming ideal B0 and B1 maps. Panel b) with a B0 map included in the optimisation, panel c) with a B1 map included in the optimisation, and panel d) with both B0 and B1 maps included. The graphs in g)–j) correspond to histograms obtained from normalised images. The upper and lower histograms belong to the areas inside and outside the ROI, or excitation and suppression, respectively. The B0 , and B1 maps are shown in e), and f ), respectively.

8

Figure 3: Scatter plot results from the water phantom using Protocol Ib; statistical comparison between the three approaches in accordance with the field variations for excitation quality. Scatter plots in a) are with a manually B0 -field shimmed reference slice as the only means to improve field homogeneity, thus assuming ideal B0 and B1 maps everywhere. Scatter plots in b) are with the B0 map taken into account, scatter plots in c) are with the B1 map taken into account, and scatter plots d) is with both the B0 and B1 maps taken into account. The red points in the scatter plots correspond to data points from the reference slice, and the blue points are from the neighbouring slices.

Figure 4: Scatter plot results from the water phantom using Protocol Ib; statistical comparison between the three approaches in accordance with the field variations for suppression quality. The letter notation a)–d) is the same as in Fig.3.

9

Table 2: Linear fit and correlation parameters from the scatter plots in Fig. 3. The neighbouring slices have been fitted individually according to the distance to the central slice. The letter notation a)–d) is the same as in Fig.3

a

Slice x=0 x = ±∆x x = ±2∆x x = ±3∆x x = ±4∆x

b

x=0 x = ±∆x x = ±2∆x x = ±3∆x x = ±4∆x

c

x=0 x = ±∆x x = ±2∆x x = ±3∆x x = ±4∆x

d

x=0 x = ±∆x x = ±2∆x x = ±3∆x x = ±4∆x

slope/shift corr slope/shift corr slope/shift corr slope/shift corr slope/shift corr slope/shift corr slope/shift corr slope/shift corr slope/shift corr slope/shift corr slope/shift corr slope/shift corr slope/shift corr slope/shift corr slope/shift corr slope/shift corr slope/shift corr slope/shift corr slope/shift corr slope/shift corr

GRAPE vs MC 0.978/0.016 0.981 0.992/0.014 0.987 0.989/0.0129 0.988 0.977/0.016 0.988 0.995/0.014 0.989 0.913/0.021 0.959 0.964/0.001 0.975 0.941/0.012 0.977 0.926/0.0128 0.977 0.937/0.015 0.979 0.939/0.025 0.966 1.004/0.003 0.975 1.003/0.004 0.977 0.958/0.023 0.978303 0.989/0.009 0.982 0.961/0.019 0.963 0.987/0.007 0.979 0.975/0.008 0.981 0.977/0.008 0.979 0.977/0.013 0.984

BFGS vs MC 1.052/0.027 0.981 0.992/0.045 0.981 0.999/0.044 0.984 1.024/0.033 0.984 1.043/0.036 0.986 1.091/0.006 0.965 0.987/0.042 0.962 1.014/0.032 0.970 1.073/-0.001 0.975 1.067/0.009 0.977 0.997/0.013 0.966 0.999/0.0126 0.975 1.0118/0.007 0.976 1.005/0.011 0.978081 1.031/-0.001 0.978 1.043/0.001 0.966 0.925/0.043 0.964 0.942/0.037 0.967 1.019/-0.002 0.974 1.011/0.013 0.977

BFGS vs GRAPE 0.913/-0.003 0.982 0.974/ -0.019 0.981 0.971/-0.022 0.984 0.938/-0.008 0.986 0.940/-0.013 0.989 0.809/0.028 0.960 0.917/-0.006 0.952 0.891/0.002 0.967 0.837/0.0267 0.972 0.857/0.019 0.979 0.911/0.027 0.966 0.978/0.004 0.974 0.966/0.0103 0.975 0.931/0.023 0.977282 0.938/0.023 0.981 0.893/0.031 0.966 1.009/-0.010 0.960 0.988/-0.005 0.968 0.926/0.026 0.971 0.942/0.015 0.981

Table 3: Quantitative histogram comparison of excitation quality in Fig. 5

f g

FWHM maximum FWHM maximum

GRAPE 0.475 0.48 0.450 0.68 10

BFGS 0.485 0.52 0.438 0.66

MC 0.470 0.48 0.431 0.64

Figure 5: The in vivo results of excitation of a square ROI with the 2D spatially selective RF pulses obtained from the three algorithms GRAPE, BFGS, and MC. a) Uncompensated images correspond to a manually B0 -field shimmed slice and assuming ideal B0 and B1 maps. b) Compensated images correspond to the RF pulse design taking the B0 and B1 maps into account. The histograms in f ) and g) correspond to areas inside and outside the ROI, or excitation and suppression, for the uncompensated and compensated images, respectively. The anatomical image, highlighting the ROI e), the B0 c), and B1 maps d) are shown as well.

11

Figure 6: Scatter plot results from the in vivo experiments using Protocol IIc; statistical comparison between the three approaches in accordance with the field variations for excitation/suppression quality. Scatter plots in a) is an excitation statistics with both the B0 and B1 maps taken into account. Scatter plots in b) is a suppression statistics with both the B0 and B1 maps taken into account. The red points in the scatter plots correspond to data points from the reference slice, and the blue points are from the neighbouring slices.

12

Table 4: Linear fit and correlation parameters from the scatter plots in Fig. 6. The neighbouring slices have been fitted individually according to the distance to the central slice.

a

Slice x=0 x = ±∆x x = ±2∆x x = ±3∆x x = ±4∆x

255

256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290

slope/shift corr slope/shift corr slope/shift corr slope/shift corr slope/shift corr

MC vs GRAPE 0.882/0.017 0.804 0.813/0.032 0.804 0.816/0.0342 0.778 0.971/0.008 0.896 0.883/0.021 0.926

MC vs BFGS 0.811/0.038 0.813 0.782/0.043 0.816 0.778/0.046 0.792 1.009/0.008 0.885 0.846/0.034 0.930

GRAPE vs BFGS 0.733/0.066 0.806 0.763/0.060 0.806 0.741/0.061 0.789 0.933/0.027 0.887 0.898/0.032 0.942

5. Discussion In this paper, we set forth to present a comparative analysis of the experimental models subjected to 2D spatially selective MRI measurements using RF pulses from three popular OC algorithms, GRAPE, BFGS and MC. The first model system, where we compared an algorithmic tolerance to field inhomogeneities was the water phantom. Using Protocol Ia we excited a slice in the central part of the phantom and analysed the quality of the measured images with visual inspection, which may give qualitative guidance. Quantitative quality of the images was evaluated via the histograms and their parameters such as the FWHM and maximal signal belonging to the excitation and suppression parts of the images. Ideally magnetisation excited with 90◦ flip-angle should manifest itself as signal in the right part of the histogram as close as possible to 1. In turn, unavoidable signal from the suppression region should be located left near 0. In the case of ideal field maps, which was assumed in the optimisation, but obviously not true in the experiments, we see in Figs. 2a,g that all approaches had more or less the same excitation fidelity. Interestingly for all cases, the suppression patterns and histograms are basically coincided and have nearly identical performances. It is clear that absence of B0 or B1 maps compensation leads to blurred excitation patterns in the images. Including the B0 map, keeping the B1 map ideal, (see Figs. 2b,h), all three algorithms substantially improved the image quality, particularly, sharpness of excitation profile (see, for example, Tab. 1 for values). Usage of the B1 map only, (Figs. 2c,i), allows one to increase the performance of each algorithms comparing to the previous case. Adding both B0 and B1 maps, (Figs. 2d,j), led MC to moderate improve the excitation profile. GRAPE and BFGS kept their results as in Fig. 2i, (check the values in Tab. 1). The water phantom data show a tendency for BFGS and MC to provide somewhat better results than GRAPE based on the more homogeneous profiles and enhanced signal amplitudes. The scatter plot data from Protocol Ib shown in Fig. 3 along with the linear regression values presented in Table 2, gives us further hints on robustness of the designed pulses in respect to B0 /B1 maps. The plots show how all image points from one scan are distributed with the same points from another scan. Thus, for one algorithm versus another, a linear trend above the diagonal, and a slope larger than 1, favours the first algorithm in the excitation case and the second algorithm in the suppression case, vice versa. Excitation-wise, assuming ideal field maps, the scatter plots in Fig. 3a do not show a huge difference between the algorithms in the central slice as well as in neighbouring slices. The experimental variations are however evident in the fitting parameters, see Tab. 2. Now including the B0 map, Fig. 3b, BFGS outperformed the other two in the central excitation slice as well as in neighbouring slices. A similar situation arose when including the B1 map only, Fig. 3c and Tab. 2. Addition of both B0 and B1 maps, Fig. 3d, showed nearly the same results, where BFGS has a slightly higher statistical efficiency, Tab. 2. Suppression-wise, we conclude that a linear regression model is not very suitable in this form of presenting the data. Thus, the Pearson correlation coefficient is not a suitable parameter of the drawn scatter plots. However, the visual inspection of scatter plots allows one to come to 13

291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342

conclusion about the character of suppression obtained by three approaches. Indeed, one can see that the BFGS excitation quality diminishes its suppression quality, statistically producing higher values outside the ROI. We can still, be guided qualitatively by the scatter plots, though. In many cases, the points distributed quite evenly around the diagonal. That coincides with the histogram observations above and its numerical characteristics. Table 2 presents additional information about the robustness and stability of the executed pulses based on the scatter plots in Fig. 3, however with a distinction between the individual neighbouring slices and the central slice. E.g. the notation x = ±2∆x means that the tabularised data are taken from the slices second closest to the central slice. We see that the statistical variation depending on the slice shift ∆x is a minimal almost in all four cases. The next test with the in vivo volunteer and head scans yielded full FoV data for a histogram analysis of excitation and suppression signals, see Fig. 5; both with “uncompensated and “compensated” pulses. Histograms based on the excitation signal showed broadened peaks corresponding to contamination by grey matter of a prevailing white matter pattern. For all three algorithms we saw that the “compensated” pulses increased the signal of the histogram peaks and decreased the FWHM values of histograms, likely, exhibiting a better performance in the narrower peaks of BFGS and MC. Histograms of the suppression signals showed an invisibility between BFGS and MC while GRAPE stood out from the other two. Robustness of “compensated” pulses in vivo is very important for clinical applications. It allows one to design a pulse for a given slice and use it for a slice in its vicinity, or it can increase the tolerance to, for example, patient movement altering the field distributions. We see that in the in vivo case the Pearson correlation coefficients in Tab. 4 are much lower comparing to the water phantom experiment due to the stronger field variations. In course of the statistical assessments for the water phantom in Tab. 2, the in vivo model also allows us to probe the partial volume effect due to the different image resolution of the field maps (slice thickness 5 mm) and the acquired images (slice thickness 2.5 mm), see Tab. 4. This effect is also insignificant, enabling faster field mapping sequences with enhanced SNR due to the larger slice thickness. The experiments with the 2D spatially selective RF pulses provided means for field inhomogeneity robustness tests. The time frame needed to develop these pulses including the field maps depends obviously on the field mapping techniques. We used a B0 mapping sequence which for one slice took 3.7 s. The B1 AFI mapping was 25 s. Handling of these maps and choosing the ROI, as we did with ITK-SNAP, took around 3 minutes. Note, that specially designed tools can improve this step significantly. For the chosen 128 × 128 spatial grid, the spiral resolution, and the other optimisation parameters, e.g., the gradient slewrate limitation, a typical optimisation time, say, with the BFGS with both field maps in the in vivo case, was 86 seconds. Similar durations exist for the GRAPE and MC algorithms. Clearly, the spatial grid and spiral resolution can be reduced significantly, if square ROIs are indeed the intention. For more advanced, say, anatomical shapes, the grid size should be chosen to reflect the details of the shape. We chose 128 × 128 for demonstration purposes. How and how far the optimisation can progress is another factor to determine and trade-off. Our software enables monitoring of the optimisation progress i.e. by the optimisation metric or numerical fidelity, but also by the visual display of the magnetisation profile obtained freely via the Bloch framework we use. Both, but especially the visual display, matter for the computation time and they can be turned off for faster optimisation. Each algorithm yielded nearly the same performance in every case, as stated, to the third decimal of the optimisation metric. The visual guidance was the second assurance that the optimisations were “well behaved”. In our opinion the time spent attempting to enhance the numerical fidelity much further can perhaps be spent on more important issues in in vivo MRI experiments. To the best of our knowledge it is still an open question how to choose optimisation parameters for each algorithm in a fast and efficient way. Typically, and as in our case, the MRI practitioner tries to apply a few parametrisations in order to find an acceptable fidelity in both excitation and suppression, depending on the application. However, we expect significant improvements by application of a mathematical framework in parameter estimations such as the Cramer-Rao approach or similar [64]. This approach can also be applied for selection of more effective imaging protocols. We hope to present related results in the future based on the Cramer-Rao approach. 14

345

Summing up, with our strategy we could in less than five minutes produce a B0 - and B1 -inhomogeneity compensating 2D spatially selective RF pulse. The total experiment time will obviously depend on the application the pulse is implemented into.

346

6. Conclusion

343 344

358

Aimed for 2D spatially selective MRI experiments, we have presented three OC algorithms (GRAPE, BFGS and MC) we see fit for in vivo applications. We tested the designed RF pulses in two subjects (a water phantom, and an in vivo head), and our statistical survey included influences of B0 and B1 inhomogeneity and compensation thereof by utilising B0 and B1 maps. We found that field maps are essential inputs in the 2D spatially selective RF pulse design in order to retrieve good excitation fidelities and obtaining images with higher tissue contrasts. Although the field is close, the candidates BFGS and MC are the recommended methods for 2D spatially selective RF pulse design with in vivo applications. They provided better image quality with respect to SNR and tissue contrast, and proved to be more robust to variations of RF and field inhomogeneities compared with the GRAPE algorithm. Our MATLAB scripts are freely available to the MR community enabling correction, modification and addition of any other features in order to promote the interest of OC in pulse sequence design.

359

Acknowledgement

347 348 349 350 351 352 353 354 355 356 357

366

We acknowledge the Deutsche Forschungsgemeinschaft (SU 192/32-1), Danish National Research Foundation (DNRF59) and the Danish Ministry of Higher Education and Science, UNIK (UNiversitetsforskningens InvesteringsKapital), MINDLab 0601-01354B. D.H.Y.T. is supported by the Marie Curie Initial Training Network Methods in Neuroimaging (MC-ITN-238593). The authors gratefully thank Prof. Torben E. Lund and Dr. Ryan Sangill for assistance with measurements, Dr. Daniel Brenner and Dr. Michael Poole for their valuable help with B0 /B1 field mapping procedures, and Prof. Ilya Kuprov for valuable discussions on exact gradients.

367

References

360 361 362 363 364 365

368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392

[1] Berstein MA, King KF, Zhou XJ. Handbook of MRI pulse sequences. Elsevier Academic Press, San Diego, 2004. [2] Rosenberg J, Maximov II, Reske M, Grinberg F, Shah NJ. Early to bed, early to rise: Diffusion tensor imaging identifies chronotype-specificity. NeuroImage 2014; 84: 428-434. [3] Wallentin M, Nielsen AH, Vuust P, Dohn A, Roepstoff A, Lund TE. Amygdala and heart rate variability responses from listening to emotionally intense parts of a story. NeuroImage 2011; 58: 963-973. [4] Kurhanewicz J, Vigneron DB, Brindle K, Chekmenev EY, Comment A, Cunningham CH, DeBerardinis RJ, Green GG, Leach MR, Rahan SS, Rizi RR, Ross BD, Warren WS, Malloy CR. Analysis of cancer metabolism by imaging hyperpolarized nuclei: prospects for translation to clinical research. Neoplasia 2011; 13: 81-87. [5] Schulte RF, Hoeske R, Peripheral nerve stimulation – optimal gradient waveform design. Magn. Reson. Med. 2014; doi: 10.1002/mrm.25440 [6] Brunner DO, Pruessmann KP. Optimal design of multiple-channel RF pulses under strict power and SAR constraints. Magn. Reson. Med. 2010; 63: 1280-1291. [7] Guerin B, Gebhardt M, Cauley S, Adalsteinsson E, Wald LL. Local specific absorption rate (SAR), global SAR, transmitter power, and excitation accuracy trade-offs in low flip-angle parallel transmit pulse design. Magn. Reson. Med. 2013; doi: 10.1002/mrm.24800 [8] Budde J, Shajan G, Scheffler K, Pohmann R. Ultra-high resolution imaging of the human brain using acquisition-weighted imaging at 9.4T. NeuroImage 2014; 86: 592-598. [9] Vinding MS, Laustsen C, Maximov II, Søgaard LV, Ardenkjær-Larsen, Nielsen NC. Dynamic nuclear polarization and optimal control spatial-selective 13 C MRI and MRS. J. Magn. Reson. 2013; 227: 57-61. [10] Van Essen DC, U˘ gurbil K, Auerbach E, Barch D, Behrens TEJ, Bucholz R, Chang A, Chen L, Corbetta M, Curtiss SW, Della Penna S, Feinberg D, Glasser MF, Harel N, Heath AC, Larson-Prior L, Marcus D, Michalareas G, Moeller S, Oostenveld R, Petersen SE, Prior F, Schalggar BL, Smith SM, Snyder AZ, Xu J, Yacoub E, WU-Minn HCP Consortium. The human connectome project: a data acquisition perspective. NeuroImage 2012; 62: 2222-2231. [11] Cloos MA, Boulant N, Luong M, Ferrand G, Giacomini E, Le Bihan D, Amadon A. kT -points: short three-dimensional tailored RF pulses for flip-angle homogenization over an extended volume. Magn. Reson. Med. 2012; 67: 72-80.

15

393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457

[12] Hsu YC, Chern IL, Zhao W, Gagoski B, Witzel T, Lin FH. Mitigate B+ 1 inhomogeneity using spatially selective radiofrequency excitation with generalized spatial encoding magnetic fields. Magn. Reson. Med. 2013; doi: 10.1002/mrm.24801 [13] Yip CY, Fessler JA, Noll DC. Advanced three-dimensional tailored RF pulse for signal recovery in T∗2 -weighted functional magnetic resonance imaging. Magn. Reson. Med. 2006; 56: 1050-1059. [14] Zhao F, Nielsen JF, Noll DC. Four dimensional spectral-spatial fat saturation pulse design. Magn. Reson. Med. 2013; doi: 10.1002/mrm.25076 [15] Zheng H, Zhao T, Qian Y, Schirda C, Ibrahim TS, Boada FE. Multi-slice parallel transmission three-dimensional tailored RF (PTX 3DTRF) pulse design for signal recovery in ultra high field functional MRI. J. Magn. Reson. 2013; 228: 37-44. [16] Zheng H, Zhao T, Qian Y, Ibrahim T, Boada F. Parallel transmission RF pulse design for eddy current correction at ultra high field. J. Magn. Reson. 2012; 221: 139-146. [17] Wu X, Vaughan JT, U˘ gurbil K, Van de Moortele PF. Parallel excitation in the human brain at 9.4 T counteracting k-space errors with RF pulse design. Magn. Reson. Med. 2010; 63: 524-529. [18] Gershenzon NI, Kobzar K, Luy B, Glaser SJ, Skinner TE. Optimal control design of excitation pulses that accommodate relaxation. J. Magn. Reson. 2007; 188: 330-336. [19] Shen J. Preserving the excitation profile of small flip angle RF pulses in the presence of rapid transverse relaxation. J. Magn. Reson. 2012; 224: 8-12. [20] Skinner TE, Gershenzon NI, Nimbalkar M, Glaser SJ. Optimal control design of band-selective excitation pulses that accommodate relaxation and RF inhomogeneity. J. Magn. Reson. 2012; 217: 53-60. [21] Nishimura DG. A multiple-pulse sequence for improved selective excitation in magnetic resonance imaging. Med. Phys. 1985; 12: 413-418. [22] Conolly S, Nishimura D, Macovski A. Optimal control solutions to the magnetic resonance selective excitation problem. IEEE Trans Med. Imaging 1986; 5: 106-115. [23] Finsterbusch J. High-resolution diffusion tensor imaging with inner field-of-view EPI. J. Magn. Reson. Imaging. 2009; 29: 987-993. [24] Jeong EK, Kim SE, Guo J, Kholmovski EG, Parker DL. High-resolution DTI with 2D interleaved multislice reduced FOV single-shot diffusion-weighted EPI (2D ss-rFOV-DWEPI). Magn. Reson. Med. 2005; 54: 1575-1579. [25] Finsterbusch J. Simultaneous functional MRI acquisition of distributed brain regions with high temporal resolution using a 2D-selective radiofrequency excitation. Magn. Reson. Med. 2014; doi: 10.1002/mrm.25143. [26] Krotov VE, Global methods in optimal control theory. Marcel Dekker, New York 1996. [27] Waren WS, Rabitz H, Dahleh M. Coherent control of quantum dynamics: the dream is alive. Science 1993; 259: 1581-1589. [28] Kehlet CT, Sivertsen AC, Bjerring M, Reiss TO, Khaneja N, Glaser SJ, Nielsen NC. Improving solid-state NMR dipolar recoupling by optimal control. J. Am. Chem. Soc. 2004; 126: 10202-10203. [29] To˘sner Z, Glaser SJ, Khaneja N, Nielsen NC. Effective Hamiltonians by optimal control: solid-state NMR double-quantum planar and isotropic dipolar recoupling. J. Chem. Phys. 2006; 125: 184502. [30] Kehlet C, Bjerring M, Sivertsen AC, Kristensen T, Enghild JJ, Glaser SJ, Khaneja N, Nielsen NC. Optimal control based NCO and NCA experiments for spectral assignment in biological solid-state NMR spectroscopy. J. Magn. Reson. 2007; 188: 216-230. [31] Nielsen AB, Bjerring M, Nielsen JT, Nielsen NC. Symmetry-based dipolar recoupling by optimal control: band-selective experiments for assignment of solid-state NMR spectra of proteins. J. Chem. Phys. 2009; 131: 025101. [32] Reiss TO, Khaneja N, Glaser SJ. Time-optimal coherence-order-selective transfer of in-phase coherence in heteronuclear IS spin systems. J. Magn. Reson. 2002; 154: 192-195. [33] Skinner TE, Reiss TO, Luy B, Khaneja N, Glaser SJ. Reducing the duration of broadband excitation pulses using optimal control with limited RF amplitude. J. Magn. Reson. 2004; 167: 68-74. [34] Khaneja N, Reiss T, Kehlet C, Schulte-Herbr¨ uggen T, Glaser SJ. Optimal control of coupled spin dynamics: design of NMR pulse sequences by gradient ascent algorithms. J. Magn. Reson. 2005; 172: 296-305. [35] Neves JL, Heitmann B, Khaneja N, Glaser SJ. Heteronuclear decoupling by optimal tracking. J. Magn. Reson. 2009; 201: 7-17. [36] Frueh DP, Ito T, Li JS, Wagner G, Glaser SJ, Khaneja N. Sensitivity enhancement in NMR of macromolecules by application of optimal control theory. J. Biomol. NMR. 2005; 32: 23-30. [37] Bjerring M, Jain S, Paaske B, Vinther JM, Nielsen NC. Designing dipolar recoupling and decoupling experiments for biological solid-state NMR using interleaved continuous wave and RF pulse irradiation. Acc. Chem. Res. 2013; 46: 20982107. [38] Serra RM, Oliveira IS. Nuclear magnetic resonance quantum information processing. Philos. Trans. A Math. Phys. Eng. Sci. 2012; 370: 4615-4619. [39] Criger B, Passante G, Park D, Laflamme R. Recent advances in nuclear magnetic resonance quantum information processing. Philos. Trans. A Math. Phys. Eng. Sci. 2012; 370: 4620-4635. [40] Floether FF, de Fouquieres P, Schirmer SG. Robust quantum gates for open systems via optimal control: Markovian versus non-Markovian dynamics. New Journal of Physics 2012, 14: 073023. [41] Grissom W, Yip CY, Zhang Z, Stenger VA, Fessler JA, Noll DC. Spatial domain method for the design of RF pulses in multicoil parallel excitation. Magn. Reson. Med. 2006; 56: 620-629. [42] Grissom WA, Yip CY, Wright SM, Fessler JA, Noll DC. Additive angle method for fast large-tip-angle RF pulse design in parallel excitation. Magn. Reson. Med. 2008; 59: 779-787. [43] Vinding MS, Maximov II, To˘sner Z, Nielsen NC. Fast numerical design of spatial-selective RF pulses in MRI using Krotov and quasi-Newton based optimal control methods. J. Chem. Phys. 2012; 137: 054203. [44] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical Recipes. The art of scientific computing. Cambridge

16

492

University Press. Cambridge 2007. [45] de Fouquieres P, Schirmer SG, Glaser SJ, Kuprov I. Second order gradient ascent pulse engineering. J. Magn. Reson. 2011; 212: 412-417. [46] Bryson AE, Ho Y. Applied optimal control:optimization, estimation, and control. Taylor & Francis Group, New York 1975. [47] Maday Y, Turinici G. New formulation of monotonically convergent quantum control algorithms. J. Chem. Phys. 2003; 118: 8191-8196. [48] Maximov II, To˘sner Z, Nielsen NC. Optimal control design of NMR and dynamic nuclear polarization experiments using monotonically convergent algorithms. J. Chem. Phys. 2008; 128: 184505. [49] Maximov II, Salomon J, Turinici G, Nielsen NC. A smoothing monotonic convergent optimal control algorithm for nuclear magnetic resonance pulse sequence design. J. Chem. Phys. 2010; 133: 084107. [50] Hogben HJ, Krzystyniak M, Charnock G, Hore PJ, Kuprov I. Spinach – a software library for simulation of spin dynamics in large spin systems. J. Magn. Reson. 2011; 208: 179-194. [51] Spinach library. http://spindynamics.org/Spinach.php [52] Rabitz HA, Hsieh MM, Rosenthal CM. Quantum optimally controlled transition landscapes. Science 2004; 303: 1998-2001. [53] Bechmann M, Clark J, Sebald A. Genetic algorithms and solid state NMR pulse sequences. J. Magn. Reson. 2013; 228: 66-75. [54] Li JS, Ruths J, Yu TY, Arthanari H, Wagner G. Optimal pulse design in quantum control: A unified computational method. PNAS 2011; 108: 1879-1884. [55] Kuprov I, Rodgers C. Derivatives of spin dynamics simulations. J. Chem. Phys. 2009; 131: 243108. [56] Reich DM, Ndong M, Koch CP. Monotonically convergent optimization in quantum control using Krotov’s method. J. Chem. Phys. 2012; 136: 104103. [57] Eitan R, Mundt M, Tannor DJ. Optimal control with accelerated convergence: combining the Krotov and quasi-Newton methods. Phys. Rev. A 2011; 83: 053426. [58] Machnes S, Sander U, Glaser SJ, de Fourquieres P, Gruslys A, Schirmer S, Schulte-Herbr¨ uggen. Comparing, optimizing, and benchmarking quantum-control algorithms in a unifying programming framework. Phys. Rev. A 2011; 84: 022305. [59] Glover GH. Simple analytic spiral k-space algorithm. Magn. Reson. Med. 1999; 42: 412-415. [60] Yushkevich PA, Piven J, Hazlett HC, Smith RG, Ho S, Gee JG, Gerig G. User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability. NeuroImage 2006; 31: 1116-1128. [61] Yarnykh VL. Actual flip-angle imaging in the pulsed steady state: A method for rapid three-dimensional mapping of the transmitted radiofrequency field. Magn. Reson. Med. 2007; 57(1): 192-200. [62] Nehrke K. On the steady-state properties of actual flip angle imaging (AFI). Magn. Reson. Med. 2009; 61(1): 84-92. [63] Blasche M, Riffel P, Lichy M. TimTX TrueShape and syngo ZOOMit: technical and practical aspects. Magnetom Flash 2012; 74-84. [64] Kay SM. Fundamentals of statistical signal processing. Estimation theory. Pretence Hall PTR. 1993.

493

———————————

458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491

17

blOCh

Graphical Abstract

GRAPE MC BFGS B0 B1

Highlights

Highlights Experimental tests of different OC algorithms Efficient workframe for in vivo MRI pulse design Computer software “blOCh” for fast OC pulse design