Theory and Generalization of Monte Carlo Models of the BOLD Signal Source John Martindale, Aneurin J. Kennerley, * David Johnston, Ying Zheng, and John E. Mayhew The dependency of the blood oxygenation level dependent (BOLD) signal on underlying hemodynamics is not well under- stood. Building a forward biophysical model of this relationship is important for the quantitative estimation of the hemodynamic changes and neural activity underlying functional magnetic res- onance imaging (fMRI) signals. We have developed a general model of the BOLD signal which can model both intra- and extravascular signals for an arbitrary tissue model across a wide range of imaging parameters. The model of the BOLD signal was instantiated as a look-up-table (LuT), and was veri- fied against concurrent fMRI and optical imaging measure- ments of activation induced hemodynamics. Magn Reson Med 59:607– 618, 2008. © 2008 Wiley-Liss, Inc. Key words: Monte Carlo simulations; biophysical model; BOLD; hemodynamics; fMRI The technique of blood oxygenation level dependent (BOLD) functional magnetic resonance imaging (fMRI) has had a great impact on our understanding of the brain in the 15 years or so since its inception, and yet the exact depen- dency of the BOLD signal on underlying hemodynamics is not well understood. Such an understanding is critical if we wish to undertake quantitative estimation of the hemo- dynamics and neural activity underlying fMRI signals. The aim of this work is to provide and test a general modeling scheme as a first step toward formalizing such an under- standing. This modeling scheme is designed to be as real- istic and general as possible. This work introduces the theory involved in constructing Monte Carlo models of the BOLD signal source, and describes the generation of a look-up-table (LuT) that can be used to predict the BOLD signal. We use concurrent BOLD and optical imaging mea- surements of activation-induced hemodynamics to test the validity and accuracy of the model. The microscopic physics of BOLD signal generation is well understood: paramagnetic deoxyhemoglobin in the blood vessels causes local field inhomogeneities, which in turn causes diffusing water protons to sustain phase shifts relative to one another and thus induces signal decay. Two approaches have been taken in attempting to understand this effect. Analytic modeling attempts to derive expres- sions for signal decay mathematically (1–3), and has pro- vided various models for asymptotic cases (such as the “static dephasing regime”). Monte Carlo methods use com- puter simulations of proton phase shifts during diffusion to calculate signal decay (4 – 6), and have provided quali- tative understanding of signal decay and validation of several assumptions. These models have certainly ad- vanced our understanding of the BOLD effect, but have not resulted in a generally applicable model: signal decay in the intravascular (IV) space is generally ignored or charac- terized as a pure phase shift; models are often limited to a specific field strength and/or echo time; models are often valid only for limiting cases. The aim of this work is to provide characterization of signal decay in both the IV and extravascular (EV) spaces across all commonly used field strengths and to give specific methods for combining mod- eling results into tissue models of arbitrary complexity. The physics underlying Monte Carlo models of the BOLD signal is simple: as protons diffuse through a mag- netically inhomogeneous medium, their magnetic spins acquire phase differences with respect to one other, caus- ing destructive interference and thus a characteristic mac- roscopic decay. This decay is, in general, complex; a point often overlooked by models which seek to characterize only the magnitude of the signal. The choice of model for generating field perturbations is all-important in these Monte Carlo models. We used three different models. The first is the “standard” model, in which vessels are treated as infinite cylinders with homogeneous susceptibility and field perturbations are calculated through analytic equa- tions (e.g., Boxerman et al. (4)). We used this model for the EV space. The second model treats red blood cells as spherical perturbers with homogeneous susceptibility, as in Boxerman et al. (7). This model is used to provide estimates of IV signal decay. The third model treats arbi- trary shapes as a collection of magnetic dipoles and uses the insight that the nonlocal magnetic dipole equation becomes local in the Fourier domain (8) to generate a three-dimensional (3D) discrete grid of field perturbations, following the method of Marques and Bowtell (9). This model gave an additional estimate of signal decay to verify our results, and allowed visualization of the EV field per- turbations predicted by the infinite cylinder and spherical perturber models. In any model of the BOLD signal, it is important to be clear which aspects of signal generation one is trying to model. The observed transverse relaxation rate is com- monly denoted as R * 2 and can be defined as a combination of intrinsic relaxation (R 2 , due to spin-spin interactions) and relaxation due to local field inhomogeneities and dif- fusion (R' 2 ). This can be expressed as, R 2 * = R 2 + R' 2 [1] Centre for Signal Processing in Neuroimaging and Systems Neuroscience (SPiNSN), Department of Psychology, University of Sheffield, Sheffield, UK. *Correspondence to: Dr. Aneurin James Kennerley, SPiNSN, University of Sheffield, Western Bank, Sheffield S10 2TP, UK. E-mail: A.J.Kennerley@shef.ac.uk Received 13 June 2007; revised 11 October 2007; accepted 20 November 2007. DOI 10.1002/mrm.21512 Published online 25 January 2008 in Wiley InterScience (www.interscience. wiley.com). Magnetic Resonance in Medicine 59:607– 618 (2008) © 2008 Wiley-Liss, Inc. 607