An Efficient Computational Approach to Characterize DSC-MRI Signals Arising from Three-Dimensional Heterogeneous Tissue Structures

The systematic investigation of susceptibility-induced contrast in MRI is important to better interpret the influence of microvascular and microcellular morphology on DSC-MRI derived perfusion data. Recently, a novel computational approach called the Finite Perturber Method (FPM), which enables the study of susceptibility-induced contrast in MRI arising from arbitrary microvascular morphologies in 3D has been developed. However, the FPM has lower efficiency in simulating water diffusion especially for complex tissues. In this work, an improved computational approach that combines the FPM with a matrix-based finite difference method (FDM), which we call the Finite Perturber the Finite Difference Method (FPFDM), has been developed in order to efficiently investigate the influence of vascular and extravascular morphological features on susceptibility-induced transverse relaxation. The current work provides a framework for better interpreting how DSC-MRI data depend on various phenomena, including contrast agent leakage in cancerous tissues and water diffusion rates. In addition, we illustrate using simulated and micro-CT extracted tissue structures the improved FPFDM along with its potential applications and limitations.


Introduction
The passage of paramagnetic or superparamagnetic contrast agents (CA) through brain tissue induces a measurable drop in T 2or T 2 * -weighted MR signal [1] that forms the basis for dynamic susceptibility contrast (DSC) MRI. When combined with appropriate kinetic models, DSC-MRI can be used to measure hemodynamic parameters quantitatively, such as blood flow, blood volume and mean transit time [2]. This imaging approach relies upon MR signal relaxation enhancement created by CAinduced susceptibility differences between tissue compartments, such as blood vessels and the surrounding extravascular space. The assessment of tumor perfusion parameters using DSC-MRI has proven to be useful for characterizing tumor grade [3][4][5][6][7][8][9] and treatment response [10][11][12][13][14]. Despite its increased use in brain tumor and stroke patients, accurate calculation of perfusion parameters using DSC-MRI relies on two assumptions: 1) a linear relationship, with a spatially uniform rate constant termed the vascular susceptibility calibration factor (k p ), exists between CA concentration and the measured transverse relaxation rate change [15]; and 2) the blood-brain barrier (BBB) is intact, so that contrast agent remains intravascular and can be treated as a nondiffusible tracer [2]. However, heterogeneous distributions of blood vessels within tissue and the dependence of susceptibility field gradients on vascular geometry may yield spatially variant k p values across tissue. Furthermore, leakage of contrast agent in tumors with BBB disruption causes additional T 1 and T 2 * shortening with subsequent distortion of DSC-MRI signal profiles [16][17][18][19][20]. Improved characterization of these potential confounding factors could shed new insights into the biophysical basis of DSC-MRI signals and direct future improvements in acquisition and post-processing strategies.
In order to better understand susceptibility-based image contrast, several theoretical [21][22][23][24][25] and computational models using fixed perturber geometry (e.g., cylinders or spheres) [25][26][27][28][29][30][31][32] have been proposed. To address the limited ability of these computational models to represent the complex vascular morphologies in both normal brain and tumors, Pathak et al introduced the Finite Perturber Method (FPM) for simulating susceptibility-based contrast for arbitrary microvessel geometries [33] and evaluating differences in k p for normal brain and tumor [34]. The FPM uses estimated magnetic field perturbations to calculate MR signal by simulating proton diffusion and phase accumulation using conventional time consuming Monte Carlo methods.
For realistic complex tissues, the MC method needs to track the diffusion of a large number of spins to capture complex structural features, which in turn can increase the computation time. As an alternative, the Bloch-Torrey partial differential equation describing the transverse magnetization can be directly solved using finite difference methods (FDM). This approach has been previously shown to improve the computational efficiency of such simulations [35,36] and used to explore water diffusion in MRI and to aid the interpretation of diffusion-weighted imaging measures and their dependence on the morphology of biological structures such as those found in tumors.
In this study, we propose to evaluate the combination of the finite pertuber and finite difference methods, termed the FPFDM, as a tool for modeling susceptibility based contrast mechanisms. Such an approach leverages the strengths of the FPM, for computing magnetic field perturbations for arbitrarily shaped structures, and the FDM, for efficiently computing the resulting MRI signal evolution. The accuracy of the FPFDM is validated by comparison to traditional Monte Carlo methods. The potential of the FPFDM to compute DSC-MRI signals arising from realistic three-dimensional cellular and vascular models as well as micro-CT based renal angiograms is demonstrated. Going forward, the FPFDM provides a useful tool with which to investigate the influence of vascular morphology, contrast agent kinetics and extravasation on DSC-MRI signals.

