Bloch‐Siegert B1+‐mapping for human cardiac 31P‐MRS at 7 Tesla

Purpose Phosphorus MR spectroscopy (31P‐MRS) is a powerful tool for investigating tissue energetics in vivo. Cardiac 31P‐MRS is typically performed using surface coils that create an inhomogeneous excitation field across the myocardium. Accurate measurements of B1+ (and hence flip angle) are necessary for quantitative analysis of 31P‐MR spectra. We demonstrate a Bloch‐Siegert B1+‐mapping method for this purpose. Theory and Methods We compare acquisition strategies for Bloch‐Siegert B1+‐mapping when there are several spectral peaks. We optimize a Bloch‐Siegert sensitizing (Fermi) pulse for cardiac 31P‐MRS at 7 Tesla (T) and apply it in a three‐dimensional (3D) chemical shift imaging sequence. We validate this in phantoms and skeletal muscle (against a dual‐TR method) and present the first cardiac 31P B1+‐maps at 7T. Results The Bloch‐Siegert method correlates strongly (Pearson's r = 0.90 and 0.84) and has bias <25 Hz compared with a multi‐TR method in phantoms and dual‐TR method in muscle. Cardiac 3D B1+‐maps were measured in five normal volunteers. B1+ maps based on phosphocreatine and alpha‐adenosine‐triphosphate correlated strongly (r = 0.62), confirming that the method is T1 insensitive. Conclusion The 3D 31P Bloch‐Siegert B1+‐mapping is consistent with reference methods in phantoms and skeletal muscle. It is the first method appropriate for 31P B1+‐mapping in the human heart at 7T. Magn Reson Med 76:1047–1058, 2016. © 2015 The Authors. Magnetic Resonance in Medicine published by Wiley Periodicals, Inc. on behalf of International Society for Magnetic Resonance in Medicine. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.


INTRODUCTION
Phosphorus magnetic resonance spectroscopy ( 31 P-MRS) provides unique insight into the metabolism of the human heart in vivo. 31 P-MRS studies have revealed the role of the creatine kinase energy shuttle-and its constituent "high energy phosphate" metabolites-in supplying energy to drive the contractile work of the heart (1,2). 31 P-MRS has been used to study a range of cardiomyopathies and to stratify risk in patient groups at 1.5 Tesla (T), 3T, and now recently with improved signal-to-noise ratios (SNRs) at 7T (3)(4)(5).
In principal, MR spectra provide a quantitative measure of metabolite concentrations in vivo. Yet, because 31 P-MR spectra typically have $10 25 x lower SNR= ffiffi t p than 1 H MRI, and 31 P-containing metabolites have T 1 relaxation times of several seconds in vivo, it is not normally practical to obtain fully relaxed 31 P-MR spectra in vivo. Using a repetition time (TR) < 5 3 T 1 causes partial saturation. So it is essential to know the flip angle in each voxel (and the metabolite T 1 ) if we are to determine accurate metabolite concentrations, or metabolite concentration ratios (6).
The phosphorus gyromagnetic ratio c31 P 517.235 MHz T -1 , is only 40% of c1 H , so there is less nutation of the magnetization for a given radiofrequency (RF) transmit field strength ðB 1 1 Þ. To mitigate both the low SNR of cardiac 31 P-MRS and the detrimental effect of the low gyromagnetic ratio on the transmit performance surface coils are often used to maximize both B 1 1 and the peak receive sensitivity (B 2 1 ). However, surface coils give large variations in B 1 1 depending on coil placement and the location of the voxel of interest in the heart. B 1 1 -insensitive adiabatic excitation pulses have been used to obtain uniform excitation flip angles despite this B 1 1 -inhomogeneity at 1.5T and 3T, thereby enabling absolute quantification of metabolite concentrations (7,8). But in the heart at 7T, the wide bandwidth of 31 P spectra, the limited peak B 1 1 and the regulatory limits on specific absorption rate (SAR) make adiabatic excitation challenging.
Approaches that compute B 1 1 using the Biot-Savart law or finite-differences time-domain simulations, or that measure B 1 1 in advance in a phantom are demanding, because dielectric effects make the spatial distribution of B 1 1 depend on each subject's anatomy (9)-and more so at 7T than at 1.5T or 3T (10). Therefore, at 7T, because it is challenging to use B 1 1 -insensitive adiabatic pulses and because calculated field maps are increasingly inaccurate, it is essential to be able to directly measure, for each subject, the spatial distribution of B 1 1 in the human heart, which is the focus of this work. There is a limited choice of existing B 1 1 -mapping methods that are compatible with cardiac 31 P-MRS at 7T. The dual-TR, flip-angle measurement method of Chmelik et al is effective in skeletal muscle, but it relies on knowledge of the metabolite T 1 values (11). This is particularly limiting for the highest SNR peak, phosphocreatine (PCr), whose apparent T 1 depends on the creatine-kinase exchange rate (12), which is known to change markedly in heart failure (13).
A second potential method is actual flip-angle imaging (AFI), which has been used for musculoskeletal 31 P-MRS at 7T (14). However, that protocol used a total TR of 4 s; a single average, 16 3 16 3 8 matrix, acquisition weighted, threedimensional (3D) chemical shift imaging (CSI) acquisition would take 40 min, leaving insufficient time for the main acquisition in a patient study.
The above methods illustrate some of the classic difficulties associated with 31 P-MRS and X-nuclear MRS in general, where long, and potentially unknown T 1 values, combined with low SNR, limit the ability to measure accurately a fine change in signal magnitude, hampering the translation of methods from 1 H-MRI and 1 H-MRS. We, therefore, introduce a novel method for B 1 1 -mapping in multivoxel spectroscopy based on the Bloch-Siegert shift (15). This approach encodes the B 1 1 measurement in the phase of the magnetisation, so it is independent of the TR and of metabolite T 1 values. 31 P-MRS, and other X-nuclear systems, typically have complex multi-peak spectra, therefore, we explore how a Bloch-Siegert module can be applied for multi-peak spectroscopy and derive expressions for the accuracy and precision of Bloch-Siegert spectroscopy B 1 1 -mapping. Finally, we optimize a protocol for cardiac 31 P Bloch-Siegert spectroscopy B 1 1 -mapping at 7T. We validate this protocol in phantoms, in human skeletal muscle and demonstrate it in the hearts of five normal volunteers at 7T.
1 -mapping uses an additional off-resonance RF pulse (a "Bloch-Siegert pulse") inserted between excitation and readout in a pulse sequence. During this Bloch-Siegert pulse, there is a transient Bloch-Siegert shift in the Larmor precession frequency (16,17). This causes an accumulated phase shift of the transverse magnetization and hence alters the phase of the spectral peak. The Bloch-Siegert pulse is placed at a large frequency offset (x RF ) cB 1 1 ) from the measured peak, to minimize additional excitation of magnetization by the Bloch-Siegert pulse ("direct excitation"), and to give a linear relationship between the Bloch-Siegert phase shift and B 1 1 (18). In this linear regime, the additional phase / BS (in radians) accumulated by a species over the duration T p of the Bloch-Siegert pulse is where t is time (in s), x RF is the frequency offset (in Hz) from the Bloch-Siegert pulse to the metabolite of interest, c is the gyromagnetic ratio (in Hz T -1 ), B 1 1 ðtÞ is the time dependent pulse amplitude of the Bloch-Siegert pulse (in T), B 1 1;peak is the maximum amplitude of the Bloch-Siegert pulse (in T), and B 2 1;normalized 5B5 Ð TP 0 B 1 1 ðtÞ 2 dt=ðB 1 1;peak Þ 2 is the normalized pulse-envelope squared-integral of the Bloch-Siegert pulse.
In the original 1 H imaging Bloch-Siegert B 1 1 -mapping method (15), two images are acquired with Fermi-shaped (19) Bloch-Siegert pulses at equal and opposite frequency offsets 6x RF relative to the water resonance. In each voxel, the images have phase / i 5/ 0 1/ BS , where / 0 is the normal phase due to the coil transmit and receive phases, B 0 inhomogeneity, sequence dead-time, flow, etc.. B 1 1 is then obtained from the phase difference between the two images Using Eq. [1], cB 1 1;peak is given by: Possible Approaches When There Are Multiple Peaks In 31 P-MRS, and X-nuclear MRS in general, multiple resonances are observed and there is often no dominant peak (unlike 1 H-MR where the water peak dominates). Therefore, for a single choice of Bloch-Siegert pulse center frequency, each peak in the spectrum will experience a different frequency offset, x RF , from the Bloch-Siegert pulse to that metabolite, and therefore, according to Eq. [1], a different Bloch-Siegert phase shift. As the peaks in a spectrum are at well-defined x RF , it is possible to measure B 1 1 from the difference in the Bloch-Siegert phase shift between multiple peaks in a single spectrum (with a single Bloch-Siegert pulse frequency), or from the totality of Bloch-Siegert phase shifts of multiple peaks in multiple spectra (with different Bloch-Siegert pulse frequencies for each acquisition). It is not immediately apparent whether incorporating signals from multiple peaks is beneficial when determining B 1 1 . And it is unclear whether it is more efficient to devote the whole measurement time to a single spectrum (with maximal SNR) using one Bloch-Siegert pulse frequency or whether it is better to acquire several spectra (with lower SNR) at several Bloch-Siegert pulse frequencies.
To investigate, we calculated the accuracy and precision of a representative selection of measurement strategies using Monte Carlo simulations and the propagation-oferrors method, applied to each strategy in turn as follows: 1. Deduce an expression that computes gB þ 1 from the measured metabolite phases, e.g., Eq. [3]. (And where relevant, generalize to an arbitrary number of metabolite peaks.) 2. Estimate experimental accuracy (bias) by Monte Carlo simulation using this expression. For each peak in a simulated spectrum, generate an array of 10 6 normally distributed phases with mean equal to the phase value from Eq. [1] (for the appropriate Bloch-Siegert-pulseto-metabolite-frequency-offset v i and B þ 1 ) and with standard deviation Df. Express the resulting accuracy as a percentage deviation from the true B þ 1 , 3. Compute experimental uncertainty DgB þ 1 by the propagation of errors (20): where @q=@i is the partial derivative of q with respect to the i th measured variable and Di is the uncertainty in the measurement of i. Express the uncertainty as the coefficient of variation in percent ðC v ¼ 100:DgB þ 1 =gB þ 1;true Þ: Each measurement strategy was considered for use in 7T cardiac 31 P-MRS. The simulation variables were set as follows to the study means in reference (21), because the hardware and acquisition method were identical to those used here. The PCr, c2; a2; and b2adenosine triphosphate (ATP) peaks were simulated at frequencies of 0Hz, -300 Hz, -900 Hz, and -1950 Hz, respectively. Single peak simulations used only the PCr peak (at 0 Hz). The phases / were calculated for a 3.5 ms Fermi pulse with cB 1 1;peak 5 277 Hz placed at x RF 562000 Hz. The phase standard deviation D/ was set equal to the Cram er-Rao Lower Bound of the phase, calculated as described in Cavassila et al (22), for the study mean in Rodgers et al (21), divided by the square-root of the number of scans per measurement (CRLBu= ffiffiffiffi ffi N p ). Each measurement strategy was also evaluated for D/ values corresponding to the mean SNR 6 standard deviation in Rodgers et al (21). The derived equations and results are given in Table 1.

