Bloch-Siegert B1+-mapping for human cardiac (31) P-MRS at 7 Tesla.

PURPOSE: Phosphorus MR spectroscopy ((31) P-MRS) is a powerful tool for investigating tissue energetics in vivo. Cardiac (31) P-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 (31) P-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 (31) P-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 (31) P 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 (31) P Bloch-Siegert B1+-mapping is consistent with reference methods in phantoms and skeletal muscle. It is the first method appropriate for (31) P B1+-mapping in the human heart at 7T. Magn Reson Med 76:1047-1058, 2016.

INTRODUCTION

Phosphorus magnetic resonance spectroscopy (31P‐MRS) provides unique insight into the metabolism of the human heart in vivo. 31P‐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. 31P‐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 31P‐MR spectra typically have x lower than 1H MRI, and 31P‐containing metabolites have T1 relaxation times of several seconds in vivo, it is not normally practical to obtain fully relaxed 31P‐MR spectra in vivo. Using a repetition time (TR) < 5 × T1 causes partial saturation. So it is essential to know the flip angle in each voxel (and the metabolite T1) if we are to determine accurate metabolite concentrations, or metabolite concentration ratios 6.

The phosphorus gyromagnetic ratio =17.235 MHz T‐1 , is only 40% of , so there is less nutation of the magnetization for a given radiofrequency (RF) transmit field strength . To mitigate both the low SNR of cardiac 31P‐MRS and the detrimental effect of the low gyromagnetic ratio on the transmit performance surface coils are often used to maximize both and the peak receive sensitivity ( ). However, surface coils give large variations in depending on coil placement and the location of the voxel of interest in the heart.

‐insensitive adiabatic excitation pulses have been used to obtain uniform excitation flip angles despite this ‐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 31P spectra, the limited peak and the regulatory limits on specific absorption rate (SAR) make adiabatic excitation challenging.

Approaches that compute using the Biot‐Savart law or finite‐differences time‐domain simulations, or that measure in advance in a phantom are demanding, because dielectric effects make the spatial distribution of 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 ‐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 in the human heart, which is the focus of this work.

