A comparison between equations describing in vivo MT: The effects of noise and sequence parameters

Quantitative models of magnetization transfer (MT) allow the estimation of physical properties of tissue which are thought to reflect myelination, and are therefore likely to be useful for clinical application. Although a model describing a two-pool system under continuous wave-saturation has been available for two decades, generalizing such a model to pulsed MT, and therefore to in vivo applications, is not straightforward, and only recently have a range of equations predicting the outcome of pulsed MT experiments been proposed. These solutions of the 2-pool model are based on differing assumptions and involve differing degrees of complexity, so their individual advantages and limitations are not always obvious. This paper is concerned with the comparison of three differing signal equations. After reviewing the theory behind each of them, their accuracy and precision is investigated using numerical simulations under variable experimental conditions such as degree of T1-weighting of the acquisition sequence and SNR, and the consistency of numerical results is tested using in vivo data. We show that while in conditions of minimal T1-weighting, high SNR, and large duty cycle the solutions of the three equations are consistent, they have a different tolerance to deviations from the basic assumptions behind their development, which should be taken into account when designing a quantitative MT protocol.


Introduction
The Magnetization Transfer (MT) effect is based on the exchange of magnetization occurring between groups of spins characterized by different molecular environments. In biological tissues, two or more ''pools'' of protons can be identified: those in free water (the free, or liquid, pool) and those bound to large molecules (referred to as restricted, semisolid, or macromolecular, pool). The latter protons are characterized by a very short transverse relaxation time (T 2 ) and therefore do not directly contribute to signal intensity in conventional magnetic resonance (MR) images. Nevertheless, it is possible to sensitise an MR experiment to the magnetic resonance characteristics of macromolecular protons by exposing the sample to radiofrequency (RF) energy several kilohertz away from the Larmor frequency. Protons in free water are relatively insensitive to such irradiation, but it can cause saturation of protons in the semisolid pool which, due to their short T 2 and correspondingly large line width, are responsive to irradiation at these frequencies. In these conditions, any exchange of magnetization between pools results in a decreased intensity of the observed MR signal.
From a quantitative MT model based on the exchange between two pools Henkelman et al. [1] derived a signal equation for the continuous wave (CW) case, in which RF irradiation of particular (constant) amplitude and several seconds duration is used to saturate the macromolecular pool. The parameters characterizing the two pools in the model are potentially interesting to measure, and they can be estimated by fitting Henkelman's equation to a set of MR measurements obtained in the presence of MT pulses with a suitable set of amplitudes x 1 and offset frequencies Df.
As CW irradiation is impractical and generally not available for in vivo imaging experiments, in vivo MTweighted MRI is generally obtained using the so-called pulsed MT acquisition, in which the long period of saturation is replaced by a much shorter irradiation pulse (typically applied just before each excitation pulse) along with intervals without irradiation (during which data is collected). For data from this type of acquisition, Henkelman's equation must be modified to allow for the short duration of the saturation pulses relative to T 1 [2]. A number of such modified signal equations for pulsed MT, all based on the same original two-pool model, have been developed [3][4][5]. While the numerical results published so far suggest reasonable consistency across the solutions predicted by these equations, no direct comparison is available, and the differing conventions and symbols used mean that evaluating discrepancies and similarities between them is not straightforward.
This paper is concerned with the comparison of three of these signal equations-two derived by Sled and Pike [3,6], plus that of Ramani et al. [5]. Ramani et al. used a CW power equivalent approximation (CWPE) [5] where the pulse is simply replaced by a CW irradiation with the mean square amplitude that would give the same power over the interval between MT pulses. By means of the CWPE approximation, Henkelman's steady state model can be straightforwardly applied to the in vivo MRI case, neglecting the imaging elements of the pulse sequence. Due its steady state nature, the implicit assumption within the equation is that the relative signal intensity in data obtained with different MT-weightings only depends on the characteristics of the MT pulse, and that T 1 and T 2 relaxations equally affect all measurements.
As Ramani's equation does not explicitly model the effects of the excitation pulses and TR, its description of the MT-weighted signal is valid only when the degree of T 1 -weighting in the acquisition sequence is minimal. As it effectively assumes that the MT pulse is applied continuously, another parameter likely to affect the accuracy of Ramani's equation is the duty cycle, i.e. the duration of the MT pulse relative to the repetition period, whose effect has never formally been investigated. Sled and Pike [6] propose an alternative equation which can be fitted directly to the measured signal. Their solution is derived by approximating the pulse sequence as a series of periods of free precession, CW irradiation and instantaneous saturation of the free pool. It has the advantage of incorporating the effect of the excitation RF pulses, and also makes it possible to account for saturation effects of the excitation. Together with this solution, the authors propose also a simpler variant which neglects free precession, thus assuming a succession of instantaneous saturation of the free pool and CW irradiation of the macromolecular pool for the total duration of the interval between pulses. Both equations presented by Sled and Pike for in vivo applications require the numerical evaluation of ordinary differential equations (when modeling the rate of saturation of the macromolecular pool with a super-Lorentzian, see next section), at least for the estimation of the effect of the MT pulse on the free pool, and they are therefore computationally intensive. Ramani's solution has the advantage of being simpler, at the price of its inability to account for the effects of the excitation pulses. A recent paper presented an evaluation of these signal equation, validated using animal data [7]. Here we first review the theory behind them and then use numerical simulations to extend the range of experimental conditions under which each can be tested (investigating how duty cycle, saturation effects of the excitation and noise affect the MT parameters fitted by each of them). We also perform a statistical comparison between MT parameters estimated using each signal equation in healthy brain tissue from in vivo data.

