Prediction of Lung Shunt Fraction for Yttrium-90 Treatment of Hepatic Tumors Using Dynamic Contrast Enhanced MRI with Quantitative Perfusion Processing

There is no noninvasive method to estimate lung shunting fraction (LSF) in patients with liver tumors undergoing Yttrium-90 (Y90) therapy. We propose to predict LSF from noninvasive dynamic contrast enhanced (DCE) MRI using perfusion quantification. Two perfusion quantification methods were used to process DCE MRI in 25 liver tumor patients: Kety’s tracer kinetic modeling with a delay-fitted global arterial input function (AIF) and quantitative transport mapping (QTM) based on the inversion of transport equation using spatial deconvolution without AIF. LSF was measured on SPECT following Tc-99m macroaggregated albumin (MAA) administration via hepatic arterial catheter. The patient cohort was partitioned into a low-risk group (LSF ≤ 10%) and a high-risk group (LSF > 10%). Results: In this patient cohort, LSF was positively correlated with QTM velocity |u| (r = 0.61, F = 14.0363, p = 0.0021), and no significant correlation was observed with Kety’s parameters, tumor volume, patient age and gender. Between the low LSF and high LSF groups, there was a significant difference for QTM |u| (0.0760 ± 0.0440 vs. 0.1822 ± 0.1225 mm/s, p = 0.0011), and Kety’s Ktrans (0.0401 ± 0.0360 vs 0.1198 ± 0.3048, p = 0.0471) and Ve (0.0900 ± 0.0307 vs. 0.1495 ± 0.0485, p = 0.0114). The area under the curve (AUC) for distinguishing between low LSF and high LSF was 0.87 for |u|, 0.80 for Ve and 0.74 for Ktrans. Noninvasive prediction of LSF is feasible from DCE MRI with QTM velocity postprocessing.


Introduction
Trans-arterial radioembolization (TARE) with yttrium-90 (Y90) microspheres has been widely adopted as a primary locoregional treatment for patients with hepatocellular carcinoma (HCC) and liver-dominant hepatic metastases [1][2][3]. Evaluation of lung shunting fraction (LSF) is a critical step in Y90 treatment planning because increased intra-hepatic arterio-venous shunting poses a risk of lung damage and treatment failure [4]. Instead of lodging into the tumor microvasculature, the Y90 microspheres may travel through these shunts and into the lungs leading to radiation pneumonitis, especially if the absorbed radiation dose to the lungs exceeds 30 Gy in a single session or 50 Gy cumulative [5]. Moreover, lung shunting may decrease the Y90 radiation dose to the hepatic tumor and reduce the treatment effect [6][7][8]. Currently, LSF estimation requires a pre-treatment invasive arteriogram using Technetium-99m macroaggregated albumin (Tc-99m MAA) injected into the hepatic artery, followed by single-photon computed emission tomography (SPECT) imaging to determine the amount of Tc-99m MAA shunted to the lungs [9,10]. This LSF estimation using Tc-99m MAA SPECT assumes similar biodistribution in human body for Tc-99m-MAA and Y90 following hepatic arterial injection [11], and requires a dedicated hepatic artery catheterization in a separate earlier hospital visit increasing the costs in time, financial expense, and risk to patients [12]. Based on a recent analysis, the cost to Medicare of a visit for mapping angiography and SPECT imaging is $14,194 in 2021 US dollars [13]. A noninvasive, cost-effective method to estimate LSF without catheterization would improve the value of Y90 therapy to patients.
We hypothesize that LSF can be estimated from pre-TARE non-invasive imaging that sensitizes hepatic tissue blood flow. Tumors with higher arterio-venous shunting are expected to have higher blood velocities, as arterio-venous shunting bypasses small, high resistance vessels such as capillaries. Arterio-venous shunting requires tumor angiogenesis and may also change the plasma space and extravascular extracellular space (EES) volume. These changes can be captured by MR perfusion-weighted imaging such as dynamic contrast enhanced (DCE) MRI. Therefore, it may be possible to predict LSF by perfusion quantification from DCE MRI.
Recently, the quantitative transport mapping (QTM) method [14,15] has been developed to estimate mass flux characterized by velocity |u| without a global arterial input function (AIF) used in traditional Kety's tracer kinetic analysis [16,17]. QTM velocity has been shown to be more accurate than Kety's parameters in validation with numerical ground truth [14] and more sensitive than Kety's parameters for differentiating benign from malignant tumors compared with biopsy [18,19]. Accordingly, we propose to investigate the feasibility of noninvasive prediction of LSF according to QTM velocity |u|, as well as Kety's parameters, derived from DCE MRI.

