Genetic algorithm- based optimization of pulse sequences

Purpose: The performance of pulse sequences in vivo can be limited by fast relaxation rates, magnetic field inhomogeneity, and nonuniform spin excitation. We describe here a method for pulse sequence optimization that uses a stochastic numerical solver that in principle is capable of finding a global optimum. The method provides a simple framework for incorporating any constraint and imple-menting arbitrarily complex cost functions. Efficient methods for simulating spin dynamics and incorporating frequency selectivity are also described. Methods: Optimized pulse sequences for polarization transfer between protons and X- nuclei and excitation pulses that eliminate J- coupling modulation were evaluated experimentally using a surface coil on phantoms, and also the detection of hyperpolarized [2- 13 C]lactate in vivo in the case of J- coupling modulation- free excitation. Results: The optimized polarization transfer pulses improved the SNR by ~50% with a more than twofold reduction in the B 1 field, and J- coupling modulation-free excitation was achieved with a more than threefold reduction in pulse length. Conclusion: This process could be used to optimize any pulse when there is a need to improve the uniformity and frequency selectivity of excitation as well as to design new pulses to steer the spin system to any desired achievable state.

get trapped in local minima, 1 and the gradient of the cost function can be difficult to calculate. Another constraint is that the time-steps be constant, 1 imposing a further restriction on the process. As a potential solution for nonconvexconstrained optimization problems, a genetic algorithm was proposed, 13 which is a stochastic solver working with a population of solutions. Convergence to a global optimum is possible with no initial guess and convexity requirements on the cost function. A further advantage is that a cost-function gradient does not need to be calculated, yielding a simple framework that can incorporate any constraint and can be used to compose arbitrarily complex cost functions. Veglia et al 2 introduced a genetic algorithm-based approach for optimized spin dynamics, but this optimizes only the phases of a fixed number of pulses and delays between the pulses in order to avoid the need for spin dynamics simulations.
Here we present a more general approach, which includes an efficient method for spin dynamics simulations and incorporates frequency selectivity. The optimization process was used here to design improved pulse sequences for transferring polarization from hyperpolarized X nuclei to spin coupled protons with the aim of increasing the sensitivity of detection of hyperpolarized 13 C and 15 N-labeled metabolites. [14][15][16] Given the reciprocity of this polarization transfer process, we used thermally polarized instead of hyperpolarized metabolites and compared the performance of the optimized pulses to the conventional insensitive nuclei enhanced by polarization transfer (INEPT) pulse sequence 17 for transferring polarization from protons to X nuclei in [1][2][3][4][5][6][7][8][9][10][11][12][13] C]lactate, [2-13 C]lactate, and [ 15 N 2 ]urea. The former shows weak J-coupling, for which relaxation losses are mainly responsible for decreased sensitivity. The latter 2 are examples of strong J-coupling where the validity of the product operator formalism, which assumes instantaneous rotations, is compromised; and for [2-13 C]lactate, the several ms long adiabatic pulses used in experiments with surface coils in vivo to improve B 1 field uniformity cannot fit into the ideal time window of the polarization transfer block. With [ 15 N 2 ] urea, we also demonstrate the feasibility of partial transfer of polarization from the hyperpolarized 15 N nucleus to the spin-coupled protons. 15 Finally, we design excitation pulses that result in the selection of a single component of the 13 C and 15 N multiplets in [2-13 C]lactate and [ 15 N 2 ]urea respectively, which eliminate J-coupling-associated image artefacts, and compare the results obtained with the solution proposed in Ref. [18].

| THEORY
The optimization treats the pulse sequence as a shaped pulse with arbitrary amplitude, phase, and duration at each pulse point on each frequency channel. The cost function contains the spin dynamics simulation and a subsequent Bloch simulation to yield the desired frequency selectivity.