Coupled Bloch equations
Assuming that the MT effect can be modeled using a liquid pool (A) and a macromolecular pool (B), the magnetization of either pool can be described by its longitudinal component ðM A z ; M B z Þ and its transverse components The exchange between pools associated with the transverse components of magnetization can be considered negligible due to the extremely short T 2 associated with the macromolecular pool [2,6]. The coupled Bloch equations for the system can thus be written as follows: In Eqs. (1)-(4), T A 2 represents the transverse relaxation time of the liquid pool, M A 0 and M B 0 are the fully relaxed values of magnetization associated with the two pools (assumed dimensionless), R A and R B are their longitudinal relaxation rates, and R is the exchange rate constant. Df represents the frequency offset of the pulse, while x 1 (t) is the time dependent amplitude of the pulse expressed in rad s À1 (i.e. the angular frequency of precession induced by the pulse). R RFB (Df, x 1 (t)) is the rate of saturation of longitudinal magnetization in pool B due to the irradiation by the amplitude defined by Df and x 1 (t), and depends on the transverse relaxation time of the macromolecular pool, Li et al. [8] show that, in brain tissue, the specra associated with macromolecular pool are better modeled by a super-Lorentzian, (with x 2 1 being the average power of the MT pulse), yielding: and we adopt this model here.
We remark here that this notation is not universal. Some authors label the A and B pools as 'F' and 'R', respectively [6], or 'f' and 'm' [9,10] and use the symbol W instead of R RFB [6]. The pseudo first-order exchange rates, RM B 0 ðA ! BÞ and RM A 0 ðB ! AÞ, are often referred to as k f (or simply k) and k r .

Sled and Pike's RP signal equation
Assuming that the pulse sequence consists of an MT pulse followed by an excitation pulse and by a period of recovery, Sled and Pike [6] decompose it into a series of periods where Eqs. (1)-(4) have exact or approximate solutions. These solutions can then be concatenated by imposing the appropriate initial conditions, leading to an expression for the measured signal which is less expensive to compute than numerically integrating the full set of differential equations. The effect of an MT pulse on the macromolecular pool is modeled as a rectangular pulse whose width is equal to the full width at half maximum (s RP ) of the curve obtained by squaring the instantaneous amplitude of the MT pulse throughout its duration, and whose amplitude is such that the pulses have equivalent average power (rectangular pulse, or RP, approximation). The effect of the pulse on the liquid pool is modeled as an instantaneous fractional saturation of the longitudinal magnetization. Such fractional saturation (S 1A ) is estimated by solving (numerically) the system of Eqs. (1), (3) and (4) when R and R A are set to 0.
In matrix form [6], considering the longitudinal components of magnetization only Instantaneous saturation of the free pool, caused by both MT and excitation pulses, is simply described by multiplying M z by the matrix S (where h is the excitation flip angle) The state of the magnetization after a period t 1 (assuming starting time = t 0 ) is given by the solution to the system of Eqs. (1) and (2) for either free precession [FP] or CW: with According to Sled and Pike's RP approximation, over the time interval T between application of MT pulses (typically the time required to excite and collect data for a single k-space line of a single image slice), M z undergoes instantaneous saturation, CW irradiation for a period s RP /2, FP for a period (T À s RP ), and CW for another s RP /2. After including all thee steps, we can impose the equality and solve for M z yielding an equation for the longitudinal components of magnetization. Recalling that the signal observed at readout is (where the repetition time TR may be PT depending of the details of the image acquisition) it is thus possible to model the MT-weighed signal. The solution to Eq. (10) can be easily computed in matrix form, and we adopt this procedure for all the following experiments.

Sled and Pike's CW signal equation
A simpler expression is presented in the same paper [6], where the effect of the MT pulse on the macromolecular pool is described by a CW irradiation of duration T. In this case, over the same period T, M z undergoes instantaneous saturation, and CW irradiation (of the restricted pool only) for a period T. Following the same procedure as described above, a more manageable analytic expression can be derived in this case: Here M A z;CW is the steady state solution obtained by Henkelman et al. [1] for CW irradiation of the restricted pool of duration T: which is equivalent to the first element of the vector M z;CW ¼ A À1 CW BM 0 . In Eq. (12), We note here that for the specific case of the 3D spoilt gradient echo acquisition described in this paper (see Section 3) T = TR.
Sled and Pike also introduce, as a useful index which is believed to correlate with myelin content, the relative size of the macromolecular pool [6] F, defined as F can be fitted directly, by substituting and previous. More details on these signal equations can be found in [3,6,11].

