the effect of diffusion on signal decay in blood oxygenation level - dependent ( bold ) 
 imaging is mainly due to the extravascular contribution of spins , especially at 
 high field strengths > 4  t . in the current study , the well - known monte carlo ( mc ) approach modeling brownian 
 diffusion of protons in a background magnetic field has been used to compute 
 extravascular ( ev ) bold signal changes . to this end , the static dipole model 
 presented previously   has been 
 extended to a dynamic model describing the sampling of phases of the individual 
 protons moving in the inhomogeneous magnetic field . earlier studies on the effect of subvoxel variations in magnetic 
 susceptibility were reported by fisel et al . . weisskoff et al . compared mc simulations with experiments with polystyrene 
 microspheres to demonstrate that enhanced relaxation can be explained 
 quantitatively for both spin - echo and gradient - echo experiments . the effect of an endogenous 
 paramagnetic agent ( deoxygenated hemoglobin ) on image contrast has been 
 addressed by several authors , for example , ogawa et al . all models 
 are based on the fact that in the vicinity of capillaries and venules , local 
 magnetic field distortions are generated by the presence of paramagnetic 
 deoxyhemoglobin in the blood . data from the models in the literature so far have mostly been 
 restricted to a magnetic field strength of 1.5  t , that is , the clinical scanner 
 field strength in the past , and mostly for gre only . however , nowadays , scanners 
 with high or ultra - high field strength for humans up to 9.4  t are available for 
 research , and ev - bold data for these field strengths both for gre and hse have not 
 yet been provided . the aim of the present work , therefore , was to investigate 
 these issues at such high magnetic field strengths . to examine the 
 contribution of extravascular spin in isolation from other factors , an 
 impenetrable vessel wall boundary for extravascular spins the 
 range of investigated susceptibility values was determined using a 
 deoxygenation content of 5% at 1.5  t as the lowest susceptibility value and up 
 to 50% at 9.4  t as the highest value . in addition , standard approaches used in 
 the literature have been evaluated as to how they influence the computed 
 relaxation rates . in particular , the choice of the time step , the diffusion step , 
 the number of angles and the influence of neighboring vessels , and the number of 
 protons were examined . our aim was to study signal decay due to the phase sampling of the 
 individual spins during their random movement . the spins in the brain 
 parenchyma are diffusing in a background magnetic field caused by deoxygenated 
 blood present in capillaries , venules , and veins . the field distribution is 
 therefore determined by vessels inside a computational volume , filled with 
 deoxygenated blood . the susceptibility creating the field distribution around 
 the blood vessel is proportional to the level of blood deoxygenation and to the 
 exterior field . proposed generalizing their results obtained using a 
 numerical model by the use of the bloch - torrey equation . we 
 briefly recapitulate their arguments in the following paragraph and thereby 
 show how their theories relate to one other . because the mc method solution must respect the bloch - torrey equations , 
 generalized scaling laws might be derived in advance to generalize the 
 numerical solutions . if the 
 length scale is made dimensionless by xy=x , we obtain 
 ( 1)ds(y , t)dt = i(y)s(y , t)+(2d)2s(y , t ) . if the time scale is made dimensionless by tt=t/ , we obtain 
 ( 2)ds(x , t)dt = i(x)s(x , t)+(d)2s(x , t ) . 
 substituting  by 1/r and  by te , the general scale independent relation 
 ( 3)dsdt=is+2s 
 with =te and =dte / r2 can be established . from this equation , ( 3 ) , ( 1 ) , and ( 2 ) follow as special cases . the minus sign in ( 3 ) indicates 
 the direction of rotation in the complex plane and therefore can be omitted . in 
 the present article , our model is similar to that of boxerman et al . . a flow chart illustrating the model is shown in figure 1 . the total signal 
 is computed from individual complex transverse magnetizations of the protons , 
 which are accumulating phases in the magnetic background field . the magnetic 
 background field is described either by an analytic formula , or by integrating 
 the individual contributions from the discretized susceptibility distribution . we used the field distribution caused by a vessel segment approximated 
 by a paramagnetic cylinder with infinite length to be the susceptibility 
 distribution for an infinite vessel , which is given by 
 ( 4)f(x)=(1r)2cos(2)sin2( )  ( r1 )  = (cos213 )  ( r1 ) , 
 where  is the angle between the cylinder and the 
 direction of the static magnetic field , r and  represent the nondimensional radial coordinate 
 relative to the vessel radius r and the azimuthal angle in the plane orthogonal 
 to the cylinder , respectively , and  is the susceptibility - induced maximal 
 frequency shift occurring at the vessel surface , given by 
 ( 5)=20hct(1y ) . 
 in ( 5 ) , hct denotes the haematocrit value , that is , the fraction of the 
 volume taken up by the red blood cells ,  is the volume susceptibility difference per 
 unit hct between fully oxygenated and fully deoxygenated blood , 0 is the static magnetic field strength in terms 
 of frequency , and y is the degree of blood oxygenation . our parameter space was given by a haematocrit value of hct = 0.4 , a 
 susceptibility of fully deoxygenated blood of =0.18  ppm ( cgs units ) . frequency thus varied 
 from 1.4454  hz ( corresponding to 1y = 5% at 1.5  t ) to 65  1.4454  hz ( 90.58  hz corresponding to 1y = 50% at 9.4  t ) in steps of 5.7816  hz . the radii varied from 3  m to 60  m in logarithmical 
 equidistant values . alternatively , one could discretize the paramagnetic cylinder 
 representing the vessel into elementary volumes , for example , of size of the 
 computational cell . in the dipole model , to each elementary volume , a dipole 
 strength pm = mv is assigned , where m is 
 the magnetization proportional to the exterior field strength m=b0 . the elementary volume then contributes to the field distribution by an 
 amount of bz=(pm / r3)(3cos21 ) . to compute the total field distortion bz diffusion was computed 
 according to 3 models which can be distinguished by the choice of the diffusion 
 step . take a fixed step ( =2d ) in each direction ( along the positive or 
 negative axis , directions chosen randomly ) . take a fixed step ( =6d ) in a randomly chosen direction . take a random step with a gaussian distribution ( =6d ) in a randomly chosen direction . in the formulae , m/ms , which is 
 typical for the cerebral cortex . as a random generator , the various mc step methods were first tested by computing the radial 
 distribution of the protons after n time steps for m protons . unless stated 
 otherwise , the time step  was chosen as 10  microseconds and te was 40  milliseconds at 
 maximum , and the number of protons was 24 . using these definitions , the routines   given by ogawa et al . the arbitrarily chosen direction must be determined from random values lying 
 inside a sphere ( except the origin ) . to determine the direction , the values have 
 been projected onto the unit sphere 's surface . the time 
 step was chosen as either 10  microseconds or 200  microseconds for the ge and se signal relaxation rate computations . initial positioning and the monte carlo steps were ruled by 
 two separate random generators from press et al . , that is , using the same 
 starting conditions , the results can be reproduced completely . computations 
 have been repeated and mean values have been calculated to reduce numerical 
 noise due to the finite - sized sample of proton spins . the model consists of the following steps : 
 place a cylinder with a given orientation relative 
 to b0 in a volume of interest;distribute protons randomly in that volume of 
 interest;compute the field at the location of each 
 proton as a superposition of the field generated by the cylinder;(optional ) add the gradient fields from the 
 mr sequence at the location of each proton ( e.g. , to simulate a cpmg echo train);advance the phase of each proton according to 
 its local field and , in case of a 180 pulse for hse , invert the direction of 
 phase accumulation;advance each proton in an mc step in an 
 arbitrary direction ( in the case of a 3d step chosen from an arbitrary position 
 on the surface of a unit sphere);if the proton has transgressed the cylinder wall , repeat step ( 6 ) with a given probability , which is determined by the vessel 
 wall porosity ( in our case : repeat step ( 6 ) always);if the proton has left the volume of 
 interest , it might reenter , depending on suitable boundary conditions . place a cylinder with a given orientation relative 
 to b0 in a volume of interest ; distribute protons randomly in that volume of 
 interest ; compute the field at the location of each 
 proton as a superposition of the field generated by the cylinder ; ( optional ) add the gradient fields from the 
 mr sequence at the location of each proton ( e.g. , to simulate a cpmg echo train ) ; advance the phase of each proton according to 
 its local field and , in case of a 180 pulse for hse , invert the direction of 
 phase accumulation ; advance each proton in an mc step in an 
 arbitrary direction ( in the case of a 3d step chosen from an arbitrary position 
 on the surface of a unit sphere ) ; if the proton has transgressed the cylinder wall , repeat step ( 6 ) with a given probability , which is determined by the vessel 
 wall porosity ( in our case : repeat step ( 6 ) always ) ; if the proton has left the volume of 
 interest , it might reenter , depending on suitable boundary conditions . the resulting signal was computed by summating 
 all ( normalized ) complex magnetizations of the ev spins . spins which did not 
 contribute ( spins inside the vessel ) were ignored . in ( 4 ) , the cosine of the angle is computed as a scalar vector 
 product . this , together with expressing the trigonometric functions by a cosine 
 function , prevents the evaluation of trigonometric series , which is highly time consuming . instead of one cylinder , we also used 5 and 9 cylinders as illustrated , 
 for example , by kennan et al . for 5 cylinders . as conditions at the limits of the computational domain , we used 
 periodical boundary conditions ( with the spins re - entering the computational 
 domain from the adjacent side of the domain ) , reflecting boundary conditions 
 ( with the spins being reflected at the limits of the domain ) , and free boundary 
 conditions ( where the spins are free to leave the domain ) . background gradients can be added to the exterior field at any time 
 during the computations , so that diffusion - sensitive sequences like cpmg can be 
 established . the occurrence of vessel direction varies with sin( ) , 
 where  is the angle between the exterior field and 
 the vessel orientation . for our computations , we used 6 , 9 , and 18 angles , 
 equally distributed between 0 and 90. the total signal from vessels with 
 varying orientation thus becomes 
 ( 6)stot=sin()s()sin( ) . the relaxation rate r2 was computed using a two - point evaluation 
 function at te = 16  milliseconds and te = 40  milliseconds assuming stot = s0exp(-ter2 ) . testing of the diffusion process is shown in figures 2and 3 : figure 2 shows the distribution of spins after a 
 fixed diffusion time td together with the theoretical distribution . both the local spin distribution as well as the 
 cpmg signal decay follows the theoretical predictions with some random noise 
 due to the limited number of protons . results for the computation of r2/r2 * for gre and hse imaging are shown 
 in figure 4 . the curves are the mean 
 values of 8 computations , each with a different random series initialization . time 
 step length has been set at 10  microseconds within these simulations . as can be seen , hse 
 relaxation rates are greatest for small radii , whereas for gre large vessel radii 
 have the highest relaxation rates . in addition , as expected , the relaxation 
 rate both for gre and hse increase with blood oxygenation and field strength . higher susceptibility for hse results in shifting the maximal relaxation rate 
 to even smaller radii . comparison of these results to computations with a time step length of 200  microseconds is shown in figure 5 . the 200-microsecond data are slightly noisier since there is no averaging 
 as in the 10-microsecond data , which have been computed 8 times for different random 
 series initialization . comparisons to a gauss3d diffusion step computation and to 
 a free boundary condition computation are also shown ( time step 200  microseconds in both ) . only the data for the free boundary computation are slightly inconsistent with 
 the others . figure 6 compares the results of the 8 
 computations with a computation using 9 or 18 angles instead of 6 ( time step 
 200  microseconds each ) . a nonnegligible deviation of the ge curve , and a slightly 
 smaller deviation also for the se curve can be seen ( in the static dephasing 
 regime : a deviation of 5% for 6 angles , i.e. , a deviation of 2.5% for 9 compared 
 to 18 angles ) . again , data from the 9-angle versus 18-angle comparison have been 
 computed only once . in figure 7 , neighbor rods are taken into account 
 and the result is compared to the result of the 8 computations with different 
 random series initialization . finally , in figure 8 , the data for different 
 numbers of spins in the computational domain can be found . the numbers of spins 
 distributed in the computational domain were 24 , 30 , 36 , 42 , and 48 ( time step 200  microseconds each ) . a change in the number of spins does not result in any significant 
 difference between the calculated relaxation rates . when compared to lower field strengths , the maxima of the spin - echo 
 decay were found to be shifted to lower radii . a set - up with too few angles as well as the use of free boundary conditions introduces errors and therefore 
 should be avoided . all other parameters ( time - step , number of protons , and 
 influence of neighboring vessels of the same orientation , choice of the mc step 
 model used ) did not influence the results in our computations . the current loop approach according to biot savart proved the use of the 
 dipole formulae for a finite voxel element to be correct . the plot of the field distribution of an infinite vessel is in agreement with the field distribution of a very long 
 discretized vessel . however , we found that results for a discretized cylinder 
 slightly deviate from the results for an infinitely long cylinder . for 
 consistency , we show only results obtained for the infinite cylinder in the 
 present work . moreover , the cpmg decay experiment showed that our diffusion 
 modeling works quite well for all mc step methods investigated . , 
 at the limit of infinite numbers of spins and time steps , the signal s( ) in ( 6 ) would be real . however , because we 
 use a finite number of proton spins , there is a residual complex part of the 
 signal . in ( 6 ) , we forced s to be real by taking its absolute value . however , as the number of spins increases , the residual part should diminish 
 and , at the limit of infinite numbers of protons , should vanish . the precursor modeling approaches according to balac et al . , bhagwandien et al . , 
 ldeke et al . , and 
 bakker et al . are all limited in 
 that they rely either on singular analytical solutions for spheres and cylinders , 
 or on complicated procedures such as boundary element methods ( bems ) , finite 
 differences ( fds ) , or finite element methods ( fems ) . the dipole model , in contrast , 
 uses arbitrarily chosen elements ( discretized to a grid ) in combination with 
 simple geometries ( spheres , finite cylinders , parallelepipeds ) . this allows the 
 design of special interventional instruments ( e.g. , a biopsy needle with 
 markers ) and the computation of the decay of the transverse magnetization in  real 
 world  geometry as in a trabecular bone model , or in a vessel network to 
 model signal decay in brain parenchyma . use of analytic field distributions remains 
 possible , thereby allowing modeling similar to boxerman et al . discretization of the susceptibility distribution can be made more 
 precise than the mean distance between protons by introducing distance factors . intra- and extravascular pools were separated by a routine which tests 
 whether a spin is intra- or extravascular . the same routine can be used to 
 detect the transgression of a boundary , and in the case of impenetrable vessel 
 walls , the diffusion step can then be rejected until a transgression - free step has 
 been attained . this is achieved by testing whether a spin was extravascular 
 before the mc step but becomes intravascular after it , or vice versa . actually , 
 vessel walls are to be regarded as impenetrable , because exchange rates are 
 much longer ( typically 500  milliseconds ) than te ( typically 100  milliseconds ) as discussed in fisel 
 et al . , and in boxerman et al . . using a rejection 
 with probability p ( drawn from a random series ) , vessel walls might be modeled 
 as partially penetrable . spins belonging to the intravascular pool are not 
 considered to contribute to the bold signal . the boundaries of the computational domain might be treated as periodical 
 boundary conditions , reflection of spins at the boundaries or by no conditions 
 at all . however , we observed that tests with such penetrable boundaries ( also 
 called  free boundary conditions  ) lead to slightly false results . let us summarize an important , specific finding revealed in figure 4 of section 3 . at high field 
 strength ( 9.4  t ) , spin - echo imaging is especially sensitive at about 3  m 
 radii . the maximum of the spin - echo decay is shifted to lower radii than with 
 lower field strengths . the 6 angles used regularly in our computations essentially yield only an 
 approximate result . however , since the effort to compute with 18 angles is 3 
 times greater , we have restricted our analysis to 6 angles . one should bear in 
 mind that this leads to a ( somehow limited and therefore minor ) error in r2/r2*. moreover , we found that placing vessels in the next 6 or 9 neighbor 
 volumes does not change the result . vessels with orientations different from 
 the orientation of the central vessel in the next neighbor volumes have not 
 been investigated and are left for future work . from 24 protons on , the 
 result does not change significantly with increasing numbers of protons . however , to obtain smoother curves ( less numerical noise due to a limited 
 number of protons ) , computations might be repeated using different random 
 initializations for the mc step procedure . computing time largely depends on 
 the hardware resources available . for the computations with te = 40  milliseconds in steps 
 of 10  microseconds , computing the complete parameter - space for one random series 
 initialization value with 24 protons took approximately 
 8 days on a pentium machine with a clock rate of 1.8  ghz . the model might be further improved by introducing bloch 's equations instead 
 of considering a constant transverse magnetization . pulses and gradients can 
 then be simulated to model a real pulse sequence with all its implications regarding 
 magnetization . the diffusion constant might also be replaced by a diffusion 
 tensor , making the treatment of problems based on diffusion anisotropy 
 possible .