| Encoding the solutions
The shaped pulse at each frequency (e.g., 13 C and 1 H) is discretized to N number of time-steps. During each time-step, every pulse has an amplitude (B 1 ), duration ( ), and phase (Φ) (Figure 1). For delays, the RF amplitude of the time-step is set to 0. The only restriction is the number of time-steps. The minimum length of the time-steps can be constrained by the time resolution of the waveform-generator and the amplitudes of the RF pulses by the performance of the RF amplifier and transmit coil. Duty cycle constraints can also be incorporated as a penalty term in the cost function. Including the length of the pulse points in the optimization variables increased the flexibility of the algorithm; however, the final pulses had to be adjusted to match the RF raster time, which was 4 µs. This entailed rounding of the timesteps to integer multiples of the raster time and a slight correction in the corresponding amplitude. This had no noticeable effect on pulse performance. Alternatively, the time-steps can be constrained during the optimization to only take values that are integer multiples of the raster time.

| Smoothness and amplitude constraints on shaped pulses
Because RF coils and amplifiers cannot produce arbitrarily high B 1 fields or produce an instantaneous change in pulse amplitude, the amplitudes of the shaped pulses were truncated to a user-defined maximum value and then projected onto a set of waveforms with prescribed upper bounds on the first and second derivatives for each iteration of the optimization. The result is a smoothed equivalent, which is the closest to the original input waveform in the function norm defined in Equation 1. This is a convex problem that can be solved efficiently with accelerated proximal gradient descent. 19 This projection does not add a noticeable increase to the run time of the cost function evaluation ( Figure 2) and gave smooth changes in spin dynamics with respect to B 1 deviations and frequency offsets. This allowed the use of very coarse B 1 and frequency offset grids and accelerated considerably solver iterations. The values of α and β were selected empirically and had no noticeable effect on the runtime.

| Cost function and spin-dynamics simulations
The spin system Hamiltonian contains interaction terms for the nuclear spins with the external magnetic field and with each other. By adding RF pulses, a general effective Hamiltonian can be simulated. In Liouville-space, the equation of motion for the density operator ( ) in the presence of relaxation can be written as 20 : (2) Parameters used in pulse sequence optimization (A). The parameters are the amplitude (B 1 ), phase (Φ i ) and the length (τ i ) of the individual pulse points. A shaped pulse is played out on the RF channel(s), which steer the initial density operator (̂ 0 ) to the target density operator (̂ tgt ). In a heteronuclear polarization transfer experiment, there are two RF channels: one for proton and one for the Xnucleus (B). Flowchart for calculation of the cost function (C). Φ i , phase; τ i , length; B 1 , amplitude F I G U R E 2 (A) Runtime of the cost function evaluation for 5 coupled spins as a function of the number of RF pulse points. The additional runtime for implementation of the smoothness constraint (B) and the Bloch simulation for enforcing frequency selectivity (C) was negligible in comparison (6-10 ms and 200-400 ms, respectively). The runtime analysis for the Bloch simulation also illustrates the speed of the proposed method for pulse sequence design when there is no spin-spin coupling where Ĥ 0 is the free-evolution Hamiltonian superoperator; Ĥ k are the RF field Hamiltonian superoperators corresponding to the different frequency channels and coordinate axes (d is the overall number of degrees of freedom, therefore for 2 RF channels d = 4); and Γ is the thermally corrected 21,22 relaxation superoperator. With a discretized shaped pulse, when the Liouvillian of the individual pulse points is timeindependent, the evolution of the density operator is given by the N propagator through the N time-steps: The density operator at the end of the pulse sequence is: The value of the cost function can be obtained by comparing the final density operator, which is the result after the propagator train, to the preferred density operator. The preferred density operator reflects the desired result, for example, that all of the available magnetization is on the proton x-axis.
The pulse sequence has to be robust to the effects of B 0 and B 1 field inhomogeneities. Therefore, the cost function is evaluated over a range of frequency offsets and B 1 amplitudes. The final fitness of the solution is the weighted average of these cost function values for which the weights account for the width of the resonance and the expected B 1 variation over this range of frequencies. These weights were assigned according to Gaussian distributions.
. In [ 15 N 2 ]uream there are 2 pairs of protons, each pair coupled to one 15 N nucleus. Therefore, the Hamiltonian of the spin system in the presence of RF irradiation in the doubly rotating frame is: The dual-tuned surface coils used enabled simultaneous pulsing on both frequency channels; therefore, the RF pulse terms for both protons and X-nuclei are present in Ĥ RF .

