Extended phase graph formalism for systems with magnetization transfer and exchange

Purpose An extended phase graph framework (EPG‐X) for modeling systems with exchange or magnetization transfer (MT) is proposed. Theory EPG‐X models coupled two‐compartment systems by describing each compartment with separate phase graphs that exchange during evolution periods. There are two variants: EPG‐X(BM) for systems governed by the Bloch‐McConnell equations, and EPG‐X(MT) for the pulsed MT formalism. For the MT case, the “bound” protons have no transverse components, so their phase graph consists of only longitudinal states. Methods The EPG‐X model was validated against steady‐state solutions and isochromat‐based simulation of gradient‐echo sequences. Three additional test cases were investigated: (i) MT effects in multislice turbo spin‐echo; (ii) variable flip angle gradient‐echo imaging of the type used for MR fingerprinting; and (iii) water exchange in multi‐echo spin‐echo T2 relaxometry. Results EPG‐X was validated successfully against isochromat based transient simulations and known steady‐state solutions. EPG‐X(MT) simulations matched in‐vivo measurements of signal attenuation in white matter in multislice turbo spin‐echo images. Magnetic resonance fingerprinting–style experiments with a bovine serum albumin (MT) phantom showed that the data were not consistent with a single‐pool model, but EPG‐X(MT) could be used to fit the data well. The EPG‐X(BM) simulations of multi‐echo spin‐echo T2 relaxometry suggest that exchange could lead to an underestimation of the myelin‐water fraction. Conclusions The EPG‐X framework can be used for modeling both steady‐state and transient signal response of systems exhibiting exchange or MT. This may be particularly beneficial for relaxometry approaches that rely on characterizing transient rather than steady‐state sequences. Magn Reson Med 80:767–779, 2018. © 2017 The Authors Magnetic Resonance in Medicine published by Wiley Periodicals, Inc. on behalf of International Society for Magnetic Resonance in Medicine. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

The EPG method is a Fourier approach to solving the Bloch equations, and therefore assumes that tissues are characterized by a single set of relaxation parameters. It is recognized, though, that a single-compartment approach fails to fully characterize complex biological tissues in many circumstances. Instead, coupled multicompartment models have been proposed. The Bloch-McConnell (BM) equations (17) are a general form for describing systems that are coupled via a general exchange process, with a modification to further describe magnetization transfer (18,19). The MT effects in particular have been shown to be strong determinants of observed signals in human tissue (e.g., ref. (20) shows that on-resonance MT effects are expected to change the signal from balanced steady-state free-precision (SSFP) sequences in brain by approximately 30% and muscle by 50%).
The EPG formalism provides a computationally efficient method for the modeling of MR sequences that also gives intuitive insight into signal formation, by isolating different pathways that can lead to echo formation. Currently, it is not possible to use EPGs to model multicompartment systems with exchange (different compartments in nonexchanging systems can simply be modeled separately and then averaged (16)). Hence, this work seeks to extend the EPG formalism to model such systems. This is increasingly relevant to emerging transient phase relaxometry approaches such as magnetic resonance fingerprinting (MRF) (21), which require numerical simulation, as opposed to more traditional steady-state methods for which analytic or closed-form solutions are often available (e.g., DESPOT or multicomponent DESPOT methods (22)).
We first outline the proposed "EPG-X" method, and then validate it against isochromat-based simulations and established steady-state solutions for gradient-echo imaging. The EPG-X calculations are then used to explore some test cases to illustrate the effect of the new approach. Test cases include multislice versus single-slice TSE, and two different relaxometry methods: multi-echo Carr-Purcell-Meiboom-Gill (CPMG) data for multicomponent T 2 estimation, and gradient-echo imaging with modulated flip angles (similar to MRF). Experimental data were collected, both in vivo and on phantoms.
THEORY For simplicity, we consider in this work sequences with equidistant timing and an unbalanced gradient in a single direction. In this case, the intravoxel magnetization distribution resulting from a sequence of RF and gradient pulses may be characterized by the gradient-induced phase c during some fixed time period Dt. An idealized voxel is defined by the interval c 2 ½Àp; p with uniform density of magnetization in this range. In the EPG representation the magnetization is represented by configuration statesF n andZ n , which, using the notation from the introductory review from Weigel (3), are defined as n¼À1Z n e inc [1] The signal at any time is the voxel average of M þ (i.e., the value ofF 0 ).