Ramani's signal equation
Henkelman's solution [1] for the CW case is obtained by solving Eqs. (1)-(4) in the steady state, i.e. setting the derivatives on the left hand side to zero. Ramani et al. [5] adopt the same approach in the pulsed MT case, simply replacing the MT pulse with a CW irradiation with the same mean square amplitude where P SAT is the mean square saturating field. In order to ease the comparison between the three signal equations, we break with the terminology of the original paper [5], where the macromolecular fraction f (with f = F/(F + 1)) was used, and instead rewrite the Henkelman-Ramani expression using F and RM B 0 to obtain where M 0 is the signal with no MT-weighting (again assuming a constant of proportionality of c = 1). Note that while in the equations presented in the previous section, M 0 is the equilibrium magnetization, in Eq. (17) it simply represents the partially recovered magnetization available prior to the application of an MT pulse.

Fitting
All the equations are written in terms of seven parameters: M 0 ; R A ; R B ; RM B 0 ; F ; T A 2 and T B 2 , but these cannot be uniquely determined [1,12]. Constraints are imposed by measuring the observed longitudinal relaxation rate of the sample, R Aobs (=1/T 1obs ) independently, linked to R A by Henkelman et al. [1] A further issue is that dependence of S(x 1 , Df) on R B is weak, making fitting of this parameter unstable. Since the estimates of the other parameters are largely insensitive to its value, R B is usually kept fixed at 1 s À1 [1,5,6]. This reduces the number of free parameters to 5, which can be estimated by fitting the equation to five or more measurements with different combinations of x 1 (t) and Df.

