A model selection framework to quantify microvascular liver function in gadoxetate‐enhanced MRI: Application to healthy liver, diseased tissue, and hepatocellular carcinoma

We introduce a novel, generalized tracer kinetic model selection framework to quantify microvascular characteristics of liver and tumor tissue in gadoxetate‐enhanced dynamic contrast‐enhanced MRI (DCE‐MRI).


| Gadoxetate-enhanced MR imaging of the liver
Gadoxetate disodium (Eovist or Primovist, Bayer Healthcare, Berlin, Germany) is a gadolinium-based hepatobiliary contrast agent taken up by hepatocytes and excreted by the biliary pathway, allowing more direct measurement of liver function than standard extracellular contrast agents. It is regularly used in the clinical assessment of chronic liver disease and cancer. [1][2][3][4] In recent meta-analyses of hepatocellular carcinoma (HCC, most common primary malignancy of the liver 5 and second leading cause of cancer mortality worldwide 6 ) detection, gadoxetate-enhanced MRI returned the highest overall per-lesion sensitivity and positive predictive value, compared with contrast-enhanced ultra-sound, CT, and MRI with extracellular contrast agents. 7,8 In addition to diagnosing HCC, relative lesion intensity to surrounding liver tissue has been shown to characterize subtypes of the tumor, indicating progression of nodules to hypervascular HCC, 9 microvascular invasion, 10 and correlating with prognostic histological biomarkers. 11 Gadoxetateenhanced imaging may also help identify and grade chronic liver disease, including portal hypertension 12 and fibrosis, 13 often co-present in patients with HCC.

DCE-MRI of the liver
Quantitative dynamic contrast enhanced (DCE) MRI, where tracer kinetic models are fitted to image time series, enables estimating parameters reflecting tissue microvascular function. In metastatic liver tumors, there is a long history of tracer kinetic modeling using standard gadolinium-based contrast agents. 14-17 A common approach assumes metastases recruit an arterial-only supply, with passive exchange of contrast agent between separately measurable plasma and interstitial spaces, enabling the fit of a 2-compartment extended-Tofts model. 18 Model parameters, and in particular, the intercompartment volume transfer constant K trans , can provide clinically useful endpoints in trials of antivascular drugs. 14,17 To account for liver parenchyma's arterial and hepatic portal vein blood supply, a dual-input single-compartment model was proposed for application outside tumors, 19 and used in patient studies showing a correlation between perfusion parameters and the severity of cirrhosis or fibrosis. 20 The singlecompartment model was generalized into a 2-compartment model suitable for parameterizing tracer kinetics in liver tissue and neuroendocrine metastatic lesions. 21,22 In an alternative approach, the dual-input single-compartment model and an extended-Tofts model were fitted to metastatic cancers, 15 and the Aikake information criterion (AIC) 23 used to select the most appropriate model at each voxel.
There have been far fewer quantitative DCE-MRI studies of HCC than metastatic lesions: applying an extended-Tofts model, decreased K trans on first scan has been shown to predict survival, 24 and the dual-input single-compartment model used to quantify HCC perfusion. 25 In the only study we are aware of applying quantitative DCE-MRI to gadoxetate enhanced HCC tumors, 26 fitting a single-input (arterial only) Tofts model showed correlations for K trans and the extracellular, extravascular volume v e to histological grade and tumor microvascular density, respectively.
Other tracer kinetic modeling studies of gadoxetate include initial attempts to estimate intracellular uptake fraction as a marker of liver function 27,28 ; however, these were limited by using models comprising a single-input blood supply, which further failed to account for active transport of gadoxetate into the intracellular space. These shortcomings were resolved in a dual-input, uptake model 29,30 that provided significantly better fit than the single-compartment model. To model gadoxetate dynamics through to the hepatobiliary phase (typically reached 20 minutes postinjection), an efflux term was added, and applied to time series acquired using an 50-minute protocol in healthy volunteers. 31 Recently, a 3-compartment model was proposed to quantify hepatic perfusion and hepatocyte function in patients with chronic liver disease, applied to a 38-minute acquisition. 13