Methods
In this section, we first describe a new approach for generating more realistic, three-dimensional tissue structures that can be used for the systematic investigation of DSC-MRI signals arising from heterogeneous tissues. We then describe the computation of the appropriate magnetic field perturbations and the associated MR signal evolution, including the influence of water diffusion, using the FPFDM.

Tissue Structures
Tissue structures consisting of cells, vessels or both were simulated in a 3D volume sampled with N 3 voxels. The motivation for exploring whether cellular properties influence DSC-MRI data originates from previous reports demonstrating that contrast agent extravasation and compartmentalization around cells can induce measurable and dynamic changes in gradient echo acquired signals [16][17][18]37]. While spheres have been used extensively for evaluating susceptibility-based contrast mechanisms they poorly represent in vivo cellular distribution and shape. In particular, packed spheres intrinsically provide no means for modeling orientation heterogeneity and are unable to achieve cellular densities that approximate those found in vivo. To overcome these limitations we explore here the use of randomly packed ellipsoids [38]. Modeling cells as ellipsoids enables the systematic investigation of several features relevant to DSC-MRI including ellipsoid orientation heterogeneity, volume, aspect ratio and higher packing fractions. For completeness, we compare results from randomly distributed spheres, closely-packed spheres on a face centered cubic (FCC) grid and randomly packed ellipsoids.
Typically, randomly oriented cylinders are used to explore susceptibility contrast mechanisms [25][26][27][28][29][30][31][32]. More recently, several groups have employed the use of microvascular angiograms in order to better model in vivo conditions [33,39,40]. In order to provide a framework that mimics in vivo conditions but also enables the systematic investigation of the influence of vascular features on DSC-MRI data we explored the use of fractal tree based vascular networks [41,42]. Starting with an initial cylindrical segment representing an arterial vessel, the vascular tree was created using bifurcation at each junction into smaller daughter segments and a target vascular volume fraction (2%) was used to terminate the fractal tree development. At each junction the diameter of each daughter vessels was calculated using Murray's law [43] and given some degree of randomness along with the branching angles to create tumor-like heterogeneous structures.
All simulated tissue structures were represented by a binary function V(x,y,z) indicating whether a given point within the simulation volume lies inside or outside the simulated compartments, i.e.: To further illustrate the versatility of the FPFDM, in addition to the simulated structures, micro-CT was used to create a threedimensional rendering of a murine kidney vasculature perfused with Microfil (MV-122, Flow Tech). Following perfusion and fixation in 10% neutral buffered formalin, the kidney was scanned in a microCT50 (Scanco Medical AG, Brüttisellen Switzerland). Cross-sectional images of the entire kidney were acquired with an isotropic voxel size of 5.0 mm using an energy of 55 kVp, 200 mA intensity, 700 msec sample time, and 1000 projections per rotation using the manufacturers 1200 mg HA/ccm beam hardening correction algorithm in a 10.24 mm field of view. Using the manufacturer's software, we assembled individual slices into a zstack and contrast-filled vessels were segmented from soft tissue by applying a threshold of 260 mg HA/ccm (determined by calibration against a hydroxyapatite phantom) and a threedimensional Gaussian noise filter with sigma 2.3 and support of 4. The resulting binary three-dimensional reconstruction of the vasculature was then subdivided into MR voxel size sections using in-house Matlab codes (Mathworks, Natick, MA) and used as an input for the FPFDM simulation.

Computation of Magnetic Field Perturbations
The magnetic field perturbations induced by susceptibility variations within the simulated volume were computed using the FPM [33]. To calculate the magnetic field shift at a given point, the FPM breaks the structure into numerous small cubic perturbers and the contribution to the magnetic field shift due to each perturber is calculated independently. The total magnetic field shift is then evaluated as the sum of the magnetic field shifts from all of the perturbers. As computational input the FPM requires a tabular listing of the vascular/cellular structure V(x,y,z), the B 0 field components, and the susceptibility difference (Dx) between tissue compartments. The magnetic field shift computed by the FPM is: where V(x,y,z) is a binary tissue structure, = represents the Fourier transform and DB cube (x,y,z) is the magnetic field change arising from a single finite perturber approximated by the magnetic field shift of an embedded sphere within the cube weighted by a volume factor 6/p, and is described by Eq. [3]: where R is one-half the side length of the cube, Dx is the susceptibility difference between tissue compartments, and r and h indicate the distance from the center of the cube and the angle with respect to B0 of the magnetic field calculation point, respectively. The accuracy of this method has been tested using simple geometries with known theoretical field perturbations [33].