Numerical simulations
In order to compare the three signal equations, and to highlight their shortcomings, we need to test their performance against data corresponding to a known set of parameters. The easiest way to obtain such data is to synthetically produce them, using numerical simulations.
We consider here the case of an MT-weighted spoiled gradient echo acquisition, where off-resonance saturation is achieved using Gaussian pulses (of duration s SAT ) applied once every TR (just prior to RF excitation), while on-resonance excitation is obtained using short 5lobe sinc pulses (in the presence of a 'slab selection' gradient).
Eqs. (1)-(4) can be solved numerically to predict the longitudinal magnetization at the end of the MT pulse (i.e. just before the excitation pulse), M A z ðs SAT Þ. Since the measured signal intensity is proportional to the transverse magnetization at readout (if we neglect T Ã 2 decay): The MT pulse is characterized by its maximum amplitude, B 1SAT, MAX , by its duration, s SAT , and by the standard deviation of the Gaussian envelope, r. The excitation pulse is characterized by its maximum amplitude, B 1EXC, MAX , by its duration, s EXC , and by its bandwidth, BW. Both pulses can be described by their equivalent on-resonance flip angle, given by the integral over pulse duration of x 1 (t).
Additionally, we assume the excitation pulse to have no effect on pool B. The measured signal intensity can be estimated by calculating the solution of Eq. (19) as M z tends towards a steady state, i.e. after solving Eqs. (1)-(4) for several TRs, until the difference between M z (TR n + s SAT ) and M z (TR n+1 + s SAT ) is less than 0.01% of M z (TR n + s SAT ). We simulate four experiments to probe the effects of different parameters on the fits. Firstly we aim to identify any systematic biases in the three signal equations with respect to variable experimental conditions, by fitting them to noise-free simulated data (Experiments 1-3); next we investigate their robustness in the presence of noise, by using a Monte Carlo approach (Experiment 4). The simulated signal is computed by using a Runge-Kutta ordinary differential equation (ODE) integrator with adaptive step-size control [13]. Spoiling of the transverse magnetization is simulated by setting the transverse components of magnetization equal to zero before the occurrence of every MT pulse. All experiments use the same set of MT parameters ðR A ; R B ; T A 2 ; T B 2 ; F ; RM B 0 Þ, chosen to be similar to values measured previously in white matter [14,15] and shown in Table 1 as 'test set'. All also use the same excitation pulse parameters (s EXC = 3.2 ms, and BW = 2.5 kHz), but other values differ: Experiment 1: In the first simulation, we investigate the accuracy of Ramani's, Sled and Pike's CW and Sled and Pike's RP signal equations as a function of duty cycle, by simulating 5 MT experiments using an excitation flip angle of 5°, TR = 30 ms, and varying s SAT between 5 and 25 ms, in steps of 5 ms. The full width at half maximum (FWHM) of the pulses varies accordingly between 2.29 (r = 0.97) and 11.45 ms (r = 4.86 ms), in steps of 2.29 ms. Each simulated set consists of 60 points, generated using only two fixed values of x 1CWPE (250.2 and 850.7 rad s À1 ) (and thus different MT flip angles for each value of s SAT ), following Sled and Pike [6,15], and 30 values of Df per flip angle. The offset frequency ranges from 400 to 30,000 Hz, sampled at regular interval on a logarithmic scale. The three equations are fitted to the synthetic datasets using the Levenberg-Marquardt method, as implemented in Numerical Recipes [13], to yield the estimated parameters. Numerical derivatives are computed where an analytical expression is unavailable. The test set of parameters (see Table 1) are used to provide the initial parameter estimates, and R Aobs is obtained by solving Eq. (18) with respect to this quantity. The same ODE integrator used to compute the simulations is used to estimate S 1A at every step for Eqs. (7), (11) and (12).
Experiment 2: Secondly, we explore the effects of saturation of the excitation on the estimated parameters, simulating the outcome of six MT experiments using regularly spaced excitation flip angles ranging from 5°to 20°and TR = 30 ms. We fix s SAT = 15 ms and keep the other parameters as in the first experiment. With the exception of the excitation flip angle, the input parameters are identical for all six cases.
In order to check that our results are not specific to the choice of the sampling scheme, we repeat Experiments 1 and 2 with an alternative sampling scheme generated using four fixed values of x 1CWPE (250.2, 450.4, 650.5 and 850.7 rad s À1 ) and 15 values of Df per flip angle.
Experiment 3: As Sled and Pike used two sequences with two different MT pulse durations and two TRs to constrain the estimate of RM 0B , we create a synthetic dataset formed by two ''sequences'' (40 MT points with TR = 30 ms, s SAT = 15 ms, h = 5°, 2 MT flip angles equal to 250°and 850°, respectively; and 20 with TR = 45 ms, s SAT = 20 ms, h = 6°, 2 MT flip angles equal to 353°and 1202°, respectively). We then compare the accuracy of the MT parameters estimated fitting the three signal equations to such a dataset and to the ''single sequence'' dataset obtained in Experiments 2 when fixing s SAT = 15 ms and h = 5°.
Experiment 4: Finally, in order to investigate the sensitivity to noise, we add complex noise with zero mean Gaussian real and imaginary parts to the dataset obtained in Experiment 2 for s SAT = 15 ms and h = 5°. We then take the modulus to obtain a noisy data sample. The standard deviation of the Gaussian noise is set to be M 0 /R, where R is the desired SNR in the unweighted image, which we vary over the interval [20,300]. (The SNR values typically observed in 3D spoiled gradient echo scans from our system, with acquisition parameters similar to those detailed below (see Section 3.2), typically range between 40 and 100 depending on the coil used, resolution, use of parallel Table 1 MT parameter values used to create synthetic data (TEST SET) and results of noise-free simulation with excitation flip angle = 5°, s SAT = 15 ms (duty cycle = 50%) using all signal equations The results obtained when using a single set of 60 points (same TR and s SAT ) are also compared obtained when using 2 sets (40 points with the same TR and s SAT as the single set, and 20 points with longer TR and s SAT ).
imaging, etc.) For each level of noise, we generate 10,000 sets of noisy independent samples and fit the three equations to each set. Look-up tables are used for the super-Lorentzian lineshape, the fractional saturation S 1A , and their derivatives, in order to speed up the computation.

In vivo data
A single subject (female, 34 years old) was scanned twice on a 1.5 T system (SIGNA Horizon Echospeed, General Electrics, Milwaukee, WI, USA) using a 3D MT-weighted fast spoiled gradient recalled-echo (SPGR) sequence [15] (TR/TE = 28/5.1 ms, Gaussian MT pulses, duration = 14.6 ms, standard deviation = 2.84 ms, bandwidth = 125 Hz, matrix = 256 · 96 · 32, FOV = 240 · 180 · 160 mm 3 , to reconstruct twenty-eight 256 by 256 voxel slices). The excitation flip angle was 5°on the first session and 15°on the second one, the interval between scans was 7 days. On each occasion a dataset of 20 MT points was obtained, using 2 MT pulse flip angles (220°a nd 820°, corresponding to a CWPE amplitude of 251.1 and 861.2 rad s À1 ), and 10 values of Df per flip angle. Df ranged between 400 and 20,000 Hz, and was stepped using a constant logarithmic interval. In addition to the MT data, three 3D SPGRs (TR = 13.1 ms, TE = 4.2 ms, same FOV and resolution as the MT sequence) were also obtained on each occasion, with three different excitation flip angles (h = 25°, 15°, 5°) in order to independently estimate the longitudinal relaxation rate of the system, R Aobs . The body coil was used for signal transmission and the manufacturer's eight-channel head coil was used for reception. The total scan time was about 45 min.
The study was approved by the Joint Research Ethics Committee of The National Hospital for Neurology and Neurosurgery and the Institute of Neurology, UCL, and the subject gave written informed consent before taking part.