| Modeling challenges
With limited previous studies, and different models applied in each study, there is no common method for applying tracer kinetic modeling to gadoxetate-enhanced images. There are several challenges to consider: contrast dynamics in tissue where hepatocytes actively transport gadoxetate will be very different to tumor cells where gadoxetate behaves more like a passively transported agent, indicating the need for (at least) 2 forms of tracer-kinetic model. Moreover, in tumors native to the liver such as HCC, there may be residual uptake of contrast agent through the organic anion polypeptide (OATP) pathway due to varying OATP function in the cells of different tumor subtypes, 1 meaning that even with well-defined tumor regions-of-interest (ROIs), it is not certain a priori which model will be appropriate in any voxel. Thus computing parameter estimates over a whole ROI risks including voxels that do not meet a model's physiological assumptions, potentially leading to unreliable estimates and misleading interpretations. This suggests voxel-wise model selection may be a more suitable method for analyzing the data. 15,32 Common co-morbidities associated with HCC provide additional complications, with patients having a range of liver | 1831 BERKS Et al.
function and signs of cirrhosis, fibrosis and splenomegaly which may further confound model assumptions if analysis is extended to the non-tumorous liver.
We investigate fitting tracer kinetic models to individual voxels of gadoxetate-enhanced MRI sequences, using models developed to explain the active transport of gadoxetate and the 2-compartment exchange model (and its derivatives including the Tofts model). 18,19,29,31,33 We introduce a novel model selection framework in which previously applied models can be derived as specific instances of generalized functional forms. In doing so, we highlight mathematical connections between the models, using these to hypothesize which will better fit different tissue types. We test these hypothesizes using extensive Monte Carlo simulations, exploring model selection criteria. We compute expectations on the precision and accuracy of derived model parameters, and show how these are affected by sequence duration, temporal resolution, and signal noise. Finally, we test our hypothesizes on 2 in vivo datasets: the healthy volunteers presented in Ref. [31] and a new dataset of 10 patients with HCC, showing links between model selection and patient co-morbidities, determining the limitations of the data and using these to make recommendations for designing further patient studies.

| Tissue models
Georgiou et al 31 presented a model to describe both the active uptake of gadoxetate contrast agent by hepatocytes and the subsequent efflux of the agent into the bile duct ( Figure 1A). The model incorporates a dual-input function to allow for arterial and hepatic portal vein blood supply, convolved with a biexponential impulse response function (IRF): The active uptake + efflux model (AUEM) is a 2-compartment biexponential model with dual arterial and portal vein plasma input, designed to model the active transport of gadoxetate contrast agent by hepatocytes. Contrast agent is assumed to be instantaneously well mixed in the extracellular volume v ecs , and then actively transported into the intracellular volume v i with uptake rate K i , which together comprise the whole voxel space. Gadoxetate is later excreted to the bile ducts at rate K ef . The model can be simplified first assuming zero efflux (K ef = 0), then zero uptake K i = 0, which reduces to the dual-input single compartment model. B, In the 2-compartment exchange model (2CXM), the interstitial v e and plasma volumes v p are assumed separable, with measurable passive-exchange between them at a rate depending on the permeability surface area product PS. It is assumed cells in these voxels do not actively uptake contrast agent, and thus there is no contribution to signal change from the intracellular space. See Sourbron and Buckley 34 for a comprehensive discussion of simplifications of the 2CXM to the other passive-exchange regime models such as the (extended-) Tofts model. The 2CXM is often considered with a single arterial input; however, our framework includes a portal vein supply (dashed box), generating a dual-input model with the same functional form as the AUEM where C a and C v are the arterial and venous input functions of plasma concentration, offset by a (minutes) and v (minutes) and f a ∈ [0, 1] controls the proportion of each input. The IRF is parameterized by: uptake rate K i (min − 1 ); efflux rate, K ef (min − 1 ); plasma flow, F p (mlmin − 1 ml − 1 ), and extracellular volume v ecs (mlml − 1 ). The model assumes v ecs comprises the plasma volume (v p ) and interstitial space v e , with rapid exchange between the 2 (ie,, endothelial permeability → ∞) that immediately reaches equilibrium, such that v ecs acts as a single compartment. 29,34 In addition, it is assumed the volume of nonhepatocyte cells in the intracellular space v i is negligible, so that v i = 1 − v ecs , and the contribution of contrast agent in bile canaliculi can be ignored. 35,36 We label this the active uptake + efflux model (AUEM), noting it extended an earlier uptake-only model described by Sourbron et al, 29 derived from the full AUEM by setting the efflux rate K ef = 0. Setting the uptake rate K i = 0, further reduces the uptake-only model to a dual-input singlecompartment model. 19 This hierarchy of nested models evokes the analysis by Sourbron and Buckley, 34 in which it is shown a biexponential 2-compartment exchange model (2CXM, Figure 1B), with physiological parameters F p (plasma flow), PS (permeability surface area product, min − 1 ), v e , and v p , reduces to the extended-Tofts (assuming F p = ∞) and Tofts model (assuming v p = 0).
In a similar vein, we note that while the 2CXM is usually considered with a single input function, the form of its IRF is identical to the AUEM, and thus if we fit a dual-input 2CXM, with no restrictions on the parameters, the modeled contrast time series is be identical to the AUEM-it is only the derivation and interpretation of the physiological parameters that differ.
We hypothesize that the AUEM (or one of its simplifications) suitably describes the dynamics of gadoxetate in healthy liver tissue, whereas a 2CXM is more likely to fit tumor voxels, where hepatocytes are not present to actively transport gadoxetate, and the plasma and interstitial spaces are separable volumes v p and v e , with a measurable bisymmetric exchange of gadoxetate between them.
This leads us to adopt the following approach: we fit as the most general model, a dual-input biexponential of the form This has 6 free parameters to optimize: the IRF functional parameters ± , ± , and physiologically meaningful parameters controlling the arterial input fraction and delay time f a and a = v (Equation (4), see also Supporting Information: Delay parameters). We then derive 4 further physiological parameters from the functional parameters, using either the active uptake or passive exchange interpretations (see Supporting Information: Model derivations for full complete intermediary equations): Thus, the 2 model regimes share 3 physiological parameters (F p , f a , and a ), differing only in the interpretation of the 2 compartment volumes and the transfer rates. Moreover, although they differ in physiological interpretation, the extracellular volume v ecs in the active-uptake regime is functionally equivalent to the vascular volume v p in the 2CXM (both being the initial compartment receiving vascular input; Equations (7) and (13)).
These interpretations impose constraints on valid ranges for the parameters-in the active-uptake regime v ecs ≤ 1 and in the passive-exchange regime, v e Hct)). However, we do not impose these constraints during model fitting. Instead, the 4 functional parameters ± , ± are simply constrained to be non-negative, and then post-fit, we analyze whether there is a valid physiological interpretation of the parameters under either regime. This allows the optimizer to explore the complete parameter space of both physiological models without requiring a prior assumption of which specific instance should be preferred.

