Introduction

Dynamic contrast-enhanced (DCE) MRI has a developing role as an imaging tool for quantifying brain tumour microvasculature and tumour response to anti-angiogenic therapy1,2,3,4,5,6. In human clinical studies, non-invasive measurement of a suitable vascular input function (VIF) is essential for deriving microvascular kinetic parameters from brain DCE-MRI. Non-invasive measurement of a VIF with high temporal resolution is especially important when first-pass bolus tracking is necessary for quantitative kinetic analysis, such as when using the extended Tofts model7,8,9,10. Use of a fixed, experimentally derived population-averaged VIF for all subjects has been previously proposed11,12,13, thereby simplifying data acquisition, but large variations can occur in the actual VIF between subjects and scan visits due to both technical (e.g. differences in injection timing and dose) and patient specific factors (e.g. cardiac output, haematocrit, caffeine intake and atherosclerosis related vessel narrowing)11,14,15. For this reason methods have also been developed that enable simultaneous measurement of plasma gadolinium-based contrast agent (GBCA) concentration changes in both the blood and tissue under study, permitting VIF to be measured on an individual patient basis11,14,16.

In DCE-MRI studies of brain tumours the ideal choice for defining an individual VIF is the feeding artery of the tumour but due to either data acquisition constraints, the small size of the feeding vessel or lack of the feeding artery within the imaging field of view (FOV) this is often not possible11,16,17,18. Large intracranial arteries such as the internal carotid artery (ICA) and middle cerebral artery (MCA) are therefore often used as a surrogate global VIF measurement11,16,19. Use of the superior sagittal sinus (SSS), a venous structure, as a surrogate global input function has also been adopted, however, in numerous cross-sectional and longitudinal DCE-MRI studies1,2,5,11,16,19,20,21,22,23,24,25,26. The presented rationales for the use of the SSS vascular output function (VOF) as a surrogate global input function were said to be an assumption that the venous and arterial GBCA concentration is equal; that the degree of dispersion between arterial and venous structures in the brain during the first-pass circulation of the GBCA bolus is minimal; and that due to its larger comparative size the SSS is less susceptible to partial volume errors (PVE) than smaller intracranial arteries such as the ICA8.

The most commonly used 3D DCE-MRI acquisition is the 3D spoiled gradient echo method, and a VIF (plasma GBCA concentration-time curve) is typically determined from GBCA related magnitude changes in signal intensity27. Using a whole-brain, low-dose high temporal resolution (LDHT) axial 3D spoiled gradient-recalled echo sequence and a magnitude-based determination method we simultaneously measured a VOF/VIF from the SSS and large intracranial arteries (ICA, MCA) in a brain tumour patient cohort. Through this we sought to compare and understand features of each VOF/VIF and evaluate the respective ability of each VOF/VIF to accurately capture interindividual changes in patient dosing and plasma GBCA concentration. Through an included test–retest study we sought to establish the respective reproducibility of parameters derived from these arterial and venous VOF/VIFs; and through comparison with resected tumour specimens in a patient cohort, evaluate the ability of kinetic parameters derived using each VOF/VIF to detect intertumoural differences in histopathological data.

Methods

Study population

Previously acquired dual temporal resolution (DTR), dual injection DCE-MRI data in three groups of patients were analysed for this study: twenty-five patients with newly diagnosed WHO grade IV glioma synonym glioblastoma (GBM); twenty-nine patients with sporadic vestibular schwannoma (VS) listed for either radiological surveillance or treatment with surgery or stereotactic radiosurgery (SRS); and twelve patients with neurofibromatosis type 2 (NF2) related VS undergoing treatment with the anti-vascular endothelial growth factor (anti-VEGF) antibody, bevacizumab (Avastin ©). Ethical approvals were in place from the National Research Ethics Service Greater Manchester North-West research ethics committee (REC references: 13/NW/0131, 13/NW/0247 and 15/NW/0429). All patients had provided informed consent for study participation and later analysis of their MRI data and all research was performed in accordance with the Declaration of Helsinki and with local guidelines and policies.

MR imaging

The 25 patients with GBM and 29 patients with sporadic VS were all imaged once on a 1.5 T scanner (Philips Achieva, Best, Netherlands). Thirteen of the included patients with GBM and twenty-two of the included patients with sporadic VS had been recruited and scanned as part of previous published studies at our institution1,2,6. The twelve patients with NF2-related VS had similarly been recruited as part of an earlier published study investigating bevacizumab (Avastin ©) related changes in DCE-MRI derived kinetic parameters in VS and these patients had been imaged twice at 1.5 T: pre-treatment (day 0) and 3 months (day 90) following bevacizumab (Avastin ©) treatment.