Method A: Dual-Acquisition þ/À
This approach is equivalent to the original Bloch-Siegert B 1 1 -mapping MRI technique (15). Two acquisitions are made with Bloch-Siegert pulses at offsets symmetrically either side of the target metabolite. This approach computes B 1 1 using signal from a single peak.

Method B: Dual-Acquisition On/Off
Single-peak (B 0 ) This approach uses one acquisition without a Bloch-Siegert pulse and one acquisition with the Bloch-Siegert pulse on one side of the resonances. B 1 1 is computed from a single peak.

Multi-peak (B 00 )
This acquisition also allows simultaneous use of data from all peaks that satisfy the linearity condition: cB 1 1;peak ( jx i j. (N.B. each peak may have a positive or negative offset.) We combine the measured B 1 1 values from each resonance peak in a maximum likelihood sense, which is equivalent to the outcome that would be achieved from least-squares fitting of the linear Eq. [1] to the simulation data. The equations generated can be expressed for any number of resonances n at offsets x i from the Bloch-Siegert pulse.

Method C: Single-Acquisition Method
Measurement of B 1 1 from a single Bloch-Siegert sensitized acquisition may be achieved for a multipeak spectrum, providing the first-order phase correction and excitation phase profile are known. For n peaks, a single acquisition is made with a Bloch-Siegert pulse placed so that the linearity condition is satisfied for all the peaks to be included in the measurement. To calculate B 1 1 Eq. [2] is applied to the phase differences between each pair of peaks and this set of B 1 1 values may be averaged in a maximum-likelihood sense.  A fuller description of the methods are given in the multi-peak approaches subsection of Theory. The values used to calculate the bias and standard deviation are taken from the study means of ref (21), the bias is calculated using Eq. [4] with gB þ A generalized closed form could not be found for the uncertainty, DgB þ 1;peak , of Method C (the single-acquisition method), however an expression may be found by expanding the equation for determining B þ 1 for a specific number of peaks and then applying Eq. [5]. We provide a Mathematica notebook in the supporting information, for the four peak case that could easily be adapted for other situations.

Bloch-Siegert Pulse Design
An ideal Bloch-Siegert pulse maximizes / BS , for some anticipated range of B 1 1 , at the same time minimizing direct excitation, and minimizing signal loss due to any additional dead-time between excitation and readout (when the transverse magnetization experiences T 2 * relaxation). Providing one remains in the linear regime (x RF ) cB 1 1 ; maximum cB 1 1 % 1200 Hz) and assuming that the SAR contribution from the excitation pulse is negligible, Eq. [6] in Duan et al (18) shows that the Bloch-Siegert phase difference is limited only by the SAR limit. Therefore, an optimized pulse must be chosen by minimizing direct excitation at the target frequency offset and by minimizing T 2 * -induced signal losses.
A commonly employed Bloch-Siegert pulse is the Fermi pulse. We use the definition in Eq. [4.14] of Bernstein et al (19): where A F is the B 1 1 -field peak amplitude, x RF is the pulse frequency offset and T 0 and a are two adjustable parameters having the dimension of time.