Image analysis
The two datasets (one from each MRI session) were processed on a Unix workstation (Sun Microsystems, Mountain View, CA, USA), as described elsewhere [15]. Briefly: the 20 MT-weighted volumes obtained with the MTweighted SPGR sequence and the three volumes obtained with the non-MT-weighted SPGR sequence were co-registered to the first MT-weighted volume using a modified [16] version of Automated Image Registration (AIR, available at http://air.bmap.ucla.edu:16080/AIR) [17]. R Aobs was estimated on a pixel-by-pixel basis by fitting the theoretical SPGR signal equation to the signal in the non-MTweighted SPGR images, as a function of the flip angle [18]. The 3 MT signal equations were fitted to the remainder of the images (as described for the synthetic data) yielding estimates of the MT parameters. Six bilateral regions of interest (ROIs), three located in white matter (frontal, temporal and internal capsule) and three in gray matter (thal-amus, putamen, caudate nucleus) were outlined on the T 1weighted images obtained from the non-MT-weighted SPGR scan (flip angle = 25°). The 12 ROI outlines were then superimposed on the MT parameter maps, yielding 6 (3 signal equations times 2 flip angles) estimates for each of the following: RM B 0 ; F ; T B 2 , and T A 2 . Paired sample T-tests were used to compare the mean estimated parameters between equation solutions, and two-sample T-tests were used to test inequalities between flip angles, considering statistically significant two-tailed p values lower than 0.01.

Results
To provide a qualitative description of the accuracy of the three signal equations, we show in Fig. 1 the MT spectra simulated using Eq. (18) and the test set in Table 1, together with those simulated using each of the three equations and the same test set, for two of the cases explored in Experiment 1 (A: duty cycle = 50%, flip angle = 5°; B: duty cycle = 17%, flip angle = 5°). All curves are normalized to the maximum intensity. As expected, the largest deviations between signal equations are observed at high power, for small offset frequencies, with Ramani's equation providing the least accurate description. The deviation between Ramani's predictions and the others are more pronounced for lower duty cycles (Fig. 1B).

Duty cycle effect
The estimates of R A ; RM B 0 ; F , and T A 2 against duty cycle are shown in Fig. 2

Saturation effect of the excitation
The MT parameter estimates obtained from noise-free simulated data at various flip angles are reported in Fig. 3    We report in Table 1 the estimated MT parameters obtained by fitting the three equations to the simulated noise-free 60 point dataset with excitation flip angle of 5°( minimum saturation effect from the excitation pulse) and duty cycle %50% (at which point duty cycle effects have reached a plateau, and the three signal equations seem to provide consistent estimates). All equations give similar results (very close to the 'test set' of parameters used in creating the simulation (see Table 1)). The two variants of Sled and Pike's solution provide more accurate estimates than Ramani's equation for all parameters, with the CW approximation yielding a more accurate value for RM B 0 , and RP approximation yielding a slightly more accurate value for F (which Ramani's equation tends to underestimate). For all the signal equations, the largest error is in the estimation of T A 2 (7.6% of true value for Sled and Pike's RP equation, 10.2% of true value for Sled and Pike's CW equation, 21.2% of true value for Ramani's equation). The results of Experiment 3 (fitting the signal equations to a combined dataset obtained with differing ''sequences'', i.e. different TR, s SAT , and flip angle) are also shown in Table 1. The estimates of T A 2 appear highly sensitive to the combination of acquisition parameters, while all other parameters are not. Interestingly, T A 2 is underestimated using both of Sled and Pike's equations when using 2 sequences (and the error becomes larger), while it is generally overestimated when using a single one.
When using the four-power scheme the main difference compared to the two-power scheme was in the estimates of RM B 0 , which were slightly lower for the two variants of Sled and Pike's solution (between 2.98 and 3.1 for the CW and between 3.3 and 3.4 for the RP approximation). Fig. 4 compares the estimates of the MT parameters using the three signal equations in the presence of noise. Overall, the estimates obtained from the 60 point noisy dataset using any of the three signal equations are characterized by similar precision, although Ramani's solution seems to provide slightly more robust results than Sled and Pike's at SNR lower than 120 for RM B 0 and T A 2 . The two equations proposed by Sled and Pike's provide very similar estimates (with similar precision at all SNR levels) for T A 2 and T B 2 . The CW variant provides more accurate values for RM B 0 . The parameter whose estimate deviates most from its true value is (for all signal equations) T A 2 , with Ramani's estimates deviating more than the others. The standard deviation associated with this parameter, however, is smaller when using Ramani's equation, at low SNR. The opposite is true (with the two variants of Sled and Pike's formulation giving very similar performances) at high SNR.

In vivo results
The MT parametric maps obtained by fitting all three equations to the 5°dataset obtained in vivo are character-ized by similar quality (examples of F maps are shown in Fig. 5). Fitting the two equations proposed by Sled and Pike to the 15°dataset was more problematic; in some voxels (most commonly in gray matter, and at the boundary between tissues) the equations gave physically meaningless parameter estimates, particularly for RM B 0 and T A 2 . It is possible that this is the result of T 1 -contrast between tissues and CSF increasing the sensitivity of partial volume effect to motion. Fitting of Eq. (17) appeared to be more robust. Table 2 shows the values (mean and standard deviation) obtained for each white and gray matter ROI using each signal equation and both flip angles.

5°data (between equations)
Although the estimates of F ; RM B 0 and T A 2 obtained with Ramani's equation were statistically different (p < 0.001, providing higher estimates T A 2 and lower estimates of F and RM B 0 ) than those obtained with either variant of Sled and Pike's equations, for all parameters the absolute difference between Ramani's and Sled and Pike's estimates was always lower or comparable to the between-voxels standard deviation (within each ROI) (see Table 2). When comparing the two solutions proposed by Sled and Pike, none of the variables differed (p values of 0.02 for RM B 0 , 0.6 for F, 0.32 for T B 2 and 0.47 for T A 2 ).

15°data (between equations)
Conversely, when comparing the parameters obtained by fitting different equations to the 15°dataset, the mean values were all significantly different (p < 0.001), with the exception of T A 2 obtained using the two variants of Sled and Pike's solution (p = 0.09). The largest differences were in the estimates of F and T A 2 obtained using Ramani's equation with respect to both Sled and Pike's equations. The absolute difference between estimates, in this case, were approximately one order of magnitude larger than the between-voxel standard deviation (within each ROI) of the same parameters.

15°data vs. 5°data
The estimates of T B 2 obtained fitting any of the three equations to the 5°dataset were significantly different from those obtained by fitting the same equation to the 15°dataset (p < 0.001). For Ramani's equation, the estimates of F and T A 2 obtained from the 5°dataset were also significantly different from those obtained from the 15°dataset (p < 0.001).

Discussion
We have shown that the estimates of RM B 0 ; F , and T B 2 obtained by fitting different equations describing the behavior of the two-pool model under conditions of pulsed MT to proton density-weighted noise-free simulated data with an MT duty cycle of approximately 50% (typical of most in vivo applications) are in good agreement and deviate only slightly from the values used to create the simulations. There are some systematic differences, however, with Sled and Pike's equations typically providing more accurate results. Overall, these findings are in line with those reported by Portnoy and Stanisz [7]. T A 2 is substantially overestimated when using all signal equations, albeit less so when using the Sled and Pike's RP approximation. The discrepancy between the estimate of this parameter and the observed T 2 of the system has been observed before [11]. As white matter is known to have multiple free water T 2 components, one of the explanations provided by researchers is that the T 2 measured by spin echo experiments and T A 2 measured by MT experiments represent different weighted averages of multiple water components. This explanation, however, does not apply to analogous results obtained in gels [1,11]. Similarly, the results of our simulations suggest that such a difference could merely result from the inadequacy of the model to estimate this parameter. It is unsurprising that, for both simulated and in vivo data, the largest quantitative difference between the results obtained when fitting different equations to the same (low flip angle) data is in the estimation of this parameter. When the effects of saturation from the excitation pulse are minimized, in fact, the main difference between Ramani's and Sled and Pike's equations is in the way the MT pulse effect on the liquid pool is modeled. It is interesting to notice that others have often reported this parameter to be underestimated [11], and our own in vivo results are consistent with such observations (Table 2). However, when fitting the MT signal equations to simulated data, we observed the opposite trend. A possible explanation for this inconsistency is the sensitivity of Levenberg-Marquardt fitting to the values used as starting points: while for synthetic data we can use the real (test set) numbers, when fitting in vivo data we can only provide a 'best guess'. In order to compensate for this effect we repeated Experiments 1 and 2 running the fitting procedure 10 times, each time perturbing the initial guesses (independently for each parameter) by a random factor from a Gaussian distribution with zero mean and standard deviation equal to 10% of the test parameter [19]. We retained and compared the set of parameters which gave the best fit out of the 10 trials in each case. This did not affect the estimates obtained using Sled and Pike's equations. T A 2 obtained using Ramani's equation was even larger in this case, with all other parameters almost unaffected (data not shown).
An alternative explanation is that we used a larger number of points for simulations than for the human brain data. The range and spacing of the sampling points is likely to affect the results of the fitting [20], and it is interesting to observe that even in the simulated data, T A 2 is underestimated when using two ''sequences'' (as recommended by Sled and Pike [6,11]). Contrary to their observations, however, the estimates of RM B 0 were virtually unaffected by this experimental parameter, and therefore our simulation experiments do not support the need for this type of acquisition. It is not clear why our result in this case should differ from that of Sled and Pike, but it is possible that this is due to our use of a different fitting routine. Because of our simulation results, we chose not to adopt the two-sequence approach for our in vivo comparison, as this would have required all the MT-weighted acquisitions to be normalized to the same maximum value, and also might introduce variable degrees of T 1 -weightings, which would further violate the assumptions underlying Ramani's equation.
Although the reason why T A 2 is so strongly affected by the use of a dual TR protocol is unclear, it has been reported by others that even when using more accurate formulations, the estimates of R and T A 2 are highly sensitive to the choice of data points [7], and this has been explained as a result of the poor sensitivity of the equations to these parameters, or to systematic errors between the model and the data, which can vary with the sampling points. The poor sensitivity of the two-pool model to changes in T A 2 partially explains the robustness of the other parameters against its misestimation. As T 2 can be determined by the use of alternative techniques, the remaining parameters are generally of greater interest. We would like to stress, however, that we are not suggesting the two-pool model is an exact representation of the distribution of macromolecular pools or of their exchange. We simply report that all of these relatively simple signal equations appear to be able to provide consistent and practical information which might be useful in a clinical/clinical research context. Although we did not directly compare the sensitivity of different signal equations to changes similar to those caused by pathology, we have attempted to provide useful guidelines to the choice of the most appropriate one for a particular experimental setting.
Although suggesting optimal sampling schemes is beyond the scope of this paper, it is clear from Fig. 1 that low offset frequencies (<1 kHz) should be avoided when using Ramani's equation. This is consistent with the set of points used by Ramani et al. in their original paper [5], and with the findings of Portnoy and Stanisz [7], and confirms that in pulsed MT experiments it is difficult to model the behavior of magnetization at low off-resonance frequencies. Our results also confirm the observation of Portnoy and Stanisz [7] that Ramani's model tends to underestimate the signal especially close to the Larmor frequency. This deviation at low frequency offsets is likely to be a consequence of the poor ability of the CWPE approximation to characterize the liquid pool. As noted by Portnoy and Stanisz [7], the cut-off of 1 kHz they empirically determined might vary with field strength and MT pulse amplitude and bandwidth.
We should also note that it is apparent from Figs. 2-4 that F and R A are strongly inversely correlated. This is a trivial consequence of the typical values of F in the human brain (of the order of 10 À1 ), for which the denominator on the right hand side of Eq. (18) is approximately equal to RM B 0 F , effectively coupling R A and F. This corresponds to conditions of rapid exchange [9].
The use of numerical simulations allows the investigation of conditions for which it is impractical to acquire in vivo data, for example very narrow or very long MT pulses. As expected, and shown by others [7], even in noise-free condition, the estimates of the MT parameters provided by Ramani's signal equation depend on the duty cycle (in a non-linear fashion), suggesting the need for pulse sequences with a duty cycle of at least 50%. The estimates obtained using Sled and Pike's equations, on the other hand are, with the exception of T A 2 , less sensitive to this parameter.
We also showed that the degree of T 1 -weighting (i.e. of saturation of the excitation) has a large effect on the MT parameters estimated by fitting Ramani's equation. This is a direct consequence of the assumptions underlying it. We therefore strongly recommend avoiding the use of Eq. (17) to fit data which deviate substantially from those assumptions. The two equations proposed by Sled and Pike, on the other hand, by removing the ''steady state'' assumption, appear to be less sensitive to this problem. The RP approximation gives estimates that are substantially insensitive to changes in T 1 -weighting, while the CW approximation becomes slightly affected by it for excitation flip angles larger than 14°(assuming a TR = 30 ms).
In this respect, the results obtained from synthetic data are in keeping with in vivo measurements. Ramani's equation's estimates of F obtained from the more heavily T 1weighted data ( Table 2, flip angle = 15°) are substantially lower than those obtained with less heavily weighted data, while estimates of T A 2 are larger. Conversely, when fitting Sled and Pike's equations to the 15°dataset, the estimates of F are slightly larger compared to those from the 5°dataset, although the difference is not statistically significant. A surprising finding from in vivo data was that estimates of T B 2 obtained from all the equations increased slightly with the excitation flip angle, while this was not observed with simulated data. This discrepancy between simulated and real data may be explained by the presence of noise, and by the use of a smaller number of points for the real data. Noise in the raw data critically affects the estimation of MT parameters, as shown by our fourth simulation (Fig. 3). Monte Carlo simulations also suggest that at high SNR (P150) Sled and Pike's estimates (at least of some parameters) are at least as precise as (and generally more precise than) Ramani's estimates, but that Ramani's equation provides more precise answers at lower SNR. We also note that at low SNR the uncertainty associated with the estimated parameters is much larger than the systematic difference between the mean parameters estimated by each signal equation. For comparison of the in vivo results with simulations, an estimate of the in vivo SNR is needed. The typical SNR for our system, measured in white matter on 5°-excitation minimally MT-weighted 3D SPGR images acquired using a 8-channel head coil and the parameters described in the paper, ranges (due to non-uniformity of the receive coil) between 60 and 100 [20]; other systems are likely to be similar. It would be interesting to compare the standard deviations obtained by Monte Carlo iterations to the variance lower bound predicted by Cramer-Rao theory [21], as well as to the standard deviation obtained from real data using the bootstrap method [22].
A further limitation to this analysis of sensitivity to noise is in the use of suboptimal schemes for both synthetic and real data. For all quantitative techniques based on model fitting, the precision and accuracy of the parameter estimates depend on the choice of the sampling points, and we have previously shown [20] that the error in parameter estimates can be reduced by factors around 2 or 3 by using optimized sampling schemes. The precision of the estimates obtained from all three equations could therefore dramatically improve by using a more suitable set of MT points; we recommend such optimization for practical applications, but this was not possible here as we needed identical schemes to allow direct comparison of the three equations.
Nevertheless, since the schemes here used were likely to be equally suboptimal for all signal equations (and are similar to those used by the authors of the original papers), we believe that our analysis should not have been unduly ''unfair'' to any of the solutions tested. It should also be noted, though, that the acquisition protocol reported by Sled and Pike [6,11] is based on a 60-point scheme to sample a single 7-mm thick slice. Here, despite the consequences on both scan time and SNR, we use only 20 points (for in vivo experiments) and smaller voxels as, for most clinical applications, such whole brain coverage and thinner slices are likely to be essential. Despite reducing the number of points by two thirds, a scan time of 35 min (required by the protocol we used in vivo) is still too long to be feasible in a clinical setting. We also estimated five parameters (RM B 0 ; F ; T B 2 ; T A 2 and M 0 ) directly from the fitting, while the original paper by Sled and Pike [11] recommended the use of an independent estimate of the relative proton density, which allows the MT curves to be normalized to one, reducing to four the number of parameters to be extracted. While this clearly would improve the precision in the estimated quantities, it also lengthens the scan time, making it again less attractive for clinical applications, where the number of acquisitions is typically restricted by time constraints. A better fit could be obtained also by iteratively repeating the fitting using the estimate of M 0 obtained from the current iteration as the starting point for the next. Preliminary results suggest that two iterations should suffice. Furthermore, simulated data suggest that when high SNR and a large number of MT points are available, fitting Sled and Pike's signal equations provides more accurate results, and therefore is preferable. Between the two variants, the main differences seems to be in sensitivity to the degree of T 1 -weighting in the acquisition sequence (as a consequence of saturation from the excitation pulse), suggesting that, in the absence of this confounding factor, the CW variant can be used without major disadvantages, given its reduced complexity. It should be noted, however, that we used a single set of MT parameters to create the simulations, without any attempt to explore the consistency of these results for a different type of tissue (e.g. gray matter). Our results, therefore, are limited to this specific case until confirmed by further experiments.
Regarding the in vivo results, it should be noted that we made no attempt to correct in vivo data for B 1 inhomogeneities (although the use of the body coil for transmission should provide a fairly uniform B 1 distribution at 1.5 T). Both the measurements of R Aobs , and the MT fitting are affected by deviations from the nominal flip angle [6,18]. However, we expect the error introduced by B 1 inhomogeneity to equally affect the three solutions, and therefore not to have major consequences on the conclusions drawn from our experiments.
In the present work we have restricted our analysis to the comparison of three signal equations derived from the two-pool model to predict signal intensity in pulsed MT experiments, without any attempt to modify them, or compensate for their limitations. Several aspects of MT modeling such as the quantification of the effects of the excitation pulses on the macromolecular pool, which is typically considered negligible [2,6], need to be addressed. It also would be interesting to investigate the dependency of these 3 equations on other sequence parameters such as TR, as this quantity controls the efficiency of magnetization transfer [23]. This analysis would be complementary to the investigation of duty cycle effects we performed, and may confirm whether the relative insensitivity of Sled and Pike's CW approximation to duty cycle is maintained for all pulse sequences.
Furthermore, it would be interesting to explore possible modifications of Ramani's equation to account for the imaging parameters of the pulse sequence (for example by incorporating an additional contribution of [1-cos(alpha)]/TR to the CW saturation rate acting on the free pool). Finally, providing reliable information about the optimal number of sampling points and their distribution would yield an additional element towards the choice of the most appropriate equation for a given application. All these areas deserve further investigation, which we hope to pursue in the future.

Conclusion
We have shown that (1) Sled and Pike's CW signal equation provides the most accurate estimates of MT parameters in the absence of noise; (2) Ramani's signal equation is sensitive to changes in the MT duty cycle, although this effect becomes stable for duty cycles P50%; (3) Ramani's equation is (as expected) not suitable for fitting T 1weighted data, and doing so leads to underestimates F and RM B 0 when the amount of T 1 -weighting is high. An error, albeit of smaller magnitude, is also introduced into Sled and Pike's CW approximation estimates when the amount of T 1 -weighting is high. The RP variant is, on the other hand, extremely robust to this effect; (4) in data with SNR typical of in vivo protocols, although the estimates of RM B 0 ; F , and T B 2 obtained from the fit of all three equations show some differences, the magnitude of the difference is smaller than the between-voxels within-ROI standard deviation, provided that T 1 -weighting of the imaging sequence is minimal; (5) Sled and Pike's equations are slightly less robust than Ramani's one in conditions of low SNR.