DCE-MRI data was acquired using a previously described DTR, dual injection technique6,16,28. Single dose macrocyclic GBCA (gadoterate meglumine; Dotarem, Guerbet S.A.) was used at a dose of 0.2 ml/kg. For VOF/VIF estimation and as the first part of this DTR technique, a low-dose fixed volume pre-bolus (either 2 or 3mls) of GBCA was administered over 1 s during acquisition of a high temporal resolution (LDHT) DCE-MRI dataset. All intravenous injections were performed using a two-cylinder power injector (MEDRAD® Spectris Solaris EP, Bayer, PA, US). The GBCA and 0.9% saline are contained within separate cylinders and the pre-bolus injection was followed by a chaser of 20 ml of 0.9% saline administered at the same rate (2 or 3 ml/s). A 3D spoiled gradient recalled echo (GRE) sequence with axial slab orientation and anterior–posterior frequency encoding was used for data acquisition and acquisition parameters for this LDHT acquisition were as follows: flip angle of 20°, TR/TE of 2.5 ms/0.696 ms, SENSE acceleration factor of 1.8, reconstructed matrix size of 96 × 96 × 22, voxel size of 2.5 × 2.5 × 6.35 mm3, pixel bandwidth of 700 Hz, frame duration (Δt) 1.0 s (n = 300)6,16,28. The minimum TE and fixed volume low GBCA dose used for the LDHT DCE series and VOF/VIF estimation was designed to avoid signal magnitude saturation and expected to produce minimal T2* and water exchange effects29,30,31.

As the second part of this DTR DCE-MRI technique, a full-dose of GBCA (dose = 0.2 ml/kg ·weight – dose of pre-bolus) was administered at the same rate (2 to 3 ml/s) as the pre-bolus (followed by a chaser of 20 ml of 0.9% saline administered at the same rate) during acquisition of a high-spatial resolution (FDHS) sequence6,28. FFT (Fast Fourier Transform) reconstruction in the z-direction was used for both the low-dose high temporal resolution (LDHT) and full-dose high spatial resolution (FDHS) acquisitions, doubling the number of slices6. Variable flip-angle (VFA; α = 2°, 8°, 15° and 20°) acquisitions were undertaken prior to both the LDHT and FDHS DCE-MRI series for baseline longitudinal relaxation rate (R10) mapping, and the spatial resolution of each VFA acquisition series was chosen to match the LDHT and FDHS DCE series respectively6.

To eliminate unsaturated flowing spins entering the imaging slab and improve the accuracy of VIF estimation27,32 a large 3D acquisition volume covering the top of the brain, the circle of Willis and the terminations of the internal carotid arteries bilaterally was used16,33. Throughout the FOV, the number of radiofrequency (RF) pulses and the gradient spoiling that spins in flowing blood received was also maximized through the use of a fast spoiled GRE sequence with short TR and phase cycling. Gradient spoilers were applied along both the read and slice/slab selection directions and phase cycling with a phase increment angle of 117° was used. This allowed for more complete dephasing of residual transverse magnetization, minimizing blood inflow-induced errors within each imaging slice27,30,34. A pulse sequence diagram of the 3D spoiled gradient recalled echo (GRE) sequence used for both the LDHT and FDHS acquisition is shown in Fig. 1A.

Figure 1
figure 1

