Unbiased Signal Equation for Quantitative Magnetization Transfer Mapping in Balanced Steady-State Free Precession MRI

Purpose: Quantitative magnetization transfer (qMT) imaging can be used to quantify the proportion of protons in a voxel attached to macromolecules. Here, we show that the original qMT balanced steady-state free precession (bSSFP) model is biased due to over-simplistic assumptions made in its derivation. Theory and Methods: We present an improved model for qMT bSSFP, which incorporates finite radio-frequency (RF) pulse effects as well as simultaneous exchange and relaxation. Further, a correction to finite RF pulse effects for sinc-shaped excitations is derived. The new model is compared to the original one in numerical simulations of the Bloch-McConnell equations and in previously acquired in-vivo data. Results: Our numerical simulations show that the original signal equation is significantly biased in typical brain tissue structures (by 7-20 %) whereas the new signal equation outperforms the original one with minimal bias (<1%). It is further shown that the bias of the original model strongly affects the acquired qMT parameters in human brain structures, with differences in the clinically relevant parameter of pool-size-ratio of up to 31 %. Particularly high biases of the original signal equation are expected in an MS lesion within diseased brain tissue (due to a low T2/T1-ratio), demanding a more accurate model for clinical applications. Conclusion: The improved model for qMT bSSFP is recommended for accurate qMT parameter mapping in healthy and diseased brain tissue structures.


Introduction
Quantitative magnetization transfer (qMT) imaging can be used to quantify the proportion of protons in a voxel attached to macromolecules. qMT has shown considerable promise for characterising myelin-related diseases, such as multiple sclerosis. Due to a high signal-to-noise ratio and short acquisition times, balanced steady state free precession (bSSFP) acquisition modules have become a popular method for quantifying MT parameters [1,2,3]. However, the derivation of the qMT bSSFP signal equation is based on two major assumptions, which limit its generality and accuracy.
Firstly, it is assumed that magnetisation relaxation and the spin exchange between the free and macromolecular pool (MT) can be modelled as independent processes. This implies that the continuous phenomenon of MT has an instantaneous effect on the magnetisation. Though the separation of exchange and relaxation simplifies the derivation of the original qMT bSSFP signal equation, this assumption does not accurately describe the physical nature of MT, as these effects happen simultaneously.
Furthermore, the originally proposed signal equation assumes an instantaneous rotation of magnetisation by the RF pulse. Bieri and Scheffler have shown [4,5] that this assumption does not accurately describe the finite nature of an RF pulse in bSSFP due to an overestimation of transverse relaxation. While this effect is negligible for short pulse durations TRF T R ≪ 1, a significant bias is introduced if that condition is not satisfied [4]. In conventional bSSFP (non-qMT), this bias can amount to 10 % (α ∼ 90 • , T2/T1 ≪ 1) [4,5]. As qMT bSSFP is based on a stepwise variation of the RF pulse duration, this condition is certainly not met in the original qMT bSSFP acquisition scheme, where the ratio TRF T R can be as high as 0.44 [1,6]. A correction to this bias has been proposed for Gaussian pulse shapes, which are, however, not commonly used in qMT bSSFP, where a sinc pulse is more typically used [1,3,2].
Here, we present an improved signal equation for qMT bSSFP, which incorporates finite pulse effects as well as simultaneous exchange and relaxation. A correction to finite RF pulse effects for sinc-shaped excitations is derived. By means of numerical simulations of the Bloch-McConnell equations, it is demonstrated that the original signal equation is significantly biased in typical brain tissue structures. Additionally, this bias is strongly dependent on the time-bandwidth (TBW) product for sinc pulses; thus, a framework to minimise this bias is presented.

