Journal of Magnetic Resonance 206 (2010) 59–67 Contents lists available at ScienceDirect Journal of Magnetic Resonance journal homepage: www.elsevier.com/locate/jmr Evaluating the accuracy and precision of a two-compartment Kärger model using Monte Carlo simulations M. Nilsson a,*, E. Alerstam b, R. Wirestam a, F. Ståhlberg a,c, S. Brockstedt d, J. Lätt e a Department of Medical Radiation Physics, Lund University, Lund, Sweden Department of Physics, Lund University, Lund, Sweden c Department of Diagnostic Radiology, Lund University, Lund, Sweden d Radiation Physics, Skane University Hospital, Lund, Sweden e MR Department, Center for Medical Imaging and Physiology, Skane University Hospital, Lund, Sweden b a r t i c l e i n f o Article history: Received 3 March 2010 Revised 27 May 2010 Available online 9 June 2010 Keywords: DW-MRI Diffusion Microstructure Monte Carlo simulations Diffusion time a b s t r a c t Specific parameters of the neuronal tissue microstructure, such as axonal diameters, membrane permeability and intracellular water fractions are assessable using diffusion MRI. These parameters are commonly estimated using analytical models, which may introduce bias in the estimated parameters due to the approximations made when deriving the models. As an alternative to using analytical models, a database of signal curves generated by fast Monte Carlo simulations can be employed. Simulated diffusion MRI measurements were generated and evaluated using the two-compartment Kärger model as well as the simulation model based on a database containing signal curves from approximately 60 000 simulations performed with different combinations of microstructural parameters. A protocol based on a pulsed gradient spin echo sequence with diffusion times of 30 and 60 ms and with gradient amplitudes obtainable with a clinical MRI scanner was employed for the investigations. When using the analytical model, a major negative bias (up to approximately 25%) in the estimated intracellular volume fraction was observed for short exchange times, while almost no bias was seen for the simulation model. In general, the simulation model improved the accuracy of the estimated parameters as compared to the analytical model, except for the exchange time parameter. Ó 2010 Elsevier Inc. All rights reserved. 1. Introduction Specific characteristics of neuronal tissue microstructure can be assessed using diffusion magnetic resonance imaging (MRI), useful for investigating disease and development of cerebral matter. Examples of such characteristics are axonal or cell diameters, cell membrane permeability and intracellular water volume fractions. Axonal diameters have been estimated using nuclear magnetic resonance (NMR) spectrometers in excised neuronal tissue [1–3]. Using clinical MRI scanners, microstructural features of phantoms have been assessed [4] and recent studies suggest that this is feasible also in vivo [5–8]. Experiments aiming at the investigation of tissue microstructure are commonly designed to measure the diffusion-weighted (DW) signal intensity for different degrees of diffusion encoding (b-values) and different diffusion times (TD), thereby probing the microstructure at different length and time scales [9]. Measurements designed to selectively probe the microstructure of white matter fibres are commonly performed in a single direction per* Corresponding author. E-mail address: markus.nilsson@med.lu.se (M. Nilsson). 1090-7807/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.jmr.2010.06.002 pendicular to the fibre structure [1–5], although models have been developed to also estimate the fibre direction from diffusion measurements in multiple directions [6]. The analytical models employed to estimate the microstructural features comprise the DW signal from different intra- and extracellular components. However, it is challenging to accurately model the signal curves from these individual components, due to, for instance, long diffusion encoding durations, exchange between components and different relaxation rates for different components. The modelling challenges are met by various assumptions, first regarding the tissue geometry and subsequently about the measurements on the geometry. For example, the tissue structure is commonly assumed to be comprised of ideal geometries, such as parallel planes, cylinders or spheres. For these ideal geometries, analytical models describing the signal curves have been derived under the assumption that the measurement produces a Gaussian phase distribution or that the diffusion encoding duration is infinitesimal [5,6,10–12]. For a tissue model based on one intracellular and one extracellular compartment, Monte Carlo simulations have been employed to investigate the validity of these measurement models. Meier et al. showed that it is possible to extract the intracellular exchange time under appropriate experimental conditions 60 M. Nilsson et al. / Journal of Magnetic Resonance 206 (2010) 59–67 [13]. Nilsson et al. confirmed this, but found that the estimated intracellular volume fraction may be biased in case of short exchange times [5]. These studies suggest that a given measurement protocol yields accurate and precise results only for certain combinations of microstructural parameters. Alexander showed that under experimental conditions obtainable using a clinical MRI scanner, axonal diameters larger than approximately 5 lm are possible to estimate accurately in a two-compartment system, using an analytical model without exchange [6]. However, in the presence of compartmental exchange, the accuracy and precision of the estimated parameters has not been established and requires further investigation. As an alternative to using analytical measurement models, microstructural parameters may be extracted using a database of signal curves obtained by using Monte Carlo simulations, similar to the approach by van Enden et al. [14]. Signal curves in the database can be considered as ‘‘fingerprints” of the microstructure and measured signal curves can be matched against these. This approach, herein denoted the simulation model, requires less assumptions regarding the measurement model than the analytical models, but relies on fast simulations in order to build a sufficiently large database without inordinate demands on computation times. The aim of this study was to investigate the accuracy and the precision when estimating microstructural parameters in a twocompartment system with exchange, using simulated diffusion MRI data. The microstructural parameters were estimated from signal curves simulated for experimental settings and noise levels obtainable with a clinical MRI scanner. The performance in estimating the microstructural parameters was compared between two different measurement models, i.e. the analytical two-compartment Kärger model and the simulation model. 2. Theory In DW magnetic resonance experiments, pulsed magnetic-field gradients cause a signal attenuation by producing a phase dispersion of the diffusing spins. The DW signal S is given by S¼ Z 1 cosð/Þf ð/Þd/ 1 ð1Þ where / is the spin phase and f(/) is the probability density of spin phases in a given volume and its Monte Carlo approximation S* is obtained according to S ¼ n 1X cosð/k Þ; /k 2 Fð/Þ n k¼1 ð2Þ where n is the number of particles. The phase /k is drawn from the distribution F(/) using simulations, according to /k ¼ c Z T t¼0 rk ðtÞ  gðtÞdt ð3Þ where c is the gyromagnetic ratio, T is the time at which the MR signal is acquired and g(t) is the magnetic-field gradient [15]. The particle trajectory, rk(t), is simulated as a random walk in a simulation geometry, defined by the tissue model and described by tissue parameters m. In the two-compartment model employed in the present study, m = {d, si, Di, De, ci}, where d is the cell diameter, si is the intracellular exchange time, Di and De are the physical intraand extracellular diffusion coefficients and ci is the intracellular volume fraction (Fig. 1). The extracellular volume fraction (ce) was defined according to ci + ce = 1. For simplicity, D = Di = De was assumed in the present study. In the one-dimensional Stejskal–Tanner pulsed gradient experiment investigated in this study, only the gradient amplitude g was Fig. 1. The two-compartment simulation geometry (tissue model) investigated in the present study is described by the diameter d, the intra- and extracelular diffusion coefficients Di and De, the intracellular volume fraction ci (where ci + ce = 1) and the intracellular exchange time si. In the present investigation, Di = De was assumed and thus four parameters described the simulation geometry. Note that ADCe is estimated in the analytical model, which is related to De according to ADCe = De/k2. The local patch illustrated in the figure was repeated infinitey, creating a hexagonal grid of simulated cells. varied in successive experiments while the direction n was fixed. Assuming rectangular gradient lobes with durations given by d and the time between their leading edges given by D, Eq. (3) is rewritten as /k ¼ cgn  Z d t¼0 rk ðtÞdt  Z Dþd t¼D rk ðtÞdt  ¼ cdgn  Drk ð4Þ where Drk is the difference between the centre-of-masses of the particle’s trajectories during the two gradient lobes. Defining q = c(2p)1dgn = qn and combining Eq. (2) with Eq. (4) gives S ðqÞ ¼ n X k¼1 cosð2pq  Drk Þ ð5Þ Signal values for an arbitrary range of q-values can be obtained from a single simulation where n samples Drk are drawn (simulated) from the F(Dr) distribution, although this leads to correlated simulation noise in the signal-versus-q curve. Discrete simulations were performed with rk discretized in steps of Dx in space and in steps of Dt(k) = Dx2/2ndD(k) in time, where nd is the number of dimensions in the simulation (two or three) and D(k) is the diffusion coefficient for the kth compartment. In each time increment, the particles were moved one step in a random direction. Note that care must be taken to ensure that the simulated gradient durations do not deviate from the specified values when the time is discretized. Particles moving from an intracellular compartment i to the extracellular compartment e were only allowed to complete the transition with a probability p(i, e), related to si according to pði; eÞ ¼ ntot nout   i  DtðiÞ si ð6Þ where ntot is the total number of possible displacements for particles inside compartment i and nout is the number of displacements that would cause a transition out of compartment i during Dt [5]. If the transition was not fulfilled, the particles were elastically reflected back to their original positions, i.e. the membranes were considered to be positioned between the discrete points where the particles were positioned. In terms of membrane permeability (P), p(i, e) is given by pði; eÞ ¼ PDx 2nd V DtðiÞ ¼ DðiÞ Dx A s i ð7Þ when assuming that P = (1/si)(V/A), where V/A is the volume-tosurface fraction of the compartment and nd = 2 [16]. When this assumption is valid, Eq. (7) is approximately equal to Eq. (6). For M. Nilsson et al. / Journal of Magnetic Resonance 206 (2010) 59–67 compartments without specified exchange times, such as the extracellular compartment, p(e, i) was given by pðe; iÞ ¼ pði; eÞ DtðeÞ DtðiÞ ð8Þ Note that this relationship assumes and implies identical particle concentrations in the different compartments. The Monte Carlo formulation in Eq. (2) shows that S* is the average of the stochastic variables cos(/). Therefore V[S*] = V[cos(/)]/n2, with V½cosð/Þ ¼ E½cos2 ð/Þ  E½cosð/Þ2 ! !2 n n 1X 1X 2  cos ð/k Þ  cosð/k Þ n k¼1 n k¼1 Z 2p 0 cos2 ð/Þ 1 1 d/ ¼ 2p 2 ð9Þ ð10Þ In order to avoid too large an influence of the simulation noise on the estimated parameters when using a database of simulated signal curves to fit the measured signal curves, the condition rs  rm pffiffiffiffiffiffiffiffiffiffi ffi should be fulfilled, where rs ¼ V½S   1=2n and rm is the noise standard deviation in the measurements. 2.1. Analytical models The employed analytical model was a two-compartment modified Kärger model, similar to models employed previously in various forms [5,7,10,13,18,19]. This model was defined using a matrix exponential according to   ci SðbÞ ¼ S0 ½ 1 1  expðb  ADC þ K  T D Þ ce ð11Þ where b = (2pq)2TD and TD is the diffusion time TD = D  d/3 (infinitesimally short gradient ramp times assumed). The apparent diffusion coefficient matrix ADC is defined according to ADC ¼  ADC i 0 0 ADC e  ð12Þ with the intracellular ADCi defined as ADC i ¼ ½d  kða; bÞ2 =2T D where k is the tortuosity factor of the extracellular space. The exchange matrix K was defined according to K¼  kie kie kei kei  ð17Þ where kie and kei are the exchange rates out of and into the intracellular compartment with kieci = keice and kie = 1/si. This model utilizes six free parameters, m = {S0, d, si, ADCe, Di, ci}, where S0 is the signal without diffusion encoding and ADCe is the apparent diffusion coefficient in the extracellular space. 3. Method where E[] and V[] are the expectation value and the variance of a stochastic variable, respectively [17]. In the asymptotic case where q ? 1 and S ? 0, /k can be considered to be drawn from a uniform distribution on [0, 2p), where the variance is given by V½cosð/Þ  61 ð13Þ where k(a, b) is calculated using the gaussian phase distribution (GPD) approximation according to The simulation kernel for simulating Drk in Eq. (5) was implemented in C, and other programming was performed in MATLAB (MathWorks, Natick, MA). The simulations employed either a single thread on the central-processing unit (CPU), an Intel Core 2 Duo @ 2.13 GHz, or the GPU, a NVIDIA GeForce 9800 GT graphics card, supporting NVIDIA’s compute unified device architecture (CUDA). The two versions were compared in a basic benchmark test of free diffusion, where the trajectories of n = 107 520 particles were simulated during a simulation time period of 200 ms in steps of D t = 2.5 ls. Random number generation was performed using the Multiply-With-Carry algorithm, in order to obtain high-quality random numbers in the parallel simulations on the GPU [21]. 3.1. Validation of the simulation framework The validity of the simulation framework was initially investigated by three categories of tests investigating free diffusion, restricted diffusion and compartmental exchange. 3.1.1. Free diffusion Simulations of free diffusion were performed where signal-versus-b curves for b-values between 0 and 2500 s/mm2 were generated for free diffusion with d = 10 ms and TD = 100 ms (Dx = Dy = 0.2 lm, D = 2.0 lm2/ms and n = 2  105). The generated signal-versus-b curves was compared with the theoretical signal curves given by S(b) = exp(bD), considering also the signal variance calculated according to Eq. (9). 3.1.2. Restricted diffusion The effects of different step lengths Dx were investigated by simulated measurements performed perpendicular to a cylinder (d = 10 lm) using Dx = Dy = d/60, d/45, d/30 and d/15 lm, with d  0 ms and TD = 500 ms (D = 1.0 lm2/ms, n = 106). Under these conditions, the intracellular signal-versus-q curve is expected to show a diffraction pattern according to the short gradient pulse approximation (SGP), given by sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 X 2aam  2 þ 2 expðaam Þ þ expðbam Þ  ð2  expðaam Þ  expðaam ÞÞ kða; bÞ ¼ a2 a3m ðam  1Þ m¼1 ð14Þ pffiffiffiffiffiffi pffiffiffiffiffiffi with am defined by J 0 ð am Þ ¼ 0, so that am are the roots of the derivative of the Bessel function of the first kind and order one [20] and ð18Þ a¼ 4dDi 2 d ; b¼ 4DDi 2 d ð15Þ The extracellular ADCe relates to De according to ADC e ¼ De =k2 ð16Þ SðqÞ ¼  2 2J1 ðpqdÞ pqd where J1 is the Bessel function of the first kind [12]. To evaluate the effect of different step lengths, the diameter in Eq. (18) was estimated from the simulated signal-versus-q curves. The effects of different diffusion encoding times d on the signal curves were investigated by comparing signal curves from simulated measurements perpendicular to a cylinder, to those expected 62 M. Nilsson et al. / Journal of Magnetic Resonance 206 (2010) 59–67 from the GPD and SGP approximations in Eq. (13) and in Eq. (18), respectively. Two separate simulations were performed with d = 15 ls and d = 15 ms for TD = 500 ms (d = 10 lm, Dx = 0.2 lm, D = 2.0 lm2/ms, n = 106). 3.1.3. Compartmental exchange The relation between si and the probability of a membrane transition p(i, e) in Eq. (6) was tested by simulations in a rectangular confinement using n = 105 intracellular particles for si = 25, 100 and 250 ms. The time that the kth particle spent in the intracellular compartment before its first membrane transition was recorded as (Ti)k and the effective exchange time s^i was estimated by fitting P ni ðtÞ ¼ n expðt=s^i Þ to nk¼1 1fðT i Þk < tg, where 1{a < b} equals unity if a < b and zero otherwise [16]. This process was repeated 100 times and s^i was determined to be significantly longer than si if more than 99% of the estimated s^i exceeded si. This procedure was repeated for si = 100 ms using different step lengths, Dx = d/60, d/45, d/30 and d/15 lm, for diffusion perpendicular to a cylinder (d = 10 lm, D = 2.0 lm2/ms, n = 105) and repeated 20 times to investigate the influence of simulation noise. 3.2. Evaluating diffusion MRI data To investigate the accuracy and precision when estimating microstructural parameters from diffusion MRI data, a database with simulated signal-verus-b curves was constructed using 60 060 different combinations of microstructural tissue parameters m = {d = 2, 3, . . ., 14 lm, D = 1.0, 1.1, . . . , 2.0 lm2/ms, ci = 25%, 27.5%, . . ., 75%, si = 50, 100, . . ., 1000 ms}. The simulations were performed using n = 107 520 particles and Dx = d/40 for TD = 30 and 60 ms, with d = 30 ms. The signal-versus-b curve for each TD was sampled by 41 different b-values linearly spaced between zero and 20 000 s/mm2, corresponding to a maximal gradient strength of approximately 100 mT/m. The microstructure parameters m were estimated for each combination m in the database, after adding noise to the simulated signal curves rm = 0.025, corresponding to SNR = 40 for S(0). Simulated experiments were constructed by adding noise to the qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi simulated signal curves according to Sn ¼ ðS þ N r Þ2 þ N 2i where Sn is the noisy signal, S is the simulated signal and Nr and Ni are gaussian noise with mean zero and standard deviation rm. The model parameters were estimated by fitting the models to the simulated experiment, by finding the signal curves in the database that produced the least squared difference to the simulated experiment, including only signal values exceeding 3rm. For this purpose, a database of analytical signal-versus-b curves was constructed for the same microstructural parameters as in the database of simulated signal curves. Since the simulation model contained only four free parameters, two of the six parameters in the analytical model were fixed to their true values, i.e. S0 and Di, in order to allow fair comparisons of the models. Moreover, specific simulations with impermeable membranes and particles only in the extracellular space were performed to estimate k in Eq. (16), in order to determine the expected value of ADCe. The estimated value of k was obtained as the average value of the simulations performed with TD = 30 and 60 ms. The distributions of the fitted parameters, obtained from fitting the models to 1000 different simulated experiments, were initially investigated for two different sets of microstructural parameters m. The accuracy and precision in the estimated parameters were then investigated by studying the bias and the coefficient of variation (CV) of the estimated parameters. The bias was defined as mean ð^ hÞ=h  1, where h is the true value of the parameter of interest and meanð^ hÞ is the mean value of that parameter estimate. The CV was defined as rð^ hÞ=h, where rð^ hÞ is the standard deviation of the parameter estimate. The bias and the CV were estimated for all combinations of microstructural parameters with D = 1.7 lm2/ms, by fitting the models to 50 different simulated experiments for each combination of microstructural parameters. As discrepancies between the signal curves obtained with the analytical model and the simulation model were found, simulations were carried out using another simulation kernel constructed to fulfill the assumptions of the analytical model. Diffusion in two freely diffusing water pools was simulated, where the ADC values of the two water pools were ADCe = D/k2, with D = 1.7 lm2/ms and k = 1.32 and ADCi as obtained from Eq. (13) with d = 10 lm. The probability of the particles to change water pools was constant at p0(i, e) = Dt/si and p0(e, i) = p0(i,e)ci/ce, with si = 50 ms and n = 2  105 particles were simulated using the same protocol as in the other simulations. Note that these simulation were physically unrealistic. 4. Results 4.1. Validation of the simulation framework The parallel GPU simulations were approximately 26 times faster than simulations on the CPU (238 s versus 9.3 s). Subsequent simulations were therefore performed using the GPU. 4.1.1. Free diffusion The simulated signal-versus-b curve was within the 95% confidence level of the theoretical signal curve for free diffusion (Fig. 2). 4.1.2. Restricted diffusion When Dx was varied between d/60 and d/15 (d = 10 lm), the estimated cylinder diameter deviated less than 0.3% from the true value. Diffraction patterns in the signal-versus-q curves were clearly observed in the simulated signal-versus-q curves when a ? 0 and b = 40 (Fig. 3), expected according to the SGP approximation for measurements performed perpendicular to a cylinder. In the second case, where the SGP approximation was violated (a = 1.2), the simulated signal curve and the curve from the GPD approximation agreed well for a signal attenuation up to approximately one magnitude. At larger signal attenuations, the GPD approximation is not valid (Fig. 3b). 4.1.3. Compartmental exchange The validity of Eq. (6) was investigated for three different intracellular exchange times si = 25, 100 and 250 ms. The Fig. 2. Signal-versus-b curve with 95% confidence interval, obtained from a simulation of free diffusion (D = 2.0 lm2/ms) using n = 2  105 particles. This result indicates a correct implementation of the diffusion encoding. M. Nilsson et al. / Journal of Magnetic Resonance 206 (2010) 59–67 63 Fig. 3. Signal-versus-q curves from a simulated measurement perpendicular to a cylinder (d = 10 lm), with all particles confined inside a cylinder. (a) When a was virtually zero (d = 15 ls), the SGP approximation was valid and the simulated signal curve showed the expected diffraction pattern. (b) For a = 1.2 (d = 15 ms), the simulated signal curve and the curve from the GPD approximation agreed well for a signal attenuation less than approximately one magnitude. estimated values were s^i ¼ 27:0 0:17 ms, 102 ± 0.65 ms and 252 ± 1.6 ms (mean ± one standard deviation), see Fig. 4a. The estimated values were significantly longer than the expected values for si = 25 and 100 ms. As Dx was varied, the estimated values varied only minimally: s^i ¼ 101:46 0:37; 101:25 0:41; 101:41 0:52; 101:72 0:46 ms for Dx = d/15, d/30, d/45 and d/60, respectively. Note that the particle distribution within the compartment became non-uniform as the simulation elapsed (Fig. 4b). 4.2. Evaluating diffusion MRI data The distributions of the estimated model parameters m, obtained using the simulation and the analytical model, are shown in Fig. 5 for two examples with different sets of microstructural parameters. The estimated distributions overlapped to a large extent for si, but differed for most of the other parameters. For example, the analytical model overestimated the diameter in the case where d = 3 lm, in contrast to the simulation model. The difference in D between the analytical and the simulation model is a matter of definitions, since D in the analytical model corresponds to ADCe, while for the simulations, it refers to the diffusion coefficient of all water (D = Di = De). However, the estimated values of ADCe were close their the expected values, as calculated from Eq. (16), with k estimated from a simulation with extracellular particles only in an identical geometry, but with impermeable membranes. For ci, both models produced the close to the correct estimate for the higher ci with longer si, while ci was underestimated in the example with short si when using analytical model. The distributions of the estimated parameters were characterised by their bias and coefficient of variation for all sets of model parameters m with D = 1.7 lm2/ms and ci P 50%. For the analytical model, the reference value of ADCe was set according to Eq. (16) with k estimated from the specific simulations with extracellular particles only. The obtained values of k, when averaged over all included entries in the database, was 1.35 ± 0.06 (mean ± one standard deviation). The bias and the CVs were averaged over all investigated ci and the resulting maps are shown in Figs. 6 and 7, respectively. Using the analytical Kärger model, the bias in d was lower than 10% when d P 4 lm and si P 500 ms, but for si = 100 ms, the bias was low only when d P 8 lm. Exchange times were accurately estimated (bias 610%) for si 6 900 ms and d 6 9 lm, but generally showed large biases otherwise. The intracellular volume fraction ci Fig. 4. (a) The percentage of particles that never passed the membrane, i.e. remained intracellular as the simulation time elapsed. Solid lines indicate the expected curve and the crosses show the results from the simulations. (b) The distribution of particles within the rectangular confinement became non-uniform as the time elapsed. The solid line shows the simulated distribution of intracellular particles and the dashed line the uniform distribution assumed in the Kärger model. 64 M. Nilsson et al. / Journal of Magnetic Resonance 206 (2010) 59–67 Fig. 5. Simulated signal-versus-b curves for TD = 30 ms (circles) and TD = 60 ms (crosses), shown in the leftmost column for two different examples in the top and bottom rows, respectively. The solid horizontal black line corresponds to 3rm, below which no signal values were included in the fitting. The distributions of the estimated parameters si, d, ADCe or D and ci are shown from left to right, estimated using the analytical model (dashed line) and the simulation model (solid line) from 1000 different simulated experiments. The estimated diffusion coefficients represent ADCe in the analytical model and D = Di = De in the simulation model. The true values are indicated by black triangles, but for ADCe it is indicated by a black circle. Note that ci was underestimated when using the analytical model in the second example. Fig. 6. The bias in percent of the true value is shown for the analytical model (a) and the simulation model (b). The x- and y-scales represent d and si, respectively. Only results simulated with Di = 1.7 lm2/ms and the bias values were averaged for ci P 50%. The true values of ADCe in the analytical model were calculated using Eq. (16) with k obtained from dedicated simulations with extracellular particles only and impermeable membranes. was negatively biased for short exchange times (si 6 200 ms) and larger diameters (d P 5 lm) and the bias in ADCe was generally low. Using the simulation model, the biases in the estimated parameters were generally similar to those obtained when using the analytical model. The bias in d was for small diameters negative as opposed to positive for the analytical model, while the bias for si was slightly higher. Notably, the bias in ci found for the analytical model nearly vanished when using the simulation model. The coefficients of variation were nearly identical for the analytical model and the simulation model. The CVs in d were below 10%, when d P 4–5 lm, but up to 50% for smaller d. The lowest CVs, generally below 10%, was found for ci. In si, the CVs were low only for si 6 200 ms, but were otherwise as high as 40%. Signal-versus-b curves from the compartment simulations were compared with the simulations of two freely diffusing water pools with a constant probability of exchange. The obtained results were M. Nilsson et al. / Journal of Magnetic Resonance 206 (2010) 59–67 65 Fig. 7. The coefficient of variation is shown for the analytical model (a) and the simulation model (b). The x-and y-scales represent d and si, respectively. Only results simulated with Di = 1.7 lm2/ms are shown and the values of CV were averaged for ci P 50%. highly different (Fig. 8). For the compartment simulations, where exchange only occurs after a successful membrane transition, a bias was induced in the parameters estimated by the analytical model (Fig. 8a). For the simulations performed with a constant probability of exchange between the water pools, the model returned almost the correct values since those simulations fulfilled the assumptions in the analytical model (Fig. 8b). 5. Discussion Signal-versus-b curves in nerve tissue are known to be nonmonoexponential, i.e. reflecting multiple water pools with different ADC values [19]. These pools are assumed to origin from one extracellular component and one or several intracellular components showing restricted diffusion [5,6,10,11,22]. If the compartments are assumed to be in the slow exchange regime, with negligible exchange occurring during the diffusion time, the bias in the estimated d is small for d P 5 lm, when an measurement protocol optimized for a clinical MRI scanner is employed [6]. Water exchange between compartments might, however, influence high b-value diffusion measurements in fixated nervous tissue as well as in vivo [5,10,23]. The modified Kärger equations employed in Eq. (11) allows for exchange to be included in analytical diffusion models [5,10], but the equations assume a constant probability over time for a compartment exchange [18]. This assumption is not valid for diffusion restricted in a compartment (Fig. 4b). As a result, the estimated parameters are biased when using the analytical measurement model (Fig. 6a). However, the parameters are in some cases biased also when using the simulation model (Fig. 6b), due to measurement noise in combination with limitations in the measurement protocol. For example, prolonging TD would most likely reduce the bias in si, for both models. A resolution limit of approximately 4 lm was observed under slow exchange conditions, similar to the results presented by Fig. 8. The plots show signal-versus-b curves simulated in two different compartments using d = 10 lm, D = 1.7 lm, ci = 60% and si = 50 ms (a) and in a specially designed simulation of water diffusing freely in two different water pools (b). In the specially designed simulation, the probability of exchange between the water pools was constant and the ADC values in the different pools were selected to match those defined by the analytical model for the microstructural parameters investigated in (a). The solid line represents the best fit of the analytical model, where the estimated parameters were (a) ci = 35%, ADCe = 0.88 lm2/ms, d = 10 lm and si = 50 ms for the compartment simulation and (b) ci = 65%, ADCe = 1.1 lm2/ms, d = 11 lm and si = 50 ms for the water pool simulation. This shows that the analytical model underestimates ci from the compartment simulation, while it returns almost the correct value for the (unphysical) simulation specially designed to fulfill the assumptions of the analytical model. 66 M. Nilsson et al. / Journal of Magnetic Resonance 206 (2010) 59–67 Alexander [6]. The cause of the bias in d is the extremely low signal attenuation of the intracellular signal, due to the nearly zero ADCi for d smaller than the resolution limit. An increase in the maximal gradient strength could lower the resolution limit by allowing shorter d. Improving the resolution limit from d1 to d2 = d1/2 could be expressed as the condition b1 ADCi(d1, a1, b1) = b2ADCi(d2, a2, b2), i.e. that the intracellular signal attenuation should be equally large at the smaller d. This condition holds if b2 = b1, a2 = a1 and b2 = b1, which requires d2 = d1/4 as well as (TD)2 = 1/4(TD)1. In order to fulfill b2 = b1, the maximal gradient amplitude must increase eightfold, i.e. g2 = 8g1. This result was described by Lätt et al. in the context of q-space diffusion MRI [24]. The observed bias in ci when using the analytical model under fast exchange conditions (Fig. 6a) might hamper interpretation of diffusion measurements. For example, Clark et al. observed a slow diffusion fraction smaller than expected if that fraction originated from intracellular water and concluded that this mismatch might be caused by a rapid exchange. Indeed, generating biexponential and diffusion time independent signal-versus-b curves using the parameters for gray matter presented by Clark et al. (cfast/slow = 71/29% and ADCfast/slow = 1.02/0.19 lm2/ms) and evaluating data using the simulation model yields ci = 51 ± 7%, d = 12 ± 1.6 lm, si = 61 ± 22 ms, using SNR = 40 and assuming D = 1.7 lm2/ms. The value of ci is still lower than the expected value of 80%, but warrants further investigations. One limitation of the present study was that possible parameter bias in relation to the true tissue-related parameters caused, for example, by an oversimplified tissue model was not evaluated. Moreover, the measurement protocol chosen in the present study represent a high b-value protocol with sequence parameters obtainable with a clinical MRI scanner. The protocol might be optimized for improved accuracy of certain parameters [6]. However, we recommend the presently described methods for detecting biases to be employed before interpreting results obtained with protocols other than the one examined in this study. The simulations and the subsequent data evaluations require a computation time of approximately a week on a desktop computer. Another limitation was that the simulation model used the same model to generate data as for estimating the parameters, in contrast to the analytical model. This could have resulted in slightly higher biases and coefficients of variation for the analytical model, although the negative bias in ci for short si clearly originated from the assumptions in the analytical model (Fig. 8). The results were also influenced by the reduction of two parameters in the analytical model. For example, instead of fixing Di, it could have been included in the fitting with ADCe calculated according to Eq. (16), knowing that De = Di and with k estimated in separate simulations. This was approach was not chosen since any deviation in D would have resulted in a deviation also in d, since these parameters are related via Eq. (13). Simulation noise, although much smaller in amplitude than the applied measurement noise, likely had a marginal influence on results from single simulations (Fig. 5). The bias and coefficient of variations maps, however, were averaged over simulations performed with different ci, resulting in less influence on the results from the simulation noise. At higher SNR of the simulated experiment, the correlated simulation noise would have a larger influence. Fast simulations were required in order to allow a large number of simulations to be performed with different microstructural parameters in a short time. This was achieved by performing parallel simulations on a discrete simulation grid using a modern GPU. Continuous particle displacements, as employed in some studies [13,25,26], may seem more realistic than discrete displacements, but other studies have confirmed that discrete simulations are feasible and of acceptable quality [15,27,31]. Discrete step sizes are also more efficient from a numerical perspective [31]. Simulations with discrete steps in all directions, as employed in the present framework, do not require the handling of multiple membrane reflections, as discussed by Hall and Alexander [31]. However, in order to allow the fractions of intra- and extracellular space to be sufficiently well specified, our design enforced shorter step lengths than required when using continuous displacements. This was a disadvantage of the present simulation framework, since longer step lengths are generally more computationally efficient [31]. The quality tests ensured that the simulation framework generated reliable results. The simulated signal curves in case of free and restricted diffusion were in agreement with the theoretically predicted results, independently of the step length within a reasonable interval (Figs. 2 and 3). However, the estimated effective exchange times were slightly longer than expected for all investigated step sizes (Fig. 4). This small bias was caused by the definition of si in the derivation of Eq. (6), since for high exchange probabilities, the concentration of particles that had never experienced a membrane transition became non-uniform as the time elapsed (Fig. 4b), thereby violating the requirements for Eq. (6) [5]. Nevertheless, this minimal bias can be neglected when si P 50 ms. The present investigation was based on the Stejskal–Tanner pulse sequence, but other gradient waveforms are easily simulated by modifying the assumed gradient waveform, for example, to simulate double-wave vector diffusion MRI [28,29]. Another advantage of using simulations as the measurement model is that it allows the use of complex geometries as the tissue model. Available analytical models for the intracellular signal curves only allows simple geometries, such as planes, cylinders and spheres as the tissue model. Similarly, the extracellular space is commonly modelled as a mono-exponential signal-verus-b curve, although it is plausible that it is non- mono-exponential for intermediate diffusion times [30]. These issues are extraneous when using simulations as the measurement model, since the signal curves are implicitly given by the simulation geometry. In conclusion, signal-versus-b curves simulated for conditions and noise levels achievable on a clinical MRI scanner were evaluated using an analytical two-compartment model as well as by matching them with signal curves produced by Monte Carlo simulations, i.e. the simulation model. Both models produced a bias in the estimated diameter for diameters smaller than approximately 4 lm. A major negative bias was found in the intracellular volume fraction for exchange times shorter than approximately 350 ms when using the analytical model. This bias nearly vanished when using the simulation model. This work could potentially improve the possibility of absolute quantification of tissue microstructural properties using diffusion MRI. The simulation framework and the required source code for this analysis are available on request as open-source software. Acknowledgments The authors are grateful to Prof. Elisabet Englund for helpful discussions. This study was supported by the Swedish Cancer Society (Grant no. CAN 2006/1272), the Swedish Research Council (Grant no. 13514), the Crafoord Foundation Lund and the Lund University Hospital Donation Funds. References [1] Y. Assaf, Y. Cohen, Assignment of the water slow-diffusing component in the central nervous system using q-space diffusion MRS: implications for fiber tract imaging, Magn. Reson. Med. 43 (2) (2000) 191–199. [2] Y. Cohen, Y. Assaf, High b-value q-space analyzed diffusion-weighted MRS and MRI in neuronal tissues – a technical review, NMR Biomed. 15 (7–8) (2002) 516–542. M. Nilsson et al. / Journal of Magnetic Resonance 206 (2010) 59–67 [3] Y. Assaf, T. Blumenfeld-Katzir, Y. Yovel, P.J. Basser, AxCaliber: a method for measuring axon diameter distribution from diffusion MRI, Magn. Reson. Med. 59 (6) (2008) 1347–1354. [4] J. Lätt, M. Nilsson, A. Rydhög, R. Wirestam, F. Ståhlberg, S. Brockstedt, Effects of restricted diffusion in a biological phantom: a q-space diffusion MRI study of asparagus stems at a 3T clinical scanner, Magn. Reson. Mater. Phys. Biol. Med. 20 (4) (2007) 213–222. [5] M. Nilsson, J. Lätt, E. Nordh, R. Wirestam, F. Ståhlberg, S. Brockstedt, On the effects of a varied diffusion time in vivo: is the diffusion in white matter restricted?, Magn Reson. Imaging 27 (2) (2009) 176–187. [6] D.C. Alexander, A general framework for experiment design in diffusion MRI and its application in measuring direct tissue-microstructure features, Magn. Reson. Med. 60 (2) (2008) 439–448. [7] J. Lätt, M. Nilsson, D. van Westen, R. Wirestam, F. Ståhlberg, S. Brockstedt, Diffusion-weighted MRI measurements on stroke patients reveal waterexchange mechanisms in sub-acute ischaemic lesions, NMR Biomed. 22 (6) (2009) 619–628. [8] D. Barazany, P. Basser, Y. Assaf, In-vivo measurement of the axon diameter distribution in the rat’s corpus callosum, in: Proceedings 16th Scientific Meeting, International Society for Magnetic Resonance in Medicine, Toronto, 2008, p. 567. [9] P.N. Sen, Time-dependent diffusion coefficient as a probe of geometry, Concepts Magn. Reson. A 23 (1) (2004) 1–21. [10] G.J. Stanisz, A. Szafer, G.A. Wright, R.M. Henkelman, An analytical model of restricted diffusion in bovine optic nerve, Magn. Reson. Med. 37 (1) (1997) 103–111. [11] Y. Assaf, P.J. Basser, Composite hindered and restricted model of diffusion (charmed) MR imaging of the human brain, Neuroimage 27 (1) (2005) 48–58. [12] W.S. Price, Pulsed-field gradient nuclear magnetic resonance as a tool for studying translational diffusion: Part 1. Basic theory, Concepts Magn. Reson. 9 (5) (1997) 299–336. [13] C. Meier, W. Dreher, D. Leibfritz, Diffusion in compartmental systems i a comparison of an analytical model with simulations, Magn. Reson. Med. 50 (3) (2003) 500–509. [14] J. van den Enden, D. Waddington, H. van Aalst, van Kralingen CG, K. Packer, Rapid determination of water droplet size distributions by PFG-NMR, J. Colloid Interface Sci. 140 (1) (1990) 105–113. [15] P. Linse, O. Söderman, The validity of the short-gradient-pulse approximation in NMR studies of restricted diffusion simulations of molecules diffusing between planes in cylinders and spheres, J. Magn. Reson. A 116 (1) (1995) 77–86. [16] D. Regan, P. Kuchel, Mean residence time of molecules diffusing in a cell bounded by a semi-permeable membrane: Monte Carlo simulations and an [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] 67 expression relating membrane transition probability to permeability, Eur. Biophys. J. 29 (2000) 221–227. B. Lindgren, Statistical Theory, Chapman & Hall, New York, 1993. J. Kärger, Der einfluß der zweibereichdiffusion auf die spinechodämpfung unter berücksichtigung der relaxation bei messungen mit der methode der gepulsten feldgradienten, Anna Physik 482 (1) (1971) 107–109. T. Niendorf, R.M. Dijkhuizen, D.G. Norris, M. van Lookeren Campagne, K. Nicolay, Biexponential diffusion attenuation in various states of brain tissue: implications for diffusion-weighted imaging, Magn. Reson. Med. 36 (6) (1996) 847–857. P. van Gelderen, D. DesPres, P.C. van Zijl, C.T. Moonen, Evaluation of restricted diffusion in cylinders: phosphocreatine in rabbit leg muscle, J. Magn. Reson. B 103 (3) (1994) 255–260. G. Marsaglia, Random number generators, J. Mod. Appl. Stat. Methods 2 (1) (2003) 2–13. A. Szafer, J. Zhong, A.W. Anderson, J.C. Gore, Diffusion-weighted imaging in tissues: theoretical models, NMR Biomed. 8 (7–8) (1995) 289–296. C.A. Clark, D. Le Bihan, Water diffusion compartmentation and anisotropy at high b values in the human brain, Magn. Reson. Med. 44 (6) (2000) 852–859. J. Lätt, M. Nilsson, C. Malmborg, H. Rosquist, R. Wirestam, F. Ståhlberg, D. Topgaard, S. Brockstedt, Accuracy of q-space related parameters in MRI: simulations and phantom measurements, IEEE Trans. Med. Imaging 26 (11) (2007) 1437–1447. S. Peled, New perspectives on the sources of white matter DTI signal, IEEE Trans. Med. Imaging 26 (11) (2007) 1448–1455. G.T. Balls, L.R. Frank, A simulation environment for diffusion weighted MR experiments in complex media, Magn. Reson. Med. 62 (3) (2009) 771–778. B. Balinov, B. Jönsson, P. Linse, O. Söderman, The NMR self-diffusion method applied to restricted diffusion simulation of echo attenuation from molecules in spheres and between planes, J. Magn. Reson. A 104 (1) (1993) 17–25. M.A. Koch, J. Finsterbusch, Compartment size estimation with double wave vector diffusion-weighted imaging, Magn. Reson. Med. 60 (1) (2008) 90–101. M.E. Komlosh, F. Horkay, R.Z. Freidlin, U. Nevo, Y. Assaf, P.J. Basser, Detection of microscopic anisotropy in gray matter and in a novel tissue phantom using double pulsed gradient spin echo MR, J. Magn. Reson. 189 (1) (2007) 38–45. E. Fieremans, Y. De Deene, S. Delputte, M.S. Ozdemir, Y. D’Asseler, J. Vlassenbroeck, K. Deblaere, E. Achten, I. Lemahieu, Simulation and experimental verification of the diffusion in an anisotropic fiber phantom, J. Magn. Reson. 190 (2) (2008) 189–199. M.G. Hall, D.C. Alexander, Convergence and parameter choice for Monte-Carlo simulations of diffusion MRI, IEEE Trans. Med. Imaging 9 (2009) 1354–1364.