Magnitude-based input function extraction. (A): Pulse sequence diagram of 3D spoiled gradient recalled echo (GRE) sequence used for the variable flip angle (VFA) and low-dose high temporal resolution (LDHT) DCE-MRI acquisition, with axial slab orientation and anterior–posterior frequency encoding. Acquisition parameters: flip angle/s of 2, 8, 15, 20°; TR/TE of 2.5 ms/0.696 ms; SENSE acceleration factor of 1.8, reconstructed matrix size of 96 × 96 × 22, voxel size of 2.5 × 2.5 × 6.35 mm3, pixel bandwidth of 700 Hz, and frame duration (Δt) 1.0 s (n = 300). Pulse of frequency encoding/readout gradient (gfrequency), phase encoding gradient (gphase) and slice select gradient (gss) are shown along with signal acquisition time (tacq). Gradient spoilers (gspoiler) were applied along both the read and slice/slab selection directions and phase cycling with a phase increment angle of 117° was also used, allowing for more complete dephasing of residual transverse magnetization, and minimizing blood inflow-induced errors within each imaging slice. (B): Sagittal (left panel) and axial (right panel) signal-intensity (SI) images of a postcontrast time frame from the LDHT DCE-MRI series for a representative patient. Sagittal view shows that the anterior-middle portion of the SSS (red arrow) is included in the FOV of the 3D acquisition. This maximizes the number of RF pulses that the flowing blood experiences before it reaches the posterior portion of the SSS, helping to reduce blood inflow-induced errors. Following manual delineation of a small rectangle ROI within the vertical/posterior SSS (blue arrow), an automatic extraction method is used to identify voxels within neighbouring axial slices of the posterior SSS (blue arrows) that display maximum enhancement area under the SI curve within 30 s of the bolus arrival time (AUC30). A mean SI-time curve is then calculated from 20 voxels with the highest AUC30 and converted to a plasma GBCA concentration–time curve Cp(t) using a literature value of blood R10 of 0.694 s-1. (C): Axial SI images of a postcontrast LDHT DCE time frame for the patient shown in panel A showing the site of MCA VIF delineation. A pair of rectangle ROIs was drawn covering the horizontal segment of the MCA bilaterally (yellow rectangles, left panel). Using the described automatic method neighbouring contiguous axial slices were searched to identify twenty voxels within the MCA bilaterally (yellow arrows, right panel) that displayed maximum AUC30. (D): Axial (left panel) and coronal (right panel) SI images of a postcontrast LDHT DCE time frame for the patient shown in panel A showing the site of ICA VIF delineation. A pair of rectangle ROIs were manually drawn on an axial image section over the ICA (just distal to the carotid syphon) bilaterally. An automatic method was then used to search and identify twenty voxels within neighbouring contiguous axial slices of the carotid syphons (red arrows, right panel) that displayed maximum AUC30. (E): Coronal GBCA concentration C(t) image (mM) derived from the difference image ((ΔS = S(t)−S(0)) between the pre- and post-contrast SI images shown above. Whilst observed ring-like artefacts are present in the coronal SI images shown above, such artefacts are not present in either the 4D GBCA concentration—time images (where the VIF/VOF are measured) shown, or pharmacokinetic parameter maps derived from these 4D concentration–time images. GBCA = Gadolinium-based contrast agent; ICA = Internal carotid artery; MCA = Middle cerebral artery; SSS = Superior sagittal sinus.

Vascular input function extraction

For VOF/VIF extraction, acquired 4D LDHT DCE-MRI datasets (voxel size 2.5 × 2.5 × 6.35 mm3) were spatially aligned with and resliced to the FDHS data (voxel size 1 × 1 × 2 mm3) using SPM35. Although native LDHT datasets can also be used for input function extraction without prior co-registration and re-slicing, use of this co-registration step improved delineation of the blood vessel from surrounding tissues and permitted greater flexibility in ROI delineation and voxel selection. Plasma GBCA concentration-time curves, Cp(t), were derived from signal intensity (SI)—time curves measured at three different sites: 1) vertical/posterior third of the SSS, VIFSSS; 2) the ICA just distal to the carotid syphon, VIFICA: and 3) horizontal segment of the middle cerebral artery (MCA), VIFMCA. The first-pass data from these GBCA concentration-time curves were then fitted using a gamma variate function with a recirculation cut-off36, i.e.

$$ {\text{C}}_{{\text{p}}} \left( t \right) \, = Q\left( {t^{r} } \right)e^{ - (t/b)} ,{\text{ where Q}},{\text{ r and b are constants}}. $$

For each extracted VOF/VIF the following semi-quantitative parameters were derived and compared: the contrast-to-noise ratio (CNR), the bolus arrival time (BAT), the bolus peak-amplitude, the bolus peak width (FWHM, full-width at half-maximum) and the bolus peak∙FWHM product (PWP), a parameter associated with the area under the bolus curve.

For measurement of VOF in the SSS, blood inflow-induced errors were reduced by including the anterior and middle third of the SSS in the field of view thereby maximizing the number of RF pulses that the flowing blood experiences before it reaches the posterior third of the SSS. Following manual delineation of a small rectangle ROI within the vertical/posterior SSS (Fig. 1B), an automatic extraction method was used to search and identify voxels within this segment of the SSS that display maximum enhancement area under the SI curve within 30 s of the bolus arrival time (AUC30)37. Through this semi-automatic extraction method voxels with maximal AUC30 and thereby less inflow and PVE were chosen for inclusion in the VOF. The above semi‐automatic extraction method was also applied for VIF extraction from the ICA and MCA. A pair of rectangle ROIs was manually drawn on an axial image section over the ICA (just distal to the carotid syphon) or MCAs bilaterally respectively (Fig. 1). Similar to VOFSSS an automatic method was then used to search and identify voxels within neighbouring contiguous axial slices that displayed maximum enhancement (AUC30).

For each VOF/VIF a mean SI-time curve was calculated from 20 voxels with the highest AUC30 and this mean SI-time curve was then converted to a plasma GBCA concentration-time curve Cp(t) using previously described methods33. Due to the difficulties of accurately measuring the pre-contrast T1 of flowing blood using standard DCE-MRI sequences and the bias introduced through in vivo experiments38,39,40, a literature value of blood R10 of 0.694 s-1 was used for the conversion41,42.

Kinetic parameter analysis

