Abstract
Diffusion Magnetic Resonance Imaging (dMRI) sensitises the MRI signal to spin motion. This includes Brownian diffusion, but also flow across intricate networks of capillaries. This effect, the intra-voxel incoherent motion (IVIM), enables microvasculature characterisation with dMRI, through metrics such as the vascular signal fraction fV or Apparent Diffusion Coefficient (ADC) D*. The IVIM metrics, while sensitive to perfusion, are in general protocol-dependent, and their interpretation can change depending on the flow regime spins experience during the dMRI measurements (e.g., diffusive vs ballistic), which is in general not known — facts that hamper their clinical utility. Innovative vascular dMRI models are needed to enable the in vivo calculation of biologically meaningful markers of capillary flow. These could have relevant applications in cancer, for instance assessing responses to anti-angiogenic therapies targeting tumor vessels. This paper tackles this need by introducing SpinFlowSim, an open-source simulator of dMRI signals arising from blood flow within pipe networks. SpinFlowSim, tailored for the laminar flow patterns in capillaries, enables the synthesis of highly-realistic microvascular dMRI signals, given networks reconstructed from histology. We showcase the simulator by generating synthetic signals for 15 networks, reconstructed from liver biopsies, and containing cancerous and non-cancerous tissue. Signals exhibit com-plex, non-mono-exponential behaviours, pointing towards the co-existence of different flow regimes within the same network, and diffusion time dependence. We also demonstrate the potential utility of SpinFlowSim by devising a strategy for microvascular property mapping informed by the synthetic signals, focussing on the quantification of blood velocity distribution moments, and of an apparent network branching index. These were estimated in silico and in vivo, in healthy volunteers and in 13 cancer patients, scanned at 1.5T. In conclusion, realistic flow simulations, as those enabled by SpinFlowSim, may play a key role in the development of the next-generation of dMRI methods for microvascular mapping, with immediate applications in oncology.
1 Introduction
In diffusion Magnetic Resonance Imaging (dMRI), water proton motion is encoded in the acquired signals through magnetic field gradients [Kiselev, 2017]. Diffusion encoding provides sensitivity not only to Brownian motion due to pure diffusion, but also to pseudo-diffusion effects arising from the incoherent flow of blood protons through intricate capillary networks [Le Bihan et al., 1986]. Flow through sets of pseudorandomly distributed capillaries leads to magnitude dMRI signal attenuation, a phenomenon known as Intra-Voxel Incoherent Motion (IVIM) effect. IVIM enables the in vivo characterisation of microvascular perfusion through dMRI [Le Bihan, 2019], relevant in a variety of diseases, as, for example, in cancer [Fokkinga et al., 2023]. Cancers feature aberrant microvasculature, whose flow patterns can differ considerably from normal tissues [Munn, 2003]. Tumour vasculature is targeted specifically by anti-angiogenic treatments, which are being used in several cancers (e.g., in liver or kidney carcinomas [Jayson et al., 2016]) and tested in combination with therapies such as immune check-point inhibitors, with promising results [Huinen et al., 2021]. The non-invasive assessment of vascular properties through dMRI can equip physicians with new tools for tumour characterisation and longitudinal assessment. It is thereby an active field of research, with studies spanning from malignancy detection to treatment response assessment [Iima et al., 2018, Perucho et al., 2021].
IVIM methods typically rely on disentangling vascular from extra-vascular tissue dMRI signals [Barbieri et al., 2016b, Barbieri et al., 2016a]. Multi-exponential models are routinely used for this purpose, providing metrics such as the vascular signal fraction fν, or the pseudo-diffusion (vascular) apparent diffusion coefficient (ADC) D*. Both fν and D* are useful indices, as they have shown value in cancer assessment [Dappa et al., 2017]. However, these metrics have limitations, since they entangle several, different characteristics of the microvascular scale into a single number, e.g., the product between the means of the blood velocity and capillary length distributions in the diffusive flow regime [Le Bihan and Turner, 1992]. More-over, they do not account for higher-order cumulants of the diffusion decay (e.g., kurtosis terms proportional to b2), and their actual numerical value can depend on the acquisition protocol in non-trivial ways [Wu and Zhang, 2019]. In practice, this makes routine IVIM metrics semi-quantitative, surrogate parameters, a fact that, together with their known high variability [Barbieri et al., 2020], hampers their practical clinical deployment.
Recently, the numerical simulation of dMRI signals within histologically-realistic voxel models is being increasingly used to inform parameter estimation [Nilsson et al., 2010, Nguyen et al., 2014, Fieremans and Lee, 2018, Buizza et al., 2021, Morelli et al., 2023]. Simulation-informed approaches increase the realism of signal models, and may thus improve the biological fidelity of dMRI parametric maps [Nedjati-Gilani et al., 2017, Palombo et al., 2019]. However, up to date dMRI simulations have been dominated by Monte Carlo Brownian random walks [Hall and Alexander, 2009, Ginsburger et al., 2019, Rafael-Patino et al., 2020, Lee et al., 2021]. Given that only a few simulation frameworks have focussed on blood flow [Van et al., 2021, Weine et al., 2024], there is an urgent need for new, histologically-meaningful, and reproducible simulation frameworks tailored for dMRI signal arising from blood flow. These could be used to inform novel numerical approaches for non-invasive microvasculature mapping based on dMRI, which could equip oncologists with biologically-meaningful vascular markers in clinical settings. The new dMRI methods could enable the characterisation of capillary flow patterns that are not captured by classical IVIM fν and D*, e.g., informing on anisotropic flow patterns, higher-order cumulants or diffusion-time dependence of the vascular signal.
With this article we aim to fill this scientific gap. We present an open-source framework for blood flow simulation within vascular networks, referred to as SpinFlowSim from here on, and demonstrate its potential to inform microvasculature property estimation in dMRI. We start by illustrating the computational engine behind SpinFlowSim, based on pipe network theory. Afterwards, we describe the synthesis of dMRI signals arising from flow within realistic vascular networks obtained from histological images of human tumours. Finally, we showcase a potential application of SpinFlowSim, by using the synthetic signals to inform microvasculature property estimation, which is demonstrated in silico and in vivo, in healthy volunteers and in cancer patients.
2 Methods
In this section we illustrate the computational engine upon which SpinFlowSim relies. Afterwards, we present the histological data used to generate realistic vascular networks, and then describe how synthetic dMRI signals were used to inform microvasculature parameter estimation in silico and in vivo. SpinFlowSim is made freely available at https://github.com/radiomicsgroup/SpinFlowSim.
2.1 Simulation framework
In SpinFlowSim we aim to reconstruct the distribution of volumetric flow rate (VFR) across the different segments of an input vascular network. The following characteristics of the vascular network are specified directly by the user:
a list of capillary segments with their radii;
the 3D coordinates of the extremities of each segment, referred to as nodes;
the inlet/outlet of the whole network;
the input VFR qin.
To obtain the VFR distribution, we solve a linear inverse problem, in which the pressure drop Δpk,n across each pair of connected nodes (k, n) is proportional to the VFR qk,n between k and n through a flow resistance coefficient Rk,n, via
The approach, valid for the laminar flow regime in micro-capillaries, has been recently proposed for capillary flow simulations [Schmid et al., 2015, Van et al., 2021].
To solve for all unknown qk,n in Eq. 1, we rely on PySpice [Salvaire, 2023], a python package for electric circuit analysis, given that our flow problem is formally equivalent to solving a passive electric circuit (electric-hydraulic analogy). Note that in a passive electric circuit, the voltage drop across a resistor is proportional to the product of the electric current through the resistor and the resistance of the element itself, i.e., it is formally equivalent to Eq. 1. In this first demonstration of SpinFlowSim, we compute the resistance between nodes k and n through a modified Hagen-Poiseuille law, as done in [Blinder et al., 2013]:
Eq. 2 models the effect of the hematocrit as well as erythrocyte-capillary wall interactions [Pries and Secomb, 2008, Blinder et al., 2013]. Above, μ is the dynamic viscosity of pure plasma [Késmárky et al., 2008] (μ=1.20 mPa s at 37ºC), rk,n is the radius of the capillary segment, and Lk,n its length.
After recovering the VFR qk,n between each pair of connected nodes, in SpinFlowSim we obtain the corresponding mean velocity
Finally, the 3D trajectory pw(t) of the generic w-th blood spin is synthesised by integrating the discrete-time system
given an initial position pw(0) = pw,0. In Eq. 4, Δt is the temporal resolution of the simulation, while νw(t) and nw(t) are the instantaneous velocity vector magnitude and direction experienced by the spin at time t. Spins’ initial positions pw,0 are seeded across the whole network, with uniform spin density in each segment. The numbers of spins assigned to each segment is proportional to its volume [Van et al., 2021]. During the integration of Eq. 4, spins reaching the termination of a capillary are assigned at random to one of the emanating branches. The probability of a spin being assigned to a specific branch is proportional to the VFR through that branch [Van et al., 2021]. More formally, once a flowing spin reaches the k-th node, the probability of it continuing its trajectory in the k → n branch emanating from k is 
. Moreover, spins reaching the network outlet continue flowing through a shifted copy of the vascular network, whose inlet position coincides exactly with the outlet itself. This ensures that no spins are lost during the simulation (periodic boundary condition). SpinFlowSim supports the visualisation of spin trajectories as a video, in order to facilitate the visual inspection of the simulation output.
Once the trajectories for W spins have been generated, we synthesise a complex-valued dMRI signal s for any input gradient wave form G(t) as [Fieremans and Lee, 2018]
given the requested total simulation duration T.
2.2 Vascular networks
We deployed SpinFlowSim on vascular networks reconstructed from 2D histological images. These consisted of biopsies obtained in patients suffering from advanced tumours and participating in an ongoing imaging study at the Vall d’Hebron Institute of Oncology (Barcelona).
The biopsied tissue, taken from liver tumours, was processed and stained. Digitised images of the stained tissue were acquired on a Hamamatsu C9600-12 optical slide scanner (resolution: 0.454 μm). For this study, we used 11 histological images, obtained from 11 patients. For each patient, we had access to either a routine hematoxylin-eosin (HE) stain (n=9) or a CD31 stain (n=2).
We drew a total of 15 2D networks. We drew networks manually, by tracing visible capillaries in non-cancerous liver parenchyma or in cancerous regions-of-interest (ROIs). Networks were drawn on approximately square ROIs, of sizes ranging from 250 to 550 μm per side. Networks were made of interconnected segments, with curved capillaries approximated by a piece-wise series of straight pipes. A characteristic radius was assigned to each segment by averaging three radius measurements, performed on at the inlet, middle point, and outlet level. For each network, we computed an approximated network size as the maximum euclidean distance between any pair of nodes. We also computed the mean and standard deviation of the capillary radii and lengths (r and L).
We generated 100 VFR distribution realisations by changing randomly the position of the network inlet/outlet 10 times, and varying the input VFR qin for each inlet/outlet pair (10 uniformly-spaced qin values in [1.5⋅10−4; 5.5⋅10−3] mm3/s), to cover plausible blood capillary velocities [Ivanov et al., 1981]. The total duration and the temporal resolution of the simulations were T = 100 ms and Δt = 0.01 ms. We characterised each realisation by computing the mean (νm) and standard deviation (νs) of the velocity distribution across capillary segments, as well an Apparent Network Branching (ANB index). ANB measures the average number of segments spins travel through during the simulation. Spearman’s correlation coefficients among all possible pairs of mean νm, νs, ANB, r and L were computed.
Finally, we synthesised illustrative dMRI signals for routine pulsed-gradient spin echo (PGSE) sequences. We probed b-values in the range [0; 150] s/mm2, and varied the gradient separation Δ in Δ = {30, 50, 70} ms, while fixing the gradient duration to δ = 20 ms. Signals were generated for two orthogonal directions within the plane containing the 2D networks, and their magnitude averaged. Signals were characterised through the computation of the vascular ADC D* via linear fitting (minimum/maximum b of 0/100 s/mm2, with 10 s/mm2 increments). We investigated whether microvascular properties are encoded in the diffusion-weighted (DW) signal by scattering D* against each of νm, νs and ANB in turn, at fixed gradient separation Δ. Linear fitting between D* and each microvascular property was performed and Spearman’s correlation coefficients were computed, in order to corroborate qualitative trends on scatter plots.
2.3 Microvascular property estimation from dMRI
We also investigated whether the synthetic signals generated with SpinFlowSim can be used to inform microvascular parameter estimation in dMRI. We hypothesised that, for a given dMRI protocol, large dictionaries of synthetic, noise-free signal arrays S = {s1, …, sM }, coupled with their corresponding vascular parameter arrays P = {p1, …, pM }, can be used to find practical numerical implementations of the forward signal model P ↦ S(P). Numerical implementations of this type could be easily incorporated in standard non-linear least square (NNLS) fitting, used routinely in dMRI, thus avoiding the need for approximated analytical signal expressions.
In the following sections, we will describe in silico analyses performed to investigate the feasibility of simulation-informed fitting. We will then describe experiments performed to demonstrate the approach in vivo, based on the acquisition of dMRI scans in healthy volunteers and cancer patients at 1.5T.
2.3.1 In silico estimation
We used SpinFlowSim to synthesise signals for 3 realistic dMRI protocols, and then analysed such signals to test whether it is possible to estimate νm, νs and ANB from noisy measurements. One of the protocols represents a rich, comprehensive PGSE acquisition, encompassing several b-values in a measurement regime with high sensitivity to IVIM effects (i.e., b smaller than approximately 100 s/mm2 [Le Bihan, 2019]), as well as multiple diffusion times. A second protocol is instead a shorter subset of the rich protocol. Finally, the third protocol matches a DW twice-refocussed spin echo (TRSE) used for in vivo imaging. Signals were generated for two orthogonal directions within the plane containing the 2D networks, and their magnitude averaged. In summary, the protocols were:
a rich PGSE protocol, referred to as “richPGSE”. It consisted of a total of 99 measurements, consisting of 9 b = 0 and 10 non-zero b-values b = {10, 20, 30, 40, 50, 60, 70, 80, 90, 100} s/mm2, each acquired for 9 unique diffusion times, corresponding to (δ, Δ) = {10, 20, 30} ms × {30, 50, 70} ms.
A second PGSE protocol, referred to simply as “PGSE”. It is a subset of the former, and describes a more realistic acquisition that could be implemented under time pressure. It encompassed 3 b = 0 and 6 diffusion-weighted (DW) measurements, namely b = {50, 100} for 3 different diffusion times. The gradient duration δ was fixed to 20 ms, while the 3 diffusion times were achieved by varying Δ as Δ = {30, 50, 70} ms.
a DW twice-refocussed spin echo (TRSE) protocol, referred to simply as “TRSE”. It matches the protocol implemented on a 1.5T Siemens Avanto system in vivo (see section 2.3.2 below). It consisted of 3 non-DW and 6 DW measurements. These were b = {50, 100}, acquired for 3 diffusion times. The gradient duration of the 4 gradient lobes (Supplementary Fig. 1) for the 3 diffusion times were: δ1 = {8.9, 13.2, 18.9} ms, δ2 = {17.6, 19.3, 21.0} ms, δ3 = {20.4, 24.8, 30.5} ms, δ4 = {6.0, 7.7, 9.5} ms. The separation between the onset of the gradient lobes (Supplementary Fig. 1) were instead: Δ1,2 = {17.4, 21.7, 27.5} ms, Δ1,4 = {63.9, 74.2, 87.5} ms.
Briefly, we performed a leave-one-out experiment. For each vascular network in turn, we used noise-free signals from 14 out of 15 substrates to learn the forward signal model (νm,νs, ANB) ↦ s(νm, νs, ANB), which we then used for estimating νm, νs and ANB on noisy signals for the remaining 15-th network (signal-to-noise ratio (SNR) at b = 0 of 20). The forward signal model (νm, νs, ANB) ↦ s(νm, νs, ANB) was learnt by interpolating the set of paired examples signals/vascular parameters with a radial basis function (RBF) regressor, so that fitting could be performed by embedding s(νm, νs, ANB) into standard maximum-likelihood NNLS routines [Panagiotaki et al., 2012]. Fitting was performed with the freely-available mri2micro_dictml.py tool, part of the bodymri-tools python repository (https://github.com/fragrussu/bodymritools). To characterise fitting performance, we generated scatter plots between ground truth and estimated νm, νs and ANB, and computed corresponding Spearman’s correlation coefficients.
2.3.2 In vivo estimation
We also investigated the feasibility of using synthetic signals from SpinFlowSim to inform microvascular property estimation in vivo, on both healthy volunteers and cancer patients. All participants were scanned after providing informed written consent, in imaging sessions approved by the Clinical Research Ethics Committee (CEIm) of the Vall d’Hebron University Hospital of Barcelona, Spain (codes: PR(AG)29/2020 and PR(IDI)109/2022)).
Healthy volunteers: data and analysis
We scanned two healthy volunteers in their thirties on a 1.5T Siemens Avanto system. The acquisition included routine anatomical imaging and a DW TRSE Echo Planar Imaging (EPI) scan, with salient parameters: resolution of 1.9 × 1.9 × mm2; slice thickness of 6 mm; b = {0, 50, 100, 400, 900, 1200, 1600} s/mm2, with each b acquired at 3 different diffusion times, with the same diffusion times as the “TRSE” protocol used simulations (see section 2.3.1 above); TE = {93, 105, 120} ms for the short, intermediate, and long diffusion time; TR = 7900 ms; trace DW imaging; NEX = 2; GRAPPA = 2; 6/8 Partial Fourier imaging; BW = 1430 Hz/pixel; acquisition of a b = 0 image at the shortest TE with reversed phase encoding.
We denoised scans with MP-PCA [Veraart et al., 2016], mitigated Gibbs ringing [Kellner et al., 2016] and corrected for motion and EPI distortions [Andersson et al., 2003]. Subsequently, we normalised the signal acquired at each TE to the b = 0 signal level at the same TE, and then estimated the vascular signal SV for b ≤ 100 s/mm2 in each voxel [Gurney-Champion et al., 2018, Wang et al., 2021], as
Above, S is the measured signal and ST is an estimate of the extra-vascular tissue signal. ST was computed by extrapolating to b ≤ 100 s/mm2 an ADC fit 
performed on signal measurements at b > 100 s/mm2.
Afterwards, we fitted (νm, νs, ANB) voxel-by-voxel, using the same fitting procedure employed in in silico experiments above, but learning the forward model P ↦ S(P) on all 1500 synthetic signals from all vascular networks. For reference, we also computed more standard IVIM metrics fV and D*, by fitting 
to the vascular signal estimated at the shortest TE, with 
. Mean and standard er-rors of νm, νs, ANB, fV and D* within manually drawn ROIs were computed. The ROIs were placed in the liver, spleen, as well as medulla and cortex of a kidney.
Cancer patients: data and analysis
Finally, we tested our simulation-informed parameter estimation on dMRI scans of 13 patients suffering from advanced solid tumours (7 females, 5 males; approximate age range 30-85 y.o.), who participated in an ongoing imaging study at the Vall d’Hebron Institute of Oncology (Barcelona, Spain). Patients were scanned on the same 1.5T Siemens Avanto system used to acquire data on healthy volunteers, and according to the same imaging protocol. dMRI scans underwent the same processing as described above, obtaining voxel-wise maps of νm, νs, ANB, fV and D*. Mean and standard deviation of such metrics within tumours were obtained, with tumours manually segmented by an expert radiologist (R.P.L).
3 Results
3.1 Vascular networks
Fig. 2 illustrates the 15 vascular networks generated in this study from liver tumour biopsies. Out of the total, 3 were segmented on non-cancerous liver parenchyma, while the remaining 12 on cancerous tissue. The 3 non-cancerous networks were drawn on liver tissue found on the histological slide, adjacent to tumour tissue (n=2 melanoma metastases; n=1 ovarian cancer metastasis). The 12 networks drawn on cancerous tissue came from primary liver hepatocellular carcinoma (HCC, n=5), or from liver metastases of colorectal cancer (CRC, n=5), endometrial cancer (n=1), and melanoma (n=1).
Outline of the proposed SpinFlowSim framework. The dashed red box indicates user-provided input information. An illustrative example of a network segmented on a biopsy with resolved volumetric flow rates for an input flow qin = 3.1 ⋅ 10−3 mm3/s and synthesised signals are shown.
Vascular networks segmented on digitised liver tumor biopsies (resolution: 0.454 μm). Each network is labelled as “Non-Cancerous” or “Cancerous”, depending on whether it was drawn on non-cancerous liver parenchyma or on tumour tissue. For the latter case, the primary cancer is also indicated (CRC stands for Colorectal Cancer, while HCC for Hepatocellular Carcinoma). The non-cancerous networks were drawn on liver tissue found on liver tumour biopsies of patients suffering from Melanoma (n=2) and Ovarian cancer (n=1).
Table 1 reports salient statistics of the vascular networks shown in Fig. 2, related to capillary radii, length, velocity distribution, and number of vascular segments sensed by flowing spins. The table shows that different network morphologies lead to different blood velocity distributions. For example, mean νm across VFR realisation can vary from as low as approximately 4 mm/s up to 25 mm/s. This range of variation is mirrored in the average number of capillaries blood travels through during the simulation (ANB metric), which varies from just over 12 up to almost 57 segments. Supplementary Fig. 2 shows distributions of νm, νs and ANB for all networks, across the 10 different inlet/outlet realisations and given an illustrative input VFR qin = 3.1⋅10−3 mm3/s. Distributions are skewed, and strong contrasts in terms of νm, νs and ANB are seen across networks (e.g., compare Net 3 with Net 4). Supplementary Fig. 3 shows correlation coefficients among all possible pairs of metrics among νm, νs, ANB, as well as mean segment length L and capillary radius r. There is a strong, positive correlation between νm and νs and ANB (0.89, 0.82), and a moderate positive correlation between νs and ANB (0.55). All of νm, νs and ANB are negatively correlated with L and r (strongest correlations between ANB and r, of -0.93; weakest for νs and L, of -0.19). Finally, L and r are positively correlated between each other (correlation of 0.68).
Fig. 3 shows examples of VFR and blood velocity fields reconstructed in two vascular networks with SpinFlowSim, along-side dMRI signals. The two networks were segmented on non-cancerous liver parenchyma of a patient suffering of melanoma (top panel, Net 6) and on metastatic CRC (bottom panel, Net 12). The figure highlights that distributions of VFRs and velocities arise across network segments, owing to their different resistance to flow. The segments with the highest VFRs do not necessarily feature the highest velocities, due to differences in terms of segment diameter. The VFR distributions lead to fast dMRI signal attenuation in both networks, with most of the signal decayed by b = 150 s/mm2. The signal decay is not mono-exponential (note the log-scale in the y-axis). Oscillatory patterns are also seen as well as some diffusion time dependence, with the dMRI signal decreasing slightly with increasing Δ at fixed b.
Examples of resolved vascular networks. The top row shows results from a vascular network segmented on a HE-stained non-cancerous liver region, found on a biopsy of a patient with metastatic melanoma (Net 6). The bottom panel shows results from a CD31-stained rectal cancer area (Net 12). From left to right, we show the vascular network, the resolved blood flow velocity field for qin = 3.1 ⋅ 10−3 mm3/s, and examples of dMRI signal decay over a range of b-values (0-150 s/mm2) and diffusion times (Δ = {30, 50, 70} ms, δ = 20 ms).
Fig. 4 reports on the relationship between D* and microvascular properties νm, νs and ANB. D* increases with increasing νm, νs and ANB. There is a strong positive correlation between D* and νm and between D* and ANB (e.g., Spearman’s r = 0.85 at Δ = 30 ms for both νm and ANB), while the correlation is moderate between D* and νs (r = 0.43 at Δ = 30 ms for νs). The dependence of D* on νm, νs and ANB is modulated by the diffusion gradient separation Δ, albeit slightly. For all three metrics, D* increases with increasing Δ.
Scatter plots relating D* to microvascular properties for different gradient separation Δ and fixed gradient duration δ = 20 ms. From left to right: relationship between D* and metric νm ; D* and metric νs ; D* and metric ANB. Spearman’s correlation coefficients between D* and each microvascular property are also reported in each plot, for any fixed Δ.
3.2 Microvascular property estimation from dMRI
3.2.1 In silico estimation
Fig. 5 reports results from in silico estimation of νm, νs and ANB from noisy vascular signals, synthesised according to protocols “TRSE”, “PGSE” and “richPGSE”. There is a moderate correlation between ground truth and estimated νm, νs and ANB values for protocols “PGSE” and “PGSE” (minimum correlation: 0.45 for νs for protocol “TRSE”; maximum correlation of 0.65 for ANB for protocol “PGSE”). Correlation is instead strong for protocol “richPGSE”, e.g., up to 0.79 for metric ANB. As an example, Supplementary Fig. 4 illustrates the complete set of synthetic signals generated for protocol “TRSE” across the 15 segmented networks. The figure highlights that the signal decay spans several order of magnitude: variations in the microvascular characteristics of the networks lead to remarkably different vascular dMRI signals.
Scatter plots relating estimated and ground truth microvascular parameters in computer simulations. First row (panels (a), (b), (c)): results for protocol “TRSE”. Second row (panels (d), (e), (f)): results for protocol “PGSE”. Third row (panels (g), (h), (i)): results for protocol “richPGSE”. From left to right: results for metric νm (panels (a), (d), (g)); results for metric νs (panels (b), (e), (h)); results for metric ANB (panels (c), (f), (i)). Spearman’s correlation coefficients between estimated and ground truth values are also reported in each plot.
3.2.2 In vivo estimation
Fig. 6 shows IVIM metrics fV and D* alongside νm, νs and ANB in the liver and spleen of healthy volunteer 1. On visual inspection, fV and D* are systematically higher in the liver than in the spleen. This contrast is mirrored by νm, νs and ANB, which are as well higher in the former organ than in the latter. Fig. 7 reports instead mean and standard errors of all metrics within several ROIs (liver, kidney medulla and cortex, and spleen), and in both healthy volunteers. Inter-organ differences are seen, as for example higher D*, νm, νs and ANB in the liver, compared to the spleen. Trends of inter-subject differences are also seen. E.g., in healthy volunteer 1, higher values of all of fV, D*, νm, νs and ANB in the kidney medulla than in the kidney cortex are seen. However, in healthy volunteer 2, D* is higher in the cortex than in the medulla, and differences between medulla and cortex among all other metrics are less marked. Intra-/inter-scanner variability is seen in all of fV, D*, νm, νs and ANB.
Microvascular maps in a representative healthy volunteer. (a): b = 0 image; (b) and (c): IVIM maps fV and D*; (d), (e) and (f): microvascular indices νm, νs and ANB. In the b = 0 image, we highlight the location of the liver and the spleen.
Bar plots reporting mean and standard error of the mean of all microvascular metrics in the different regions-of-interest (ROIs) of both healthy volunteers. (a): trends for metric fV ; (b): trends for metric D*; (c): trends for metric νm; (d): trends for metric νs; trends for metric ANB.
Fig. 8 shows representative microvascular maps in cancer. The figure refers to a colorectal cancer adrenal gland metastasis. Both IVIM metrics fV and D* as well as microvascular νm, νs and ANB show intra-tumour contrasts. For example, we observe a core of lower fV and D* as compared to the outer ring of the tumour. This spatial trend is mirrored by metrics νm, νs and ANB: the lower νm, νs and ANB point towards slower and less variable blood velocity in the core of the tumour, and predict blood to travel through fewer vessel segments, as compared to the outer ring. Overall, νm, νs and ANB exhibit similar contrasts among each other, but certain differences are also seen e.g., voxels with high νm that do not necessarily show the highest ANB values.
: Parametric maps obtained in an adrenal gland metastasis of a 61 y.o. male patient, suffering from advanced rectal cancer (patient 8, scanned on a 1.5T Siemens Avanto system with a DW TRSE sequence). Top row: b = 400 s/mm2 image and IVIM maps fV and D*. Bottom row: microvascular parameters obtained via simulation-informed model fitting, namely: mean blood velocity νs, blood velocity standard deviation νs, and Apparent Network Branching ANB.
Table 2 reports mean and standard deviation of all vascular metrics within tumours. The metrics highlight inter-tumour differences in terms of vascularisation, as seen in dMRI. For example, breast cancer metastases feature the highest vascular signal fraction fV among all tumours. Conversely, the highest D* is seen in a lung cancer adrenal gland metastasis (patient 11), which also features the highest νm, νs and ANB across the whole cohort. The lowest D* is instead seen in liver metastasis of rectal cancer (patient 8), a trend that is mirrored by νm and ANB, which in this case are the lowest across all tumours.
Summary of patients’ demographics and key clinical data (primary cancer type, location of the imaged tumours, patients’ sex and approximate age, expressed in 5-years ranges) and trends of microvascular metrics fV, D*, νm, νs and ANB within the manually-segmented tumours (mean and standard deviation). For sex, F indicated female, while M male. Routine IVIM metrics fV and D* were obtained at fixed diffusion time, on the dMRI scan with the shortest TE.
Supplementary Fig. 5 reports Spearman’s correlation coefficients between all possible pairs of vascular metrics, as obtained across the 13 cancer patients. IVIM D* is significantly, positively correlated with all of νm, νs and ANB (rs = 0.55, p = 0.049 with νm; rs = 0.57, p = 0.044 with νs; rs = 0.64, p = 0.019 with ANB). No significant correlations are instead seen between fV and any of νm, νs and ANB. D* and fV are weakly, negatively correlated between each other (rs = -0.25, p = 0.394).
4 Discussion
4.1 Summary and key findings
This work presents SpinFlowSim, an open-source simulator of blood flow based on pipe network analysis. The simulation framework, tailored for the laminar flow regime at the micro-capillary level, enables the synthesis of DW signals for any desired input dMRI gradient waverform. We demonstrate Spin-FlowSim on 15 microvascular networks, reconstructed from biopsies in a variety of liver cancers and in non-cancerous liver parenchyma. These allowed us to simulate micro-perfusion IVIM signals for realistic dMRI protocols, in the low b regime. The signals exhibit complex, non-mono-exponential behaviour, pointing towards the co-existence of spin pools experiencing different flow regimes. Key microvascular properties, such as moments of the blood velocity distribution the local network branching sensed by flowing spins, are shown to be encoded in the dMRI signal decay. This fact enabled the design of a strategy for the practical estimation of these properties from any input dMRI signal measurements, given simulations of the corresponding acquisition protocol. We showcase the approach in silico and on in vivo scans of healthy volunteers and cancer patients, obtaining patterns of microvascular metrics that are plau-sible with the known anatomy and cancer pathophysiology.
4.2 Simulation framework
Our simulator relies on a well-established computational approach for laminar flow characterisation in capillaries. This links the pressure drop across a capillary to the VFR passing through it, via a flow resistance proportionality factor [Schmid et al., 2015, Van et al., 2021]. In this study, as a first proof-of-concept, we borrowed an empirical expression for this resistance from [Blinder et al., 2013], and used the freely available PySpice [Salvaire, 2023] package to convert the VFR estimation problem into the analysis of a passive electric circuit. Our strategy, computationally efficient, retrieves the VFR distribution across all segments of a vascular network. These are used to estimate the mean blood velocity in each capillary and, finally, the trajectories of flowing spins, by numerical integration of the velocity field over time. By superimposing arbitrary dMRI gradient wave forms to spins flowing in networks reconstructed from histology, our framework enables the synthesis of realistic IVIM signals, without making assumptions on the specific flow regime in which measurements take place (e.g., diffusive/ballistic [Kennan et al., 1994, Scott et al., 2021]). Our approach offers a practical way to characterise the salient characteristics of micro-capillary perfusion, and its relationship to dMRI. It may therefore play a key role in the development of innovative dMRI methods for vascular characterisation with un-precedented specificity to physiology, urgently needed for non-invasive cancer characterisation.
4.3 Vascular networks
We studied HE and CD31-stained histological images from liver tumour biopsies, obtained from cancer patients suffering from advanced solid tumours. From these data, we manually reconstructed 15 2D vascular networks, on which we simulated blood flow by varying the input VFR and the inlet/outlet positions. We characterised the networks in terms of the underlying blood velocity distribution (νm and νs parameters), and by introducing a metric quantifying the average number of capillary branches spins flow through, referred to as ANB. Over-all, our simulations generated a total of 1500 network realisations, which provide insight into microvascular blood perfusion. The most important observation is that the networks exhibit blood velocity distributions that can vary considerably from each other, with mean velocity νm ranging from approximately 4 to 25 mm/s. This is exemplified, for example, by Net 11 and 12 in Table 1, which feature a mean νm of 19 and 5 mm/s, despite exhibiting a similar mean capillary length of circa 60 μm. This finding suggests that, for the typical diffusion times that can be probed in clinical settings (15-65 ms approximately), spins in the vascular compartment likely experience flow regimes that can vary considerably from subject to subject. On the one hand, this implies that hypothesising a specific regime in IVIM modelling (e.g., diffusive versus ballistic [Kennan et al., 1994, Scott et al., 2021]), may not suffice to capture the full complexity of blood micro-perfusion in real-world cohorts. On the other hand, this also implies that remarkably different patterns of vascular dMRI signals are to be expected, even for short, clinically-feasible IVIM dMRI protocols, as exemplified by two examples in Fig. 3. Our simulated signals exhibit complex patterns as a function of the b-value and the diffusion time, e.g., fast decay, typical of the diffusive regime, as well as oscillatory behaviours, as instead expected in the ballistic regime (note that the PGSE signal for a set of uniformly distributed straight capillaries, characterised by the same blood velocity ν, is s = sinc(γ ν G δ Δ) [Scott et al., 2021]). Moreover, they also feature a clearly non-mono-exponential behaviour as a function of b, pointing again towards the co-existence of different flowing spin pools within the same network, potentially characterised by different flow regimes. All in all, these results suggest that numerical approaches such as SpinFlowSim may lead to more accurate characterisation of unexplored properties ofvascular dMRI signals — e.g., concerning flow anisotropy or apparent pseudo-diffusion and kurtosis tensors, as illustrated in Supplementary Fig. 6 for the apparent pseudo-diffusion tensor in an exemplificative case —, ultimately opening up new opportunities for the development of more specific biomarkers of micro-perfusion.
4.4 Microvascular property estimation
We also investigated whether it is possible to use the synthetic signals generated through SpinFlowsim to inform the non-invasive estimation of microvascular properties. For this purpose, we interpolated the full set of paired synthetic signals and microvascular parameters, obtaining numerical forward models that can be fitted through standard NNLS approaches. We specifically investigated the feasibility of estimating νm, νs and ANB, since the analysis of simulated signal suggests that even standard PGSE dMRI protocols carry some sensitivity to these properties (note that strong correlations between the vascular ADC D* and νm and ANB, or the moderate correlation between D* and νs, visualised in Fig. 4).
Firstly, we studied νm, νs and ANB estimation on noisy in silico signals. We considered 3 protocols: two were based on PGSE, with one of these very rich in terms of b-values and diffusion times, and another shorter and clinically feasible. One additional protocol was instead based on DW TRSE, matching that of our in vivo data. All protocols point towards the feasibility of estimating νm, νs and ANB: we observed strong correlations between ground truth and estimated νm and ANB, and moderate correlations for νs. As expected, performances were the highest for the richest protocol, with correlations as high as 0.79 for the ANB metric, yet still acceptable for the shorter protocols (e.g., correlation of 0.63 for ANB and the TRSE protocol). These promising results, obtained without requiring any explicit analytical modelling of the signal, highlight the potential utility of simulation-informed microvascular property estimation, motivating its testing in vivo.
Following in silico experiments, we moved on and tested whether νm, νs and ANB can also be estimated in vivo. For this purpose, we analysed dMRI scans acquired according to the TRSE protocol on two healthy volunteers and in 13 cancer patients. We fitted νm, νs and ANB alongside standard IVIM fV and D*, and assessed trends qualitatively in several organs in the healthy volunteers, and in the patients’ tumours.
In healthy volunteers, all metrics show high level of variability on visual inspection, which is confirmed by cross-organ trends in Fig. 7. The variability, qualitatively comparable between fV /D* and νm, νs and ANB, is in line with the well-known challenge of estimating microvascular property accurately with dMRI [Barbieri et al., 2020]. This finding suggests that more robust parameter estimation procedures may be needed than those used here (e.g., Bayesian fitting or deep learning [Barbieri et al., 2016a, Barbieri et al., 2020]), for the effective deployment of simulation-informed fitting in clinical settings. However, despite the variability, metrics show trends that are compatible with known physiology. For example, in healthy volunteers the liver shows much higher fV, D*, νm, νs and ANB than in the spleen. This finding is plausible considering that the liver is a highly vascularised organ, a blood reservoir receiving approximately 25 % of the cardiac output, despite representing only 2.5 % of the body weight [Lautt, 2010].
We also observe higher νm, νs and ANB in the kidney medulla than in the cortex, a finding that may be reflecting their different vascularisation. Regarding kidneys, we do not observe a clear trend in terms of cortex-medulla differences in standard IVIM fV and D* (e.g., fV is higher in the medulla than in the cortex for both healthy volunteers, while D* is in one case higher, and in the other lower). This is in line with recent studies, which have found high variability and strong inter-subject/inter-machine differences of kidney IVIM [Barbieri et al., 2016a,Ljimani et al., 2018, Stabinska et al., 2023].
Finally, we also demonstrated the feasibility of simulation-informed microvascular quantification in a pilot cohort of 13 cancer patients suffering from advanced solid tumours. While this demonstration only represents a first, exploratory proof-of-concept, it serves to highlight that contrasts seen in νm, νs and ANB are physiologically plausible, and consistent with patterns seen on fV and D*. For example, reduced νm, νs and ANB is seen in areas of low fV and D* compatible with reduced perfusion, expected in the tumour core [Karsch-Bluman et al., 2019,Herman et al., 2011], exemplified by Fig. 8. In vivo, νm, νs and ANB are positively correlated among each other, and they correlate moderately to strongly to IVIM D*. These correlations agree with the correlations observed in simulations (compare Supplementary Fig. 3 and 5), and may indicate that νm, νs and ANB, while providing complementary information to each other, are sensitive to similar, characteristics of the network morphology. For example, the strong correlation between νm and νs, indicating that higher variability in blood velocity has to be expected as the mean velocity increases, may be a signature of heteroscedasticity of the blood velocity distribution across capillaries.
All in all, our in vivo results demonstrate the feasibility of simulation-informed microvascular mapping in dMRI. While further confirmation and more detailed metric characterisation is required in future studies, realistic flow simulations informed by histology may increase the accuracy of dMRI microvascular signal models. Ultimately, these may provide innovative, biologically-specific indices of micro-perfusion, urgently sought for the non-invasive evaluation of cancer neo-angiogenesis, vascular heterogeneity and in treatment during the design of anti-angiogenic drugs.
4.5 Methodological considerations and limitations
In this article, we show the potential utility of flow simulations to inform dMRI signal modelling and analysis. We provide a first demonstration, based on a simplified simulation framework as a preliminary proof-of-concept. For example, we rely on an empirical expression for the resistance to flow across a capillary, borrowed from a model of cortical perfusion in the mouse primary sensory cortex [Blinder et al., 2013]. While this model accounts for salient features of blood flow resistance in capillaries (e.g., the effect of an average hematocrit and erythrocyte-wall interactions), a more realistic characterisation of the capillary resistance would be obtained by accounting for the Fåhræus-Lindqvist’s, the Fåhræus’ and the phase separation effects [Schmid et al., 2015, Van et al., 2021]. This would have required the simulation of the propagation of actual erythrocytes through the network, until a steady-state is reached, so that a per-capillary hematocrit (and hence, effective blood viscosity) can be calculated. Here we did not simulate erythrocyte propagation, being this computationally demanding. Nevertheless, we acknowledge that it would enable more realistic representations of flow patterns within micro-capillary networks. We plan to include erythrocyte flow in future work, and also extend SpinFlowSim to account, for example, for oscillatory pressure patterns and vessel deformation, and for fluid exchange between capillaries and the interstitial space.
Furthermore, for this first demonstration, we simulated vascular dMRI signals on 2D capillary networks. While Spin-FlowSim is designed to work with generic 3D networks, here we focussed on 2D representations due to the availability of 2D data (i.e., HE and CD31-stained biopsies). We accounted for this by averaging synthetic dMRI signals generated for two, orthogonal, in-plane gradient directions. However, in future we plan to increase the fidelity of our flow simulations by reconstructing 3D networks.
Related to the point above, the vascular networks reconstructed from histology for this article were obtained at the capillary level. Therefore, our synthetic signals may not be representative of larger vessels, including smaller feeding arterioles and small veins or venules. This implies that maps of νm, νs and ANB from our approach has to be taken with care in presence of larger vessels. In future, we plan to expand our vascular signal dictionary to include realisations of larger vessels, and thus improve the generalisability of our simulation-informed fitting.
Another point to acknowledge is that in this study we focussed on the characterisation of vascular dMRI signals, and devised a simulation-informed fitting procedure requiring pure vascular signals as input. For this reason, the analysis of in vivo signals required disentangling vascular from extra-vascular tissue signals, since low b measurements include contributions from both. This was achieved by extrapolating and ADC fit performed on b-values with negligible vascular signal contribution, and thus required identifying a b-value threshold. An approach of this type, i.e., splitting the vascular-tissue signal characterisation in two steps, is sometimes referred to as segmented IVIM fitting [Gurney-Champion et al., 2018,Wang et al., 2021]. While segmented fitting reduces the variability of vascular metrics, since it avoids the challenging, joint estimation of vascular and tissue properties [Barbieri et al., 2020], it may lead to biases in fV estimates, since fV may depend, at least slightly, on the b-value threshold. In future, we plan to improve the simulation-informed fitting performed here, by employing more advanced estimation techniques for the joint computation of vascular and tissue properties.
Lastly, we acknowledge that the results reported here should be confirmed by future studies. These would require the acquisition of data from additional healthy volunteers and from larger patient cohorts, and should include diffusion images from different MRI scanners and from more advanced dMRI protocols. Here, we used a simple acquisition scheme in the low b regime (b = 0 and b = {50, 100} s/mm2 at 3 diffusion times). However, the accurate characterisation of the complex signal patterns arising from microvasculature would likely require denser b samplings. Similarly, higher image quality and increased sensitivity to micro-perfusion could also be achieved, for example, by improving the robustness of the dMRI acquisition with cardiac/respiratory gating, or by employing flow-compensated [Wetscherek et al., 2015] gradient wave forms, or advanced b-tensor encoding [Nilsson et al., 2021].
5 Conclusions
SpinFlowSim, our open-source, freely-available python simulator of blood micro-perfusion in capillaries, enables the synthesis and characterisation of realistic microvascular dMRI signals. Perfusion simulations in vascular networks reconstructed from histology may inform the non-invasive, numerical estimation of innovative microvascular properties through dMRI, whose feasibility is demonstrated herein in vivo in healthy subjects and in cancer patients.
Data Availability
SpinFlowSim is made freely available as a GitHub repository at the permanent address: https://github.com/radiomicsgroup/SpinFlowSim. The repository includes the 15 vascular networks presented in this study that can be used to generate synthetic signals to inform model fitting. The code for simulation-informed fitting is freely available as part of BodyMRITools at the permanent address: https://github.com/fragrussu/bodymritools (script mri2micro dictml.py). The in vivo human data cannot be made freely available at this stage due to ethical restrictions.
Competing interests
This study has received support by AstraZeneca. K.B. was a researcher at VHIO (Barcelona, Spain), and is now an employee of AstraZeneca (Barcelona, Spain). AstraZeneca was not involved in the acquisition and analysis of the data, interpretation of the results, or the decision to submit this article for publication in its current form.
Data and code availability
SpinFlowSim is made freely available as a GitHub repository at the permanent address: https://github.com/radiomicsgroup/SpinFlowSim. The repository includes the 15 vascular networks presented in this study that can be used to generate synthetic signals to inform model fitting. The code for simulation-informed fitting is freely available as part of BodyMRITools at the permanent address: https://github.com/fragrussu/bodymritools (script mri2micro_dictml.py). The in vivo human data cannot be made freely available at this stage due to ethical restrictions.
Credit author statement
Conceptualization: AV, FG, RPL, DSN, EF, MP. Data curation: AV, FG, RPL, KB, AG, GS, PN, RT, EG. Formal analysis: AV, FG. Funding acquisition: RPL, EG, PN, RT, FG. Investigation: AV, FG, KB, GS, SS, MV. Methodology: AV, FG, AG, RPL, MP, DSN, EF. Project administration: AV, FG, KB, RPL, EG, PN, RT. Resources: AV, FG, RPL, AG, EG, PN, RT, ME, NR. Software: AV, FG, AG. Supervision: FG, RPL, RSL, DSN, EF, MP. Validation: AV, FG. Visualization: AV. Writing – original draft: AV, FG. Writing – review and editing: all authors.
Supplementary Material
Illustration of the PGSE (panel (a)) and DW-TRSE sequences used in simulations and implemented for actual in vivo data acquisition.
Histograms describing νm, νs and ANB distributions obtained across the 10 realisations of each vascular network for an input volumetric flow rate of qin = 3.1 ⋅ 10−3 mm3/s.
Spearman’s correlation coefficients, in form of a correlation matrix, for parameters νm (mean blood velocity), νs (standard deviation of blood velocities), ANB (apparent network branching), L (mean segment length) and r (mean radius length) from the 15 vascular networks segmented for this study. Values of νm, νs, ANB, L and r used to compute the matrix are reported in Table 1 (mean values).
Synthetic signals generated according to in silico protocol “TRSE”, and used for learning numerical forward signal models.
Spearman’s correlation coefficients, in form of a correlation matrix, for parameters νm (mean blood velocity), νs (standard deviation of blood velocities), ANB (apparent network branching), fV (IVIM vascular signal fraction) and D* (vascular ADC) from the 13 cancer patients included in this study. Values of νm, νs, ANB, fV and D* used to compute the matrix are reported in Table 2 (mean values).
Scatter plots illustrating the relationship between routine IVIM D* and the mean diffusivity of the apparent diffusion tensor 
for spins flowing in the ballistic regime [Kennan et al., 1994]. Above, t is the diffusion time (t ≈ Δ for PGSE in the narrow-pulse regime), while i = 1, 2, 3 and j = 1, 2, 3 iterate over the x, y and z components of the capillary velocity vectors, and ⟨⋅ ⟩c; is the ensemble average over the capillary set. The figure was obtained by pooling together data points for all 100 realisations of the 15 vascular networks drawn for this study, with qin = 3.1 ⋅ 10−3 mm3/s and δ = 0.1 ms. From top to bottom, clock-wise: (a): Δ = 1 ms; (b): Δ = 10 ms; (c): Δ = 30 ms; (d): Δ = 50 ms; (e): Δ = 70 ms; (f): Δ = 90 ms. The mean diffusivity of Di,j, by construction, increases with increasing diffusion time, a fact that is replicated by the D* metric obtained directly from our synthetic signals. However, there is a relatively weak-to-moderate correlation between D* and the apparent diffusion tensor mean diffusivity, suggesting that complex features of the vascular signal (which influence D* estimation) may not be captured by simple analytical expressions for diffusion cumulants derived in specific flow regimes.
Acknowledgements
VHIO would like to acknowledge: the State Agency for Research (Agencia Estatal de Investigación) for the financial support as a Center of Excellence Severo Ochoa (CEX2020-001024-S/AEI/10.13039/501100011033), the Cellex Foundation for providing research facilities and equipment and the CERCA Programme from the Generalitat de Catalunya for their support on this research. This research has been supported by PREdICT, sponsored by AstraZeneca. This study has been cofunded by the European Regional Development Fund/European Social Fund ‘A way to make Europe’ (to R.P.L.), and by the Comprehensive Program of Cancer Immunotherapy and Immunology (CAIMI), funded by the Banco Bilbao Vizcaya Argentaria Foundation Foundation (FBBVA, grant 89/2017). R.P.L is supported by the “la Caixa” Foundation CaixaResearch Advanced Oncology Research Program, the Prostate Cancer Foundation (18YOUN19), a CRIS Foundation Talent Award (TALENT19-05), the FERO Foundation through the XVIII Fero Fellowship for Oncological Research, the Instituto de Salud Carlos III-Investigación en Salud (PI18/01395 and PI21/01019), the Asociación Española Contra el Cancer (AECC) (PRYCO211023SERR) and the Generalitat de Catalunya Agency for Management of University and Research Grants of Catalonia (AGAUR) (2023PROD00178). The project that gave rise to these results received the support of a fellowship from “la Caixa” Foundation (ID 100010434). The fellowship code is “LCF/BQ/PR22/11920010” (funding F.G, A.V., and A.G.). This research has received support from the Beatriu de Pinós Postdoctoral Program from the Secretariat of Universities and Research of the Department of Business and Knowledge of the Government of Catalonia, and the support from the Marie Sklodowska-Curie COFUND program (BP3, contract number 801370; reference 2019 BP 00182) of the H2020 program (to K.B.). M.P. is supported by the UKRI Future Leaders Fellowship MR/T020296/2. A.G. is supported by a Severo Ochoa PhD fellowship (PRE2022-102586). The authors are thankful to the Vall d’Hebron Radiology department and to the ASCIRES CETIR clinical team for their assistance with MRI acquisitions, and to the GE/Siemens clinical scientists for their support with diffusion sequence characterisation.
Footnotes
↵** joint senior (last) authors
References