METHODS
All experiments used a Magnetom 7T scanner (Siemens, Erlangen, Germany). Localizers were acquired with a 10 cm 1 H loop coil (Rapid Biomedical, Rimpar, Germany). This was replaced by a T/R switch and preamplifier module (Virtumed LLC, MN, USA) connected to a 10-cm 31 P transmit-receive loop coil (21) for the 31 P acquisitions. The 31 P coil was tuned and matched using an RF sweeper (Morris Instruments Inc, Ottawa, Canada) for each subject. All subjects in this study were recruited in a manner approved by the local Research Ethics Committee.

Fermi Pulse Optimization
The suitability of a range of Fermi pulses was examined using a Bloch simulation grid search. T P was varied from 1 ms to 10 ms in steps of 0.5 ms, T 0 was varied from 0:1T P to 0:5T P in steps of 0:05T P and a was varied from 0:01T P to 0:3T P in steps of 0:001T P . Spin evolution was simulated in Matlab (Mathworks, Natick, MA) at cB 1 1 values ranging from 100 Hz to 1000 Hz in steps of 100 Hz. The deviation from the ideal flat stop-band profile was quantified as The minimum DM xy of the shortest pulse which maintained a DM xy < 0:01 across all simulated peak B 1 1 values had T P 5 3:5 ms; T 0 50:875 ms; and a50:224 ms. These parameters were used for all experiments presented below.
To assist in the implementation of the Bloch-Siegert method and pulse optimization step, we provide the Matlab implementation of this optimization in the file "FermiPulseOptimisation.m" in the Supporting Information, which is available online.