To evaluate the effect of VIF approach on kinetic parameter estimates, high spatial resolution (1 × 1 × 2 mm3) voxelwise maps of the microvascular kinetic parameters Ktrans (transfer constant), vp (fractional plasma volume) and ve (the fractional volume of extravascular extracellular space or EES) were derived using the extended Tofts model (ETM)43 and the previously described LEGATOS (LEvel and rescale the Gadolinium contrast concentrations curves of high-temporal TO high Spatial DCE-MRI) method6. Pre-surgery DCE-MRI datasets from 15 sporadic VS with available comparative tissue histology were chosen as test group for this analysis and for each VOF/VIF (SSS, ICA, MCA) separate kinetic parameter estimates were derived43. In addition to derivation of Ktrans, vp and ve within each tumour voxel through the ETM, the voxelwise intracellular fraction (vi) was also estimated through the relationship vi = 1−ve−vp.

The LEGATOS method for deriving high-spatial resolution kinetic parameter maps from DTR, dual-injection DCE-MRI data has been previously described6. In key step 1 of this method, errors through temporal jitter uncertainty are reduced through construction of a merged DTR 4D GBCA concentration volume containing a high temporal (HT) resolution ‘arterial’ phase followed by a later low temporal but high spatial (HS) resolution ‘parenchymal’ phase6,44. In key step II the high temporal but low spatial resolution arterial phase of each pixel concentration curve is then re-scaled using the LEGATOS method and a derived pixelwise calibration ratio, to increase the spatial resolution of derived kinetic parameter maps. For the LEGATOS method, a combined VOF/VIF is adopted. The whole Cp (t) from the LDHT-derived VOF/VIF is concatenated with the dose-calibrated late part of the Cp (t) measured from the FDHS-derived VOF/VIF, and is used for kinetic analysis of the LEGATOS-generated 4D high spatiotemporal resolution GBCA concentration volume6.

The BAT for each tissue voxel is calculated as part of each fitting procedure and the Cp(t) measured from each VOF/VIF time-shifted and aligned with the BAT of each tissue voxel contrast agent concentration-time curve6. As part of the fitting procedure and to assess the discrepancy between the original data and the derived curve a map of scaled fitting error (SFE) was also generated, with voxels displaying an SFE value > 50% being excluded from the tumour statistics6,45. For all patients, the SFE and derived kinetic parameter maps, both before and after exclusion of voxels with SFE > 50%, were visually inspected to confirm the acceptance of using SFE > 50% for outlier tumour voxel exclusion6.

Tissue analysis

For the 15 resected sporadic VS, previously obtained tissue metrics were compared against derived kinetic parameter estimates using each VOF/VIF6. Collected paraffin blocks from each case were cut into serial 5-µm tissue sections and assessed for cell density (haematoxylin and eosin, H&E), vascular permeability (fibrinogen) and microvessel surface area (CD31) using immunoperoxidase immunohistochemistry and established protocols1,2,6. Ethical approval was obtained for tissue analyses (REC reference 15/NW/0429 and 19/NS/0167) and detailed protocols are described in prior publications1,2,6.

Statistical analysis

The SPSS statistical software package (version 25, IBM Corp.) and Stata version 11 were used for all statistical tests. Extracted semi-quantitative parameters were compared across each vessel (SSS, ICA and MCA) using a repeated-measures ANOVA with Greenhouse–Geisser correction for non-sphericity. Post hoc analysis of pairwise comparisons between different VOF/VIF locations was performed using the Bonferroni method. Due to the design of a fixed volume injection approach, the dose of pre-bolus slightly varied with the patient body mass, whilst keeping the same length of the bolus (1 s) for all subjects. This allowed the sensitivity of different VOF/VIFs to small variations in GBCA dose to be assessed. For each vessel region (SSS, ICA, MCA) the relationship between extracted features and administered GBCA dose (mmol/kg) was assessed using scatterplots and correlation analysis. Correlation analysis was also used to evaluate the relationship of the GBCA bolus arrival time delay between each arterial (ICA, MCA) VIF and VOFSSS, and differences in bolus peak-amplitude, bolus peak width (FWHM) and the bolus PWP between the SSS and either ICA or MCA. In particular, the correlation of the BAT delay with either the absolute difference in each semiquantitative parameter (bolus peak-amplitude, bolus peak FWHM, bolus PWP) or the ratio of each parameter between the SSS and ICA/MCA (e.g., PeakSSS/PeakICA/MCA ratio) was assessed.

The intra-subject repeatability of each VOF/VIF semiquantitative parameter was assessed across the twelve patients with NF2 related VS who were imaged twice using the test–retest coefficient of variation (CoV). The CoV is the standard deviation, σ, across all measurements for each subject, divided by the mean, μ, for that subject. For a group of N subjects the global test–retest CoV is defined as \(\sqrt{\sum {(\sigma /\mu )}^{2}/N}\)46,47. As a supporting measure of repeatability the average measures intraclass correlation coefficient (ICC) of each VOF/VIF semiquantitative parameter across the two visits was also calculated using an absolute-agreement, 2-way mixed-effects model48. The inter-tumour correlation between DCE-MRI derived parameter estimates (Ktrans, vp, ve and vi) and tissue-derived metrics (H&E cell density, CD31% microvessel surface area, fibrinogen optical density) for the 15 resected sporadic VS are reported as Pearson’s product moment correlation coefficient (r).