| Approximation used in the spin dynamics simulations
Simulation of the spin dynamics scales exponentially with the number of spins and quickly becomes computationally intensive and often unfeasible. Numerous approaches have been introduced to tackle this problem. [23][24][25][26][27][28] These build mostly on the inherent symmetries of the spin system and reduce the matrix size by means of state space restriction methods. Simulation of the 5 spins in [1-13 C] lactate yields a Liouvillian matrix with dimensions of 1024 × 1024. Exponentiation of such a matrix is computationally expensive. Repeating it at every pulse sequence element for every member of the genetic algorithm population in every iteration makes the optimization unfeasible, even when Krylov propagation 29 is used. State space reduction methods are designed to reduce significantly the dimensionality but are more efficient for larger spin systems. Instead, the following approximation was used: The error of this approximation for noncommuting operators has an upper bound, 30 and the discrepancy is small provided that the relaxation rates≪ effective rotation field. For the optimization process, close approximation is not required if the monotonicity of the cost function is preserved. This means that the location of the minima are preserved, and a minimizer of the transformed cost function is also a minimizer of the original cost function. Splitting of the Liouvillian is especially useful for hyperpolarized experiments, for which the thermal correction in the relaxation superoperator becomes redundant, 14,15,31 yielding a diagonal matrix when the phenomenological relaxation model is used. Therefore, the approximate spin dynamics can be simulated in Hilbert space, which reduces the matrix size to 32 × 32 and, in the uncorrelated random fields relaxation regime, the exponentiation of the relaxation matrix can be performed by using a matrix composed of the exponentials of the entries. The situation is similar for thermally polarized samples with long T 1 times, for which signal recovery during the RF pulse sequence can be neglected. This approximation reduced the computational time by more than a hundredfold. This means that the relaxation superperator had diagonal entries containing only the relaxation rates. Higher-order spin terms were assumed to relax with the sum of the individual relaxation rates. For example, the time required for 1 iteration of the solver for the [1-13 C]lactate spin system decreased from > 500 s to ~2-5 s ( Figure 2A). The polarization transfer sequences required approximately 1000 solver iterations, resulting in a total runtime of approximately 1 h. The J-coupling artefact-free excitation pulse required approximately 10 min. The convergence criterion was met when the decrease in the cost value was < 1% over the course of 50 solver iterations. If a specific problem needs higher performance, or to test optimality, the solver can be restarted; the output solution can be used as 1 member of the initializer population, which is otherwise a completely random sample. In off-line optimization-based design processes, it is often acceptable to allow several days of run time. 6,32

| Validation of the simulations
The spin dynamics in the cost function of the optimization and the predicted transfer values were simulated with in-house MatLab (version 2020a; Mathworks, Natick, MA) scripts. The peak transfer values were confirmed using SpinDynamica (www.spind ynami ca.soton.ac.uk) in Wolfram Mathematica (version 11; Wolfram Research, Inc, Champaign, IL). When relaxation losses were estimated, the PhenomenologicalRelaxationSuperoperator function was used with the indicated T 1 and T 2 relaxation times.

| Achieving frequency selectivity
Following the spin dynamics simulation, a subsequent Bloch simulation was inserted to yield the desired frequency selectivity. This is different from the robustness to frequency offset, which aims to yield the same action on the resonance of interest over a range of frequency offsets, whereas the frequency selectivity aims to prevent unwanted excitation of other resonances. This selectivity is achieved by prescribing frequency regions by means of a weight function (W) with ideally 0 excitation, and the deviation from this ideal excitation profile is used to calculate a penalty term (Equation 8)

| Adiabatic INEPT pulse sequence
The adiabatic INEPT pulse sequence used for proton to carbon polarization transfer in the [1-13 C]lactate phantom is described in Ref. [14]. The transfer block employed 4 ms 90° sinc pulses and 10 ms hyperbolic secant (HS8) adiabatic refocusing pulses. The delay between the center of the first and second 90° pulses was 121.95 ms, and the delay between the center of the second 90° pulse and the start of acquisition was 47.78 ms.

