Simulation of Intravoxel Incoherent Perfusion Signal Using a Realistic Capillary Network of a Mouse Brain
Valerie Phi van, Franca Schmid, Georg Spinner, Sebastian Kozerke, Christian Federau
SSimulation of Intravoxel Incoherent Perfusion Signal Using a Realistic Capillary Network of a Mouse Brain
Valerie Phi van , Franca Schmid , Georg Spinner , Sebastian Kozerke , Christian Federau University Hospital Zürich, Zürich, Switzerland Institute of Pharmacology and Toxicology, University of Zurich, Zürich, Switzerland Institute of Fluid Dynamics, ETH Zurich, Zurich, Switzerland Institute for Biomedical Engineering, ETH and University of Zürich, Zürich, Switzerland *Corresponding Author:
Christian Federau University and ETH Zürich Institute for Biomedical Engineering Gloriastrasse 35 8092 Zürich Switzerland Email: [email protected] urpose To simulate the intravoxel incoherent perfusion magnetic resonance magnitude signal from the motion of blood particles in three realistic vascular network graphs from a mouse brain. Methods In three networks generated from the cortex of a mouse scanned by two-photon laser microscopy, blood flow in each vessel was simulated using Poiseuille’s law. The trajectories, flow speeds and phases acquired by a fixed number of simulated blood particles during a Stejskal-Tanner monopolar pulse gradient scheme were computed. The resulting magnitude signal as a function of b-value was obtained by integrating all phases and the pseudo-diffusion coefficient D* was estimated by fitting an exponential signal decay. To better understand the anatomical source of the IVIM perfusion signal, the above was repeated by restricting the simulation to various types of vessels. Results The characteristics of the three microvascular networks were respectively: vessel lengths [mean±std. dev.]: 67.2±53.6µm, 59.8±46.2µm, and 64.5±50.9µm; diameters: 6.0±3.5µm, network 2: 5.7±3.6µm, and network 3: 6.1±3.7µm; simulated blood velocity: 0.9±1.7µm/ms, 1.4±2.5µm/ms and 0.7±2.1µm/ms. Exponential fitting of the simulated signal decay as a function of b-value resulted in the following D* [10 -3 mm /s]: 31.7, 40.4 and 33.4. The signal decay for low b-values was the largest in the larger vessels, but the smaller vessels and the capillaries accounted more to the total volume of the networks. Conclusion This simulation improves the theoretical understanding of the IVIM perfusion estimation method by directly linking the MR IVIM perfusion signal to an ultra-high resolution measurement of the microvascular network and a realistic blood flow simulation. NTRODUCTION
Measuring specifically microvascular perfusion is of particular clinical interest, because microvascular dysfunction plays an important role in many important diseases, such as vascular dementia, lacunar infarcts, or diabetes. For example microvascular complications in patients with diabetes are the cause of blindness, renal failure non-traumatic amputations, and a powerful predictors of cardiovascular complications. Microvascular perfusion can be assessed with Intravoxel Incoherent Motion (IVIM) magnetic resonance (MR) perfusion imaging without injection of intravenous contrast agent. The method has seen a regain of research interest in the last decade, both in the brain and in the body, and many recent experimental and clinical studies have demonstrated the feasibility to measure microvascular perfusion with IVIM. For example, IVIM perfusion has been shown to be altered early in patients with diabetes and the dependence on the local microvascular perfusion has been used to assess the quality of collateral blood flow in the context of stroke or to detect vasospasm at the microvascular level after aneurysm rupture. But the microscopic anatomical origin of the Intravoxel Incoherent Motion (IVIM) perfusion signal is not well understood. It is assumed to arise from motion of the blood inside the microvasculature, but an ultimate experimental validation of this assumption is lacking. In addition, the exact relationships between the IVIM perfusion parameters, the micro-vascular network structure, and the blood flow, are not clearly established. The purpose of this study was to simulate blood flow in three realistic microvascular networks, which were imaged at ultrahigh resolution, and to derive the expected IVIM magnetic resonance magnitude signal and from that, the expected IVIM pseudo-diffusion coefficient D*. he microvascular networks, embedded in a tissue volume of 1mm , were obtained by two-photon laser scanning microscopy and contained on average 15’000 vessels. Realistic pressure boundary conditions at all in- and outflows were assigned, and Poiseuille’s law was used to compute flow rate in each vessel. The trajectories of a fixed number of blood particles were then computed for a defined time interval, and used to calculate the phase acquired by each blood proton during standard Stejskal-Tanner monopolar pulse gradient scheme as used in a typical IVIM experiment. The sum of all phases resulted in the simulated MR magnitude signal as a function of the b-value, and the IVIM pseudo-diffusion coefficient D* was obtained by fitting an exponential signal decay. Finally, to better understand the anatomical source of the IVIM effect, we decomposed the signal by restricting the simulation to the various types of vessels: artery, arterioles, capillary bed, venules and veins.
METHODS
Microvascular Networks
The three microvascular networks were acquired by two-photon laser scanning microscopy . Each network is embedded in a tissue volume of approximately 1 mm and contains on average 15’000 vessels. An in-depth description of the pre-processing steps for the blood flow simulations is available here. In brief, the microvascular networks are stored in a graph format, i.e. they consist of a set of nodes (bifurcations), which are connected by edges (vessels). The tortuosity of the vessels is stored in coordinates as edge attributes. Further edge attributes are the effective vessel diameter, the vessel length and the vessel type.
Decomposition into vessel types
The vessel type was assigned by following individual penetrating vessels from the cortical surface and by applying a diameter criterion: If two consecutive vessel branches had a diameter < 7µm, the second branch and the following were considered to be part of the capillary bed. To obtain a capillary diameter distribution consistent with literature data, a histogram-based up-scaling approach was applied, which was based on a beta distribution with a mean of 4 µm and a standard deviation of 1 µm. A three dimensional representation of the microvascular networks under investigation is depicted in
Figure 1 . Blood Flow Simulation
In the following, a brief summary of the key aspects of the blood flow model is provided. A more detailed description is available here. The blood flow model is based on the small Reynolds number in the microvasculature, which allows the use of Poiseuille’s law to compute the flow rate between vertices, i.e. 𝑞 " = 𝜋𝑑 " (𝑝 " − 𝑝 )128𝜇 𝐿 " where 𝑞 " is the flow rate in the vessel connecting node 𝑖 and 𝑗 , 𝑑 " and 𝐿 " are the vessel diameter and length, respectively. The variables 𝑝 " and 𝑝 are the pressure at node 𝑖 and 𝑗 respectively, and 𝜇 is the effective viscosity within the vessel. The effective viscosity is a function of the vessel diameter and the local red blood cell (RBC) density, i.e. the hematocrit, and is computed based on empirical equations. To compute the local hematocrit, we track ndividual RBCs through the microvascular network. Here, we account for Fahraeus-Lindqvist effect, Fahraeus effect and phase separation.
In non-capillaries, the fractional RBC flow at divergent bifurcations is obtained from empirical equations. At capillary bifurcations, the RBCs follow the path of the largest pressure force. By assigning pressure boundary conditions and a physiological inflow hematocrit of 0.3 a well-posed problem is obtained, which together with the continuity equation can be transformed into a linear system of equations. By solving the linear system of equations, the pressure value at each vertex can be computed. Based on the pressure field, the blood flow rate, the flow velocity and the flow direction can be obtained for each vessel. Subsequently, the RBCs are propagated for the time step 𝛥𝑡 and the vessel resistances are updated based on the novel RBC distribution. This procedure is repeated until a statistical steady state is obtained. Due to the fluctuating RBC distribution the flow field is non-steady and hence, time-averaged simulation results are used for further analysis. Particle Tracking
The trajectories of a fixed number of 20’000 particles were computed along the network architecture, using the simulated characteristic flow speed and pressure gradient for a defined time interval, which is the diffusion time of the Stejskal-Tanner monopolar pulse gradient scheme (see below). The particle track was created as follow: 1.)
The particle was placed in a vessel chosen randomly. The starting position inside this vessel itself was chosen randomly as well, according to the relative volume of the vessel segment compared to the volume of the full vessel. .)
The velocity of the particle in this vessel segment was calculated by the flow within the vessel and its diameter. 3.)
The particle followed the tortuosity of the given vessel according to the tortuosity coordinates with the given velocity until it reaches a branching point, at which it entered a new vessel with lower pressure gradient according to the relative flow of all vessels at this branching point. 4.)
The steps 2 and 3 were repeated until the defined time interval was reached. 5.)
If the particle left the network during the defined time interval, the particle was not considered in the signal simulation.
IVIM MR Magnitude Signal Simulation
A Stejskal-Tanner monopolar pulse gradient scheme was used for the simulation, with following parameters : b-values: 0-800 s/mm ; fixed diffusion time: Δ = 100 ms; fixed gradient time: δ = 50 ms. The used gradient strength was calculated for 100 b-values between 0-800 s/mm using following formula: 𝑏 = 𝑒 :; − 1 , with a = 0.0621661 and x = [0,1,2....100] giving a higher b-value density near b=0 s/mm and a lower density for the larger b-values, near b = 800 s/mm . The gradient ramps were neglected. The track of each particle 𝑛 was divided into parts 𝑖 in each vessel segment with constant velocity and constant gradient 𝑔 . The accumulated phase acquired by this particle was 𝜙 C,E = 𝛾 G 𝑟(𝑡) ∙ 𝑔𝑑𝑡 J KLM J N = 𝛾 O P𝑟 Q," ∙ 𝑔 ∙ R𝑡 ",0CS − 𝑡 ",Q
T + 12 𝑣 " ∙ 𝑔 ∙ (𝑡 ",0CSW − 𝑡 ",QW )X " here 𝛾 is the gyromagnetic ratio of a proton, 𝑟 the position, 𝑣 the velocity, and 𝑡 the time. The signal amplitude over all 𝑁 particles was 𝑆 = 1𝑁 [O 𝑒 "\ L C [ The gradient strengths of the Stejskal-Tanner monopolar standard sequence were calculated according to following equation: |𝑔| = √𝑏𝛾Δ ‘ Wa √𝐹 where 𝐹 is 12 for the Stejskal-Tanner sequence and Δ is the diffusion time. IVIM MR Signal Fit
The trace of the MR signal was calculated with an applied gradient in x-, y and z-direction. A mono-exponential fit of the trace of the signal was used to obtain the pseudo-diffusion coefficient D*:
𝑆 = 𝑒 cde ∗ Restriction of the Simulation to the Various Anatomical Vessel Types
To better understand the anatomical source of the IVIM perfusion signal, the above was repeated by restricting the simulation to the various types of vessels: capillary, arterioles, venules and pial vessels. This was done by ignoring all of the particles outside the vessel under consideration, as well as by ignoring all of the particles leaving the vessel type under consideration in the signal simulation. tatistics
Statistical analysis was performed using MATLAB version 2015a (MathWorks, Natick, MA). Statistical significance was assessed using Student’s t-test with a significance level α of 0.05.
RESULTS
Network Characteristics : All three networks differed in their structures. Mean capillary segment lengths of the extracted capillary graphs was for network 1: [mean ± std. dev.] 67.2 ± 53.6 µm, network 2: 59.8 ± 46.2 µm, and network 3: 64.5 ± 50.9 µm. The mean of the diameter distribution was for network 1: 6.0 ± 3.5 µm, network 2: 5.7 ± 3.6 µm, and network 3: 6.1 ± 3.7 µm (
Figure 2 ). Most capillary branches in all three networks were connected to 3 branches, with a maximum of 6 connections (
Figure 3 ). Blood Flow Simulation
The mean simulated blood speed was 0.9 ± 1.7 µm/ms in network 1, 1.4 ± 2.5 µm/ms in network 2 and 0.7 ± 2.1 µm/ms in network 3 (
Figure 2 ). Arteries showed, as expected, significantly higher blood velocities compared to capillaries (12.5 ± 10.2 versus 0.7 ± 1.1 µm/ms, p < 0.05), see
Figure 3 . R Signal Simulation
The IVIM MR signal decay as a function of b-value was faster than mono-exponential for all three capillary networks (Figure 4) . Exponential fitting resulted in the following IVIM pseudo diffusion coefficient D*: 31.7 x 10 -3 mm /s (network 1), 40.4 x 10 -3 mm /s (network 2) and 33.4 x 10 -3 mm /s (network 3). Dependence of the IVIM MR Signal on the Vessel Type
The signal decay was strongly dependent on the type of vessels considered (Figure 5, Table 1) . Slowest signal decay for low b-values was seen in the capillary network (D* range for all three networks: 6.3 - 9.8 x 10 -3 mm /s), a moderate decay for the descending arterioles (21.7 – 40.4 x 10 -3 mm /s) and veins (61.9 – 219.3 x 10 -3 mm /s), and a steep signal decay for the pial artery (172.1 – 833.3 x 10 -3 mm /s). The contribution from the capillaries to the total volume of the network was the largest of all vessel types. DISCUSSION
We simulated blood motion in three realistic microvascular networks obtained by two-photon laser microscopy in the mouse brain and computed the effect of this motion on the MR IVIM perfusion signal for b-values between 0 and 800 s/mm . Our findings are in good agreement with the assumption of a microvascular source of the IVIM perfusion signal, although our simulation derived pseudo-diffusion coefficient D* values were in the upper range compared ith in-vivo measured signal. The pseudo-diffusion coefficient D* has been reported between 7 x 10 -3 mm /s and 17 x 10 -3 mm /s for young healthy adults, but D* of 31 x 10 -3 mm /s have been reported in the gray matter after functional activation, and between 28.49 and 73.96 x 10 -3 mm /s in region of interests drawn in pathological hyperperfused lesions. Interestingly, we found a non-monoexponential behavior of the simulated signal as a function of b-value, with a steep decline of the IVIM signal at very low b-value. This correlates well with our in-vivo experience, and steep declines at low b-value were already reported, see for example Fig. 6 (i) and Fig. 7 (n) in. This could also be linked with the recent observation that a two pool model seems to better describe the IVIM cerebral perfusion in the rat. The question arises why such a steeper decline is not always observed in in-vivo measurements. This might possibly be due to dephasing effects by the imaging gradients or noise causing an underestimation of the real signal at b = 0 s/mm . The largest part of the relative volume of the network consisted of the capillaries. The signal decay was strongly dependent on the specific type of vessel (capillary, arterioles/venules or pial vessels) and we found that an important portion of the IVIM perfusion signal is coming from the arteries and veins inside the voxel. Of importance, only the effects arising from motion inside a specific vessel type was studied here, the effects arising from the blood particles moving between the vessel types were not considered. Further, one should note that the sensitivity with respect to boundary conditions in the blood flow simulation is the highest in the pial vessels. However, while this might affect the absolute velocity values, the impact of the pressure boundary conditions on the velocity ratio between the different vessels types is likely very small. In summary, our findings suggest that the IVIM perfusion signal does arise from all omponents of the microvasculature, not only the capillary bed as suggested by early theoretical assumptions. In conclusion, this simulation improves the theoretical understanding of the IVIM method, by directly linking the MR IVIM signal to ultra-high resolution measurements of the microvascular network and realistic blood flow simulation. The simulated pseudo-diffusion coefficient D* was found in the upper range of corresponding in-vivo measurements. In light of our results, not only the capillary bed, but all anatomical components of the microvascular network contribute to the MR IVIM perfusion signal.
Network 1 Network 2 Network 3
D* Relative Volume D* Relative Volume D* Relative Volume
Whole network
Artery
Vein
Arteriole
Venule
Capillary
Table 1 : Decomposition of the source of IVIM perfusion signal as function of the type of vessel, with D* [10 -3 mm /s] as obtained from a mono-exponential fit for the simulated MR signal decay for b-values (0-800 s/mm ) and the relative volume of the whole network volume [%] of each vessel type. This table shows that an important portion of the IVIM perfusion signal originates from the arteries and veins inside the voxel. Figures Figure 1. Three-dimensional visualization of the three microvascular networks. The vessels are colored based on the vessel diameters. The vessel segments are represented as tortuous tubes with constant diameter. For illustrative purposes the tube radii are enlarged. Figure 2. Velocity, vessel length, vessel diameter and connectivity distributions for all three networks. Figure 3. Capillaries and arteries velocity, vessel length, vessel diameter and connectivity distributions for all three networks. Figure 4. Phase and MR magnitude signal simulation for all three networks. Left : simulated phase distribution for b-values of 10 ( black curve ) and 800 s/mm ( red curve ). Middle : simulated diffusion-weighted magnetic resonance magnitude as function of b-value. The fitted pseudo-diffusion coefficient D* were 31.7 x 10 -3 mm /s for network 1, 40.4 x 10 -3 mm /s mm /s for network 2 and 33.4 x 10 -3 mm /s mm /s for network 3, respectively. Right: simulated diffusion-weighted magnetic resonance magnitude as function of b-value in log-scale, showing the non-exponential behavior of the signal decay at low b-value. The phase is given in radian and the b-value in s/mm . Figure 5. Decomposition of the source of the diffusion-weighted magnitude signal from specific vessel types for network 2: For low b-values, signal decay was the fastest in the artery network, while the signal decay from the arterioles and venules was slower, and the capillary had the slowest signal decay. AKNOWLEDGEMENTS
Christian Federau was supported by the Swiss National Science Foundation (Grant No PZ00P3_173952 and CRSK-3_190697). Franca Schmid received funding from the European Union’s Horizon 2020 Framework Program for Research and Innovation (Specific Grant Agreement No. 720270 [Human Brain Project SGA1] and Specific Grant Agreement No. 785907 [Human Brain Project SGA2]) and the Forschungskredit of the University of Zurich (grant no. FK-19-045).