Results

Compared to large intracranial arteries (ICA, MCA) VOFSSS provides a global surrogate input function with higher CNR, higher peak and higher peak∙FWHM product (PWP)

In Fig. 2, representative plasma GBCA concentration curves Cp(t) derived from the SSS, the ICA, and the horizontal segment of the middle cerebral artery MCA following a bolus injection of 0.016 mmol/kg of GBCA are shown. As shown in Table 1, across all patient datasets the CNR of VOFSSS (mean CNR = 197.2 ± 95.7) was significantly higher than either VIFICA (mean CNR = 51.3 ± 25.9, p < 0.001) or VIFMCA (52.9 ± 22.7, p < 0.001, repeated measures ANOVA). There was no difference in CNR between VIFICA and VIFMCA (p > 0.05, repeated measures ANOVA, Table 1). Compared to either the VIFICA or VIFMCA, VOFSSS displayed significantly longer BAT (p < 0.001), higher peak (p < 0.001), larger bolus PWP (p < 0.001) and a non-significantly narrower FWHM (p > 0.05, repeated measures ANOVA, Table 1). There was no significant difference in either BAT (p > 0.05) or FWHM (p > 0.05) between VIFICA and VIFMCA, but VIFMCA displayed a higher peak (p < 0.001) and higher PWP (p < 0.001) than VIFICA.

Figure 2
figure 2

Typical plasma GBCA concentration–time curves Cp(t) (top row) derived from the vertical segment of the SSS (A); the ICA (B); and horizontal segment of the middle cerebral artery MCA (C) following a bolus injection of 2mls of GBCA (0.016 mmol/kg) at a rate of 2 ml/s. The first-pass data are fitted using a gamma variate function, which excludes contrast-agent bolus recirculation (bottom row). Images obtained from a patient with a sporadic VS imaged at 1.5 T. GBCA = Gadolinium-based contrast agent; ICA = Internal carotid artery; MCA = Middle cerebral artery; SSS = Superior sagittal sinus.

Table 1 Comparison of semi-quantitative VOF/VIF features extracted from superior sagittal sinus (SSS), internal carotid artery syphon (ICA) and middle cerebral artery (MCA).

In Fig. 3 a comparison is shown between an individual patient derived VOFSSS and the population VIF measured in the descending aorta by Parker et al13. In keeping with the observed lack of bolus widening (bolus FWHM) between the large intracranial arteries and the SSS, after converting the patient derived VOFSSS to a full-dose input function by summing several time-shifted low-dose GBCA concentration–time curves49, there was a close resemblance between the summed SSS defined VOF and the Parker population VIF13.

Figure 3
figure 3

Comparison between an individual patient derived vascular output function from the SSS (VOFSSS) and the population VIF measured in the descending aorta by Parker et al. VOF extracted from the vertical segment of the superior sagittal sinus (SSS) in the same patient as shown in Fig. 2. Six low-dose (0.016 mmol/kg of GBCA) VIFs (blue solid curves) derived from the SSS were time-shifted and summed to generate a full-dose (0.1 mmol/kg) input function (red dashed curve). The shape and amplitude of this summed SSS derived VIF shows high similarity with Parker's full-dose VIF (black solid curve), in keeping with the observed lack of bolus widening (bolus FWHM) between the large arteries and the SSS.

VOFSSS demonstrates a greater sensitivity to interindividual changes in plasma GBCA concentration compared to arterial approaches

Across all sixty-six patients, body weight adjusted GBCA dose for the pre-bolus injection and LDHT acquisition varied from 0.0091 to 0.027 mmol/kg with a mean injected dose of 0.016 mmol/kg. As shown in Fig. 4 both bolus peak (mM) and bolus PWP (mMs) correlated significantly (p ≤ 0.002) with body weight adjusted GBCA dose with the strongest correlation observed for VOFSSS peak (r = 0.70, p < 0.001) and VOFSSS bolus PWP (r = 0.65, p < 0.001). No significant correlation was observed between body weight adjusted GBCA dose and bolus peak width (FWHM) for any derived VOF/VIF (p > 0.05, Pearson’s correlation co-efficient). As shown in Table 2, no significant correlation (p > 0.05) was observed between GBCA bolus PWP ratio (PWPSSS/PWPICA or PWPSSS/PWPMCA) and the BAT difference between SSS and ICA or between SSS and MCA respectively. A longer time interval between BATICA and BATSSS correlated with an increase in the FWHMSSS/FWHMICA ratio (r = 0.35, p = 0.004) and was non-significantly associated with a lower PeakSSS/PeakICA ratio (r = − 0.18, p = 0.15).