Computation of MR Signal
Estimation of the MR signal relaxation induced by the inhomogeneous field gradients requires simulation of proton diffusion. To track the Brownian motion of thousands of protons over a large number of time steps and calculate their phase accumulation, a Monte Carlo (MC) simulation is frequently used [25][26][27][28][29][30][31]33]. The MC method is potentially time consuming for complex tissue structures because in order to accurately calculate the phase distribution it must track a large number of spins that encounter tissue boundaries during their random walks. An alternative approach is to directly solve the Bloch-Torrey partial differential equation using the FDM [35]. The FDM discretizes the tissue sample to a spatial grid and updates the magnetization at each grid point over a series of time steps. To increase the computational efficiency and eliminate edge effects encountered with traditional FD methods we previously developed a matrixbased FDM with a revised periodic boundary condition [36]. For tissue structures sampled with N 3 voxels the discretized solution of the Bloch-Torrey equation for transverse magnetization (M) using the matrix-based FDM is described by: A is an N 3 6N 3 diffusion transition matrix containing the tissue structural information given in terms of the jump probabilities (probability that a spin starts at one grid point and diffuses to another grid point after a time interval Dt), I is an identity matrix with the same size of A, and : represents element-by-element vector multiplication. The W(t) term is a 16N 3 vector containing the phase accumulation and relaxation for each voxel at each time step and is given by: where c is the proton gyromagnetic ratio (267.5610 6 rad s 21 T 21 ), Dt is the simulation time step, the subscript k indicates a spatial index, DB k (t) is the field perturbation at point k , and T 2,k is the transverse relaxation time at location k within the simulation grid. When a GE sequence is used T 2,k represent the intrinsic tissue T 2 * , and in the case of SE it represents the intrinsic T 2 . In general, the jump probability from simulation grid a to b is described by: where D a?b is the diffusion coefficient if a and b are within the same compartment, Dx a?b is the distance between simulation grid a and b, P m is the permeability of the membrane when a and b are in different compartments, c f is the free concentration of water. The explicit form of the 1D transition matrix can be found in [36]. The MR signal normalized to the initial magnetization M 0 is estimated by summing the magnetizations over all grid points at a particular t and is given by: The associated spin echo and gradient echo change in transverse relaxation rates are calculated at a particular echo time t = TE using: For spin echo imaging, the phase was inverted at t = TE/2. This model is designed to handle cases where the three tissue compartments within a voxel can have different intrinsic transverse relaxation times. By updating T 2,k in equation 5, for each grid point at each simulation time step, the total transverse relaxation, which includes the microscopic and mesoscopic relaxation effects, can also be calculated. The decay of signal from large static perturbers is known not to be exponential (e.g. diffusion in a static linear field gradient) but a simple exponential fit is a good approximation for realistic cases, and other functions can be easily fit. All simulations were performed in the Matlab environment (Mathworks, Natick, MA) running on Intel Core 2 Duo at 2.66 GHz and 4 GB of RAM processors. For clarity, the computational steps involved to obtain the final results are illustrated in Figure 1

Contrast Agent Kinetics
To demonstrate the potential of the FPFDM to simulate DSC-MRI signals arising from the dynamic passage of contrast agent through the vascular and extravascular spaces, such as would occur in brain tumors with a breakdown of the blood brain barrier, we used tissue structures composed of ellipsoids packed around fractal based vascular network. Concentration time curves were sampled using 150 time points for a total of 9 minutes. The arterial input function (AIF) was generated as previously described [44]. The plasma and extravascular extracellular concentration time curves were computed using the pharmacokinetic two compartmental model described by Brix et al [45]. The relevant input physiological, pulse sequence and physical parameters (e.g. blood flow, blood volume, contrast agent transfer coefficient, T 1 , T 2 , etc) were selected from values measured in previous MRI, PET and CT brain tumor studies as previously described [16]. To investigate the influence of extravascular features on DSC-MRI, the signal is computed for two cellular structures with a similar cell volume fraction (,60%) but different ellipsoid radii (5 and 15 mm). The ellipsoids were packed around a fixed vascular tree with a 4% volume fraction.