Refined Balanced SSFP Signal Equation
In this section, a new qMT bSSFP signal equation is derived allowing for simultaneous magnetization exchange and relaxation, and correcting for the instantaneous rotation by the RF pulse. To model the magnetisation dynamics a single bSSFP acquisition cycle of duration T R is considered, that is repeated until steady-state is reached. Each cycle can be split into two epochs: I. Excitation by the RF pulse II. Free relaxation (including spin information exchange between pools) To derive the magnetisation at steady-state, each epoch can be modelled independently and subsequently unified by the steady-state condition.
Analogous to the original derivation [1], the excitation by the RF pulse (epoch I.) is initially assumed to be instantaneous T RF → 0 (correction follows below). Thus the magnetisation state is instantly rotated at t = nT R, n ∈ N 0 , which is described by the following formalism M − (n) before rotation via the RF pulse M + (n) after rotation via the RF pulse This convention was established by Freeman [7] and is commonly used in bSSFP. The RF pulse leads to a rotation of the free-pool magnetisation around the x-axis and can therefore be modelled via the rotation matrix R x (α), representing a clockwise rotation in the x-plane with angle α for an anticlockwise polarized RF field. Simultaneously, the pulse saturates the macromolecular pool, which can be modelled using the mean saturation rate W (∆ → 0) . Thus, the operator, representing the action of the pulse on the magnetisation satisfying the relation The mean saturation rate W (∆) used in this derivation is equivalent to the one proposed in the work by Gloor and is described in detail elsewhere [1].
where t ′ = t − n · T R with n ∈ N 0 is the time of the n th acquisition cycle and the magnetisation is M = and R 2f are the longitudinal and transversal relaxation rates of the free pool, respectively, k mf and k f m are the exchange rates from macromolecular to free pool and free to macromolecular pool, respectively, and the pool-size-ratio F = M 0m /M 0f describes the ratio of the free pool M 0f and the macromolecular pool M 0m . Within the range n · T R < t < (n + 1) · T R, Equation 3 results in a first order linear inhomogeneous matrix ODE as the relaxation and exchange matrix ξ 1 (t ′ ) = ξ 1 is time-independent. The solution to this first order linear inhomogeneous matrix ODE is given by For repeated iterations of the pulse (n → ∞), the magnetisation reaches a dynamic steady-state, satisfying the condition where the rotation matrix R z := R z (Φ = 180 • ) represents the change in sign of the flip angle after each iteration The magnetisation during one pulse cycle (epochs I. and II.), can be modelled by combining Equations 2 and 5. This allows to relate the magnetisation before the (n + 1) th pulse to the magnetization before the n th pulse This equation can be solved under the dynamic steady-state condition (Equation 6) for n → ∞ with the resulting in the solution at steady-state The operator R x (α, t) represents an instant rotation of the magnetisation in the free pool. This assumption is commonly used in MRI, but leads to an overestimation of transverse relaxation during excitation [4,5]. Throughout the rotation caused by a pulse of finite pulse duration, the magnetisation spends a period when it has parallel alignment with the static magnetic field, i.e. its equilibrium orientation. This reduces the transverse relaxation, which can be accounted for by the correction suggested by Bieri for the one-pool where the hard pulse time equivalent T RF E is a pulse-shape-dependent constant. While this constant has previously been derived for Gaussian pulse shapes, here we present a solution for sinc pulse shapes (details in Appendix), as these are commonly used in qMT bSSFP for sinc pulse (proof in Appendix (15) Here, Si denotes the sine integral defined in Equation A.7. The correction accounts for the overestimation of transverse relaxation during excitation and therefore considers the finite pulse duration; the derivation to the correction given by Equation 13 can be found at [4].
As the finite pulse duration correction only affects the transverse magnetisation, which is generally neglected in the macromolecular pool, the two-pool model can be corrected by transforming R 2 within the As M xf is decoupled from the other components, the coupled equations can be reduced to This leads to the corrected steady-state solution for bSSFP in matrix notation, taking into account finite pulse duration effects and concurrent magnetization exchange and relaxation Note that while Equation 17 describes the magnetisation as a whole, experimentally only the transverse component of the free pool M yf is measured.

Numerical Studies
To validate the analytical solution, simulations were performed by numerically solving the Bloch-McConnell Equations for typical brain tissue parameters (Table 1).
Similar to the originally proposed acquisition, the flip angle α and the pulse duration T RF have been varied while setting all other acquisition parameters constant. As suggested in the original paper by Gloor [1], an on-resonance (∆ω = 0) sinc pulse shape has been chosen for excitation (Equation A.2).
The effect of the sinc pulse on the free pool magnetisation, given by ω 1 , has been simulated based on where the amplitude of each pulse has been calculated according to and the half-width of the central lobe t 0 is related to T RF and T BW according to Equation A.3. Due to its superior performance in tissue [8], a super-Lorentzian lineshape has been chosen for absorption according to for the n th iteration n ∈ {1, 2, ..., N } and Note that the super-Lorentzian absorption lineshape has a singularity at g m (∆ω = 0). Analogous to previous studies [9,1], the absorption lineshape has been extrapolated from 1 kHz to the asymptotic limit, resulting in g m (∆ω → 0) = 1.4 × 10 −5 s for which a constant T 2m = 12 µs has been assumed.

In-Vivo Studies
In addition to the simulations, the performance of the refined signal equation in comparison to the original one of Gloor et al. [1] was investigated in previously acquired human brain data. The images used in

Quantitative MT Parameter Analyis
In order to determine the qMT parameters in each voxel, the qMT bSSFP signal equation was fitted to the acquired steady-state magnetisation by means of a non-linear least-squares fit. Both refined and original qMT bSSFP signal equations are dependent on five parameters: F , k mf , T 1f , T 2f and T 1m . However, the additional acquisition of T 1f allows that parameter to be fixed and T 1m can be set equal to T 1f due to its insensitivity to the magnetisation [1,12,13]. Thus, the remaining parameters F , k mf and T 2f were fitted on a voxel-by-voxel basis within the ranges 0.01 ≤ F ≤ 30 %, 0.0001 ≤ k mf ≤ 100 s −1 and 0.01 ≤ T 2f ≤ 0.2 s using the starting pointsF = 10 %,k mf = 30 s −1 andT 2f = 0.04 s. M 0f has been set to one as has been done previously [14]. All computations were performed in Matlab (MathWorks, Natick, MA) and code has been partially taken from the qMRlab toolbox [10]. A non-parametric Wilcoxon signed-rank test was used for statistical testing.

Results
Numerical Studies Figure 1 shows the original (red) and refined (blue) qMT bSSFP signal equations, along with the numerically simulated data (black dots) for a standard acquisition scheme for typical brain tissue parameters ( Table 1).
The simulation has been performed for a standard acquisition scheme of varied flip angles (left) and pulse durations (right). Taking the numerical simulation as the mathematical ground truth, Figure 1

In-Vivo Studies
The resulting qMT parameter maps of a voxelwise least-squares fit on the human data are shown in Figure   2. The analysis was performed using the original and the refined signal equations. The mean values of all fitted qMT parameters within regions of interest (ROIs) in grey and white matter are listed in Table 3. Finite Pulse Duration Correction for Sinc-Shape In the case of a sinc pulse shape, the correction factor is as follows

Discussion
The simulations have shown that the original qMT bSSFP signal equation is biased by the assumptions made in its derivation (firstly separation of exchange and relaxation and secondly instantaneous rotation of the RF pulse). This bias has been seen to be tissue dependant, amounting to deviations of up to 7 % and 11 % in white and grey matter of healthy brain tissue and exceeding 20 % in an MS lesion. The tissue dependence is expected, as the bias linearly depends on the relaxation time ratio T2 T1 (Equation 13), which varies amongst different tissue types. Further, the bias has been shown to increase at higher pulse durations T RF . This reflects the fact that while the assumption of an instantaneous rotation by the RF pulse might be sufficient for short pulse durations, it is increasingly violated at longer T RF . In qMT bSSFP, this bias is particularly strong, as the acquisition involves long T RF relative to T R [1]. The bias is passed on to the qMT parameters, as they are determined by fitting the signal equation to the acquired data.
To address this, the suggested refined signal equation for qMT bSSFP (Equation 17) has been derived, accounting for the assumptions made in the original model, and describes the simulated data with minimal bias (<1 %).
The comparison of original and refined signal equations in-vivo shows significant differences between the resulting qMT parameters (24-31 % for pool-size-ratio, 0-21 % for exchange rate and 9-13 % for transversal relaxation time). In agreement with the simulation results, the difference between qMT parameters, determined by the original and refined model, is significantly greater in white matter compared to grey matter in in-vivo brain tissue data. This is in agreement with the theoretically predicted T2 T1 dependency of the bias (Equation 13).
In previous studies of different qMT modalities, pool-size-ratios in the range of 10-16 % and 3-8 % have been reported in white and grey matter structures, respectively [15,16,17,18,19]. The pool-sizeratios, determined by the original model in this work, exceed the previously reported range in white matter (18.4 ± 1.2 %, 19.3 ± 0.8 %) and approach the upper limit in grey matter (8.6 ± 1.2 %, 7.9 ± 0.6 %). In contrast, the refined model estimates of the pool-size-ratio are in good agreement with the findings in other studies in both white matter structures (12.7 ± 0.8 %, 14.1 ± 0.5 %) and grey matter structures (6.5 ± 1.1 %, 5.9 ± 0.7 %). This indicates that the refined model outperforms the original one not only in simulation, but also in-vivo. This conclusion is further supported by the significantly lower residual sum of squares (RSS) found in the fits of the refined model compared to the original one. The wide range of previously reported exchange rates in different qMT methodologies 10-40 s −1 includes the results in both original and refined signal equations in this work.
The pool-size-ratios determined by qMT bSSFP in [1] are 13-16 % and 6-7 % for white and grey matter structures, respectively. Although they fall at the upper end of previous findings, they are lower than the values found with the original model in this work. The reason for reduced biases in the original findings [1] can be explained by means of the pulse shape analysis, established in Section . While the parameters of the original findings have been acquired with a T BW = 2.7, in this work a T BW = 2 has been used.
Their respective hard pulse equivalent durations, which correlate with the bias ∆R 2 ∝ T RF E , differ by TRF E (T BW =2.7) = 0.58 . This implies a reduction of the bias in the original acquisition scheme for a T BW = 2.7 and explains why bias is reduced in the original publication [1]. While the bias is only reduced and not removed, much higher biases are expected for a T BW ≤ 2.5. Alternatively, the refined signal equation allows for a general solution with accurate parameter estimation for a wide range of T BW .
Additionally, the derived Equation 15 predicts that the bias oscillates for varying T BW and even approaches zero for a T BW ∈ {4, 8, 12, ...}. The physical explanation to the oscillation lies in the sinc-shape specific side lobes. These side lobes cause a temporary increased deflection of magnetisation from the equilibrium alignment, for which transverse relaxation is underestimated. The underestimation, induced by the negative side lobes, counterbalances the overestimation, resulting from the main lobe. Therefore, the bias in quantitative bSSFP methods (qMT and non-qMT) can be removed by choosing an appropriate time-bandwidth product without using the correction given by Equation 13. This might be be useful for applications, where the correction is inaccurate due to strong magnetic field inhomogeneities, such as is the case at high magnetic field strengths.
Recent work by Wood et al. [20] has demonstrated that the PLANET method [21,22,23] for phase-cycled bSSFP can be applied to qMT at higher field strengths to derive qMT parameter estimates free from banding artefacts. However, Wood et al. [20] utilised the signal model from Gloor et al. [1], which translated into increased errors in their parameter estimation, particularly in white matter. We hypothesize that combining the method from Wood et al. [20] with our methods here would provide increased accuracy, leading to a method which can produce qMT parameter estimates quickly over all clinical field strengths. However, this is beyond the scope of this paper, and is left for future work. T 1f -ratio), demanding a more accurate model for clinical applications. Therefore, the refined signal equation is recommended for accurate qMT parameter estimation in healthy and diseased brain tissue, especially when using a T BW ≤ 2.5.