Figure 4
figure 4

Relationship between extracted VOF/VIF features and administered GBCA dose for VOFSSS, VIFMCA and VIFICA. Top row: Scatter plots of bolus-peak (mM) vs administered GBCA dose (mmol/kg). Bottom row: Scatter plots of GBCA bolus peak∙FWHM product (PWP, mMs) vs administered GBCA dose. Pearson’s product moment correlation coefficient (r) shown with associated p value. FWHM = full-width at half-maximum; ICA = Internal carotid artery; MCA = Middle cerebral artery; PWP = bolus peak∙FWHM product; SSS = Superior sagittal sinus.

Table 2 Relationship between the GBCA bolus arrival time delay and differences in semiquantitative parameters between the arterial (ICA, MCA) VIFs and venous (SSS) VOF.

VOFSSS demonstrates greater repeatability compared to arterial VIF approaches

Table 3 shows the intra-subject variability and repeatability of each semiquantitative VOF/VIF parameter across the twelve patients with NF2-related VS. In the case of VOFSSS, global CoV values for BAT, bolus peak, bolus FWHM and bolus PWP were 3.98%, 17.0%, 16.8% and 12.4% respectively. Except for BAT, global CoV values for VIFSSS were lower than the corresponding global CoV values for VIFICA or VIFMCA. Across all semi-quantitative parameters extracted from the VOFSSS repeatability was good to excellent (ICC = 0.717–0.888)50,51.

Table 3 Repeatability of semi-quantitative parameters extracted from VOFSSS and arterial (ICA, MCA) VIFs.

Kinetic parameters obtained using a SSS derived VOF permitted detection of intertumoural differences in histopathological data

In Fig. 5 the inter-tumour correlation between LEGATOS derived kinetic parameter estimates and tissue metrics for each VOF/VIF are shown. There was a significant correlation between cell density and mean tumour vi when using VOFSSS (r = 0.54, p = 0.04, Fig. 5A). No such correlation was seen, however, when using either VIFICA or VIFMCA. In many tumours there was overestimation of ve when using arterial VIF (6/15 and 4/15 VS had ve > 0.7 when using VIFICA and VIFMCA respectively), and such high EES fractions were not evident on collected tissue from these tumours (Fig. 6). For VOFSSS and VIFICA a significant positive correlation was seen between CD31% microvessel surface area and mean tumour vp (p < 0.05), but this correlation was strongest for vp maps derived using VOFSSS (r = 0.85, p < 0.001, Fig. 5B). Estimates of both vp and Ktrans when using an arterial VIF were higher than VOFSSS derived estimates. Across all VOF/VIF approaches, there was a correlation of Ktrans with both CD31% microvessel surface area (p < 0.05) and perivascular leak as measured through fibrinogen optical density (p < 0.05).

Figure 5
figure 5

Scatterplot comparison of histopathological data with kinetic parameter estimates (vi, vp, Ktrans) derived using different VOF/VIF approaches. (A): Intertumour scatterplot comparison of mean tumour intracellular fraction (vi, no units) estimates against mean H&E cell density (nuclei/ x20HPF). vi estimates derived using VOFSSS (top row), VIFICA (middle row) and VIFMCA (bottom row) shown. (B) Intertumour scatterplot analysis of mean vascular fraction (vp, no units) against mean CD31% microvessel surface area (SA). vp estimates derived using VOFSSS (top row), VIFICA (middle row) and VIFMCA (bottom row) shown. (C) Intertumour scatterplot analysis of mean tumour Ktrans (min-1) against mean fibrinogen optical density (OD). Ktrans estimates derived using VOFSSS (top row), VIFICA (middle row) and VIFMCA (bottom row) shown.

Figure 6
figure 6

Parameter maps and tissue sections from a large sporadic VS demonstrates overestimation of kinetic parameters (ve, vp, Ktrans) with use of VIFICA. (A): Representative LEGATOS derived kinetic parameter maps from a patient with a large right sided sporadic VS. From left to right: parametric ve map; parametric vp map and parametric Ktrans map. Note the comparatively higher ve vp and Ktrans estimates within the tumour (arrow) and nasal mucosa (*) when using the ICA derived vascular input function (VIFICA). (B): Representative tissue sections from the tumour shown in panel A. From top: Haematoxylin and eosin-stained section (HE-x20HPF); CD31 immunostain for microvessels (brown; immunoperoxidase –x20HPF); fibrinogen immunostain for perivascular leak (brown; immunoperoxidase –x20HPF). In keeping with the VOFSSS kinetic parameter estimates, the tumour is moderately cellular (mean tumour vi ~ 0.35) with some regions of moderate/high microvessel density (CD31) and perivascular leak (fibrinogen). The VIFICA and VIFMCA kinetic parameter estimates, however, overestimate both the size of the tumour extravascular-extracellular space (ve) and the degree of tumour vascularity (vp).

Discussion