Point-Source Phantom Validation
The Bloch-Siegert effect was demonstrated in a saline-filled (18 L of 73 mM NaCl) phantom containing a single, small cube (2 3 2 3 2 cm 3 ), of 1.8 M K 2 HPO 4 , which gives one singlet signal that originates only from the cube. The optimized Fermi pulse (T P 53:5 ms; T 0 50:875 ms; a50:224 msÞ was inserted between the excitation RF pulse and the 3D phase-encoding gradients of a chemical shift imaging (CSI) sequence (Fig. 1a) (21,23). The Fermi pulse offset was swept from -10 kHz to 110 kHz.
Acquisition parameters were: 4 3 4 3 4 CSI matrix, 80 3 80 3 80 mm 3 field of view, 2048 spectral points, 6 kHz bandwidth, 200 ls hard excitation pulse with a transmit voltage, i.e., the RMS voltage at the "coil plug" on the patient table, of 270 V. The voxel covering the 2 3 2 3 2 cm 3 phosphate cube was selected for analysis. Spectra were fitted using a custom Matlab implementation of the advanced method for accurate, robust and efficient spectral fitting (AMARES) algorithm (24), with prior knowledge specifying a single unconstrained peak (25).
A cB 1 1;peak value was measured using a series of fully relaxed, nonlocalized free induction decays (FIDs) with a 4-ms hard pulse transmitted at 10-200 V. Under these conditions, where T R > 5T 1 , the amplitude of the observed FIDs will follow a sina relationship. The FIDs were fitted in Matlab, and then a sinusoid was fitted to the complex peak amplitudes with two adjustable parameters: the maximum signal amplitude, and the B 1 1 -per-volt scaling factor that relates the (known) applied transmit voltage to the observed signal amplitudes (i.e. flip angle). cB 1 1;peak may be calculated as the product of the scaling factor and the Bloch-Siegert pulse transmit voltage and subsequently used to predict the phase response for the Bloch-Siegert pulse.

Uniform Phantom Validation
Bloch-Siegert B 1 1 -mapping was validated against a multi-TR magnitude reference method in a uniform phantom. The phantom was a 120 3 270 3 270 mm 3 box containing 0.04 M K 2 PO 4(aq) , giving rise to a singlet signal from throughout the whole volume of the phantom. The T 1 of the phantom was separately determined to be 8.57 s. The acquisition parameters were: 300 ls hard pulse excitation transmitted at 270V, 6 kHz bandwidth, 2048 spectral points, 150 3 320 3 320 mm 3 field of view, acquisition weighting, 16 3 8 3 8 resolution, the first dimension perpendicular to the plane of the coil. The multi-TR method was adapted from a previously published dual-TR method (11), acquiring spectra at eight TRs (0.5, 1, 2, 3, 4, 6, 8, 10 s), to ensure functional sensitivity to a broad range of flip angles simultaneously. The number of averages at k 5 0 was adjusted for each acquisition, to achieve similar high SNR for a "best possible" validation; acquisitions times were between 80 and 120 min each (giving a total of 12.5 h run overnight). The partial saturation equation, [8] was fitted to the magnitude data to determine the flip angle, cB 1 1 was calculated from the relationship h5 180 :T P :cB 1 1 =500Hz where T P is the duration of the hard excitation pulse.
Bloch-Siegert acquisitions used the same parameters, with the optimized Fermi pulse (T P 5 3:5ms; T 0 50:875 ms; a50:224Þ inserted as in Figure 1a. Two acquisitions were made with the pulse placed at 12 kHz and -2 kHz offsets. TR was 500 ms, with 13 averages at k 5 0, resulting in a total duration of 2 3 14 min. To avoid artifacts due to phase-wrap present in the raw phase maps, phase differences were calculated using complex division and the Matlab function atan2() that returns phase between 2p and p, then cB 1 1 was computed using Eq. [3]. Voxels with D/ < 0 (indicating either very low cB 1 1 or the presence of uncorrected wrap artefacts in the phase difference map, i.e., true D/ > 180 ), voxels outside the phantom, or voxels with a CRLB / > 20 (indicating low SNR) were excluded.

In Vivo Validation (Thigh)
Full details of the acquisition strategy analysis are given in "Simulation results" below. Method A (dual-acquisition 1/2) showed highest precision and accuracy and was used for all measurements in this work.
Bloch-Siegert B 1 1 -mapping was compared with a published dual-TR flip-angle mapping method (11)  In the dual-TR comparison, excitation was centered on a-ATP (rather than PCr), to avoid complications due to the creatine-kinase mediated exchange of PCr and c-ATP (6). A literature value of T 1;a2ATP 5 1.8 s was used to calculate the flip-angle according to the published method (26). Three acquisitions were made with TRs of 2.2 s, 0.73 s, and 0.37 s, using 10, 20, and 43 averages at the center of k-space for acquisition times of 25, 15, and 15 min, respectively, as recommended in Chmelik et al (11).
Two Bloch-Siegert acquisitions were made with the optimized Fermi pulse (T P 53:5 ms; T 0 50:875 ms; a5 0:224Þ placed at 12 kHz and -2 kHz relative to PCr; excitation was centered at PCr. TR was 500 ms, with 25 averages at k 5 0, resulting in a total duration of 2 3 12.5 min. Fitting was performed with AMARES in Matlab using prior knowledge specifying 10 Lorentzian peaks (a,b,c-ATP multiplet components, PCr, phosphodiesters [PDE], and inorganic phosphate [P i ]), with fixed amplitude ratios and scalar couplings for each multiplet. cB 1 1 was calculated as described in "Uniform phantom validation" above. Voxels were excluded from the comparison if the voxel was centered outside the thigh, D/ < 0 or CRLB / >10 .  1 -mapping sequence. The Bloch-Siegert sensitizing pulse is inserted between the excitation pulse and the readout module. In this work, the readout module consists of 3D phase encoding gradients, acquisition of a free induction decay, spoiler gradients and a final delay to produce the desired TR. The sequence is adapted from that in Figure 1 of reference (21). b: Illustration of the real part of 31 P spectra acquired for each of the measurement strategies detailed in the Theory Section and corresponding to a row in Table 1. The peaks used for analysis in each case are shown in blue. The position of the Bloch-Siegert pulses are shown by red arrows. The phase accumulated by each peak is proportional to the inverse of its frequency offset from the Bloch-Siegert pulse (Eq. [1]). Note that when the Bloch-Siegert pulse is placed at -2 kHz, the b-ATP peak is almost entirely saturated. Fitting was performed with AMARES in Matlab using prior knowledge specifying 11 Lorentzian peaks (a,b,c-ATP multiplet components, PCr, PDE, and the two peaks of 2,3-diphosphoglyceric acid [DPG]), with fixed amplitude ratios and scalar couplings for each multiplet. cB 1 1 was calculated as described in "Uniform phantom validation" above. Voxels were excluded if D/ < 0 or CRLB / >15. Phase wrap, caused by the upper output limit of Eq. [9] was unwrapped manually.

Cardiac Scans
In three subjects, the scans were repeated with the central frequency adjusted to a-ATP to allow comparison of cB 1 1 values obtained using different peaks.

Simulation Results
The precision and accuracy of the single-peak measurement strategies (A and B 0 ) are shown in Figure 2. Methods A and B 0 show <10% error and <10% uncertainty for 1000 Hz < x RF < 3000 Hz and cB 1 1 > 150 Hz. Method B 0 has both precision and accuracy ffiffiffi 2 p worse than Method A . These calculations do not take into account the effects of direct excitation, which is shown in the Bloch simulations. The error caused by the direct excitation is generally larger than the underlying error arising from D/, emphasizing the importance of effective pulse design.
The analysis of the multi-peak methods (B 00 and C) results in smoothly varying functions (Fig. 3). Increasing c B 1 1;peak always improves the precision and accuracy. In Figure 3a (Method B 00 ) minimizing jx i j, maximizes / BS and improves the precision and accuracy, while adding additional peaks (x 2;3;... Þ with a larger jx RF j has only a small effect. In Figure 3c, Method C shows lower accuracy and higher bias than the other methods when all the offsets have the same sign. However, if the Bloch-Siegert pulse is placed symmetrically between two peaks, we calculate an identical profile to Method A with a ffiffiffi 2 p improvement in both accuracy and precision for a matched scan duration, although this will only be possible in situations where the inter-peak separation is sufficient to avoid direct excitation. Table 1 summarizes the methods considered for our 7T cardiac 31 P-MRS application. All the measurement strategies show reasonable accuracy (<10% error) and most show an acceptable precision. Method A (dualacquisition 1/2) shows better precision and accuracy than any other method (0.2% error and 10% standard deviation). Methods B 00 and C, were simulated with only positive offsets from the Bloch-Siegert pulse for our specific cardiac application. The greatest peak separation (PCr to b-ATP) in 7T 31 P-MRS is 1950 Hz, which is smaller than the separation required to maintain acceptable levels of direct excitation using the optimized 3.5 ms Fermi pulse. In Bloch simulation, using the optimized 3.5 ms pulse, a separation of 1950 Hz (x 1 52975 Hz; x 2 5975 Hz) gives a B 1 1 error of 19%.

Pulse Design
The Fermi pulse selection criteria were to minimize the pulse duration whilst minimizing the direct excitation of the magnetization in the interval 1500Hz < x RF < 2500Hz. The minimum of the direct excitation metric, DM xy (Eq. [7]), fell below 0.01% for pulses with T p ! 3:5ms (Fig. 4d). Thus the criteria, to select the shortest pulse whilst also minimizing direct excitation, were satisfied by selecting the T 0 and a parameters corresponding to the DM xy minimum of the parameter range with T P 53.5 ms. The resulting optimized Fermi pulse had T P 53:5 ms; T 0 50:875 ms; a50:224, and was used for all experiments and simulations in this manuscript. Figure 4 shows a range of possible Fermi pulses, the degree of direct excitation, and the B 1 1 error from each.

Point-Source Phantom Validation
The phosphate peak phase and magnitude in the single voxel phantom are displayed in Figure 5a. The predicted Bloch-Siegert phase, calculated using the B 1 1 separately measured by nonlocalized, fully relaxed FIDs, and Eq. [1], and the directly measured Bloch-Siegert phase of the single phosphate peak agree within 20% for jx RF j > 1000 Hz and within 8% for 1000 Hz jx RF j 3000 Hz. As jx RF jincreases from 1000 Hz the CRLB / is constant but the Bloch-Siegert phase / BS decreases, resulting in a higher relative uncertainty. When jx RF j < 1000 Hz direct excitation of the phosphate peak by the Fermi pulse is observed, which causes the magnitude to drop substantially. For j x RF j > 1000 Hz the magnitude varies by <1.5%, indicating that direct excitation has been minimized. Figure 5b confirms that cB 1 1 calculated from the Bloch-Siegert phases closely approaches that given by the Multi-FA reference method when 500 Hz < jx RF j < 4000 Hz.

Uniform Phantom Validation
Bloch-Siegert B 1 1 maps across the four central CSI slices covering the uniform phantom are shown in Figure 6a. The dynamic range of the Bloch-Siegert cB 1 1 was 0 Hz to 873 Hz (corresponding to a p phase shift). Figures 6b,c show the per-voxel correlation and Bland-Altman plots of the Bloch-Siegert method and the multi-TR method. A strong correlation (Pearson's r 5 0.90) was observed. The Bland-Altman plot shows the Bloch-Siegert method measured cB 1 1 15 Hz lower than the reference method on average, with 95% confidence intervals of 6154 Hz. The bias and confidence intervals are equivalent to -3.1% and 31.8%, respectively, of the cB 1 1 value (484 Hz) of a centrally located voxel in the phantom. compared with the Dual-TR method, with 95% confidence intervals of 6211 Hz. The bias and confidence intervals are equivalent to -6.9% and 33.4%, respectively, of the cB 1 1 value (631 Hz) of a centrally located voxel in the skeletal muscle.

Cardiac Scans
Five mid-septal cardiac cB 1 1 maps are shown overlaid on localizers in Figure 8. The Bloch-Siegert cB 1 1 dynamic range was 0 Hz to 1231 Hz. The short-axis localizers were manually segmented to extract voxels in the heart and skeletal muscle in the three CSI slices corresponding to the mid, apical and basal segments of the myocardium. In the mid-septal slices 39.6 6 16.0 (mean 6 SD) voxels were selected. After masking by CRLB / 30.6 613.7 voxels remained. The B 1 1 maps are smoothly varying from 100-1200 Hz without visible boundaries between skeletal muscle, myocardium and the ventricular blood pools. B 1 1 drops-off with increasing distance from the coil. Data for the adjacent apical and basal slices, after masking by CRLB / , had 28 6 10.1 and 21.6 6 15.6 voxels suitable for analysis, respectively, which were used in the subsequent comparisons.
Bloch-Siegert measurement using a-ATP instead of PCr is compared in Figures 9a,b. The mean Pearson's correlation was 0.62. The Bland-Altman plot shows a bias of 22.3 Hz with 95% confidence intervals of 6561 Hz. This result confirms the T 1 -independence found in Sacolick et al (15) and is an indication that the B 1 1 values measured in this work are correct.

DISCUSSION
In this study, we proposed a Bloch-Siegert spectroscopy B 1 1 -mapping technique suitable for 31 P-MRS in the human heart at 7T. It is straightforward to register the B 1 1 maps with other 31 P scans using the same 3D CSI readout.

T 1 Insensitivity
Bloch-Siegert methods are known to be T 1 -insensitive (15). We verified this for our 31 P-MRS method by comparing B 1 maps derived from PCr (T 1 5 3.09s) and a-ATP (T 1 5 1.39 s) in three normal volunteers. This is an advantage over multi-TR and dual-angle methods, which rely on metabolite T 1 values to calculate B 1 1 . In normal volunteers, one could use cardiac 31 P T 1 values from the literature in a multi-TR or dual-angle method (21). Yet, in studies of patients with cardiac disease, or in protocols involving exercise or pharmacological stress, these T 1 values may be different, as has been reported for other nuclei (27). Recording 3D-localized myocardial T 1 values at 7T took 45 min in Rodgers et al (21)-too long to be a calibration step in patients.
The problem of T 1 -changes in disease is particularly acute for the highest SNR peak: PCr. Chemical exchange through the creatine-kinase shuttle dominates the apparent T 1 s of PCr and c-ATP, (6) and the forward rate constant for this process changes markedly in disease (28). Therefore, T 1 -sensitive B 1 1 -mapping methods should only be applied to nonexchanging peaks such as a-ATP or b-ATP, but these peaks have lower SNR than PCr. In contrast, exchange mediated phase effects during the short (3.5 ms) Fermi pulse are predicted to be <0.1% by Bloch-McConnell simulation (not shown). Consequently, Bloch-Siegert methods can be applied to the highest SNR peak: PCr.
The Bloch-Siegert shift is independent of the excitation pulse or TR of the host sequence (15). Therefore, a Bloch-Siegert B 1 1 -mapping sequence can be run with the Ernst flip-angle (29) and a short TR to maximize its SNR efficiency (i.e., SNR/ ffiffi t p Þ. Bloch-Siegert 1 H MRI B 1 1 - 8. a: Bloch-Siegert B þ 1 maps from slices of an 16 Â 8 Â 8 CSI grid in a mid-short axis plane of five healthy volunteer subjects. The maps are overlaid on 1 H localizers registered to the CSI grid. b,c: Example magnitude and real spectra from a mid-interventricular septal voxel, marked by a black "x" in the lowest B þ 1 map in a. The maps in a are calculated from the phase of the PCr peak. The phase difference between the PCr peak of the spectrum with a 2 kHz offset (blue) on the Fermi pulse and the -2 kHz offset (red) is 86.9 corresponding to a value of gB þ 1 ¼ 605 Hz. mapping was shown to be more SNR-efficient than dualangle B 1 1 -mapping (30). Theory predicts a similar SNRefficiency advantage with 31 P-MRS.
We, therefore, believe that Bloch-Siegert B 1 1 -mapping approaches are the only known viable methods for human cardiac 31 P B 1 1 -mapping at 7T.

Bloch-Siegert Sampling Strategy
The framework to assess spectroscopic Bloch-Siegert measurements established that the dual-acquisition approach is optimal for cardiac 31 P-MRS, because there is insufficient frequency separation to place the Fermi pulse between PCr and b-ATP without creating artefacts due to direct excitation. In a situation where the peak separation is much wider (4000 Hz) or where T 2 * is sufficiently long to permit longer, narrower bandwidth Bloch-Siegert sensitizing pulses, placing the pulse between two peaks will give an improvement in accuracy and precision of ffiffiffi 2 p over the analogous Method A measurement, or equivalently, it would achieve the same results but in half the total scan time. Furthermore, the single-acquisition Method C could still be attractive for time limited acquisitions, e.g., hyperpolarized 13 C-MRS.

In Vivo B 1 Maps
We believe these are the first 31 P-MRS B 1 1 -maps made in the human heart at 7T. Therefore, there was no goldstandard method to compare against in the heart at 7T. Instead, we validated our method in phantoms and in the leg, where it was possible to obtain reference data with a published method. In both validation experiments, the Bloch-Siegert methods accurately matched the chosen reference methods across the whole dynamic range. The final cardiac maps present a physically reasonable range 100-1000 Hz of cB 1 1 (5-60 lT B 1 1 ) for the experimental setup; B 1 1 decreases smoothly with increasing distance from the surface coil, and repeated scans in three volunteers on PCr and then on a-ATP showed excellent reproducibility. Therefore, we deem that the implementation of a Bloch-Siegert B 1 1 -mapping sequence has been successful.

Limitations
The Fermi pulse optimization presented here is valid for the range of B 1 1 and SAR/ðB 1 1 Þ 2 appropriate for our coil, and for T 2 values typical in the human heart. This optimization ought to be repeated for studies that have markedly different values for these parameters. This is because there is a trade-off between using short, high-amplitude Fermi pulses which reduce the dead-time and minimize T 2induced signal losses but which give higher SAR and need to be placed at a greater offset frequency x RF for acceptable levels of direct excitation, and which, therefore, produce less Bloch-Siegert phase shift for a given B 1 1 (and hence lead to lower B 1 1 precision). Or conversely, long, lowamplitude Fermi pulses which suffer greater T 2 -induced signal losses but produce lower SAR and may be placed at a smaller offset frequency x RF for greater B 1 1 precision. The 3D-CSI sequence localizes signals to comparatively large voxels. For a surface coil, there will be spins in each voxel that experience appreciably different B 1 1 , and, therefore, experience different Bloch-Siegert phase shifts. This causes intravoxel signal cancellation and, therefore, reduces the SNR efficiency of the Bloch-Siegert method. Furthermore, because the Bloch-Siegert phase depends on ðB 1 1 Þ 2 , the measured voxel B 1 1 value will be a weighted average of B 1 1 within the voxel, not the mean voxel B 1 1 value. If a surface receive coil is used, then B 2 1 will also vary across each voxel. B 2 1 inhomogeneity will weight B 1 1 values measured by any method. To confirm that the additional B 1 1 weighting due to intravoxel difference in the Bloch-Siegert phase shift is small, we simulated 512 3 spin isochromats for a voxel 10 cm from the coil using the CSI protocol and RF coil described above. The Bloch-Siegert B 1 1 deviated from the B 2 1 -weighted mean of the B 1 1 values in this voxel by <7%. The drop in SNR due to phase cancelation was <10% for cB 1 1 < 500 Hz but increased to 50 % at 1000 Hz. Note that this simulation is separate to the simulations described in the Theory section.
We note that during the preparation of this manuscript, a single voxel spectroscopy Bloch-Siegert B 1 1 measurement technique, based on Point RESolved Spectroscopy (PRESS) has been published (31). The PRESS Bloch-Siegert method allowed highly reproducible measurements of B 1 1 , using the 1 H water signal, in the hearts of six volunteers. The focus of this manuscript has been to implement these promising Bloch-Siegert methods in a multivoxel X-nuclear sequence, addressing the challenges B 1 1 -mapping in in the presence of low SNR and multiple significant peaks.  9. a: Three volunteer comparison of consecutive Bloch-Siegert B þ 1 maps, collected with the Fermi pulse placed symmetrically around the PCr peak and then subsequently the a-ATP peak in three of the subjects in the study. b: Bland-Altman plot of the metabolite comparison. The means and 95% confidence intervals are calculated from all points shown.