Patients and MRI Protocols
This retrospective study was Institutional Review Board approved and HIPAA compliant. Twenty-five patients (age range: 42 to 87 years; 8 female and 17 male) were diagnosed with HCC as determined by Liver Imaging Reporting and Data System (LIRADS) [20]. All patients underwent trans-arterial radioembolization therapy utilizing Y90 resin or glass microspheres after diagnosis with the dose determined by either the body surface area (BSA) or medical internal radiation dose (MIRD) method. LSF was calculated from a single-photon computed emission tomography (SPECT) scan after Tc-99m-macroaggregated albumin (MAA) administration in the target hepatic artery [21,22]. Before Y90 treatment, the patients were scanned on a 1.5 T scanner (MAGNETOM, Siemens Medical Solutions, Erlangen, Germany) with 3D CAIPI-VIBE sequence before, during and after gadolinium contrast agent (Gd) injection [23]. Gd was injected with a dose of 0.1 mmol/kg (gadobutrol; Bayer Healthcare Pharmaceuticals Inc., Whippany, NJ, USA) at 2.0 mL/s, followed by a 20 mL saline flush using an MR compatible injector. The scanning parameters were as follows: In plane resolution 0.84 mm, slice thickness 3 mm, temporal resolution 5 s, repetition time 4.68 ms, echo time 2.39 ms, flip angle 10 • .

Quantitative Perfusion Processing of DCE MRI
We assumed a linear relationship between tracer concentration and relative enhancement on DCE MRI [24]. The tracer concentration was then processed using quantitative transport mapping (QTM) and traditional kinetic modeling method. In QTM, the tracer concentration profile is modeled by a transport equation [14,25]: where c(r, t) is the tracer concentration scalar field, r = r x , r y , r z the spatial location in a volume of size N x , N y , N z , ∂ t the time derivative, ∇ = ∂ x , ∂ y , ∂ z the spatial gradient operator, and time index t ∈ {1, 2, . . . N t − 1} the time index with N t as the number of time frames. u(r) = (u x (r), u y (r), u z (r)) is an average velocity vector field [25]. Both time derivative and gradient operator are implemented as finite differences in a discretized 4D spacetime-resolved image space. Equation (1) is a linear equation system for velocity that is solved as an optimization problem with L1 total variation regularization as in a recent QTM study with the regularization parameters λ = 10 −4 chosen according to the L-curve method [14]: In traditional Kety's tracer kinetics (also known as extended Tofts' model), the tracer concentration profile is modeled by [26]: where c a (t) is the global AIF, K trans is volume transfer constant, V e (r) is the volume fraction of extravascular extracellular space (EES), and τ is traveling delay. Equation (3) is a linear equation system for K trans and k ep = K trans V e , but is nonlinear to τ. A voxel-wise non-linear least squares method is used to solve for kinetic parameters and traveling delay τ of AIF with the regularization parameters λ = µ = 10 −3 chosen according to the L-curve method [27][28][29]: (4) For each case, an experienced radiologist drew the AIF in the tumor-feeding artery and manually segmented regions of interest (ROI) consisting of the tumor targeted by Y90 radioembolization. |u|, K trans , V e and τ were averaged over these ROIs and used for further statistical analysis. All reconstructions were performed on a computer using an Intel i7-8700K 6-core CPU with 64GB memory.

Data Analysis
The patients were separated into 2 groups according to SPECT measured LSF values: a low-risk group (LSF ≤ 10%) and a high-risk group (LSF> 10%) [30,31]. Tumors were manually segmented on post-Gd MRI to define lesion region of interest (ROI) and to measure tumor volume. Only the lesions that went through Y90 treatment were included in further study, and the average tumor size was 18.6 ± 15.7 cm 3 . Statistical analysis was performed on 6 parameters: ROI averaged |u|, K trans , V e , τ, as well as tumor volume and patient age. A Spearman correlation test was performed to test the relationship between these parameters and lung shunt fraction. The Mann-Whitney U test was performed comparing these values between the low-risk and the high-risk groups. p-values at or below 0.05 were considered to indicate statistical significance. A receiver operating characteristic curve (ROC) analysis was performed to investigate the risk prediction performance of all parameters, and the optimal threshold was calculated by maximizing sensitivity plus specificity. Statistical analysis was performed using the R Statistical Software (Foundation for Statistical Computing, Vienna, Austria).