extravascular signal decay rate r2 or r2 as a function of blood oxygenation , geometry , and field strength was calculated using a monte carlo ( mc ) algorithm for a wider parameter 
                             range than hitherto by others . 
 the relaxation rates of gradient - recalled - echo ( gre ) and hahn - spin - echo 
                             ( hse ) imaging in the presence of blood vessels ( ranging from capillaries to veins ) have been computed 
                             for a wide range of field strengths up to 9.4 t and 50% blood deoxygenation . 
 the maximum hse decay 
                             was found to be shifted to lower radii in higher compared to lower field strengths . for gre , however , the 
                             relaxation rate was greatest for large vessels at any field strength . 
 in addition , assessments of computational 
                             reliability have been carried out by investigating the influence of the time step , the monte carlo step procedure , 
                             boundary conditions , the number of angles between the vessel and the exterior field b0 , the influence of 
                             neighboring vessels having the same orientation as the central vessel , and the number of proton spins . 
 the results were compared with those obtained from a field distribution of the vessel computed by an analytic 
                             formula describing the field distribution of an ideal object ( an infinitely long cylinder ) . 
 it was found that the 
                             time step is not critical for values equal to or lower than 200 microseconds . 
 the choice of the mc step procedure 
                             ( three - dimensional gaussian diffusion , constant one- or three - dimensional diffusion step ) also failed to 
                             influence the results significantly ; in contrast , the free boundary conditions , as well as taking too few angles 
                             into account , did introduce errors . 
 next neighbor vessels with the same orientation as the main vessel did 
                             not contribute significantly to signal decay . 
 the total number of particles simulated was also found to play 
                             a minor role in computing r2/ r2.