Compared to large intracranial arteries (ICA, MCA), an arrival-corrected VOF derived from the SSS (VOFSSS) can provide a superior surrogate global input function for kinetic parameter analysis in brain DCE-MRI, demonstrating higher CNR, higher peak and a greater sensitivity to interindividual changes in plasma GBCA concentration. Through an included test–retest study we demonstrate that semi-quantitative parameters derived from VOFSSS display greater repeatability than arterial based approaches. Furthermore, through comparison with matched tissue datasets in patients with resected sporadic VS we demonstrate that microvascular kinetic parameters obtained using an arrival-corrected VOFSSS permit evaluation of intertumoural differences in microvessel surface area and cell density, a feature not seen with large artery based VIFs.

Previous studies comparing arterial and venous based VOF/VIFs have reported similar results to our findings, and several factors can be hypothesized to contribute to the observed higher CNR and higher peak measured within VOFSSS compared to large arteries19,32. Keil et al. demonstrated that arterial VIF showed lower peak Cp values and more often led to implausibly high ve values in subsequent pharmacokinetic parameter fitting compared to a SSS derived VOF19. The larger cross-sectional area of the SSS relative to both the ICA and MCA may serve to reduce PVE and contribute to the higher CNR, higher peak GBCA bolus measured32,39. Reduced ‘in-flow’ effects secondary to lower flow velocities within the SSS may also contribute to the higher CNR observed25,52,53,54. To maximise T1 weighting, a 3D T1W GRE sequence with the shortest TE was used for the DCE-MRI acquisition. Such GRE sequences are, however, prone to ‘in-flow’ related enhancement, which can lead to a significant attenuation of contrast enhancement by GBCA and subsequent reduction of measured VIF CNR27,30. Within our study, measures were taken to prevent these ‘in-flow’ artefacts such as RF phase cycling, use of a large acquisition volume covering the whole brain, and orientation of the read gradient of the 3D slab parallel to the vessels. It is, however, possible that such ‘in-flow’ related artefacts still occurred39 and because the contrast between flowing and static tissue depends directly upon the flow velocity53,54, the higher velocity arteries are much more prone to increased ‘in-flow’ effects and resulting lower CNR than lower velocity venous structures such as the SSS25,39,52.

An additional factor that may have contributed to the observed decreased CNR and peak within arterial VIFs relative to the SSS is differences in inter- and intravoxel velocity dispersion within these vessels39. Compared to large intracranial arteries, blood flow within the SSS is thought to be less turbulent and more laminar, with a parabolic velocity profile55,56. Assuming laminar flow and a parabolic velocity profile, then the velocity v of a given proton spin across a circular vessel, is a function of distance to the isocentre of the flow, i.e., v(r) = vmax (1−r2/a2), where a is the radius of the vessel, r is the distance between the position of the flowing spin and the centre of the vessel, v(r) is the velocity at r, and vmax is the maximum velocity at the centre of the vessel. Whereas stationary spins are completely rephased during data acquisition, moving spins retain a degree of phase shift that is a function of both the first moment of the gradient M1 and the flow velocity v. The phase increases linearly with the flow velocity but there is no effect on the spatial location of the signal by 2D Fourier transform reconstruction because the velocity is set constant from view to view. The location would, however, be affected if acceleration in the vessel is heterogenous and in vessels where there is significant intra- and intervoxel velocity heterogeneity, phase dispersion and loss of CNR can occur. Moreover, changes of the parabolic velocity profile can influence the width of first-pass of the contrast bolus. In vivo blood flow measurement is considerably more complex, blood is non-Newtonian for example and the flow profile depends on the geometry of the vessels and the sampling point of the pulsed flow, i.e., the time points in the cardiac cycle. The field gradients for spatial encoding and the spoiler gradients used in our acquisition were not specifically designed for ‘flow compensation’ and extra spin dephasing caused by spins moving during the actual signal encoding may have occurred. Such flow-induced loss of phase coherence increases with higher velocity blood flow and may in addition to PVE and residual ‘in-flow’ effects explain the lower CNR and lower peak measured within arterial VIFs compared to VOFSSS.

Within this study a rapid bolus injection protocol was used, with a compact low-dose bolus of GBCA injected over 1 s. Following injection into an upper limb vein, the GBCA bolus is expected to undergo a degree of dispersion and widening as it passes through first the pulmonary and then cerebral circulation. This dispersion process is commonly modelled as a mathematical convolution with the vascular system of the brain and a simple yet useful model for the dispersion process might contain a single artery that feeds a set of parallel pathways (an arteriole, a blood-tissue exchange unit, and a venule) that drain into a single vein57,58,59,60,61. Data from our study, however, suggests that when first-pass low-dose GBCA bolus uptake curves are measured within the distal SSS such dispersion is minimal, with no significant difference in peak width (FWHM) between VOFSSS and the VIF measured from feeding larger arteries such as the ICA and MCA. This absence of significance bolus dispersion between the ICA and SSS has also been reported in other studies and several factors may underlie this unexpected observation8,39. A large contributing factor is likely inaccuracies in the measurement of the VIF within the smaller arterial vessels due to increased PVE, increased vessel inflow effects and increased velocity dispersion of flowing blood relative to the SSS. Such effects likely contribute to the increased CNR seen within the VOFSSS but also counterbalance the broadening effect by dispersion8.