| BINEPT pulse sequence
The pulse sequence was implemented as described in Ref. [33]. The transfer block starts with 2 ms BIR4 pulses segmented into 0.5 ms-1 ms-0.5 ms subpulses. The 3 subpulses were separated by two 2.78 ms delays. The segmented pulse was followed by an evolution period of 2.78 ms, with a 2 ms 180° BIR4 pulse in the middle to refocus anti-phase magnetization.

| J-coupling modulation-free excitation of thermally polarized [2-13 C] lactate and [ 15 N 2 ]urea
Forty-five degree excitation pulses were designed to eliminate J-coupling modulation in [2-13 C]lactate ( 13 C doublet peak) and [ 15 N 2 ]urea ( 15 N triplet peak). Spectra were acquired into 1024 complex points with a spectral width of 4006 Hz. The spectra were compared to those obtained using a 100 μs hard pulse with matching flip angle.

| Imaging a [2-13 C]lactate phantom
Single-shot 3D non-Cartesian images were acquired using a 1-cm diameter 13 C transmit-receive surface coil (Agilent, Palo Alto, CA). Images were acquired using a multi-spin echo pulse sequence. 32 The FOV was 3.2 cm with a nominal image resolution of 2 mm. The excitation pulse was either a 2 ms 45° BIR4 pulse or the optimized pulse designed for elimination of J-coupling modulation. The 7 uniformly spaced spin echoes were generated using 12 ms hyperbolic secant (HS8) adiabatic refocusing pulses without slice selection, with a TE of 25.4 ms. The refocusing pulses were phase-compensated. 35 Readout band-width was set to 250 kHz. The sampled k-space was gridded to a twofold oversampled 32 × 32 × 32 3D Cartesian grid using a Kaiser-Bessel function with density pre-and postcompensation, followed by 3D inverse Fourier-transform and deapodisation. The same pulse sequence with a 2 ms 90° BIR4 excitation pulse was used for T 2 measurements with the readout gradients switched off. At each of the 7 spin echoes, 3157 complex points were acquired; and after Fouriertransform, the 7 peaks were fitted to a monoexponential decay function.

| Dynamic acquisition of J-coupling modulation-free spectra in vivo
A series of 13

| Polarization transfer from proton to carbon in [1-13 C]lactate
Polarization is transferred from the methyl protons to the 13 C nucleus (J = 4.1Hz). However, the methine proton is also coupled to both the methyl protons (J = 7 Hz) and to the carboxyl carbon (J = 3.2 Hz). The simulation was insensitive to the values of the 1 H and 13 C T 1 relaxation times; therefore, published values were used: T 13C 1 = 50 s (at 7 Tesla), 36 T 1H 1 = 1.73 s (at 4.7 Tesla). 37 T 2 relaxation times were measured using a multi-spin echo sequence, 32 giving values of T 13C 2 = 1.625 s and T 1H 2 = 463 ms. The target state for the optimization was to place the 13 C magnetization in-phase in the x-y plane immediately after polarization transfer. The resultant spectrum was compared to spectra obtained without polarization transfer and using the pulse sequence described in Ref. [14], which employed adiabatic refocusing pulses in a refocused INEPT experiment ( Figure 3B). The SNR of the optimized experiment was 1.41 ± 0.19 times higher compared to detection without polarization transfer and 1.29 ± 0.17 times higher than in the INEPT experiment. The simulated SNR improvement of the optimized sequence relative to detection without polarization transfer was 1.42.