Validation of FPFDM
For validation, FPFDM and Monte Carlo based MRI signals were computed and compared for models consisting of randomly (6) orientated cylinders and packed spheres. The dependence of gradient-echo (DR 2 * ) and spin-echo (DR 2 ) relaxivity on perturber (vessel) size has previously been characterized using Monte Carlo techniques [27]. To replicate these findings we created 10 different structures composed of approximately 40 randomly distributed cylinders for each vessel radius between 1 mm and 100 mm, each with total cylinder volume fraction equal to 2% of the simulation cube. Using the previously reported simulation parameters [27,33] (susceptibility difference Dx = 10 27 , cylinder volume fraction (V p ) = 2%, B 0 = 1.5T, water diffusion coefficient D = 10 25 cm 2 /s, simulation time step Dt = 0.2 ms, GE TE = 60 ms and SE TE = 100 ms), we computed the vessel size dependence of DR 2 * and DR 2 averaged over all cylinder arrangements. The computed DR 2 * and DR 2 values show negligible changes as the number of averaged structures increases beyond 10. As shown in Fig. 2a, there is excellent agreement between the FPFDM results and those obtained in the Monte Carlo-based comparison studies, which used analytical expressions [27] and FPM [33] for field perturbation calculations.
To compare the computational efficiency of the FPFDM with that of the MC method, we computed DR 2 * values using both techniques. For each technique DR 2 * values were computed for 18 radii using a TE = 60 ms and Dt = 0.2 ms. The computation time for the FPFDM was approximately 140 seconds per structure. Using 1000 randomly distributed spins, the computation time for the MC method was approximately 220 seconds per structure. Table 1 summarizes the simulation parameters used in the MC and FPFDM along with the respective computational times to generate DR 2 * values for 18 cylinder radii. To further validate the accuracy of the FPFDM we also computed DR 2 * for simulated 3D cellular models consisting of packed spheres. Two packing conditions were considered: randomly distributed spheres and sphere packing on FCC gird. For each model, the sphere size was fixed at 9 mm radius corresponding to an approximate pertuber size where the SE relaxivity peaks and the GE relaxivity reaches plateau [27]. The DR 2 * dependence on cell (sphere) volume fraction for the FPFDM was compared to that for the MC method [27] using similar simulation parameters. The MC method was carried out on a different computer system using the approach described previously [27]. Fig. 2b compares the volume fraction dependence of DR 2 * for each of the two sphere packing techniques computed by both  3 . The FPFDM results were obtained by averaging the MR signal for 5 different sphere distributions for each packing and cell volume fraction using a simulation grid size of 128 3 . In contrast, the MC method involves tracking 15,000 random walks for each cell volume fraction and the redistribution of the spheres prior to each random walk. The FPFDM results are in excellent agreement to those produced from the MC methodology.
To investigate the convergence of the FPFDM for randomly distributed structures such as those used above, DR 2 * values obtained from [27] for vessel sizes of 10 mm and 15 mm were compared to the FPFDM results as a function of the number of structures used for averaging. Fig. 2c shows the percentage difference between the MC and FPFDM derived DR 2 * values. For both vessel sizes the computed FPFDM DR 2 * values converge to the corresponding reported values [27,33] to within 7% with only five structure averages. This percentage difference decreases to 0.8% as the number of averaged structures increases to 30.

Modeling the Effects of Contrast Agent Extravastion on DSC-MRI
To demonstrate the potential of the FPFDM for investigating the complex susceptibility effects that occur when contrast agent extravasates and compartmentalizes around cells, we computed the volume fraction dependence of DR 2 * and DR 2 for 3D cellular models consisting of packed spheres or ellipsoids (9 mm radius). To model contrast agent leakage effects relevant contrast agent levels corresponding to a Dx value that is half the peak value of single dose Gd-DTPA injection was assumed. Fig. 3a illustrates the packed ellipsoid model with a representative 2D slice through the computed magnetic field perturbations. Fig. 3b and Fig. 3c demonstrates the influence of packing arrangements on the computed DR 2 * and DR 2 values for Dx = 5610 28 , B 0 = 1.5T,  The highly ordered FCC packing of spheres resulted in the smallest relaxivity, reflecting the more homogeneous magnetic field perturbations and proton phase distributions. Randomly distributed spheres yielded slightly greater relaxivities with a nonlinear relationship with packing fraction. Finally, the packed ellipsoids, which better approximate cell shape in vivo, enable higher random non-overlapping packing fractions (.65%), are less ordered and also yielded a non-linear relationship between relaxivity and cell volume fraction. For all cell volume fractions, the DR 2 * and DR 2 values associated with the ellipsoid-based structures were greater in magnitude than those found with spheres.