| Nested models
In addition to the full biexponential model, we fit a hierarchy of progressively simplified models, setting first + = 0, then + = 0. The nested models are defined by the form of their impulse responses, taking 4, 3, and 2 parameters, respectively, labeled as the I 4 , I 3 , and I 2 models. When we derive physiological parameters in either the active-uptake or passive-exchange regimes from a nested model, we obtain parameters that match one of the reduced forms of the AUEM or 2CXM previously outlined in Tissue models. We present these specific forms below.

| I 3 model:
In the active-uptake regime, when Bpos = 0 from Equations (6)-(9), we see that F p , v ecs , and K i are well-defined and positive, while K ef = 0. Thus, the I 3 model is functionally equivalent to the uptake-only model. 29 In the passive-exchange regime, v ecs is ill-defined, while F p , PS, and v p are positive. This is described as the 2-compartment uptake model in Ref. [33], and may be appropriate in tissues where backflow of the tracer into the plasma space is negligible, where there are poorly mixed compartments, or extraction/convection/ diffusion into additional compartments. Moreover, from Equations (8) and (11), with + = 0 we have thus given the prior observation that v ecs = v p , with a simple relabeling, the 2 regimes generate equivalent physiological parameters in the I 3 form.

| I 2 model:
Fixing + = 0, sets K i = 0; as a simplification of the AUEM, this matches the dual-input single-compartment model 19 (although now there is no active transfer of contrast agent). In the passiveexchange regime, the I 2 functional form is most commonly interpreted as a standard Tofts model, 18 with negligible vascular component (v p = 0). However, as noted in Ref. [34], there are several other interpretations matching the same analytical form, one of which, with negligible exchange (PS = 0, v p > 0, v e undefined), is equivalent to the single-compartment model derived from the AUEM. Thus again, subject to a simple relabeling, the 2 regimes produce identical physiological parameters.

| Model optimization and selection
Given the framework outline above, there are 2 aspects to model selection: choosing a suitable analytical form, and choosing the best interpretation of parameters given this form.

Analytical form
An analytical form can be selected examining residuals between the signal-derived and modeled concentration time series and determining which model best fits the datachoosing the model with minimum sum-of-squares errors (SSE) or applying model selection criteria such as AIC (23) to take into account model complexity.
Our results suggest that the percentage of voxels in a region attributed to each model form using AIC may itself be a useful marker. We label these as