Extension to Two Compartments
Consider a two-compartment system, arbitrarily labeled as a and b with thermal equilibrium magnetizations M a 0 ¼ ð1 À f ÞM 0 and M b 0 ¼ fM 0 , respectively. M 0 is the total magnetization and f is the fraction in compartment b, which is conventionally assumed to be smaller. The BM equations governing evolution of this system in the absence of RF pulses may be written as z T correspond to transverse and longitudinal magnetization in each compartment. The other matrices are [3a] [3c] In these expressions, k a is the exchange rate from compartment a to b, which is related to the reverse exchange rate k b via k a M a 0 ¼ k b M b 0 to preserve balance at thermal equilibrium. R 1,a is the longitudinal relaxation rate for compartment a (i.e., 1/T 1,a ); and d b is the frequency offset (Hz) for compartment b relative to compartment a.
An EPG description of the two-compartment model must include "states" that correspond to magnetization from a and b. Taking the Fourier transforms of Equation [2] and writing in terms of F n ¼ ½F where superscripts a and b indicate the compartment, and it is understood that the full expression for the intravoxel magnetization distribution consists of sums over n as in Equation [1]. The solution to Equation [4] is F n ðt þ DtÞ ¼ exp ðK T þ XÞDt F n ðtÞ. We may take advantage of the fact that matrices K T and X commute to re-express the matrix exponential as a product of terms expðK T DtÞexpð XDtÞ and define operators as W expðXDtÞ [6] N T expðK T DtÞ [7] such that F n ðt þ DtÞ ¼ WN T F n ðtÞ. Defining the dephasing during time Dt as c ¼ Àv z Dt, we identify W as the familiar "shift" operator, which increments the index of transverse states as with the standard EPG algorithm. Equation [5] for Z n is homogeneous for n 6 ¼ 0, but inhomogeneous for n ¼ 0. The solutions for the two regimes are Z n ðt þ DtÞ ¼ N L Z n ðtÞ ðn 6 ¼ 0Þ [8] Z 0 ðt þ DtÞ ¼ N L Z 0 ðtÞ þ ðN L À IÞK À1 L C ðn ¼ 0Þ [9] N L expðK L DtÞ [10] This is in accordance with the standard EPG framework in which longitudinal recovery occurs only in the n ¼ 0 state. The relaxation-exchange operators N T and N L replace the relaxation operators from the standard EPG approach. Relaxation and exchange are treated simultaneously using these combined operators; it is not generally possible to separate these processes, because the corresponding components of the K matrices do not commute. The operators do not have a simple analytic form; instead, the matrix exponentials are evaluated numerically. The form of the K matrices means that exchange couples only states of the same "type" and "order" (e.g., F a n $F b n andZ a n $Z b n ).

Solution for RF Pulses
In the classic EPG framework, RF pulses mix together all states for a given order n via transition operator (3) as follows: where a and f are the RF pulse flip angle and phase. In the two-compartment case, the overall transition matrix to apply to the full system ½F a nF Ãa ÀnZ a nF b nF Ãb ÀnZ b n T is simply where a a is the flip angle for compartment a, and so on. This paper only considers scenarios in which the flip angle and phase are the same for both compartments. For situations in which d b is large compared with the RF pulse bandwidth, these could be different.

Magnetization Transfer
A different formulation is generally used when describing MT effects in tissues with a "semisolid" component (18). In this case, compartment b is often referred to as the "bound" or "restricted" pool and is assumed to represent highly immobile protons whose T 2 is very short (in the order of 10 ms). In this case we assume that compartment b has no transverse magnetization; hence, statesF b n and F Ãb Àn are dropped such that F n ¼ ½F a nF Ãa Àn T . In this formulation, F n are treated exactly as in the classic EPG case (i.e., subject to T 2 relaxation and shifts caused by gradients), and d b is not defined. The coupled longitudinal states Z n still evolve as per Equations [8] and [9].
The effect of RF pulses on compartment a is to rotate the magnetization as previously described. However, for compartment b (the "bound pool"), RF pulses act so as to directly saturate the longitudinal component with saturation rate W ðv z Þ, which for pulsed saturation is defined as (19) where B 1 (t) is the RF pulse waveform and t rf is its duration. This is a function of off-resonance frequency v z , because it depends on the absorption lineshape Gðv z Þ. Different candidate lineshapes have been proposed for modeling semisolids in biological tissues with Gaussian (19) and super-Lorentzian shapes (23) used primarily. The overall RF transition matrix is therefore

Summary of Proposed Theory
To summarize, we have introduced extensions to the EPG formalism to account for multicompartment systems with exchange (EPG-X). There are two variants: one for systems governed by the BM equations, and one for the variant of BM often used for MT, in which one compartment has negligible transverse magnetization. Both effectively use two coupled EPG calculations: one for each compartment, although the MT variant uses a second compartment with longitudinal components only. These are summarized diagrammatically in Figure 1.