| Polarization transfer from proton to carbon in [2-13 C]lactate
Polarization is transferred from the methine proton to the 13 C nucleus (J = 140 Hz). However, the methyl protons are also coupled to both the methine proton (J = 7 Hz) F I G U R E 3 (A) Optimized pulses for transfer of 1 H polarization to 13 C in [1-13 C]lactate. (B) Comparison of the spectrum acquired with the optimized sequence to the spectra acquired with the adiabatic insensitive nuclei enhanced by polarization transfer (INEPT) sequence described in Ref. [14] and by detection without polarization transfer (direct 13 C detection). The latter spectra are shifted by ±1000 Hz. (C) Optimized pulses for transfer of 1 H polarization to 13  and to the C2 carbon (J = 4.2 Hz). Because use of the exact relaxation time values is not critical in cases in which the J-coupling constant is large, the simulation was run with the 1 H and 13 C T 2 , and the 1 H T 1 relaxation times used for [1-13 C]lactate. The 13 C T 1 (T 13C 1 ) was set to 7.2 s. 38 The target for the optimization was to place the 13 C magnetization in-phase in the x-y plane immediately after polarization transfer. The resultant spectrum was compared with the spectrum obtained by detection without polarization transfer and with the conventional INEPT sequence, which used 1-ms hard pulses ( Figure 3D). Adiabatic INEPT sequences could not be used due to the relatively long adiabatic pulses (2-10 ms) when compared to the short interpulse delays in the polarization transfer sequence, which are the result of the large coupling constant between the methine proton and C2 carbon. The SNR of the optimized experiment was 3.81 ± 0.83 times higher compared to detection without polarization transfer and 1.91 ± 0.67 higher compared to the INEPT sequence. The theoretical maximum for the transfer of proton polarization to 13 C is 1H ∕ 13C = 3.97, when neglecting the other much weaker couplings and relaxation. The insensitivity of the 13 C signal to variation in the 13 C and 1 H B 0 and B 1 fields in the optimized experiment was demonstrated by sweeping the 13 C and 1 H center frequencies and changing the transmitter gains. The results obtained were compared with those obtained using the hard pulse INEPT sequence (Figure 4).

| Polarization transfer from proton to nitrogen in [ 15 N 2 ]urea
Polarization is transferred from the 2 pairs of magnetically equivalent protons to each of the two 15 15 N triplet with the optimized sequence was 1.66 ± 0.20 and 9.96 ± 1.83 times higher than that obtained using the B 1 -insensitive heteronuclear adiabatic polarization transfer (BINEPT) pulse sequence and detection without polarization transfer, respectively ( Figure 3F). The optimized sequence better preserved the multiplet structure, which may be lost in the BINEPT sequence due to antiphase magnetization resulting from the relatively long adiabatic pulses. The optimized sequence required a more than twofold lower peak B 1 compared to the BINEPT sequence.

| Polarization transfer from hyperpolarized nitrogen to protons in [ 15 N 2 ]urea
In order to allow serial acquisition of spectra, polarization has to be transferred from the hyperpolarized 15  Simulations showed ΔB 0 and ΔB 1 robustness comparable to the IRRUPT pulse sequence, 15 which is a partial transfer version of the BINEPT sequence ( Figure 5A). Signal was observable for up to 24 s (12 repetitions) even though injection of the hyperpolarized urea solution significantly degraded B 0 homogeneity ( Figure 5C). The required B 1 amplitude for this optimized sequence was approximately twofold lower than that needed for the IRRUPT sequence.