Modeling the Effects of Vascular Tree Heterogeneity on DSC-MRI
To illustrate the potential of the FPFDM for modeling the complex geometries of the microvasculature, we used fractal-based branching networks as input to the FPFDM. Fig. 4 illustrates the effect of branching angle heterogeneity (Dh) on the concentration dependence of DR 2 and DR 2 * for typical DSC-MRI contrast agent concentrations. For these simulations we generated three different vascular networks within a 1 mm 3 volume containing 128 3 voxels. Fig. 4a-4c shows sample vascular trees with homogenous rotation (w) angle and bifurcation index (a), which measures the relative diameter of daughter branches at each branching node, with increasing branching heterogeneity (h). The model for normal vasculature is shown in Fig. 4a, with branching angles ranging from 25u-40u. To represent the tortuous and chaotically organized morphology of tumor vessels, the range of branching angle heterogeneity were increased to 25u-80u (Fig. 4b) and 25u-140u (Fig. 4c). Fig. 4d shows three orthogonal 2D slices through the body center of the magnetic field perturbations computed using the FPM for the vascular structure in Fig. 4c. Fig. 4e and 4f plot the concentration dependence of DR 2 * and DR 2 for the three hs  The computed DR 2 * and DR 2 dependence on cell volume fraction and packing arrangement. For all packing arrangements, the relaxivity increases and then decreases with cell volume fraction. Ellipsoid packing yields greater relaxivity than spheres. DR 2 exhibits qualitatively similar behavior to DR 2 * yet with a reduced magnitude. doi:10.1371/journal.pone.0084764.g003 radii ranging from 12 mm to 80 mm with a 2% target vascular volume fraction (v p ). The computed relaxation rates were averaged over five different orientations for each simulated vascular network. Using the slope of the DR 2 * dose response curves, k p values ranging from 100-295 (mM-sec) 21 were obtained. For this low vascular volume fraction, the k p values for the more tumor-like vascular trees (higher Dh) were higher than those in normal trees, up to three fold for this simplified simulation. Similar dependency on branching angle with a reduced susceptibility effect was observed for DR 2 dose response curves.