Physiological interpretation of parameters
We assign voxels to the active-uptake or passive-exchange regimes based on the physiologically plausible limits on compartmental volume parameters derived from the I 4 model, defining and the percentage of all valid voxels in a region-of-interest assigned to the active-uptake regime as % a (see Supporting Information: Selecting physiological regime for details).

| Monte Carlo simulations
Given model parameters, a set of time points and corresponding arterial and portal vein input functions (AIF/PIF), we can use Equations (5)-(13) to simulate dynamic time series. We generated multiple 1,000 sample datasets, each defined by their model form  i , temporal resolution , post-contrast duration d and noise distribution . We used a single fixed ground truth parameter set for each of 4 model forms (units omitted for clarity): For each form, physiological parameters were converted to functional parameters ( ± , ± ) using the active-uptake ( 1 - 3 ) and passive-exchange ( 4 ) regimes. The functional parameters were then used to compute modeled concentration time series using all combinations of: • Temporal resolution = 3.8 seconds (≈ patient dataset), 6.0 seconds (≈ volunteers dataset), 12.0, 30.0 seconds. • Post-bolus duration d = 6, 12, 18 minutes.
In each case, a population AIF 37 and derived PIF (Supporting Information: Input functions, 15 ) were sampled at the generated times.
Each concentration time series was converted to simulated signals with baseline S(0) = 100 (Supporting Information: Computing signal from concentration), and replicated 1,000 times with randomly sampled Rician noise distributions = Rice(S(t), ) for = 10, 20, 30, 40. These signal data were then converted back to concentration time series.
In addition, using the group average AIF and derived PIF from the patient dataset, we generated "real" noise datasets with temporally varying noise distribution r , by (1) randomly sampling 1,000 voxels from across the patient dataset; (2) subtracting the fitted I 4 model concentration from the signal-derived concentration; and (3) adding the resulting residuals to the simulated concentration time series for each model form  1 - 4 . A second experiment using varying parameters randomly sampled from measured distributions of real data is described in Supporting Information.

| Model fitting
The I 2 , I 3 , and I 4 models were fitted to each Monte Carlo sample, and AIC used to compute % AIC 2 , % AIC 3 , and % AIC 4 for each dataset. Physiological parameter sets in the activeuptake and passive-exchange regimes were derived for each fitted model, and the I 4 − volume criterion used to compute % a . In addition, for each dataset, parameter distributions were generated using (a) the AIC-selected model at each sample; and (b) the model with minimum SSE.