| Pulses immune to the effects of Jcoupling in [2-13 C]lactate and [ 15 N 2 ]urea
An optimized 13 C pulse ( Figure 6A) designed to yield a 45° excitation of the upfield peak in the [2-13 C]lactate doublet ( Figure 6B) and minimal excitation of the pyruvate and pyruvate-hydrate resonances ( Figure 6C) took only 9.76 ms. The algorithm finds the shortest pulse that minimizes signal loss but can still perform selection of the singlet component over the prescribed range of resonance frequency offsets. The pulse, which is applied at the center frequency of the 13 C doublet, steers the 13 C magnetization to a final state, where there is an equal mixture of in-phase and anti-phase 13 C magnetization. A better selectivity was observed between the frequency offsets of interest at +J/2 = 70 Hz and −J/2 = −70 Hz when compared with a 32 ms hard pulse and a 10 ms sinc pulse applied at the resonance frequency of 1 of the components of the doublet (+J/2 = 70 Hz), as described in Ref. [18] (Figure 7). A similar excitation pulse was designed for J-coupling artefact-free imaging of the triplet resonance from [ 15 N 2 ] urea. The analytical derivation is complicated even for a system of 2 coupled spins 18 and does not hold for a system of 3 or more spins. For the numerical optimization approach, it is sufficient simply to define the desired state of the spin system immediately after excitation as that containing a singlet peak. This is a mixture of in-phase (e.g., Î x ,Î y etc.) and certain anti-phase (e.g., Î xŜz , Î yŜz etc.) magnetization terms so that the other components of the multiplet cancel. This is trivial to find and constitutes the input for the optimization. For example, the desired state is ̂ = a x I x + a y I y + 2(a x I x + a y I y )S CH z for the 13 C doublet case and ̂ = b x I x + b y I y − 4 b x I x + b y I y S 1z S 2z for the 15 N triplet. The coefficients a x , a y , b x , b y are arbitrary real numbers such that ̂ obeys density matrix normalization rules. They also indicate that the final phase of the spectrum relative to the receiver is unimportant, but it needs to be same for the in-phase and anti-phase terms so that the respective components of the multiplets cancel. It is important to note here that the antiphase terms relax according to the sum of the relaxation rates of the individual spin operator components. The optimization used the T 1 and T 2 values for [2-13 C]lactate and [ 15 N 2 ]urea given above. The optimized excitation pulse for [ 15 N 2 ]urea and the resultant singlet spectrum are shown in Figure 6D,E, respectively. The pulse was longer, at 13.84 ms, reflecting the smaller J-coupling and therefore slower spin evolution. The 13 C pulse was used for single-shot 3D imaging of the lactate phantom using a multi-spin echo readout. 32 The image acquired with the optimized excitation pulse showed only readout point spread function-associated artefacts, whereas the image acquired using a BIR4 excitation pulse introduced blurring and image distortion originating from J-coupling dependent modulation of the 13 C signals (Figure 8). A series of 13 C spectra acquired from the rat brain following intravenous injection of hyperpolarized [2-13 C]pyruvate using the same optimized pulse contained only the upfield peak of the [2-13 C]lactate doublet, which was observable for up to 40 s. The spectra showed no detectable pyruvate or pyruvate-hydrate signal ( Figure 9).