Acknowledgements
The Wellcome Centre for Integrative Neuroimaging is supported by core funding from the Wellcome Trust (203139/Z/16/Z). We also thank the Dunhill Medical Trust and the NIHR Oxford Biomedical Research Centre for support (PJ). AKS acknowledges support from the Whitaker International Program and St.
Hilda's College at the University of Oxford. FMB acknowledges support from the Cusanuswerk Scholarship.

Appendix: Finite Pulse Duration Correction for Sinc-Shape
As shown by Bieri [4,5], the overestimation of transverse relaxation in bSSFP can be corrected by means of the substitution of R 2 →R 2 according to Equation 13. This correction has been derived under the assumption of excitation by a constant RF magnetic field (i.e. a hard pulse). By means of the hard pulse equivalent duration T RF E , the correction can be transferred to different pulse shapes where < B > is the mean B 1 amplitude. This relation has previously been solved for a Gaussian pulse shape, resulting in T RF E = 1.2 TRF T R [4]. In this section, Ansatz A.1 is solved for a sinc pulse shape.
Consider a sinc pulse of form where t 0 is the half-width of the central lobe, A is the amplitude of the pulse and t is the time. N L and N R are the numbers of zero-crossings of the sinc pulse to the left and right of the central peak, respectively (if symmetric: N L = N R ). The full width at half maximum (FWHM) of a sinc pulse can be approximated by FWHM = ∆f ≈ 1 t0 [24] and the time bandwidth product is defined as with the number of total zero-crossings N .
As shown by Bieri [4], Ansatz A.1 leads to a relation between T RF E and magnetisation M xy where M − xy is the magnetisation before excitation and M xy + is the time average magnetisation during excitation. By calculating the pulse shape dependent M xy + , Equation A.4 can be used to find the hard pulse equivalent duration T RF E . For sufficiently small flip angles, the magnetisation can be approximated where the the relations of flip angle α(t) = t −TRF /2 ω 1 (t ′ )dt ′ and ω 1 (t) = γB 1 (t), have been used.
Using the symmetry property of the sinc function sinc(x) = sinc(−x), the flip angle can be expressed as where the substitution πt t0 = θ has been used. The integral can be solved using the sine integral definition leading to a relation between flip angle, pulse duration and half-bandwidth Calculating the time average magnetisation during excitation, and exploiting the symmetry around t = 0, results in where the definition of flip angle α and the symmetric nature of the sinc function have been used. This simplifies Equation A.9 to Further substitutions and the use of the integral in Equation A.7 leads to The integral of the sine integral Si t is found by means of partial integration          Occipital GM 7.9 ± 0.6 5.9 ± 0.7 24.0 ± 3.7 25.2 ± 8.6 32 ± 6 29 ± 5