Modeling First-pass DSC-MRI Data in Brain Tumors
To demonstrate the feasibility of using the FPFDM as a tool to simulate DSC-MRI signals in the presence of contrast agent extravasation, we used two sample tissue structures composed of 60% cells and 4% blood vessels. The two tissue structures were constructed by packing ellipsoids with radii of 5 mm and 15 mm around a fixed fractal-based vascular network. Fig. 5a shows a sample 3D volume rendering of such a tissue structure, which contains ellipsoids of 15 mm average radius and vascular network with vessel size ranging from 5 mm to 45 mm. The simulated vascular (C p ) and extravascular (C e ) contrast agent concentration time curves are shown in (Fig. 5b). Fig. 5c shows a representative 2D slice through the computed magnetic field perturbations at a particular time. The standard deviation of the field perturbation (std DB) for all simulated time points is shown in (Fig. 5d). The simulated C p and C e curves along with model tissue structure, in Fig. 5, were used as an input to compute the dynamic DSC-MRI signal. Fig. 6 shows the GE post-contrast to pre-contrast DSC-MRI signal ratio time curves (S/S 0 ), both in the presence (K Trans = 0.2 min 21 ) and absence (K Trans = 0 min 21 ) of contrast agent extravasation. Fig. 6a-c show the time series for the tissue structure composed of ellipsoids with a 5 mm mean radii, at precontrast longitudinal relaxation times values of T 10 = 500 ms, T 10 = 1000 ms and T 10 = 1500 ms, respectively. The corresponding time series for the tissue structures modeled with 15 mm cellular radii are shown in (Fig. 6d-6f). The following input parameters were used to compute the DSC-MRI signal: B 0 = 3T, D = 1.3610 25 cm 2 /s, Dt = 0.2 ms, TE = 50 ms, TR = 1500 ms, flip angle a = 90u, pre-contrast transverse relaxation time T 20 * = 50 ms. The CA T 1 and T 2 relaxivity values, r 1 and r 2 , were set to 3.9 and 5.3 mM 21 s 21 , respectively [47]. The compartmental membrane water permeability values were set at P m = 0, to model restricted water diffusion. For a fixed cell volume fraction, the simulated time series demonstrate a marked cell size dependency. In general, for both tissue structures, as T 10 increases from 500 ms to 1500 ms the influence of T 1 leakage effects becomes more substantial, as indicated by the increased signal recovery. For the small cell size structure, the T 1 leakage effects eventually result in a signal overshot from baseline (e.g. Fig. 6c). However, the structure constructed with larger cell sizes is dominated by T 2 * leakage effects (as apparent from the low signal recovery well after the CA's first pass) even at T 10 = 1500 ms (Fig. 6f). The simulation time to compute the signal for 150 time points, for 3 T 10 values, 2 contrast agent leakage conditions and 2 tissue structures was approximately 240 mins.

Application to Image-based Tissue Structures
To further illustrate the versatility of the FPFDM, micro-CT images of perfused mouse kidney vasculature (1,242 slices with 142861012 matrix size and 5 mm isotropic voxels) were used as source data for multiple 1 mm 3 vascular models with 200 3 finite cubic perturbers, each 5 mm in size. Fig. 7 shows the entire extracted kidney vascular structure, with sample MR voxel-sized sub-structures and their respective vascular volume fractions. Fig. 8a and 8b shows the FPFDM derived SE and GE k p values obtained from the slope of the DR 2 and DR 2 * dose response curves, respectively. These results are normalized to the vascular fractional volumes and were computed using B 0 = 4.7T, D = 10 25 cm 2 /s, Dt = 0.2 ms, GE TE = 40 ms, SE TE = 80 ms, and a clinically relevant range of Dx values ranging from 0 to 9.4610 28 , corresponding to a Gd-DTPA concentration ranging from 0 to 3.5 mM. In general, the SE and GE k p values were highest for low vascular volume fractions and tended to decrease as the vascular volume fraction increased, with SE and GE k p values ranging from 3.6-27.8 and 53.8-174.3 (mM-sec) 21 , respectively.

Discussion
The FPFDM is a novel efficient computational tool combining features of the FP and FD techniques for calculating susceptibilityinduced relaxivity changes for realistic simulated or imaging-based 3D vascular and cellular geometries that might be observed in vivo. The FPM can compute the induced magnetic fields around arbitrary microvasculature structures without necessitating any assumptions about the underlying vessel geometry [33]. Although the Fast Fourier transform (FFT) improves the computational efficiency of the FPM for computing magnetic field perturbations, the application of MC techniques for tracking proton diffusion through the tissue in order to derive the resulting relaxivity change reduces its computational efficiency. The replacement of the MC component of FPM with matrix-based FDM can increase the computational efficiency by computing the evolution of the discretized magnetization simultaneously [36]. The transition matrix A is either invariant or requires only partial updating for most tissues under consideration, further increasing the computing efficiency of matrix-based FDM that also benefit from optimized MATLAB packages for computations involving large matrices [36].
A Gaussian diffusion kernel convolution can also be used to model CA and water diffusion [39,40,48]. This approach is computationally more efficient than MC approaches, but limited to unrestricted water diffusion. Although non-Gaussian diffusion, a consequence of tissue structure that creates diffusion barriers and compartments, could be modeled by adding a kurtosis term to the kernel, it is not clear how this will affect the slower diffusion process observed in the restricted CA diffusion model [40]. Modeling restricted water diffusion using the MC method [33] or the Gaussian diffusion kernel approach [39,40,48] requires either incorporation of elastic collisions at membrane boundaries or neglecting proton diffusion steps that involve membrane crossing. Unlike the case of unrestricted water diffusion, using these later methods to model restricted water diffusion and/or water diffusion in complex tissue with different compartmental diffusion coefficients will require additional computations, thereby increasing the overall processing time. Given compartmental diffusion coefficients and membrane permeability values, the FPFDM can be used to model restricted water and CA diffusion and water exchange across compartments. For the FPFDM, including these additional structural features requires the computation of multiple versions of the diffusion transition matrix, A. Since A can be determined at the start of the simulation, a library of diffusion transition matrices, for a range of tissue structures, can be established to increase the computing efficiency. For example, computing a dynamic signal for the same structure only requires loading the transition matrix corresponding to the structure once from the library of diffusion transition matrices.
We validated the FPFDM in two ways. First, we replicated the vessel size dependence of DR 2 * and DR 2 (Fig. 2a) using identical simulation parameters to previously described MC and FP techniques [27,33]. Next, we found excellent agreement for relaxivity from packed spheres across a range of packing densities and packing strategies using traditional MC technique versus FPFDM (Fig. 2b). The agreement between MC and FPFDM converges as the number of structures included in the average for the FPFDM increases (Fig. 2c). Unlike MC simulation, which tracks a large number of particles in the simulation or, equivalently, runs the same simulation many times to obtain an accurate average result, the FPFDM converges to the average result with only a few simulation runs. This behavior can be explained in the following way. In the MC simulation, a population of particles distributes in the whole system and the particles that encounter membranes within the complex tissue are only a small portion of all the particles such that the echo signal does not contain sufficient enough information about the tissue features that restrict diffusion. Hence, to solve this problem, more particles are considered in the simulation or, equivalently, the same simulation is run many times to obtain an accurate average result. In contrast, the FDFDM determines the diffusion transition matrix at the start of the simulation, which already contains the tissue structural information and results in a faster convergence of the average signal.
For a simplistic structure containing randomly oriented cylinders with a total of 18 different radii, the FPFDM, as compared to MC, reduced the computation time to calculate DR 2 * values from 220 s to 140 s. For complex tissue structures, and under conditions of restricted water diffusion, the increase in computational efficiency afforded by the FPFDM will improve even further. In such cases, the MC method requires a larger number of spins and additional computation steps in order to converge and capture sufficient information about the tissue structure [49]. In contrast, for these more complex structures, the FPFDM does not require additional computing time and is not limited by restricted water diffusion [36].
The FPFDM has the potential for modeling nonstandard geometries that may better simulate cells and microvasculature in vivo. We computed relaxivities for simulated 3D cellular models consisting of packed spheres and ellipsoids (Fig. 3), and found greater relaxivity for packed ellipsoids over all volume fractions compared to the sphere packing. This suggests that the additional degree of freedom in spatial orientation for ellipsoids increases field perturbation heterogeneity. The greater orientational heterogeneity and packing density afforded by ellipsoids compared to spheres would seem to make ellipsoids better suited for simulating susceptibility-based contrast mechanisms of cellular structures.
Although simulations in this study are based on a simple 2compartment model, at the expense of computational time, the same approach to model water diffusion and exchange can be used to model CA diffusion and transport across compartments, yielding a more realistic heterogeneous CA distribution within a voxel. This can be achieved by updating CA concentration for each voxel at each simulation time step using a CA diffusion transition matrix (A CA ) calculated using appropriate CA diffusion coefficients and permeability across membranes from literature [50,51]. For the purposes of this study, we assumed that the contrast agent equilibrates within each compartment over the time it takes to acquire each DSC-MRI image (1-2 seconds). Such an assumption is traditionally employed and consistent with current DSC-MRI analysis techniques.
We also computed the contrast agent concentration dependence of transverse relaxation rates for vascular trees. Traditionally, randomly oriented cylinders were used to investigate the influence of vascular properties (e.g. vessel size, vessel volume fraction) on relaxation rates. Fractal-based vascular trees better approximate the microvascular network in vivo, but this complex geometry with variable vessel rotation, size distribution, branching angles, and diameters of daughter vessels is very difficult to model and require high resolution to achieve structural details. Our results demonstrate the feasibility of using FPFDM for complex geometries, and suggest that although the generally accepted linear relationship between relaxation rate and CA concentration is reasonable, the proportionality constant k p depends upon the microvascular geometry, a finding that is consistent with previous studies [27,33,52]. The higher relaxation rate for vascular structures with greater range of branching angles is most likely due to the greater heterogeneity of vessel branch orientation with respect to B 0 and their larger space occupancy which impacts the frequency offset of a larger volume compared to narrow branching angles that pack vessels in a small region. These preliminary computational results show marked k p heterogeneity across vascular networks, suggesting that further work is needed to better characterize the influence of vascular heterogeneity on DSC-MRI derived perfusion parameters in brain tumors. More systematic studies are underway for wide range of morphological, physiological and pulse sequence parameters to investigate these findings at the higher vascular volume fractions more typically encountered in DSC-MRI experiments, such as in grey matter or brain tumors.
The systematic evaluation of fractal-based vascular trees using the FPFDM could shed new insights into the relationship between DSC-MRI relaxation rate and vascular geometry. Furthermore, the use of fractal trees enables the application of well-established flow models [53][54][55] such that contrast agent kinetics and the associated DSC-MRI time series can be considered for each vascular structure. This would enable a more rigorous investigation of DSC-MRI-based voxel-wise measures of vessel size, transit time and flow distributions and oxygen extraction. Realistic 3D vascular and flow models could then be expanded to incorporate the extravasation of contrast agent and its subsequent diffusion around cells in the extravascular space. Such expansions would create a powerful framework with which to investigate DSC-MRI and susceptibility-based imaging methodologies in brain, tumors and other organs of the body.
The FPFDM also provides a computationally reasonable approach for simulating DSC-MRI derived transverse relaxation rates both in the presence and absence of CA extravasation, and restricted water diffusion induced by membrane permeability (Fig. 5 and Fig. 6). The results shown in Fig. 6 demonstrate that contrast agent leakage can lead competing T 1 and T 2 * effects as the CA traverses the extravascular extracellular space. For a given T 10 , the structure with smaller sized cells exhibited higher signal intensity recoveries as compared to that with larger sized cells. The compartmentalization of CA around the larger cells creates stronger magnetic field perturbations and greater relaxation rate changes (T 2 * effects). In general, as T 10 increases, T 1 leakage effects will be more pronounced and may dominate any T 2 * leakage effects, as is the case for the smaller-sized cells. In such cases, the characteristic signal overshoot may be observed (Fig. 6c). For the tissue structure with larger perturber sizes, the signal intensity exhibits less recovery due to the presence of substantial T 2 * leakage effects ( Fig. 6d-6f). As shown in Fig. 3b, cell density may also influence the shape of DSC-MRI signals, with the magnitude of T 2 * leakage effects decreasing (and T 1 leakage effects increasing) as the cell density increases. Consequently, DSC-MRI data from tumors with tightly packed, smaller-sized cells would likely present with pronounced T 1 leakage effects (e.g. signal overshoot). Given the clinical importance of DSC-MRI signal recovery characteristics to help differentiate among tumor types [37,56], a systematic in silico study of DSC-MRI signal recovery and its dependence on physiological, pulse sequence and physical parameters is currently under investigation.
Prior studies have shown the potential and value of incorporating image-based vascular structures into susceptibility simulations [39]. Similar to these previous studies we sought to demonstrate the versatility of the FPFDM by determining the dose-response of relaxation rates for vascular structures derived from ex vivo micro-CT scans of perfused kidney vasculature. The dose-response curves from MRI voxel-sized regions of the kidney vasculature were used to determine the distribution of vascular susceptibility calibration factors, k p , within the kidney. For vascular volume fractions up to 30%, k p values were very heterogeneous (Fig. 8), with decreased heterogeneity for vascular volume fractions greater than 5%. The k p decreased over vascular volume fractions between 0 and 5% with a slower decrease above 5%, consistent with a previous study in rodent brain that found grey matter k p to be nearly twice that of tumor [34]. It should be noted that the kidney microvascular structure presented in this study is limited by the spatial resolution of the micro-CT data. With a 5 mm, resolution individual capillaries could not be resolved and capillary dense regions, such as in the glomeruli, present as a single large perturber. The differentiation and inclusion of these capillaries will likely influence the overall k p heterogeneity across voxels for both SE and GE computations. For the purposes of this study, this example illustrates the ability of the FPFDM to explore susceptibility contrast in tissue structures acquired using ex vivo imaging modalities. As the FPFDM only requires that structures consist of a digital format it could accept structural input from any imaging modality (e.g. optical, CT, electron microscopy, MRI).
One of the limitations of the FPM is the use of FFT to calculate the spatial convolution of the vascular structure with the finite perturber magnetic field perturbation. As demonstrated in [33] the resulting field perturbation is equivalent to the field perturbation from a periodic array of the tissue structure under consideration. Although realistic tissue structures extend beyond the boundary of the simulation space, which introduces a ''boundary problem'', we used zero-padding of the tissue structure to avoid additional field perturbation at the boundaries from the neighboring array. The padding size to eliminate boundary field effects depends on the perturber size and the tissue structure. Here we used a zero-pad size of one-tenth of the simulation box, since the field perturbation changes we observed by using higher zero-pad sizes were negligible.
The FPM was designed to compute the magnetic field changes from a single finite perturber convolved with a digitized tissue structure array, and hence this approach cannot be used for arbitrary magnetic susceptibility distributions. While methods capable of computing arbitrary susceptibility distributions are more comprehensive and should be explored [57,58], it is typically assumed that contrast agent instantaneously distributes within each tissue compartment (e.g. intravascular and extravascular extracellular space) at each imaging time point. Accordingly, the FPFDM is a practical approach to compute field perturbations  Fig. 7. SE k p values ranged from 3.6-27.8 (mM-sec) 21 , and GE k p values ranged from 53.8-174.3 (mM-sec) 21 . Above 5% volume fraction, the GE k p values were relatively constant with a mean value of 103.3(mM-sec) 21 . doi:10.1371/journal.pone.0084764.g008 An Efficient Method to Model DSC-MRI Signals arising from tissue structure with only a few susceptibility compartments, such as the intravascular, intracellular, and extravascular extracellular spaces.
The sampling of tissue structures at higher resolution increases the computational accuracy of the FPM but it comes at the expense of computational time. Such increases in resolution would also add to computational time needed to compute the MR signal using the FDM. This is particularly true if a need arises to reduce the simulation time step (Dt) due to increased resolution or decreased perturber size (Dx), in order to satisfy the constraint that the jump probability (see Eq. 6) should be #1/6. This is because when the number of spins leaving a given node exceeds the number that was present, the FDM becomes unstable [35]. With the parallel high-performance computing techniques we previously developed [36], we are exploring ways to increase the computational efficiency of the FPFDM at higher resolutions so that we can more accurately characterize fine tissue microstructure across a broader range of structural dimensions (e.g. a few microns up to a hundreds of microns).

Conclusion
The FPFDM is an alternative computational tool for efficiently modeling susceptibility induced MR signal relaxation from complex perturber geometries. In general, the proposed FPFDM could be used to investigate the influence of realistic tissue microstructure on any susceptibility based contrast mechanism such as vessel size imaging, BOLD contrast, single cell imaging, and quantitative susceptibility mapping. Currently, the proposed method is being utilized to assess the influence of geometrical, morphological and physiological parameters of microvessels and cells on susceptibility induced MR relaxation rate changes. Such studies should shed new insights into DSC-MRI contrast mechanisms and enable the systematic evaluation of how acquisition and analysis methods influence the measurement of reliable perfusion parameters in brain and tumor tissue.