| DISCUSSION
Metabolic imaging using hyperpolarized 13 C-labeled substrates has provided a new tool for investigating tissue metabolism in vivo [39][40][41][42] and has translated to the clinic with studies in oncology [43][44][45][46][47] and cardiology. [48][49][50][51] Despite the enormous gain in sensitivity, 31 there has, nevertheless, been attempts to further improve sensitivity by transferring polarization from nuclei such as 13 C and 15 N, in which the polarization is long-lived, to protons. A reverse INEPT sequence has been used for 1 H imaging of [1-13 C]lactate in vivo following injection of hyperpolarized [1-13 C]pyruvate. 14,16 However, the use of adiabatic pulses, which allowed its implementation in vivo using a surface coil, required a compromise to be made between the frequency selectivity of these pulses and their robustness to B 1 inhomogeneity. Polarization transfer efficiency can also be reduced by longer adiabatic pulses, which is particularly severe when there is strong J coupling and the polarization transfer element of the pulse sequence is consequently much shorter. 15 The high-energy deposition and peak B 1 amplitude of adiabatic pulses is also a problem for clinical translation.
The pulse sequences generated using the optimization approach described here showed improved performance in transferring polarization from protons to X nuclei These experiments should show a similarly improved performance in transferring polarization from hyperpolarized X nuclei to protons when compared to a reverse INEPT experiment. The INEPT transfer efficiency decreases due to field imperfections, finite pulse lengths, and relaxation. The latter is an issue with small coupling constants, for example, the 4.1 Hz coupling between the C1 carbon and methyl protons in [1-13 C]lactate. Nevertheless, proton detection of hyperpolarized [1-13 C]lactate still improved the SNR. 14 However, to achieve an acceptable specific absorption rate (SAR) for proton detection of [1-13 C] lactate on a 3 Tesla clinical scanner would require a minimum TR of 10 s, which would result in a considerable loss of dynamic information. The SAR of the optimized pulse sequence was 2.3 times lower. Furthermore, the F I G U R E 7 (A) Simulated frequency profile of the 9.76 ms optimized RF pulse designed to give J-coupling artefact-free excitation of [2-13 C]lactate when the pulse is applied exactly on-resonance at the center frequency of the 13 C doublet. The profile is compared to that of a +J/2 = 70 Hz off-resonant 32 ms hard pulse and a 10 ms sinc pulse. (B) The same profiles measured experimentally. (C) 13 C spectra acquired from [2-13 C]lactate with the optimized pulse at different offset values of the doublet center frequency. (D) 13 C spectra acquired from [2-13 C] lactate with the 10-ms sinc pulse applied at +70 Hz from different offset values of the doublet center frequency hyperbolic secant (HS8) refocusing pulses used in this sequence also have a high adiabatic threshold, which may not be achievable on a clinical scanner. The same problem exists for BIR4 pulses. Both BIR4 and HS8 are relatively long, which significantly reduces polarization transfer efficiency when the coupling constants are large. When compared to a polarization transfer sequence that allows serial partial transfers of polarization (IRRUPT), 15 the  13 C images acquired with the pulse sequence described in Ref. [32] using an adiabatic 45° BIR4 excitation pulse. (C), (F) The same images acquired using the optimized excitation pulse shown in Figure 6A F I G U R E 9 (A) 13 C spectra acquired from a mouse brain following injection of hyperpolarized [2-13 C]pyruvate. The inset displays the sum of the acquired spectra (solid line) and the doublet component (dotted line) that would have been obtained without the optimized pulse, which is based on the linewidth and frequency offset of the observed singlet. The excitation pulse was designed to give J-modulationfree imaging of [2-13 C]lactate. (B) Time-course of the singlet peak intensity optimized sequence showed a similar performance but required a much lower peak B 1 field, and therefore had a lower SAR. For example, the maximum B 1 field for a clinical coil is usually less than 1 G for 13 C, whereas a typical BIR4 pulse used in the IRRUPT sequence has an adiabatic threshold of several G. Also, the BIR4 pulse, and therefore the IRRUPT sequence, are not frequency-selective.
Removal of 1 component of the [2-13 C]lactate doublet, either by using a frequency-selective pulse on the other component 18 or the optimized pulse described here applied at the center frequency of the doublet, is achieved in both cases through phase cancelation. However, the 9.76-ms optimized pulse shows a comparable bandwidth to the 10-ms sinc pulse and a much wider bandwidth than 32-ms hard pulse at the frequency of the observed component but, importantly, much less signal from the other component than that observed with either the hard or sinc pulses. Effective removal of 1 of the components of the doublet using the optimized pulse was demonstrated in the phantom images and in the spectra acquired in vivo.
A limitation of the optimization process is the slow speed of convergence, typical of genetic algorithms. In combination with the computational complexity of the spin dynamics simulation in the cost function, the optimization process is limited to small spin systems. However, the hyperpolarized substrates that are of interest in MR metabolic imaging experiments are often small, frequently no more than 5 coupled spins. 52,53 Furthermore, if the coupling of interest is significantly stronger than the others, the further spin-spin pairs can be neglected. Even in the case of larger molecules, the coupling networks are usually short-range, and the scale of simulation can be reduced by means of state-space restriction techniques. The waveform smoothing technique and the approximation applied in the spin dynamics simulation (Equation 8) also markedly decreased the required computational time. A limitation of the approximation (Equation 8) is that accounting for cross-relaxation and other nondiagonal terms in the relaxation model is not straightforward. Nevertheless, the proposed method could have general applications in studies in which a hyperpolarized X nucleus is monitored 40,54-60 and the detection via spin-coupled protons is expected to result in an enhanced SNR. The optimization process could also be employed in the design of excitation and refocusing pulses with prescribed frequency profiles, both with 6 and without 61,62 the presence of J-coupling modulation. The full Bloch simulation in the cost function also allows the optimal design of selective excitation pulses for UTE MRI. [63][64][65] A potential extension would be to include gradient channels in the optimization. This would enable, for example, the optimized design of spectral-spatial pulses, water-suppression sequences, and diffusion-encoding blocks and it is expected to be computationally tractable. However, incorporating gradients into the spin dynamics simulation, for example, for coherence selection, is anticipated to increase significantly the runtime of the cost function evaluation.

| CONCLUSION
In conclusion, the proposed optimization method provides a simple framework for the design of pulse sequences that are robust to the effects of B 1 and B 0 inhomogeneity, typical for in vivo experiments, and reduce energy deposition and therefore could facilitate clinical translation of preclinical imaging sequences.