| Healthy volunteers
The healthy volunteer dataset comprised of 8 volunteers (5 male, aged 18-29 years; 3 female aged 18-22 years; mean age 23 years) and were imaged with contrast injection at 2 visits on a Philips 1.5 T Achieva MRI scanner. DCE-MRI was performed (see Supporting Information: In vivo data for protocol details), with contrast agent gadoxetate disodium (Primovist/Eovist, Bayer) was administered at the 20th time point (120 seconds) using a power injector (dose 0.025 mmolkg −1 ; 3 mLs −1 flow rate, flushed with 20 mL of saline at the same rate).
These data were collected for a previous study (approved by the university's ethics committee; all subjects gave informed consent prior to imaging) by Georgiou et al, with full acquisition details and parameter estimates for the activeuptake and efflux model fitted to ROI-averaged time series presented in. 31 For the experiments presented here, to enable efficient motion correction and voxel-by-voxel model fitting, the volumes were cropped to the margins of the liver, and only the first 264 time points used.

Patients with HCC
The second dataset collected for this study comprised of 10 patients with HCC (8 males and 2 females, aged 46-80 years; median 63 years) who were yet to begin treatment and had at least 1 lesion that could be measured in one dimension, according to the Response Evaluation Criteria in Solid Tumors (RECIST, see Supporting Information for eligibility criteria, tumor properties, and known co-morbidities). The study received Institutional Board Approval. Informed consent was obtained from all patients.
Each patient underwent a free-breathing coronal 3D FLASH DCE-MRI protocol on a Siemens 1.5 T Avanto system (see Supporting Information: In vivo data for protocol details). Gadoxetate was administered at the 8th time point ( 30 seconds) using a power injector (dose 0.025 mmolkg −1 ; 3 mLs -1 flow rate).

| Model fitting and selection
Each dataset were motion corrected, annotated, and used to compute vascular input function (Supporting Information: Data preparation). The I 2 , I 3 , and I 4 models were fitted to the DCE time series of all voxels in the liver volume for each patient (Supporting Information: Model optimization). For healthy volunteers, the model framework was applied, testing 3 durations of sequence: 6 minutes post-contrast (≈ patient dataset; n t = 80), 12 minutes (n t = 138), and 24 minutes (n t = 264). Due to the computational cost of fitting models to longer time series, the models were fitted to 1000 randomly selected voxels across the liver volume, and for a spatial visualization, to all voxels in the central slice of the liver.
Model parameters were converted to physiological parameters in both regimes, and the I 4 − volume criterion used to assign each voxel as active uptake or passive exchange. In addition, both AIC and minimum SSE were used to assign a model form at each voxel, with the percentage of AIC selections of each form used to analyze trends in the data, and minimum SSE used to select a single parameter set at each voxel. Parameters were summarized for each ROI taking the median over all voxels assigned to each regime.

| Effect of temporal resolution
Increasing temporal resolution led to improved parameter estimates ( Figure 3 shows K i estimation for model form  1 , see Supporting Information Figures S14-S37 for all other parameters), and better ability to select the ground-truth model form (Figure 2, Supporting Information Figures S7-S9). At low resolution ( = 30 seconds), with the exception of the single-compartment model form, parameter estimation is generally unreliable (at long duration, it may be possible to estimate v ecs and K i in active-uptake forms, and v e in the passive-exchange form, although precision is significantly worse than at higher temporal resolutions). As a result, at all but the lowest noise levels, the I 2 model is selected by AIC score in the majority of samples, for all model forms (Supporting Information Figure S9).
As reduces from 3.8 to 12.0 seconds, some parameters appear more sensitive to the loss of resolution. For example, f a , which is largely characterized by the shape of the contrast time series in a short period after bolus injection, is affected more than parameters such as v ecs and K i (Figure 3, Supporting Information Figures S14-S19).

| Effect of duration
At short duration (d = 6 minutes), I 3 is the AIC-preferred functional form for both active-uptake forms at all but very low noise levels (Figure 2, Supporting Information Figure S7). Parameter predictions of F p , v ecs and K i are consistent in both models (Supporting Information Figure S14, S17, S18). I 4 predictions of K ef lack precision but are not biased ( Figure S19). As the duration increases (d = 12 minutes), it becomes possible to resolve K ef , which in turn reduces the ability of the I 3 model to properly represent the latter time points of the time series, leading to the increasing preference for the I 4 model in AIC selection ( Figure 2) and the tendency for the I 3 model to underestimate K i (Supporting Information Figure S18). Despite this, at noise levels observed in the real datasets, the I 3 model is still the AIC-preferred model in approximately half the samples, leading to the overall sample median of both K i and K ef being underestimates of the ground truth if AIC selection is used to generate final parameter F I G U R E 2 Analytical form and physiological regime model selection percentages for the fixed Monte Carlo datasets, with temporal resolution = 3.8 seconds (see Supporting Information for other resolutions). Ground truth model forms are shown in columns: active uptake and efflux (AUEM); active-uptake only (AUM); single-compartment (no uptake, SCM); 2-compartment exchange (2CXM). At 6 minutes duration, solid colored circles show results for the "real" noise datasets (see main text). SNR was computed by dividing the mean signal by the root-meansquared error of the noisy signal from the ground truth signal estimates. This trend continues in longer durations, so that for time series d = 24 minutes, the I 4 is preferred in the large majority of samples (>90% at real noise levels).
Similar effects are seen for the 2-compartment exchange ground truth form. At 6 minutes, while the I 4 model is correctly selected with negligible noise, approximately half of samples at real noise levels are assigned to the I 3 model (Figure 2), producing strongly biased estimations of F p (Supporting Information Figure S20) and (unidirectional) rate constant PS (Supporting Information Figure S23), as well as leaving v e undefined. This is largely resolved at 12 minutes, where the I 4 model is consistently selected for the large majority of samples, and % a ≈ 0 at all noise levels. Table 1 shows model selection percentages and parameter estimates for the volunteers dataset. At all durations, nearly all voxels were categorized as active-uptake (95.5% at 6 minutes, 97.8% at 12 minutes, and 98.3% at 24 minutes, Figure 4A), thus we only include active-uptake regime parameters in further analysis.

Model selection
At all durations, few voxels were best matched to the I 2 model in AIC selection (0.7%, 0.3%, and 0.2%, respectively). At 6 minutes, most voxels attributed to the I 3 model over the I 4 model (82.3% vs 16.4%), with model selection shifting to the I 4 model as duration increased (27.4% vs 71.9% at 12 minutes; 4.4% vs 94.8% at 24 minutes). These trends were spatially consistent ( Figure 5A-C), and were observed in all 8 subjects, suggesting a consistent pattern of behavior in healthy liver tissue. This supports the hypothesis that at durations before efflux significantly effects gadoxetate concentration an uptake only model is sufficient to fit the data. However as duration increases the additional efflux term is necessary. 13,31 Parameter estimates If, as suggested by the Monte Carlo simulations, estimates for F p , v ecs , K i , and K ef are the most accurate using the I 4 model over the longest duration time series, we can observe how well estimations from the less complex models and/or shorter post-contrast durations correspond with these. I 4 estimations of K ef correspond well at 12 minutes (linear correlation coefficient = 0.91, P <. 001, F I G U R E 3 Estimation of K i for the AUEM in the fixed Monte Carlo datasets at a single noise level (most closely matching measured real noise). Dashed black lines show parameter ground truth (K i = 0.07); circles show median estimation over the 1,000 samples in each dataset using each IRF form, AIC, and minimum SSE, with vertical lines extending from the 25th to 75th percentiles. For temporal resolution = 3.8s and duration d = 6 mins, the cyan diamond and black lines show median and interquartile range for parameters estimation using minimum SSE on the dataset with directly sampled "real" noise added (see Supporting Information Figures S14-S37 for other parameters, model forms and noise levels)

BERKS Et al.
slope m = 1.09, offset c = 0.00) but not at 6 minutes ( = 0.45, P =. 08, m = 0.6, c = 0.01); I 4 estimates of F p , v ecs and K i were accurate at both 6 and 12 minutes (all > 0.98, P <. 001, m ≈ 1). I 3 estimates of F p , v ecs and K i correspond best at 6 minutes-at 12 and 24 minutes, estimations become inaccurate because the model is insufficiently flexible to fit the plateauing concentration of gadoxetate as efflux starts to take effect, and thus the final fit reflects a compromise between fitting the early and later parts of the time series. F p and v ecs are poorly estimated by the singlecompartment I 2 model at all durations.
The arterial fraction f a proved harder to fit in some volunteer visits than either the Monte Carlo or patient data, and resulted in consistent over-estimations (ie,, f a ≈ 1 for the majority of voxels) in 3 DCE sequences. These data had lower SNR than the rest of the dataset, most likely due to unresolved motion artefacts. The lower temporal resolution of the volunteer data may also be a factor. Analysis suggests f a overestimation results in F p under-estimation; however, other parameters are largely unaffected. Parameter medians varied more between visits than between subjects, although across the group no parameter median was significantly different between visits.  Table 2 gives per patient data. Figures 4B,C and 5D,E show the spatial consistency of model selections and the extent to which they correspond with (independently demarcated) tumor borders.  Parameter estimates for the group are shown in Table 3. For the parameters shared in both model interpretations (F p , f a , a ), mean plasma flow F p was similar in liver and tumor (cohort median 0.45 vs 0.44 mlmin − 1 ml − 1 ); arterial fraction f a was significantly higher in tumors (0.37 vs 0.18, Willcoxon signed rank test P =. 011), and arterial delay time significantly shorter (0.07 minutes vs 0.14 minutes, P = . 027). Figure 6 depicts waterfall plots of active-uptake regime physiological parameter estimates in the liver in comparison with healthy volunteers. There were no significant differences in K i , f a , and a ; however, the patient group has significantly lower F p (Wilcoxon ranked sum test P = .002) and higher v ecs (P = .003).

Pathophysiology of non-tumorous liver
In patient data, the pattern of model selection in nontumorous tissue broadly matched the healthy volunteer cohort. However, there was a general reduction in gadoxetate uptake as seen in the average concentration time-series plots ( Figure 7D,E), reflected in reduced estimates of F p and v i = 1 − v ecs , higher values for % AIC 2 and lower values for % a . This suggests reduced perfusion and a lower volume of hepatocytes in liver tissue, which may indicate reduced function associated with co-morbidities to HCC.
This deviation from the healthy pattern was particularly strong in 2 patients, both of whom had an abnormally low % a . One of these (patient 8, % a = 58.1) also had an abnormally high proportion of I 2 voxels (% AIC 2 = 33.2); SNR and model residuals were in the mid-range of the cohort, suggesting this was not caused by poor model fitting, but is instead consistent with restricted hepatocyte uptake of gadoxetate. This accords with the clinical observation that patient 8 has chronic hemochromatosis, resulting in a severely cirrhotic liver.
The second (patient 6, % a = 65.7) had higher than average I 4 voxels (% AIC 4 = 29.8), and an uptake rate in the nontumorous liver K i = 0.015 minutes − 1 less than a third of any other patient or healthy volunteer. Again, SNR was in range, suggesting the results were a genuine reflection of abnormal liver tissue function, and is consistent with complications observed in the subsequent treatment of patient 6, in which the location and size of the tumor caused a partial obstruction of the biliary tree.

Gadoxetate uptake within HCC
While it was expected tumor voxels would be more likely than non-tumorous liver voxels to be assigned to the passiveexchange regime, there are significant variations between patients and the active-transport regime is still assigned to a significant number of tumor voxels (34.6% across the cohort). In some cases, these may correspond to misalignment in ROI boundaries, or as a result of the dynamic duration being too short to properly resolve the I 4 model form.
However, an analysis of parameters for active-uptake voxels in the tumor suggests they are systematically different from those in non-tumorous tissue (Table 3). Specifically, they have increased arterial and reduced venous supply (ie, higher f a , 0.36 vs 0.17), and significantly reduced uptake parameter K i (0.019 vs 0.064).
Plausible interpretations of these observations are that such voxels are composed of cells exhibiting residual OATP function, allowing uptake of gadoxetate through the OATP pathway, or that there is a partial volume effect between T A B L E 3 HCC patient data: median (IQR) a values for each physiological parameter of the active-uptake and passive-exchange regimes  tumor cells and functioning liver tissue (or a combination of both). In either interpretation, if further clinical or histological information were available, a per tumor analysis of this effect may help understand clinically relevant physiological differences.

| DISCUSSION
Analysis of the in vivo data conforms to our initial hypotheses and reflects the results of the Monte Carlo simulations: dynamic time series in healthy liver tissue are best fitted by models that account for active transport of gadoxetate by hepatocytes from a well-mixed arterial and venous blood supply. In comparison, tumors tend to exhibit passive-exchange of contrast agent, with a more rapid and arterially dominated blood supply. However, there are also significant differences within the non-tumorous liver-where chronic disease such as cirrhosis and fibrosis may reduce (or even eliminate) portal blood supply or inhibit hepatoctye function-and within tumor sub-types. These results highlight the necessity of fitting multiple models to characterize tissue across the liver volume. Relying on a single preselected model interpretation risks fitting a model that does not adequately match the data and thus produces unreliable parameter estimates. Instead, fitting models in the functional space, and then post-fit, choosing the most appropriate physiological interpretation at each voxel provides a more robust analysis, capable of dealing with variations between individual patients and between/within heterogeneous tissues. Moreover, the data, and in particular the varying ways in which patients 6 and 8 differ from the cohort norm, suggest analyzing both the appropriate physiological interpretation and the selected level of model complexity provides clinically relevant information.
Analyzing model functional form helps clarify the relationship between different physiological interpretations and their associated assumptions of tissue microstructure. From a practical perspective, it is much more efficient to extract multiple physiological derivations from a single optimization, as the former has negligible cost compared with the computationally intensive process of fitting a complex non-linear model to each voxel. The functional forms in our framework do not allow direct measurement of sinusoidal backflux. This could be incorporated by extending the toplevel 2-compartment models to a 3-compartment model, 13 although this would significantly increase computation and may require longer acquisitions to provide sufficient data points to unambiguously resolve all model parameters.
In common with many previous studies of DCE-MRI in the liver, we did not take steps to suppress liver fat in these data. This may lead to errors in T 1 estimation and that may subsequently bias the model parameters. The relaxivity of gadoxetate should also be carefully considered. This is reported as 6.9 mmol −1 s −1 plasma at 1.5 T 38 ; however, a subsequent study measured values approximately twice as high in liver tissue (14.0 mmol −1 s −1 , 39 ). We used the former for converting signal to concentration for AIF generation, and the latter for model fitting in the liver (see Supporting Information), while accepting the limitation that patient disease may also effect gadoxetate relaxivity (eg, due to reduced albumin).
With only 10 patients imaged once using a 6-minute postcontrast acquisition, there are limits to the conclusions that can be drawn from the HCC data. A larger cohort, longer DCE protocol, with repeat baseline imaging would enable us to assess the repeatability of the analysis, while follow-on imaging would allow testing of the clinical utility of specific parameters as biomarkers in tracking disease progression or treatment response. However, our results correlate with the clinical presentation of the patients. Moreover, the spatial consistency of model selection and parameter maps suggests that the framework is suitable for quantifying the microvascular characteristics of both non-tumorous liver tissue and HCCs on a per voxel basis, despite the challenging data quality resulting from a free-breathing protocol. Thus we propose that while HCC patients are a study group of clinical significance, there is sufficient evidence to repeat the analysis described here in a larger cohort featuring longitudinal imaging.

| Recommendations for future studies
The in vivo data, together with the Monte Carlo simulations, suggest using a DCE protocol of at least 12 minutes postcontrast when using gadoxetate as a contrast agent. This allows estimation of efflux (K ef ) in the active-uptake regime, and significantly improved estimation of the extracellular, extravascular volume fraction v e in the 2CXM, as well as improving the differentiation between the 2 regimes in potentially ambiguous tissue (eg,, in HCC, or compromised non-tumorous liver).
If only shorter acquisitions are possible, as with the patient data in this study, it must be accepted that estimations of K ef will lack precision, and there will be greater ambiguity in determining the appropriate physiological regime, and thus care should be taken in any conclusions drawn from the data. Where longer protocols are used, it is necessary to fit a full biexponential (I 4 ) model form (or alternative model with sufficient complexity) to describe both the initial uptake and subsequent plateauing of gadoxetate concentration. A simpler mono-exponential (I 3 ) form may appear to fit the time series (and may even by preferred under AIC selection in a significant proportion of voxels), but parameter estimates derived from the I 3 model will likely be biased.
A corollary of these observations is that if using a nested model approach, it is better to use the model with minimum SSE (rather than model selection criteria such as AIC or an F-test) from which to estimate physiological parameters. Despite this, there are still potential advantages in fitting the nested I 3 and I 2 models and applying AIC selection. A consistent pattern in model selection was observed in non-tumorous liver tissue across the healthy volunteers and in patients that did not have clinical observations of significant liver co-morbidity. This suggests that deviations from this pattern are a sign that either there is insufficient data quality to fit a full I 4 model to individual voxels, or that there may be clinically relevant reduction in liver function to assess. Such observations may be noticed more easily by a simple check on the model selection percentages than trying to interpret a more complicated distribution of multi-variate physiological parameters. Although there is a computational cost to fitting the additional models, optimization search space scales combinatorially with additional parameters, and thus the I 2 and I 3 models do not add significant overhead to fitting the I 4 model.
As with the limitations of AIC analysis, it is not possible to state definitively the most appropriate physiological interpretation at all voxels, and it is inevitable some will be miscategorized by the framework. Nevertheless we propose it is beneficial to compute summary statistics over only those voxels assigned to each regime, and present these with the regime's percentage composition. The alternative, preselecting a model for each region, and then averaging over all voxels, will almost certainly include more voxels that do not meet the physiological assumptions of the model, with greater risk of producing invalid parameter summary statistics.
The difficulty in fitting f a to individual voxels in some healthy volunteer datasets, together with the results of varying temporal resolution in the fixed Monte Carlo simulations (Figure 3), suggest a temporal resolution of at least 6 seconds (and ideally higher) is necessary to adequately separate the arterial and venous components of vascular input to the liver. If accurate estimations of f a (and a ) are not required, a temporal resolution of 12 seconds appears sufficient to estimate all other parameters in the framework (above 12 seconds and estimations of several other parameters become unreliable).

| CONCLUSIONS
Although the individual tracer kinetic models used in this study are not new, we believe their presentation heregrouped by functional form, cross-matched with vascular input and then physiologically interpreted post-fitting on a per-voxel basis-is novel and reveals a general framework within which to understand the family of models. Extensive Monte Carlo simulations, and the application of the framework to real datasets show the value in this approach, and that attempting to model the dynamics of gadoxetate in HCC using a single model, or with fixed prior assumptions about which model will be appropriate in a specific region, may not adequately cope with the varying physiology of individual patients, and ultimately will produce unreliable tracerkinetic parameter estimates.
The unique physiology of HCC, together with varying levels of associated chronic disease in the surrounding liver, makes these tumors a challenging target for quantitative DCE imaging. We have shown that using gadoxetate as a contrast agent and applying a framework of tracer kinetic models, it is possible to investigate and quantify different aspects of tumor microvascular and general liver function in a free-breathing DCE protocol tolerated by patients. However, careful consideration should be given to sequence length and temporal resolution, to maximize the ability of model selection to differentiate between potentially ambiguous tissue types and increase the robustness of parameter estimation.