There is a limited choice of existing ‐mapping methods that are compatible with cardiac 31P‐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 T1 values 11. This is particularly limiting for the highest SNR peak, phosphocreatine (PCr), whose apparent T1 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 31P‐MRS at 7T 14. However, that protocol used a total TR of 4 s; a single average, 16 × 16 × 8 matrix, acquisition weighted, three‐dimensional (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 31P‐MRS and X‐nuclear MRS in general, where long, and potentially unknown T1 values, combined with low SNR, limit the ability to measure accurately a fine change in signal magnitude, hampering the translation of methods from 1H‐MRI and 1H‐MRS. We, therefore, introduce a novel method for ‐mapping in multivoxel spectroscopy based on the Bloch‐Siegert shift 15. This approach encodes the measurement in the phase of the magnetisation, so it is independent of the TR and of metabolite T1 values.

31P‐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 ‐mapping.

Finally, we optimize a protocol for cardiac 31P Bloch‐Siegert spectroscopy ‐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.

THEORY

Bloch‐Siegert Mapping

Bloch‐Siegert ‐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 ( ) 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 18. In this linear regime, the additional phase (in radians) accumulated by a species over the duration T p of the Bloch‐Siegert pulse is where t is time (in s), is the frequency offset (in Hz) from the Bloch‐Siegert pulse to the metabolite of interest, is the gyromagnetic ratio (in Hz T‐1), is the time dependent pulse amplitude of the Bloch‐Siegert pulse (in T), is the maximum amplitude of the Bloch‐Siegert pulse (in T), and is the normalized pulse‐envelope squared‐integral of the Bloch‐Siegert pulse.

In the original 1H imaging Bloch‐Siegert ‐mapping method 15, two images are acquired with Fermi‐shaped 19 Bloch‐Siegert pulses at equal and opposite frequency offsets relative to the water resonance. In each voxel, the images have phase , where is the normal phase due to the coil transmit and receive phases, B0 inhomogeneity, sequence dead‐time, flow, etc.. is then obtained from the phase difference between the two images Using Eq. (1), is given by:

Possible Approaches When There Are Multiple Peaks

In 31P‐MRS, and X‐nuclear MRS in general, multiple resonances are observed and there is often no dominant peak (unlike 1H‐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, , 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 , it is possible to measure 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 . 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‐of‐errors method, applied to each strategy in turn as follows:

Deduce an expression that computes from the measured metabolite phases, e.g., Eq. (3). (And where relevant, generalize to an arbitrary number of metabolite peaks.)

Estimate experimental accuracy (bias) by Monte Carlo simulation using this expression. For each peak in a simulated spectrum, generate an array of 106 normally distributed phases with mean equal to the phase value from Eq. (1) (for the appropriate Bloch‐Siegert‐pulse‐to‐metabolite–frequency‐offset and ) and with standard deviation . Express the resulting accuracy as a percentage deviation from the true ,

Compute experimental uncertainty by the propagation of errors 20: where is the partial derivative of with respect to the measured variable and is the uncertainty in the measurement of i. Express the uncertainty as the coefficient of variation in percent

Each measurement strategy was considered for use in 7T cardiac 31P‐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, adenosine 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 = 277 Hz placed at . The phase standard deviation was set equal to the Cramér‐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 ( ). Each measurement strategy was also evaluated for values corresponding to the mean SNR ± standard deviation in Rodgers et al 21. The derived equations and results are given in Table 1.

Table 1

Comparison of Possible Bloch‐Siegert Spectroscopy Measurement Strategies Applied to 31P‐MRSa

Measurement strategy	Peaks (31P equivalent)	Scan 1	Scan 2	B1+ Determination equation γB1,peak+	Propagation of error ΔγB1,peak+	Bias ± SD (Hz) SNR = 18 (μ‐1σ)	Bias ± SD (Hz) SNR = 31 (μ)	Bias ± SD (Hz) SNR = 44 (μ+1σ)
A: Dual‐acquisition +/‐	1 (PCr)	3.5ms Fermi +2000Hz	3.5ms Fermi −2000Hz	(ϕ_2−ϕ_1)ω_RF2πB˜	ω(Δϕ+2+Δϕ−2)8π(ϕ_+−ϕ_−)B˜	2.21 ± 53.7	0.52 ± 27.9	0.28 ± 21.5
B′: Dual‐ acquisition on/off	1 (PCr)	3.5ms Fermi +2000Hz	No Bloch‐Siegert pulse	(ϕ_on−ϕ_off)ω_RFπB˜	ω(Δϕon2+Δϕoff2)4π(ϕ_on−ϕ_off)B˜	8.19 ± 76.0	2.44 ± 39.5	1.32 ± 30.3
B″: Dual‐ acquisition on/off	4 (PCr, γ‐ATP, α‐ATP, β‐ATP)	3.5ms Fermi +2000Hz	No Bloch‐Siegert pulse	∑i=1nπ^−12ωi−2(ϕion−ϕioff)ω_iB˜∑i=1nωi−2	∑i=1n(Δϕion)^2+(Δϕioff)^24πB˜(∑i=1nωi−1)^2(ϕion+ϕioff)ω_i	9.03 ± 70.1	3.97 ± 36.5	2.27 ± 28.0
C: Single‐ acquisition	4 (PCr, γ‐ATP, α‐ATP, β‐ATP)	3.5ms Fermi +2000Hz	No scan	∑i=2,j=1n,i−1ϕ_j−ϕ_iB˜(ωj−1−ωi−1)^2πωj−1−ωi−1∑i=2,j=1n,i−1(ωj−1−ωi−1)^2	b	8.67 ± 153.8	7.34 ± 83.1	5.58 ± 63.9

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 .

A generalized closed form could not be found for the uncertainty, , of Method C (the single‐acquisition method), however an expression may be found by expanding the equation for determining 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.

Method A: Dual‐Acquisition +/−

This approach is equivalent to the original Bloch‐Siegert ‐mapping MRI technique 15. Two acquisitions are made with Bloch‐Siegert pulses at offsets symmetrically either side of the target metabolite. This approach computes using signal from a single peak.

Method B: Dual‐Acquisition On/Off

Single‐peak (B′)

This approach uses one acquisition without a Bloch‐Siegert pulse and one acquisition with the Bloch‐Siegert pulse on one side of the resonances. is computed from a single peak.

Multi‐peak (B″)

This acquisition also allows simultaneous use of data from all peaks that satisfy the linearity condition: . (N.B. each peak may have a positive or negative offset.) We combine the measured 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 at offsets from the Bloch‐Siegert pulse.

Method C: Single‐Acquisition Method

Measurement of 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 Eq. (2) is applied to the phase differences between each pair of peaks and this set of values may be averaged in a maximum‐likelihood sense.

Bloch‐Siegert Pulse Design

An ideal Bloch‐Siegert pulse maximizes , for some anticipated range of , 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 T2 * relaxation). Providing one remains in the linear regime ( ) 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 T2 *‐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 is the ‐field peak amplitude, is the pulse frequency offset and and 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 1H 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 31P transmit–receive loop coil 21 for the 31P acquisitions. The 31P 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. was varied from 1 ms to 10 ms in steps of 0.5 ms, was varied from to in steps of and was varied from to in steps of . Spin evolution was simulated in Matlab (Mathworks, Natick, MA) at 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 of the shortest pulse which maintained a across all simulated peak values had . 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 × 2 × 2 cm3), of 1.8 M K2HPO4, which gives one singlet signal that originates only from the cube. The optimized Fermi pulse ( 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 +10 kHz.

Figure 1

a: Pulse sequence timing diagram for a Bloch‐Siegert spectroscopy ‐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 31P 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 β‐ATP peak is almost entirely saturated.

Acquisition parameters were: 4 × 4 × 4 CSI matrix, 80 × 80 × 80 mm3 field of view, 2048 spectral points, 6 kHz bandwidth, 200 μs 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 × 2 × 2 cm3 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 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 , the amplitude of the observed FIDs will follow a 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 ‐per‐volt scaling factor that relates the (known) applied transmit voltage to the observed signal amplitudes (i.e. flip angle). 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 ‐mapping was validated against a multi‐TR magnitude reference method in a uniform phantom. The phantom was a 120 × 270 × 270 mm3 box containing 0.04 M K2PO4(aq), giving rise to a singlet signal from throughout the whole volume of the phantom. The T1 of the phantom was separately determined to be 8.57 s. The acquisition parameters were: 300 μs hard pulse excitation transmitted at 270V, 6 kHz bandwidth, 2048 spectral points, 150 × 320 × 320 mm3 field of view, acquisition weighting, 16 × 8 × 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 = 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, was fitted to the magnitude data to determine the flip angle, was calculated from the relationship where is the duration of the hard excitation pulse.

Bloch‐Siegert acquisitions used the same parameters, with the optimized Fermi pulse ( inserted as in Figure 1a. Two acquisitions were made with the pulse placed at +2 kHz and ‐2 kHz offsets. TR was 500 ms, with 13 averages at k = 0, resulting in a total duration of 2 × 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 −π and π, then was computed using Eq. (3). Voxels with (indicating either very low or the presence of uncorrected wrap artefacts in the phase difference map, i.e., true ), voxels outside the phantom, or voxels with a (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 +/−) showed highest precision and accuracy and was used for all measurements in this work.

Bloch‐Siegert ‐mapping was compared with a published dual‐TR flip‐angle mapping method 11 in a healthy volunteer's quadriceps (male, 28 years). Acquisition parameters were: 300 μs hard pulse excitation transmitted at 270 V, 200 × 200 × 200 mm3 field of view, 16 × 8 × 8 resolution (16 perpendicular to the coil), acquisition weighting, 6 kHz bandwidth and 2048 spectral samples.

In the dual‐TR comparison, excitation was centered on α‐ATP (rather than PCr), to avoid complications due to the creatine‐kinase mediated exchange of PCr and γ‐ATP 6. A literature value of = 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 ( placed at +2 kHz and ‐2 kHz relative to PCr; excitation was centered at PCr. TR was 500 ms, with 25 averages at k = 0, resulting in a total duration of 2 × 12.5 min.

Fitting was performed with AMARES in Matlab using prior knowledge specifying 10 Lorentzian peaks (α,β,γ‐ATP multiplet components, PCr, phosphodiesters [PDE], and inorganic phosphate [Pi]), with fixed amplitude ratios and scalar couplings for each multiplet. was calculated as described in “Uniform phantom validation” above. Voxels were excluded from the comparison if the voxel was centered outside the thigh, or CRLBϕ>10°.

Cardiac Scans

Cardiac 31P Bloch‐Siegert ‐mapping was demonstrated in five healthy subjects (male, 25–46 years, 75–85 kg, 1.79–1.93 m). The scan protocol followed Rodgers et al 21, with subjects positioned head‐first supine.

Two Bloch‐Siegert acquisitions were made with the optimized Fermi pulse ( placed at +2 kHz and ‐2 kHz from PCr, and excitation by a 200 μs hard pulse centered on PCr. The excitation pulse voltage was set to achieve a flip angle of 30° at the mid‐interventricular septum [the Ernst angle for a and a PCr 21]. Acquisition parameters were: 240 × 240 × 200 mm3 field of view, 16 × 8 × 8 resolution (16 perpendicular to the coil), acquisition weighting with 20 averages at k = 0, 6 kHz bandwidth and 2048 spectral samples. Total scan time was 2 × 10.5 min.

Fitting was performed with AMARES in Matlab using prior knowledge specifying 11 Lorentzian peaks (α,β,γ‐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. was calculated as described in “Uniform phantom validation” above. Voxels were excluded if or CRLBϕ>15°. Phase wrap, caused by the upper output limit of Eq. (9) was unwrapped manually.

In three subjects, the scans were repeated with the central frequency adjusted to α‐ATP to allow comparison of values obtained using different peaks.

RESULTS

Simulation Results

The precision and accuracy of the single‐peak measurement strategies (A and B′) are shown in Figure 2. Methods A and B′ show <10% error and <10% uncertainty for and . Method B′ has both precision and accuracy 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 , emphasizing the importance of effective pulse design.

Figure 2

Precision and accuracy of the single‐peak acquisition strategies [Method A (blue line) and Method B′ (red line)] computed from the analytical expressions in Table 1. A numerical Bloch simulation of both methods, with the optimized Fermi pulse, is shown for comparison (dashed lines). Note that the differences apparent in the Bloch simulations in panel (a) and (b) arise due to direct excitation and correspond to the oscillations in fractional error in Figure 4c. The magnitude of the bias introduced scales with the simulation (277Hz). The standard deviation calculated analytically is not affected by direct excitation. a: The fractional error, as defined by Eq. (4) of the simulated mean from the true (277Hz) as a function of Fermi pulse offset from the central frequency . b: The fractional error from as a function of . c,d: Standard deviation of as a function of pulse offset and .

Figure 4

a: Example Fermi pulse envelopes, with fixed duration , but varying the shaping parameters and . b: Bloch simulation of five Fermi pulses showing the transverse magnetisation immediately after the pulse. = 1000 Hz. The orange line 1 is the pulse chosen for the experimental section of this work. c: Percentage error in calculation of as a result of the direct excitation caused by the Fermi pulse. d: Percentage deviation of the transverse magnetisation from 0 as defined by Eq. (7) for a 3.5‐ms Fermi pulse as a function of the shaping parameters and . The “x” mark the locations of the pulses used to generate the lines in panels a–c.

The analysis of the multi‐peak methods (B″ and C) results in smoothly varying functions (Fig. 3). Increasing always improves the precision and accuracy. In Figure 3a (Method B″) minimizing , maximizes and improves the precision and accuracy, while adding additional peaks ( with a larger 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 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.

Figure 3

Fractional error (Eq. (4)), “Bias”, and standard deviation, “Precision”, of the multipeak acquisition strategies for a two peak spectrum. a: Method B″. b: Method C. The error and standard deviation are plotted as functions of the individual peak offsets from the Bloch‐Siegert sensitizing pulse, and ; and as the second peak position (with fixed at 2 kHz) and . The calculations do not take into account any direct excitation caused by the Bloch‐Siegert sensitizing (Fermi) pulse, therefore, the optimum separation will be a compromise between minimizing excitation and minimizing error and standard deviation. Method C (multipeak single‐acquisition) shows a improvement over Method A when the Bloch‐Siegert pulse can be placed symmetrically between the peaks. Method C has worse performance than Method B″ if the peaks are on the same side of the Bloch‐Siegert pulse, while Method B″ is independent of the sign of .

Table 1 summarizes the methods considered for our 7T cardiac 31P‐MRS application. All the measurement strategies show reasonable accuracy (<10% error) and most show an acceptable precision. Method A (dual‐acquisition +/−) shows better precision and accuracy than any other method (0.2% error and 10% standard deviation). Methods B″ and C, were simulated with only positive offsets from the Bloch‐Siegert pulse for our specific cardiac application. The greatest peak separation (PCr to β‐ATP) in 7T 31P‐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 ( ) gives a 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 . The minimum of the direct excitation metric, (Eq. (7)), fell below 0.01% for pulses with (Fig. 4d). Thus the criteria, to select the shortest pulse whilst also minimizing direct excitation, were satisfied by selecting the parameters corresponding to the minimum of the parameter range with =3.5 ms. The resulting optimized Fermi pulse had , 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 error from each.

The phosphate peak phase and magnitude in the single voxel phantom are displayed in Figure 5a. The predicted Bloch‐Siegert phase, calculated using the 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 and within 8% for . As increases from 1000 Hz the CRLBϕ is constant but the Bloch‐Siegert phase decreases, resulting in a higher relative uncertainty. When direct excitation of the phosphate peak by the Fermi pulse is observed, which causes the magnitude to drop substantially. For the magnitude varies by <1.5%, indicating that direct excitation has been minimized.

Figure 5

a: Fitted phase (red “x”) and amplitude (blue “x”) of a single peak in the presence of a Fermi pulse, as a function of Fermi‐pulse offset. The amplitude shows saturation as the narrow passband of the Fermi pulse overlaps the single phosphate peak. The phase closely follows the phase predicted by Eq. (2) using a separate measurement (solid black). The measured amplitudes are closely approximated by a Bloch simulation of the pulse sequence (magenta). The inset shows a transverse 1H image of the phantom, comprising a small cube of phosphate suspended in a tank of brine. b: calculated using the pairs of symmetric offsets from the same experiment. The values calculated by the Bloch‐Siegert method match those from a nonlocalized multi‐flip‐angle method in the range 500–4000 Hz.

Figure 5b confirms that calculated from the Bloch‐Siegert phases closely approaches that given by the Multi‐FA reference method when .

Bloch‐Siegert maps across the four central CSI slices covering the uniform phantom are shown in Figure 6a. The dynamic range of the Bloch‐Siegert was 0 Hz to 873 Hz (corresponding to a π 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 = 0.90) was observed. The Bland‐Altman plot shows the Bloch‐Siegert method measured 15 Hz lower than the reference method on average, with 95% confidence intervals of ±154 Hz. The bias and confidence intervals are equivalent to ‐3.1% and 31.8%, respectively, of the value (484 Hz) of a centrally located voxel in the phantom.

Figure 6

a: Bloch‐Siegert maps, overlaid on 1H localizer images, across four central slices of the 16 × 8 × 8 CSI grid of the uniform phantom. Each rectangle indicates the measured value in a single ideal voxel. The maps are masked for CRLBϕ, phase wrap and by the limits of the phantom. b: Scatter plot of the per‐voxel fitted multi‐TR validation method against the value measured by the Bloch‐Siegert method. The color indicates the CRLBϕ on the fitted phase difference and Pearson's correlation coefficient r = 0.90. c: Bland‐Altman plot (average of two methods versus difference) of the validation and Bloch‐Siegert methods. The color indicates the CRLBϕ on the fitted phase difference and lines show the mean difference and the ±95% confidence intervals.

Figure 7 shows a Bloch‐Siegert map across a central CSI slice covering a volunteer's quadriceps and example spectra from a centrally located voxel. The Bloch‐Siegert dynamic range was 0 Hz to 873 Hz. Figures 6e,f show the per‐voxel correlation and Bland‐Altman plots of the Bloch‐Siegert and Dual‐TR methods 11. Good correlation is observed (Pearson's r = 0.84), despite the 37.5% scan time reduction from the reference to Bloch‐Siegert method. The Bland‐Altman plot shows a ‐25 Hz lower mean from the Bloch‐Siegert method compared with the Dual‐TR method, with 95% confidence intervals of ±211 Hz. The bias and confidence intervals are equivalent to ‐6.9% and 33.4%, respectively, of the value (631 Hz) of a centrally located voxel in the skeletal muscle.

Figure 7

a: Bloch‐Siegert map from a central slice of an 16 × 8 × 8 CSI grid in a transverse plane of a healthy volunteer's quadriceps. The map is overlaid on a 1H localizer registered to the CSI grid. b,c: Example magnitude and real spectra from a central voxel 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 55.6° corresponding to a value of = 484 Hz.

Five mid‐septal cardiac maps are shown overlaid on localizers in Figure 8. The Bloch‐Siegert 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 ± 16.0 (mean ± SD) voxels were selected. After masking by CRLBϕ 30.6 ±13.7 voxels remained. The maps are smoothly varying from 100–1200 Hz without visible boundaries between skeletal muscle, myocardium and the ventricular blood pools. drops‐off with increasing distance from the coil. Data for the adjacent apical and basal slices, after masking by CRLBϕ, had 28 ± 10.1 and 21.6 ± 15.6 voxels suitable for analysis, respectively, which were used in the subsequent comparisons.

Figure 8

a: Bloch‐Siegert 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 1H 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 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 = 605 Hz.

Bloch‐Siegert measurement using α‐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 ±561 Hz. This result confirms the T1‐independence found in Sacolick et al 15 and is an indication that the values measured in this work are correct.

Figure 9

a: Three volunteer comparison of consecutive Bloch‐Siegert maps, collected with the Fermi pulse placed symmetrically around the PCr peak and then subsequently the α‐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.

DISCUSSION

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

T1 Insensitivity

Bloch‐Siegert methods are known to be T1‐insensitive 15. We verified this for our 31P‐MRS method by comparing maps derived from PCr (T 1 = 3.09s) and α‐ATP (T 1 = 1.39 s) in three normal volunteers. This is an advantage over multi‐TR and dual‐angle methods, which rely on metabolite values to calculate .

In normal volunteers, one could use cardiac 31P T1 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 T1 values may be different, as has been reported for other nuclei 27. Recording 3D‐localized myocardial T1 values at 7T took 45 min in Rodgers et al 21—too long to be a calibration step in patients.

The problem of T1‐changes in disease is particularly acute for the highest SNR peak: PCr. Chemical exchange through the creatine‐kinase shuttle dominates the apparent T1s of PCr and γ‐ATP, 6 and the forward rate constant for this process changes markedly in disease 28. Therefore, T1‐sensitive ‐mapping methods should only be applied to nonexchanging peaks such as α‐ATP or β‐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 ‐mapping sequence can be run with the Ernst flip‐angle 29 and a short TR to maximize its SNR efficiency (i.e., SNR/ . Bloch‐Siegert 1H MRI ‐mapping was shown to be more SNR‐efficient than dual‐angle ‐mapping 30. Theory predicts a similar SNR‐efficiency advantage with 31P‐MRS.

We, therefore, believe that Bloch‐Siegert ‐mapping approaches are the only known viable methods for human cardiac 31P ‐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 31P‐MRS, because there is insufficient frequency separation to place the Fermi pulse between PCr and β‐ATP without creating artefacts due to direct excitation. In a situation where the peak separation is much wider (4000 Hz) or where T2* 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 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 13C‐MRS.

In Vivo B1 Maps

We believe these are the first 31P‐MRS ‐maps made in the human heart at 7T. Therefore, there was no gold‐standard 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 (5–60 μT ) for the experimental setup; decreases smoothly with increasing distance from the surface coil, and repeated scans in three volunteers on PCr and then on α‐ATP showed excellent reproducibility. Therefore, we deem that the implementation of a Bloch‐Siegert ‐mapping sequence has been successful.

Limitations

The Fermi pulse optimization presented here is valid for the range of and SAR/ appropriate for our coil, and for T2 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 T2‐induced signal losses but which give higher SAR and need to be placed at a greater offset frequency ωRF for acceptable levels of direct excitation, and which, therefore, produce less Bloch‐Siegert phase shift for a given (and hence lead to lower precision). Or conversely, long, low‐amplitude Fermi pulses which suffer greater T2‐induced signal losses but produce lower SAR and may be placed at a smaller offset frequency ωRF for greater 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 , 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 , the measured voxel value will be a weighted average of within the voxel, not the mean voxel value. If a surface receive coil is used, then will also vary across each voxel. inhomogeneity will weight values measured by any method. To confirm that the additional weighting due to intravoxel difference in the Bloch‐Siegert phase shift is small, we simulated 5123 spin isochromats for a voxel 10 cm from the coil using the CSI protocol and RF coil described above. The Bloch‐Siegert deviated from the ‐weighted mean of the values in this voxel by <7%. The drop in SNR due to phase cancelation was for but increased to 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 measurement technique, based on Point RESolved Spectroscopy (PRESS) has been published 31. The PRESS Bloch‐Siegert method allowed highly reproducible measurements of , using the 1H 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 ‐mapping in in the presence of low SNR and multiple significant peaks.

CONCLUSIONS

We have implemented a Bloch‐Siegert spectroscopy ‐mapping technique, which is the first method capable of measuring per‐subject 31P maps in the human heart at 7T in a clinically acceptable time. We have demonstrated the optimal measurement strategy for Bloch‐Siegert spectroscopic ‐mapping and have numerically optimized a Fermi pulse for our application. Our method has been validated and shown to be accurate in phantoms and in skeletal muscle. It has been demonstrated successfully in the heart in a study on five normal volunteers.

Supporting Material S1. SuppInfo_Table1FootNote_MathermaticaNotebook.nb – Original Mathematica notebook referred to in the footnote of Table 1.

Click here for additional data file.

Supporting Material S2. SuppInfo_Table1FootNote_MathermaticaNotebook.pdf – PDF rendering of SuppInfo_Table1FootNote_MathermaticaNotebook.nb

Click here for additional data file.

Supporting Material S3. SuppInfo_Table1FootNote_WolframCDF.cdf – .cdf file conversion of SuppInfo_Table1FootNote_MathermaticaNotebook.nb. This may be viewed using the free software CDF Player available from Wolfram. The software is available here: http://www.wolfram.com/cdf‐player/

Click here for additional data file.

Supporting Material S4. FermiPulseOptimisation.m – Matlab file containing the pulse optimisation process detailed in this work.

Click here for additional data file.

Supporting Material S5. Bloch_CTR_Hz.m – Auxiliary Matlab function for S4.

Supporting Material S6. Bloch_CTR_Hz.c – C code called by FermiPulseOptimisation.m

Click here for additional data file.