METHODS
The EPG-X framework could in principle be used to simulate the response of a two-compartment system for any pulse sequence. The theory is illustrated by simulating four separate scenarios. The model tissue parameters used are outlined in Table 1; the subset used in each simulation is indicated in the text. In all models, compartment b is the smaller one. The myelin-water exchange model is motivated by literature on exchange between water trapped within myelin layers (compartment b) and intra-and extra-axonal water (24), but uses rounded values rather than taking numbers directly from a single source. In this model the myelin-water mean residence time (t b Þ is given by (the parameters in Table 1 make t b ¼ 100 ms consistent with (22), although estimates vary significantly) (24). There is evidence in literature for small nonzero d b in this system (25) hence some simulations used d b 6 ¼ 0.
The white matter and caudate nucleus MT models are taken from measurements made at 1.5 T by Gloor et al ( Table 1, ref. (26)). The super-Lorentzian lineshape function has been used for Gðv z Þ; as in (26), the function was extrapolated between 61 kHz by fitting a spline to avoid the singularity at zero frequency. Systems with MT effects have been shown to exhibit bi-exponential evolution of longitudinal magnetization (27,28) with relaxation rates given by the two eigenvalues of K L . The T 1 observed from typical inversion recovery measurements (denoted as T obs 1 ) is derived from the smaller (less negative) eigenvalue (27) as : [15] For the white matter MT model, T obs 1 ¼ 779 ms. All numerical simulations and analyses were performed using MATLAB R2015a (The MathWorks, Natick, MA, USA). A fully functional implementation is available to download at http://www.github.com/mriphysics/ EPG-X . Code and experimental data for generation of all results presented in this paper are included (hash 4e8ba3d was the version at time of submission).

Test 1: Steady State and Transient Behavior of Gradient-Echo Sequences
Steady-state expressions for spoiled gradient-echo (SPGR) and balanced steady-state free precession (bSSFP) sequences have been derived for BM and MT models elsewhere in literature (22,26,29,30). The Appendix gives these solutions in terms of the matrices already introduced in the Theory section.
To validate the transient behavior, isochromat-based simulations were also performed by solving the BM equations directly for a range of dephasing angles c, then integrating over the interval c 2 ½Àp; p. Simulations were performed using varying numbers of isochromats (N iso ) spaced evenly over range ½Àp; p.
In this test, the EPG-X(MT) simulations used the white matter MT model and EPG-X(BM) simulations used the myelin-water exchange model (Table 1). For comparison, single-component (classic EPG) calculations were performed using T 1 ¼ 799 ms and T 2 ¼ 45 ms. Sequences with repetition time (TR) of 5 ms and flip angle (a) of 10 were simulated for 5 Â T 1 , after which a steady state was assumed to have formed. For SPGR, simulations were repeated for different values of the quadratic RF spoiling phase increment F 0 from 0 to 180 . W was computed by assuming hard pulses with maximum amplitude 13.5 mT.
Equation [1] implies that M þ ðcÞ and M z ðcÞ may be obtained from the EPG predictions by performing an The RF pulses directly saturate these states (Eq. [14]), and they exchange directly with their equivalents in compartment a. For clarity, in all diagrams, relaxation and exchange effects are only depicted for the periods following the first two RF pulses. Although depicted by separate arrows, relaxation and exchange processes are governed by single combined operations. inverse fast Fourier transform (FFT) over "order" parameter n (31). For bSSFP, because net gradient area is zero, the familiar off-resonance sensitivity of the bSSFP method may be found by applying iFFT to EPG predictions. These were compared with the steady-state solutions of Equation [A2].

Test 2: Multicomponent T 2 Analysis of CPMG Data
The multicomponent analysis of multi-echo CPMG spinecho data is a well-known method for estimation of myelin fraction in white matter (32). If perfect 180 refocusing pulses are assumed, the multi-echo signal can be analyzed using nonnegative least squares (NNLS) fitting to an exponential model (33). It is recognized that B þ 1 inhomogeneity will introduce other stimulated echoes that make the data deviate from this simple model (16).
To explore any additional effects from exchange, we used EPG-X(BM) simulations with the myelin-water exchange model (Table 1) to simulate multi-echo data for a range of exchange rates, B þ 1 scaling factors, and offset frequencies d b . In each case, the simulated data were analyzed using the classic NNLS fitting approach, from which the estimated small pool fraction was taken as the area of smaller peak in the T 2 spectrum.
Simulations used 50 echoes with 5 ms spacing. k a was varied from 0 to 2.5 s À1 (for f ¼ 0.2, this corresponds to infinite t b (i.e., no exchange) down to 80 ms). B þ 1 scaling factors from 0.75 to 1.25 were included. d b in the range of 6128 Hz (6 1 ppm at 3 T) was investigated with fixed B þ 1 scaling factor of 1.0. The NNLS was performed using the MATLAB function lsqnonneg.

Test 3: MT in Transient Gradient-Echo Sequences
Magnetization transfer effects have been shown to strongly affect the SPGR signal in the steady state (27,29). Transient gradient-echo sequences with variable flip angle, often following inversion pulses, have been used for MRF (21,34); hence, we simulated a simple example of such a sequence to predict transient behavior. The sequence used an adiabatic inversion pulse followed by a series of 256 low flipangle RF pulses whose amplitude was varied sinusoidally (shown in the Results section); some pulse amplitudes were zero to allow for magnetization recovery. The RF pulse energies were 433 ms mT 2 for the inversion and 54:3 Â a 2 ms mT 2 for the small flip-angle pulses (a is the flip angle, rad). A constant TR of 12 ms was used with constant gradient area in each TR period, even those with zero flip angles. The sequence was simulated with EPG-X(MT) using the white matter MT model (Table 1). Both bSSFP and SPGR were simulated using the same timing.
For consistent comparison, single-compartment EPG was simulated using T obs 1 and T 2,a , as these are the parameters that would be measured using standard inversion recovery and spin-echo methods.

Experimental Measurements
Experiments were performed on a Philips Achieva 3T MRI scanner (Best, Netherlands). Two phantoms were made from sample tubes: water doped with 0.1 mM MnCl 2 (expected to have no MT effect), and crosslinked bovine serum albumin (BSA; Sigma-Aldrich, Dorset, United Kingdom), which has been suggested as a good material for replicating MT in human tissue (35). A 10% BSA solution (by weight) was prepared in distilled water, before treatment with glutaraldehyde (Sigma-Aldrich) as described in (35).
The samples were imaged using the SPGR sequence described previously (TR ¼ 12 ms, echo time ¼ 2.9 ms), with frequency encoding aligned with the longitudinal axis of the tubes and phase-encoding switched off to directly record echo amplitudes. Experiments were repeated with RF spoiling phase increment F 0 set to 150 (default) and 117 . Single-compartment relaxation times and apparent diffusion coefficient (D) were measured for each phantom, using inversion-recovery TSE (T obs 1 ), multi-echo spin echo (T 2 ), and diffusion-weighted spin echo (D). For the water phantom T obs 1 ¼ 899 6 5 ms, T 2 ¼ 92 6 2 ms, and D ¼ 2.35 6 0.18 Â 10 À3 mm 2 s À1 ; for BSA T obs 1 ¼ 1290 6 25 ms, T 2 ¼ 92 6 5 ms, and D ¼ 1.92 6 0.29 Â 10 À3 mm 2 s À1 . Separate long TR multi-flip-angle measurements were made to precisely measure the effective B 1 amplitude and M 0 (including receiver coil scaling) in these phantoms; these measurements were used to match the measured echo amplitudes as closely as possible to EPG simulations.
As will be shown later, it was found that diffusion effects needed to be taken into account to accurately model the SPGR sequence. Diffusion can readily be accounted for in the EPG framework (12). We experimented with extending this to multicompartment models by applying the same treatment independently to each compartment; validity of this approach is discussed later. For the results presented in this paper, diffusion effects were only included for those relating to Test 3.

Test 4: MT Effects in Multislice TSE Imaging
Multislice TSE is sensitive to MT effects (36,37), as from the point of view of a given slice location, the acquisition of the other slices may be viewed as repeated offresonant irradiation, leading to attenuation of signals from some tissues in multislice compared with singleslice acquisitions. To investigate this effect, we acquired data from a single healthy adult male volunteer (age 25; written consent was obtained before enrollment) using a Philips 1.5T Ingenia MRI system. A series of multislice TSE images were acquired using refocusing flip angles 180 and 120 and from one to 15 slices (odd numbers only). All acquisitions used 25 echoes, interecho spacing ¼ 7.7 ms, and TR ¼ 5 s. Acquired resolution was 1.5 Â 1.5 mm 2 , slice thickness ¼ 2 mm, and slice spacing ¼ 2.78 kHz for the 180 sequences, and 3.13 kHz for the 120 sequences (RF pulse details are given in Table 2). A gap of 4 mm was used to avoid cross-talk and default "odd-even" slice ordering was used. Images were analyzed by drawing regions of interest in white matter, caudate nucleus, and cerebrospinal fluid.
The sequences were modeled using EPG-X(MT) for both white matter and caudate nucleus MT models listed in Table 1. Off-resonant excitation of other slices can be modeled by trains of pulses with zero flip angle for compartment a, but with saturation still applying to compartment b. The simulations were done from the point of view of the slice at the center of the group, as this was present in all acquisitions; these were run for three TR periods to ensure equilibrium was reached. Figure 2a shows the approach to steady state for SPGR with F 0 ¼ 117 for standard EPG and the proposed variants, compared with the ideal spoiling steady-state values predicted by Equation [A1]. All three curves approach the ideal spoiling steady state (arrows). Figure  2b shows the steady-state values reached by the EPG methods for a range of F 0 compared with the ideal spoiling prediction (as expected, strong variation with F 0 is seen). Supporting Figure S1 explores this behavior further for the EPG-X(BM) simulation by varying d b from 0 to 256 Hz (2 ppm at 3 T). For the values of F 0 that result in good spoiling, not much variation with d b is seen; however, for the "spike" values such as 0 or 120 , an oscillatory dependence on d b is observed.

Test 1: Comparison With Existing Steady-State Gradient-Echo Solutions
Supporting Figure S2 compares the predictions of the transient behavior for varying numbers of isochromats (N iso ) with EPG and EPG-X. Large discrepancies are seen for smaller N iso ; however, for N iso greater than or equal to the number of RF pulses, the two types of simulation agree exactly (differences % 10 À15 ), in line with (31). Figure 3 compares the EPG and direct steady-state solutions for bSSFP for all models; in these simulations, the myelin-water exchange model was used twice, for d b ¼ 0 and 12.8 Hz (0.1 ppm at 3 T). The profile for d b 6 ¼ 0 (Fig.  3d) is markedly asymmetric. Each model produces quite different steady-state behavior; however, in each case the agreement between EPG-X and direct calculation is excellent, as is the agreement with the on-resonance analytic expression for the MT case (red asterisk, ref. (26)).
Test 2: Multi-echo CPMG Relaxometry Figure 4a shows example echo amplitudes from the simulated multi-echo CPMG data; Figure 4b has the corresponding T 2 spectrum from the NNLS analysis. There are two peaks corresponding to the two compartments (T 2,a ¼ 100 ms and T 2,b ¼ 20 ms), and the fraction is estimated by taking the ratio of the peak areas (shaded in Fig. 4b); the NNLS estimated parameters are denoted aŝ f andT 2;b . Figures 4c and 4d show that althoughT 2;b varies strongly with B 1 scaling,f is more strongly dependent on the exchange rate k a . For k a ¼ 2 s À1 (i.e., t b ¼ 100 ms), Figure 4c indicatesf ¼ 0.133, a 33% underestimate. Figure 4e shows thatf also varies weakly with d b .

Test 3: MRF-Style Transient Gradient Echo
The variable flip angle profile used is illustrated in Figure 5a. The figure also illustrates the expected signals (Fig. 5b) and evolution of longitudinal magnetization (unmodulatedZ 0 states (Fig. 5c)) in the white matter model for the SPGR sequence, comparing EPG-X with single-pool EPG using T obs 1 . The signal profiles are different, particularly immediately after the inversion (at the beginning of the sequence) when the magnetization in the MT system recovers more quickly. Supporting Figure  S3 shows an equivalent result for bSSFP. Figure 6 compares experimentally the obtained SPGR data with standard EPG predictions; the data and model are not fitted together, and all necessary parameters and scaling coefficients were measured in calibration experiments. For the MnCl 2 phantom, reasonable agreement is obtained using EPG; however, this is improved by including diffusion effects (blue arrows). Note also how the signal profiles from the two different RF spoiling phase increments F 0 are quite different. For BSA, the match to single-compartment EPG is slightly improved by adding diffusion effects; however, there are systematic differences particularly immediately after the inversion, similar to those observed in Figure 5 (yellow arrows). We hypothesized that these are caused by MT. The experimentally obtained data for BSA were fitted to EPG-X(MT) predictions using least-squares minimization, optimizing over f, T 1,a , T 1,b , k a , and G(0). Fitting used the MATLAB function fmincon including a nonlinear constraint enforcing consistency between estimated parameters and measured T obs 1 (via Eq. [15]). Data for F 0 ¼ 150 8 and F 0 ¼ 117 8 were fit simultaneously. T 2,a , D, and the overall scaling constant were fixed at the measured values. Diffusion was implemented in EPG-X using the same approach as for standard EPG; the validity of this approach will be discussed later. Figure 7 shows the fit that could be obtained; there is very good agreement (root mean square deviation of 1.3%) using k a ¼ 6.2 s À1 , T 1,a ¼ 1763 ms, T 1,b ¼ 363 ms, G(0) ¼ 29.4 ms, f ¼ 0.100. This combination would yield T obs 1 ¼ 1283 ms, which is consistent with the inversion recovery measurement (1290 6 25 ms). As a control, the single-compartment EPG model was also fitted to the data by varying T 1 . In this case, the best fit was T 1 ¼ 1215 ms, but significant residuals remain ( Fig. 7; root mean square deviation of 6.9%). Figure 8 shows single-slice and multislice TSE images from a healthy volunteer, using 180 refocusing pulses. As expected, white matter signal clearly decreases as the number of slices increases. The region-of-interest analysis (plots) compares the image data (error bars) with EPG-X(MT) predictions (red triangles) for white matter, caudate nucleus, and cerebrospinal fluid. Data were scaled to compare with EPG-X predictions by normalizing the mean signal over all experiments for each region of interest to the equivalent mean value from the EPG-X predictions. The data and predictions match to within experimental error in all cases. The EPG-X model correctly predicts the trends seen in white matter and caudate nucleus (the former having a stronger effect) and the effect of changing from 180 to 120 refocusing pulses. As expected, no change was observed in the cerebrospinal fluid signal (no MT effect), indicating that direct cross-talk between slices is negligible.

DISCUSSION
This work has introduced a general framework for the modeling of processes governed by the BM equations or their modified form for MT, using EPG. Essentially, the different compartments are described by separate phase graphs, which exchange with each other during evolution periods. For the MT case, the bound protons have no transverse components, so their phase graph contains only longitudinal states. The new model, referred to as the EPG-X, has been validated by comparing directly with isochromat-based simulations, and against existing steady-state solutions. In addition to numerical validation, the EPG-X(MT) model was found to agree well with in vivo measurements of signal attenuation in multislice TSE, and provided a plausible fit to transient gradientecho signal measurements using crosslinked bovine serum albumin (an MT phantom).
No fundamentally new biophysical model has been proposed in this work; rather, we have shown how existing methods can be incorporated into the EPG framework. As Supporting Figure S2 shows, EPG-X methods produce identical results to direct isochromat-based integration of BM equations, as long as the latter approach is sufficiently sampled (31). The EPG method is wellestablished as a computationally efficient and intuitive approach to modeling MRI sequences, particularly those with multiple echo pathways such as TSE (see (3) for a detailed discussion). The EPG-X extensions seek to retain these advantages for the modeling of more complex multicompartment systems. Calculations retain the same structure as the original EPG approach, allowing the use of established methods for truncating and hence accelerating EPG-X calculations (10,38). In addition to efficiency, the major advantage of the EPG approach is that it allows for intuitive analysis of sequences in which one or more echo pathways are isolated for measurement (TSE sequences and the many variants such as "hyperechoes" (6) are the classic example, MRF using an echo splitting technique (39) is a more recent application). In all of these cases, the proposed EPG-X framework would allow for the effect of exchange or MT to be considered, and this may have a significant effect on predictions.
The examples in this paper focused on two different types of models, both of which are commonly used for quantitative MRI in the brain. The myelin-water exchange model is used for steady-state approaches like  (Table  1) and single-pool comparison. The EPG-X signals are normalized by (1-f) to account for compartment b being "invisible." MT leads to different behavior, particularly immediately after the inversion pulse. c:Z 0 from the EPG and EPG-X models. For EPG-X, the level of saturation of compartment b changes dynamically, leading to altered dynamics of the observed signal when compared with a single-compartment model.
FIG. 6. Experimental SPGR data compared with the single-compartment EPG model. No fitting was performed (relaxation times, B 1 scaling, and receiver/M 0 scaling factors were experimentally measured). When diffusion is not included (top two rows), the match to EPG is not perfect (blue arrows). For the MnCl 2 phantom, once diffusion is included (bottom two rows), the match is very good (normalized root mean square error (NRMSE) $2%). For the bovine serum albumin (BSA) phantom, there remain systematic differences (yellow arrows), suggesting that the single-compartment model is not sufficient. Moreover, the observed signal profiles are quite different for the two different values of F 0 , as predicted by the EPG model. , the k a , T 1,a , T 1,b , f, and G(0) were varied; the best-fit parameters were k a ¼ 6.2 s À1 , T 1,a ¼ 1763 ms, T 1,b ¼ 363 ms, G(0) ¼ 29.4 ms, f ¼ 0.100, which would yield T obs 1 ¼ 1283 ms. For single-compartment EPG, the T 1 was varied; the best fit was 1215 ms. The NRMSE before fit ¼ 9.8%, single-compartment EPG fit ¼ 6.9%, EPG-X ¼ 1.3%. multicomponent DESPOT (22) and multi-echo CPMG relaxometry (24), although for the latter, exchange between compartments is usually neglected. Typically, two significant peaks are observed in white-matter T 2 spectra from multi-echo CPMG data (32), with the faster relaxing component (T 2 $ 20 ms) identified as myelin water; the fractional size of this component is used to estimate the myelin-water fraction. More recent work has found large variability in apparent myelin-water fraction from tracts with similar myelin content in rat spinal cord (40), with differing levels of exchange suspected as a potential reason. Results from Test 2 (Fig. 4) suggest that exchange can lead to significant underestimation of the myelin-water fraction, and hence might be able to explain these experimental observations. Therefore, EPG-X could potentially be used for incorporating estimation of exchange parameters into multi-echo relaxometry.
We also investigated the effect of compartment b frequency offset d b ; this is relevant to white-matter imaging in particular, as many studies have shown that there is a nonzero mean frequency offset resulting from local susceptibility effects, which is also orientation dependent (25,41). Supporting Figure S1 shows that the effect on SPGR is relatively minor; however, surprisingly, Figure  4e shows that d b does interact with NNLS analysis and alters the estimated myelin-water fraction for CPMG data. For the very small frequency shifts expected (<0.1 ppm (41)), however, this would lead to only small biases. Results from Test 1 concerning bSSFP (Fig. 3d) suggest that a frequency shift of this size would lead to an asymmetric off-resonance profile. This result is very similar to measurements made in white matter by Miller et al (42); the proposed framework would allow for further investigation, including exchange. Similarly, recent work on chemical exchange saturation transfer (CEST) (47) has proposed using observed asymmetries in the bSSFP frequency profile as an alternative detection mechanism to more standard z-spectrum saturation measurement. The latter approach to CEST measurement usually requires long saturation pulses; EPG is not useful for simulation of this interaction, as it is designed to characterize sequences rather than individual RF pulses. However the "sequence level" CEST detection method proposed in (47) can be simulated with EPG-X(BM), and further work could use the framework to investigate CEST effects for other pulsed sequences. The upper row of plots shows data for 180 refocusing pulse sequence; the lower row is for 120 refocusing pulses. The signal from white matter is strongly attenuated, particularly with the 180 pulse sequence. Caudate nucleus shows similar, but slightly less severe, attenuation. The EPG-X(MT) predictions match the experimental data within error. As expected, cerebrospinal fluid shows no attenuation, indicating no direct slice cross-talk.
The second type of model considered in this work focused on MT effects. For MRF, a simulation using white-matter parameters ( Fig. 5 and Supporting Fig. S3) shows that MT will lead to altered behavior, which is particularly pronounced after inversion because of biexponential recovery, which has been observed in human white matter (e.g., (28)). Measurements on doped water and BSA phantoms (Fig. 6) also suggest that although the former can be well modeled using standard EPG (good agreement is obtained with no fitting to the data), a large residual is observed in the period after the inversion pulse in the BSA. This type of deviation might be a source of bias in MRF; however, the manner in which this would affect the final result remains to be investigated. The fact that the observed "fingerprint" profiles are affected by MT could also imply that there is potential for quantitative MT characterization with this type of sequence. To demonstrate this, we used nonlinear fitting to estimate the MT specific parameters (Fig. 7) from the BSA data. The fitted parameter values are consistent with the inversion recovery-measured T obs 1 . Furthermore, the estimated f ¼ 0.1 agrees with the fact that the phantom was prepared with 10% BSA by weight. Contrastingly, fitting the single-pool EPG model gave a best-fit T 1 that is inconsistent with the inversion recovery, and still produces large residuals. The implication is that these data are not fully characterized by a single T 1 value, which is consistent with observations made by others (28,29). Our results indicate that the EPG-X model can explain the observed data more accurately, producing numbers that are consistent with independent inversion-recovery measurements. Future work could explore the use of EPG-X as the basis for generation of dictionaries for MT measurement via MRF.
Although most of the preceding exchange-related effects are quite subtle, Test 4 showed that significant effects attributable to MT (signal attenuation in multislice TSE) can also be accurately predicted by EPG-X (Fig. 8). Previous work found that a semi-empirical model (37) could predict this signal loss, but required experimental measurement of sequence and tissue-dependent coefficients to do so. The present work made accurate predictions using only literature tissue parameters (Table 1, ref. (26)) and scanner-reported sequence parameters. Hence, the EPG-X method could be useful for pulse-sequence parameter optimization that accounts for MT effects; these are not insignificant, for example more than a 30% reduction in white-matter signal was observed when moving from single slice to a 15-slice multislice acquisition.
Multicompartment models seek to explain complex underlying biological systems, but the choice of model is key. White matter is a particularly complex tissue; in this paper, we used "exchange" (BM) and MT models that are both relevant to white matter. This reflects the range of models that currently exist in the literature; the choice of model depends on the type of sequence being modeled as well as the tissue. For multi-echo T 2 relaxometry, the BM model is most relevant, as it is seeks to measure multiple compartments with appreciable T 2 . To evaluate the effect of off-resonant saturation in multislice imaging, the MT model is most relevant. More complex mixed models that include multiple BM and MT compartments (43,44) have also been proposed; Liu et al proposed a general framework for such systems (45). In this work we focused on two compartment models, but in principle the same arguments used in (45) can be used to extend the EPG-X formalism to more compartments as well.

Implicit Assumptions and Limitations
The EPG-X model assumes that the underlying magnetization in both compartments forms a spatial distribution at a subvoxel-length scale. Exchange couples Fourier configurations in one compartment with the same configuration in the other compartment, which is equivalent to assuming that exchange interactions couple these distributions locally (i.e., the magnetization in compartment a at one location couples with compartment b at the same location, but not adjacent locations). Isochromat-based modeling methods (e.g., (45,46)) make the same assumption in the spatial domain. This is physically reasonable, as both chemical exchange and MT occur at the level of individual molecules, and compartmental exchange occurs over diffusion-length scales, smaller than required for modeling the subvoxel magnetization distribution at spatial resolutions relevant to MRI.
Our data showed that diffusion effects must be accounted for to accurately match observed signals to EPG models. As far as we are aware, there is no commonly adopted equivalent model for the multicompartment case, as we are effectively combining the Bloch-McConnell and Bloch-Torrey equations. For the measured BSA data (Figs. 6 and 7), we took the most basic approach, which was to treat diffusion effects independently for both compartments but with the same diffusion coefficient for each. It might be expected that actually diffusion coefficients would be quite different for each compartment, and this could readily be achieved within the same framework. Further work is needed to identify the most appropriate biophysical model.
Finally, this work has been presented in terms of the "regular time increment" version of the EPG framework, which allows configuration states to be considered in terms of integer indexes only. Sequences with variable gradient directions and/or nonuniform timing are described instead using a continuous Fourier transform (see (3) for a detailed discussion). The theory put forward in this paper would generalize readily to this approach, as none of the exchange-related operators are explicit functions of space.

CONCLUSIONS
Extensions to the EPG framework to systems governed by the BM equations and modified forms for pulsed MT have been proposed. The new formalism named EPG-X may be used to efficiently model a wide range of pulse sequences, and results indicate that for steady-state sequences EPG-X gives equivalent predictions to commonly used solutions. The EPG-X model could prove to be useful for quantitative imaging, particularly for nonsteady-state sequences in which accurate modeling of the transient response is necessary.

ACKNOWLEDGMENT
The authors thank Dr. Ana Baburamani for the preparation of the phantoms. This work was supported by the

Spoiled Gradient Echo
For an "ideally spoiled" sequence, we consider only a steady state formed by longitudinal components. The measured signal is given by where time increment Dt ¼ TR is used in the definition of N L . R represents the excitation and sampling of longitudinal magnetization, and is defined differently for the BM and MT scenarios:

Balanced SSFP
For bSSFP, a coherent steady state among all components is considered. For the BM case written out for M ¼ ½M a þ M a À M b þ M b À M a z M b z T , the steady-state signal is given by D represents a rotation of 180 about the z-axis for both pools (accounting for the phase alternation in bSSFP), and R ¼ ½1 0 1 0 0 0. H is defined in the same way as T (Eq. [12]), except that the rows and columns are reordered to group terms for transverse components first, and longitudinal components second.
For the MT case, the same expression is used, except the rows and columns corresponding to M b þ and M b À are deleted in all matrices, as are all coupling terms relating to M a þ and M a À : The system is a 4 Â 4 matrix, the modified RF pulse matrix (Eq. [14]) is used, C ¼ ½0 0 R 1;a M a 0 R 1;b M b 0 T , and R ¼ ½1 0 0 0:

SUPPORTING INFORMATION
Additional Supporting Information may be found in the online version of this article.  Fig. 2b). As resonance offset d b is varied, no change in signal is observed for values of U 0 that yield good spoiling (i.e., flat signal that is close to the desired ideal spoiling steady-state value). However, for the spike values such as U 0 5120 or 90 , the steadystate signal oscillates as a function of d b . Fig. S2. The EPG and EPG-X predictions plotted against isochromat predictions for SPGR approach to steady state with repetition time 5 5 ms, flip angle 5 10 . Isochromat ensemble simulations were repeated with increasing numbers of isochromats N iso (from 10 to 1000); signal prediction comes from averaging the transverse magnetization M 1 over the whole ensemble. a, c, e: The EPG and EPG-X predictions (solid black line) are compared with isochromat simulations using differing N iso. Each different colored line is a different N iso ; most prominent are blue 5 10, rust 5 30, and yellow 5 50. b, d, f: Root mean square deviation between EPG and isochromat simulation for each case. The root mean square deviation drops as N iso is increased and suddenly falls to approximately 10 215 for N iso ! 200 (number of RF pulses). At this point, EPG and isochromat predictions are effectively identical. Fig. S3. a: Predicted signal as a function of pulse number and w for balanced steady-state free precision. b: Profiles for w50 and w5 p 2 = (see dotted lines in (a)). The oscillations in the w50 case are caused by the "stepped" nature of the changing flip angles (see Fig. 5a). The effect of MT alters the signal dynamically, particularly after the inversion as with SPGR (Fig. 5). The difference between EPG and EPG-X also changes as a function of off-resonance parameter w. c: Z 0 profiles; for EPG-X, the saturation of compartment b varies dynamically.