Previous works have demonstrated that incorrect estimation of VIF peak and the bolus arrival time can have significant impact on pharmacokinetic parameter fitting and accuracy62,63,64,65. In particular the bolus arrival time in the SSS is 4–5 s later than the true input arrival time and a correction must be made for the delayed arrival of the contrast agent within the venous system by incorporating voxel-wise BAT estimation in DCE-MRI analysis66,67. Within this study use of a DTR acquisition incorporating high temporal sampling during pre-bolus injection was adopted, improving voxel-by-voxel estimation of BAT delay, and reducing fitting errors induced by uncertainty in time alignment of the VOF/VIF and tissue uptake curves. Alongside measurement effects, physiological factors may also play a role in increasing CNR and reducing the observed dispersion seen within the SSS relative to arterial VIF sources. The SSS serves as the principal draining vein for the cerebral cortex and it’s comparatively straight course, absence of valves, and wide diameter relative to large intracranial arteries means that laminar blood flow characteristics within it resemble that of similar sized arteries22,24,25,68. Indeed a similarity between the measured VOF within the SSS and Parker’s population averaged VIF has been reported in both this and other studies13,19. Human histological studies and microsphere studies in rhesus monkeys have demonstrated that within normal brain under physiological conditions, a small percentage of blood from the arteries passes through small pial arteriovenous shunts (12.5 µm or more in diameter) directly into cerebral veins and the SSS69,70,71 Direct shunting of blood into cerebral veins has also been demonstrated angiographically within supratentorial GBM and VS, although in the latter case the source of principal venous drainage is the transverse and sigmoid sinuses rather than the SSS itself72,73,74,75,76,77,78. Although in theory the presence of such shunts may serve to reduce dispersion of the GBCA bolus within the venous system, the predicted low volume of blood passing into the SSS via these direct channels makes them unlikely to be a significant factor in driving the observed VOFSSS characteristics.

The presented study is one of the largest comparisons of different arterial and venous based input functions for brain DCE-MRI and to our knowledge one of the first studies to evaluate the sensitivity of different VOF/VIF approaches to interindividual differences in plasma GBCA concentration and histopathological data. A limitation of the present study though is that characterization of plasma GBCA concentration curves, Cp(t), was limited to semi-quantitative parameters such as the bolus peak-amplitude, bolus width and GBCA bolus PWP. Future studies should seek to undertake more detailed and sophisticated shape analysis of the whole length of the VOF/VIF. In this study, a semi-automatic extraction method was employed to measure each VOF/VIF, in which the twenty vessel voxels with the highest AUC30 are selected and averaged together to create the final input curve 18,37,79. This is under the assumption that these voxels demonstrate less inflow effects and PVE. Although it is possible that such an approach may overestimate the true value of VOFSSS, our demonstration that the VOFSSS showed an initial peak height and first-pass bolus shape very close to the ‘gold standard’ Parker VIF model13,80,81, and that the variation in the VOFSSS bolus peak was strongly correlated with dose variation, suggests that such overestimation was minimal and that peak estimates are not dominated by bolus shape distortion. Our demonstration that kinetic parameter estimates, obtained using VOFSSS correlated well with tissue derived measures of microvessel surface area and cell density, in contrast to arterial VIFs that showed lower peak values and overestimation of ve, further supports the robustness of input function surrogate measurements from the SSS and the absence of significant peak overestimation through this semi-automatic voxel selection method. Larger studies incorporating matched imaging-tissue cohorts in a range of different tumours should, however, be undertaken to better evaluate the effect of VOF/VIF location and voxel extraction method on kinetic parameter accuracy.

Conclusion

Accurate derivation of a vascular input function (VIF) is essential for quantitative kinetic analysis of brain DCE-MRI data. In this in vivo patient study, we compared VIFs extracted from either the internal carotid artery and its branches with an arrival-corrected vascular output function derived from the superior sagittal sinus (VOFSSS). We demonstrated that compared to large intracranial arteries VOFSSS can provide a superior surrogate global VIF, with lower noise, higher repeatability, and greater sensitivity to interindividual changes in plasma GBCA concentration. Through comparison with matched histopathological data, we furthermore demonstrate that microvascular parameters obtained using a SSS derived VOF permitted evaluation of intertumoural differences in microvessel surface area and cell density. These results support the use of venous sinus-based approaches for input function extraction and pharmacokinetic parameter mapping in brain DCE-MRI.