Discussion
These preliminary data from 25 patients undergoing Y90 radioembolization of HCC demonstrate a significant correlation between QTM velocity |u| derived from DCE MRI and lung shunt fraction (LSF) derived from SPECT post-Tc-99m-MAA administration, suggesting the potential for using noninvasive DCE MRI to predict LSF for distinguishing low lung shunting risk subjects (LSF ≤ 10%) from high lung shunting risk subjects (LSF > 10%). This result encourages the development of noninvasive imaging, including DCE MRI with QTM velocity |u| postprocessing, for accurate estimation of LSF to eliminate the cost and risk of catheterization in planning trans-arterial radioembolization with Y90 treatment.
Lung shunting is caused by the vasculature remodeling with tumor progression, creating high flow arterio-venous shunts that bypass high resistance capillaries [32,33]. The changes in tumoral vascularity may increase blood flow, velocity, plasma space and EES space volume, which are reflected on DCE MRI images and quantified by transport physics. In this study, we compared QTM |u| with parameters from traditional kinetic modeling (also known as extended Toft's model) [26]. In the patient cohort, QTM velocity |u| but not Kety's parameters demonstrated significant correlation with LSF. Although Kety's method has not been applied to lung shunting fraction estimation in previous studies, this result is consistent with the hypothesis that increased artery-vein connections bypassing capillaries (shunts) increases the mean liver blood velocity, and also consistent with previous reports showing QTM improves upon Kety's method by replacing a global arterial input function in Kety's model with the local mass flux gradient in QTM [14,18,19]. In addition to a significantly higher velocity |u| observed in high LSF group, we also observed an increase in Kety's parameters K trans and V e in high LSF group, which may reflect a higher tissue exchange rate and EES space in tumors with abnormal vasculature compared to normal tissue [34,35]. These observations should be evaluated in further studies, especially using histopathology for validation.
The LSF prediction accuracy may be further improved by increasing the spacetime resolution of DCE MRI [36]. Other MRI sequences that reflect vasculature properties such as arterial spin labeling, 3D/4D phase contrast MRA and susceptibility weighted imaging may also be incorporated into noninvasive LSF prediction. The artery-vein shunting process may change the arterial and venous blood volume and the average vessel size in a voxel. Velocity-sensitive pulses such as multiphase balanced steady-state free precession (bSSFP) pulse train and velocity-selective RF pulse train can be used to directly measure arterial and venous blood volume [37][38][39]. The average vessel size in the voxel can be estimated by measuring the R2 and R2* changes after intravascular superparamagnetic contrast agent injection [40][41][42]. These sequences may allow detailed microvasculature modeling and computational fluid dynamics prediction of the passage of Y90 microspheres through arteriovenous shunting [14].
Accurate LSF prediction may require detailed image features and microcirculation information in addition to simple lesion ROI values used in this study. For example, texture analysis or radiomic features of the transport quantity maps, deep learning and larger data sets may help construct the LSF prediction model. We combined QTM velocity ROI values with Kety's parameter ROI values but the combination failed to improve the LSF prediction accuracy, which may suggest that QTM already contains all the information in Kety's parameters as well as better information than Kety's parameters, consistent with previous publications [14,18,19]. Historically, Tc-99m-MAA has been a well-developed radiotracer for pulmonary vascular imaging, especially for studying pulmonary embolism since 1960s [43] and has been conveniently adopted for estimating LSF since the beginning of radioembolization practice [44][45][46]. Fundamentally, the use of LSF assumes the good agreement between pretreatment Tc-99m-MAA distribution and final Y90-microsphere distribution. However, this assumption is highly problematic [47] since there may be differences in catheter position between the two hospital visits, physiologic variances in hepatic blood flow, and size and morphology differences between Tc-99m-MAA particles and Y90-microspheres. Accordingly, future development should focus on predicting biodistribution of Y90 microspheres based on maps of transport quantities such as using computational fluid dynamics (CFD) modeling [48,49]. This CFD modeling would estimate both Y90 microsphere distribution in the targeted liver region and Y90 microsphere scape into the lung.
There are several limitations in this study that need to be addressed in the future work. Firstly, the sample size is relatively small (n = 25), thus no significantly difference in AUC between QTM parameters and Kety's parameters was observed. More cases will more thoroughly evaluate the LSF prediction performance of our proposed model. Secondly, our DCE MRI images contain only five temporal points and cover only 5 s of the arterial phase after Gd injection. A longer scan with higher spacetime resolution acquired with motion compensation [50,51] may help to increase the reconstruction accuracy of quantitative transport quantities. Thirdly, the assumption of linear relationship between tracer concentration and relative enhancement of DCE MRI may be problematic, which can be improved using quantitative susceptibility mapping [52,53].
In conclusion, noninvasive LSF prediction from DCE MRI with quantitative transport mapping postprocessing is feasible. Future work should include enlarging the dataset and exploring the combination of DCE MRI and other MR pulse sequences to improve the LSF model performance for predicting Y90 biodistribution.

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki, and approved by the Institutional Review Board (or Ethics Committee) of Weill Cornell Medicine (#0909010639, approved 10/21/2022).

Informed Consent Statement:
Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.