Computational Modeling of Interstitial Fluid Pressure and Velocity in Head and Neck Cancer Based on Dynamic Contrast-Enhanced Magnetic Resonance Imaging: Feasibility Analysis.

We developed and tested the feasibility of computational fluid modeling (CFM) based on dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) for quantitative estimation of interstitial fluid pressure (IFP) and velocity (IFV) in patients with head and neck (HN) cancer with locoregional lymph node metastases. Twenty-two patients with HN cancer, with 38 lymph nodes, underwent pretreatment standard MRI, including DCE-MRI, on a 3-Tesla scanner. CFM simulation was performed with the finite element method in COMSOL Multiphysics software. The model consisted of a partial differential equation (PDE) module to generate 3D parametric IFP and IFV maps, using the Darcy equation and K t r a n s values (min-1, estimated from the extended Tofts model) to reflect fluid influx into tissue from the capillary microvasculature. The Spearman correlation (ρ) was calculated between total tumor volumes and CFM estimates of mean tumor IFP and IFV. CFM-estimated tumor IFP and IFV mean ± standard deviation for the neck nodal metastases were 1.73 ± 0.39 (kPa) and 1.82 ± 0.9 × (10-7 m/s), respectively. High IFP estimates corresponds to very low IFV throughout the tumor core, but IFV rises rapidly near the tumor boundary where the drop in IFP is precipitous. A significant correlation was found between pretreatment total tumor volume and CFM estimates of mean tumor IFP (ρ = 0.50, P = 0.004). Future studies can validate these initial findings in larger patients with HN cancer cohorts using CFM of the tumor in concert with DCE characterization, which holds promise in radiation oncology and drug-therapy clinical trials.

Introduction

In solid tumors, interstitial fluid pressure (IFP) is elevated owing to increased permeability of abnormally formed tumor blood vessels and a lack of functional lymphatic removal pathways (1, 2). Elevated IFP serves to nullify the hydrostatic pressure differential between vasculature and the interstitium, hampering extravasation of drug into the tumor (3). Previous studies have shown that elevated IFP is known to have serious implications for the effective delivery of anticancer drugs (4–6).

Elevated tumor IFP has been invasively assessed in previous clinical and preclinical studies and is not routinely measured in clinical practice (7–10). The commonly accepted gold standard procedure for acute measurement of IFP is the wick-in-needle (WIN) method, which was established by Fadnes et al. (11) in rats. A previous report (12) attests its performance under different levels of hydration in dogs. In clinical settings, Gutmann et al. (8) observed high IFP values in head and neck (HN) tumors using the WIN technique (0.53–4.39 kPa [4–33 mmHg]). Use of these direct methods is limited to measuring IFP in locations where the tumor can be accessed easily. The development of a noninvasive method to estimate IFP as a surrogate biomarker in HN tumors, based on dynamic contrast-enhanced (DCE) MRI, can be a significant step forward if it can be incorporated into the early assessment of tumor response to therapy.

Multimodal therapy in various combinations of surgery, radiation, or chemotherapy is used in the management of HN cancers, depending on the primary site and TNM stage (13–15). For nasopharyngeal carcinoma (NPC), treatment options include neoadjuvant chemotherapy (NAC), as it has high incidence of distant metastatic disease (16). Concurrent chemoradiation is the main approach for nonsurgical treatment of locoregionally advanced HN cancers (17–19). RECIST v1.1 is the commonly used guideline tool to assess patient response to treatment and is limited by usage of bidimensional lesion measurements on anatomical imaging (20). It may be insightful to include positron emission tomography (PET) metabolic criteria (PERCIST) (21), which requires an additional fluorodeoxyglucose (18F-FDG) PET scan. MRI can provide both anatomical and functional results, which can be used to model tumor IFP estimates (22, 23).

DCE-MRI uses low-molecular-weight gadolinium-based contrast agents (CAs) to measure the influx of CA across the capillary wall (24). Pharmacokinetic modeling of the DCE data provides physiologically relevant biomarkers of perfusion and permeability for clinical endpoints such as prognosis and prediction of response in HN cancers (25, 26). For estimates of IFP in clinical HN settings, there is a need to develop computational fluid modeling (CFM) methods that use noninvasive techniques, such as DCE-MRI, and may help to understand the mechanisms of drug or radiation dose delivery in relation to tumor IFP (27), of particular relevance in regions not readily accessible by WIN.

CFM studies based on DCE-MRI, although in early stages of development, have introduced a noninvasive approach toward assessing IFP in preclinical settings (23, 28, 29) and in clinical studies on brain tumors (22, 30). The method implemented in the aforementioned studies all use , influx of CA, as an input in the Starling equation of fluid exchange accounting for both osmotic and hydrostatic pressure, and capillary hydraulic conductivity (23, 28, 31). The CFM assumes a porous medium to simulate interstitial fluid transport through tumor tissue extracellular space (pores). To the best of our knowledge, application of CFM to explore IFP in clinical HN cancers is yet unexplored, where the aforementioned modeling method has the potential to help treating physicians determine optimal radiation treatment dose for tumor control and/or assess early treatment response of novel targeted therapies, which are unmet challenges to date.

The purpose of the present clinical study is to develop and test the feasibility of quantitative CFM estimation of IFP and IFV from DCE-MRI in patient with HN cancer with locoregional lymph node (LN) metastases.

Methods

Patient Studies

The institutional review board approved and granted a waiver of informed consent for this retrospective clinical study, which was compliant with the Health Insurance Portability and Accountability Act. A total of n = 22 patients (male, 21; female, 1) with HN squamous cell carcinoma were enrolled between June 2014 and October 2015. Twenty patients were human papillomavirus–positive (HPV+) and 2 were HPV-negative (HPV−); all had neck nodal metastases with histologically proven squamous cell carcinoma. Patient characteristics are given in Table 1. Figure 1 shows the workflow pipeline steps for the CFM presented in this study.

Table 1.

Patient Characteristics

Characteristic	Value (%)
Number of patients	22
Number of LN	38
Mean age (years)	59
Range (years)	44−69
Sex
Male	21 (95)
Female	1 (5)
Location of Primary Tumor
Base of tongue	15 (68)
Oropharynx	5 (23)
Unknown	2 (9)
HPV Status
Positive	21 (95)
Negative	1 (5)

Figure 1.

Workflow for magnetic resonance imaging (MRI)-based computational fluid modeling (CFM) simulations. The patient undergoes a magnetic resonance (MR) examination with dynamic contrast-enhanced imaging. The images are contoured by a neuroradiologist on all slices to properly demarcate the 3-dimensional tumor structure. The images are then processed through extended Tofts model (ETM) using MRI-QAMPER software. A representative patient's individualized arterial input function (AIF) and tissue signals are plotted here. The ETM map is generated and incorporated into the simulation, along with the 3-dimensional mesh of the tumor region of interest (ROI). Finally, the CFM solves the dynamical equation to generate estimates of interstitial fluid pressure (IFP) in the tumor mesh domain.

DCE-MRI Data Acquisition

MRI studies were performed using a neurovascular phased-array coil on a 3-Tesla (T) scanner (Ingenia; Philips Healthcare; The Netherlands). The standard magnetic resonance (MR) acquisition parameters were as follows: multiplanar (axial, coronal, and sagittal) T2-weighted (T2w), fat-suppressed, fast spin-echo images (repetition time [TR] = 4000 milliseconds; echo time [TE] = 80 milliseconds; number of averages [NA] = 2; matrix = 256 × 256; slice thickness = 5 mm; field of view [FOV] = 20–24 cm), and multiplanar T1-weighted (T1w) images (TR = 600 milliseconds; TE = 8 milliseconds; NA = 2; slice thickness = 5 mm; matrix = 256 × 256; FOV = 20–24 cm). Standard T1w and T2w imaging were followed by multiple flip angle T1w imaging (to measure precontrast T1 [ie, T10]), and dynamic T1w imaging.

Precontrast T1w images were acquired using a fast spoiled gradient recalled (SPGR) echo sequence with acquisition parameters as follows: TR = 7 milliseconds; TE = 2.7 milliseconds; flip angles θ = 5°, 15°, 30°, and single-excitation; NA = 1, pixel bandwidth = 250 Hz/pixel, FOV = 20–24 cm2, 256 × 128 matrix zero-filled to 256 × 256 during image reconstruction.

Dynamic acquisition was performed with parameters identical to precontrast T1w imaging using a flip angle of 15° as mentioned above. Antecubital vein catheters delivered a bolus of 0.1 mmol/kg gadolinium-based clinically approved contrast at 2 mL/s, followed by saline flush using an MR-compatible programmable power injector (Spectris; Medrad, Indianola, PA) after acquiring the first 5–6 precontrast images. The entire node was covered contiguously with 5-mm-thick slices, zero gap, yielding 8–10 slices and 40 or 60 phases were acquired with ≤8-second temporal resolution. Total acquisition time was ≤8 minutes.

DCE-MRI Pharmacokinetic Analysis

The signal intensity in a voxel, measured from a SPGR T1w acquisition, is given by the following equation (32): where S(t) is the voxel signal intensity at time t, M0 is the equilibrium magnetization of the protons, θ is the flip angle, TR is repetition time and TE is echo time. R1(t) (R1 = 1/T1) and R2*(t) (R2* = 1/T*2) are the time courses of longitudinal and transverse relaxation rates, respectively.

For the R1(t) calculation, TE ≪ T2*, and so approximating , then Equation [1] becomes the following equation (33):

The solution for the relaxation rate R1 (t) follows from Equation [2] and is given by:

In the limit of fast water exchange, the change in water proton relaxation rate (ie, ΔR1 [s−1]) owing to CA relaxivity, r1 = 4.0 [(mM)−1s−1], is linearly related to the tissue CA concentration, Ct (mM) as follows: where R10 is the precontrast longitudinal relaxation rate.

The extended Tofts model (ETM) is based on a 2-compartment model (vascular space and extravascular extracelluar space [EES]). The ETM expression for modeling Ct(t) is given by the following equation (24): where, (min−1) is the volume transfer constant of CA, Cp(t) is the time-course of plasma CA concentration, which is called arterial input function (AIF). The rate constant describing CA back-flux into the vascular space from the EES is defined as kep ≡ /ve, where ve and vp are the volume fractions of the EES and blood plasma, respectively.

Individual AIF for each patient was determined from an algorithm that uses cross correlation to rate arterial voxel time series in comparison with an ideal AIF signal. Signal from arterial voxels with strongest similarity to ideal AIF was averaged to be used for each patient AIF (34). Regions of interest (ROIs) on tumors were contoured manually by an experienced neuroradiologist on a late phase of the T1w DCE images using ITK-SNAP (35). The extent of necrosis was evaluated by the same neuroradiologist on T2w and postcontrast T1w MRI, reflecting tumor heterogeneity. The total tumor volume was calculated on T2w images for the 38 neck LN metastases of the 22 patients.

T10 values were calculated from multiple flip angle precontrast T1w images (36). The time course of tissue concentration data, , was fitted (Equation [5]) using a nonlinear curve fitting technique that minimizes the sum of squared errors (SSEs) between model fit and data. The fitting procedure produces parametric maps of , ve, and vp, and parameter estimation values were bounded in the fitting routine as follows: ∈ [0, 5] (min−1), ve, and vp ∈ [0, 1]. Automation of AIF findings, quantitative pharmacokinetic analysis of DCE-MRI data, and generation of parametric maps were performed using in-house MRI-QAMPER software (Quantitative Analysis Multi-Parametric Evaluation Routines) approved by NCI Quantitative Imaging Network for Technical Benchmark (December 2019) written in MATLAB (The MathWorks, Inc., Natick, MA), as detailed in Paudyal et al. (37).

CFM Mathematical Model

The fluid mechanics of the system are given by the Navier–Stokes hydrodynamic equation. The mass-balance describing the continuity of flow (1/s) in and out of interstitial space with source and sink conditions gives a generalized description of fluid velocity from forces acting on, within, and between mass elements (38).

Where is the fluid velocity vector (m/s), is the (interstitial) fluid pressure (Pa), is the permeability of the porous membrane (m2, describing the packing of the extracellular matrix), is the interstitial fluid density (kg/m3), is interstitial fluid viscosity (Pa*s), additional forces acting on the fluid volume element may be represented as (N), and I and T are the notation for the identity matrix and the transpose of the operator, respectively.

Fluid movement through space is approximated with low Reynolds number flow (39) and modeled under assumption of steady-state velocity. Neglecting both friction within fluid and exchange of momentum between fluid and solid phases, the case of an incompressible fluid () in the Navier–Stokes equation (Equation [6]) simplifies to Darcy's law where interstitial fluid velocity (IFV) is related to the gradient in IFP ():

This Darcy velocity describes bulk fluid movement within the interstitial space and expresses velocity as the gradient of instantaneous pressure, proportional to the local hydraulic conductivity KH≡:

Fluid source-sink assumptions are as follows: fluid can enter the interstitium via the vascular compartment (); fluid is drained from the system via lymphatic pathways () (1).

Blood flow transport across capillary walls was regulated according to the Starling law (1, 31):

Where is the hydraulic conductivity of the capillary wall, or vessel permeability, S/V is microvascular surface area per unit volume, is the blood pressure in the microvessel, is IFP, is osmotic pressure in microvasculature, is osmotic pressure in interstitial space, and is the osmotic reflection coefficient.

Lymphatic drainage of fluid is assumed to be available in normal tissue outside of the metastatic LN (40, 41), dependent on the difference between local interstitial pressure,and pressure of the lymphatic system, and the lymphatic clearance rate, :

In the subsequent sensitivity analysis, we observe the change in IFP profile if we create a simulated sink, or removal, term in the tumor.

Combining Equations [9], [10], and [11], the continuity equation (becomes):

Numerical values for physical parameters were based on the default selections from prior literature, listed with references in Table 2.

Table 2.

Tissue and Vascular Parameters Used in Simulations

Parameter	Unit	Description (# References)	Value
Lp	m Pa−1 s−1	Vessel permeability (1, 23)	2 × 10−11 (tumor)3 × 10−12 (normal)
LpLSL/V	Pa−1 s−1	Lymphatic filtration coefficient (23)	1 × 10−7
K	m2 Pa−1 s−1	Hydraulic conductivity (23)	1.9 × 10−12 (tumor)3.8 × 10−13 (normal)
S/V	m−1	Microvascular surface area per unit volume (1, 47)	2 × 104 (tumor)7 × 103 (normal)
pV	Pa	Microvascular pressure (59)	2,300
πi	Pa	Osmotic pressure in interstitial space (1, 46, 47)	3,230 (tumor)1,330 (normal)
πV	Pa	Osmotic pressure in microvasculature (46)	2,670
σT	Unitless	Average osmotic reflection coefficient for plasma (1, 50)	0.82 (tumor)0.91 (normal)

Computational Fluid Modeling

As mentioned above, the influx of fluid through the capillary wall to the interstitial space is described by the continuity equation, consisting of the source and sink terms. The ETM-derived is incorporated into the Starling equation (Equation [11]), and scaled by the mean tumor value, <> to account for the heterogeneous bulk vascular influx (23, 28). The final expression for the continuity equation (Equation [12]) in terms of the dependent variable interstitial pressure, :

The continuity partial differential equation (PDE) (Equation [13]) was implemented using the COMSOL CFM simulation PDE module. Solving (Equation [13]) provides the basis for estimation of pi and 3D parametric maps of IFP and IFV.

The physiological 3D mesh model was generated from each patient's T1w DCE tumor images. The contoured tumor ROIs in the model were converted to binary masks and extended by 10 pixels to represent normal tissue around tumor. The ROIs for tumor with normal surrounding tissue were resliced to be 1 mm3 isotropic in MATLAB using the NIfTI toolbox (42), and converted to stereolithography (STL) file format. STL files were imported into the simulation software and interpreted as boundary meshes for the model. In the present study, a simplified geometry was used to represent a complex LN tumor structure. On imaging, intranodal tumor necrosis is observed, which is a mix of tumor necrosis, keratin pooling, fibrous tissue, edema, viable tumor cells, and possibly hemorrhage (43). Usually, a localized group of nodes is present in an expected nodal draining area for a specific primary tumor (43). Once extracapsular extension has occurred, the LN tumor can extend to invade surrounding tissue (43, 44).

ETM-estimated maps were resliced to matching space to ensure accurate spatial interpolation (45), and converted from units of min−1 to s−1 to carry out calculations in standard SI units; the values were interpreted in COMSOL as a scalar field over the mesh domain. The field values were normalized by the mean value in the tumor ROI (23). Numerical values for physical constants in normal and tumor tissue were respectively assigned to the appropriate regions of the 3D STL domain mesh, as given in Table 2. These physiological values are taken from those previously cited in literature (1, 23, 46, 47). To generate IFP maps, a stationary solution of the PDE in Equation [13] was computed on the 3D extended domain ROI, and a no-flux condition was implemented at mesh boundary edges.

Simulation was conducted using the general coefficient form PDE module in a commercial multiphysics software package (COMSOL Inc., Stockholm, Sweden) which uses finite element method to solve PDE equations.

IFP and IFV Sensitivity Analysis

A sensitivity analysis was performed by parameter variation in the CFM models to test the influence on IFP and IFV estimates as detailed in a preclinical setting (23). The parameter values in the CFM were sequentially varied in separate simulation runs: vessel permeability (in normal and tumor tissue, Lp), lymphatic clearance (LpLSL/V in normal tissue and tumor), and ratio of hydraulic conductivity (allowing KH,tumor to vary relative to KH,normal). The range of high and low values was chosen to scale 2× and 0.5× relative to the baseline value (Table 2). Values were assigned ([LpLSL/V]t/[LpLSL/V]n,baseline = 0.125 and 0.5) to study the effect of fluid clearance from the tumor region (either via lymphatics or vascular reabsorption).

Statistical Analysis

Spearman correlation coefficient (ρ) was calculated between the total tumor volume and CFM estimates of IFP and IFV. The results were reported as mean ± standard deviation. In addition, heterogeneity measures from the histogram were calculated as skewness and kurtosis. A P-value < .05 was considered statistically significant.

Results

All 22 patients with HN cancer underwent MRI pretreatment (pre-TX), and 38 neck LN metastases were analyzed. Out of the 22 patients, 10 patients had ≥2 metastatic nodes. Figure 2 shows the acquired MR images of a representative patient (patient #1) with a right lateral LN (total tumor volume = 14.26 cm3).

Figure 2.

T1-weighted postcontrast image of patient #1 with right lateral neck nodal metastasis (yellow outline) (A); T1-weighted (T1w) postcontrast image with overlaid ROI for analysis (green denotes normal tissue, red denotes viable tumor tissue, blue denotes necrotic core) (B); map of tumor ROI (C); preview of the imported geometry mesh of the tumor (blue) for patient #1 for CFM (D).

The CFM estimation maps of IFP and IFV are shown in Figure 3 for representative patients #1, 2, 3 with and 4 with HN cancer. The presence of cystic necrosis for the nodes in patients #1 and 4 was evidenced by hypointense signal in postcontrast T1w images and hyperintensity in the T2w images. In contrast, the nodes of patients #2 and 3 did not exhibit necrosis on imaging. The subtle differences in IFP and IFV profiles are evident for the necrotic tumors.

Figure 3.

Visualization of IFP and interstitial fluid volume (IFV) maps as estimated by simulation in patient #1, 2, 3, and 4. The profiles were calculated along the lines drawn on the corresponding maps. Row 1: T1w postcontrast image (inset: T2w view of tumor) (V1 = 14.26 cm3, V2 = 25.20 cm3, V3 = 31.27 cm3, V4 = 28.80 cm3); Row 2: Estimated IFP map of tumor and surrounding normal tissue (mean tumor pressure = 1.73 kPa, = 2.05 kPa, = 2.09 kPa, = 1.77 kPa); Row 3: IFP profile along a vertical bisector line (tumor boundary denoted by dashed line); Row 4: Estimated IFV map of tumor and surrounding normal tissue (mean tumor velocity = 1.60 × 10−7m/s, = 1.80 × 10−7m/s, = 1.67 × 10−7m/s, = 1.51 × 10−7m/s); Row 5: IFV profile along vertical bisector line, with tumor boundary (dash).

In Figure 3, mean IFP (pmean) for the cystic node in patient #1 was 1.73 kPa with maximum IFP at tumor core pmax = 2.23 kPa; the IFP profile demonstrated slowly changing pressure with little curvature throughout the tumor core. The resulting IFV profile for patient #1 (IFV gradient of IFP, by the Darcy law) results in negligible velocity throughout the tumor core. The estimated pmean and pmax for the metastatic node in patient #2 was 2.05 kPa, and 2.64 kPa, respectively. Values of pmean and pmax for patients #3 and 4 were 2.09, 1.65 kPa (pmean) and 2.67, 2.04 kPa (pmax), respectively. In the IFP profiles of patients #2 and 3 with non-necrotic tumors, a moderate curvature is seen, and the resulting IFV is more variable and elevated in the tumor core.

Table 3 summarizes the mean(± standard deviation) values for total tumor volume (cm3), (min−1), IFP (kPa), and IFV (m/s) obtained from the LN metastases. The mean total tumor volume for the 38 nodes was V = 15.83 ±10.80 cm3 ranging from 1.36 to 40.02 cm3.

Table 3.

Summary of DCE-MRI Ktrans- and CFM-estimated IFP and IFV from 22 Patients with Head and Neck Squamous Cell Carcinoma

Parameter (Unit)	Value (Mean ± SD)
Total Tumor Volume (cm3)	15.83 ± 10.80
Ktrans (min−1)	0.02 ± 0.02
Interstitial Fluid Pressure(IFP [kPa])	1.73 ± 0.39
IFP Skewness	−0.48
IFP Kurtosis	3.10
Interstitial Fluid Velocity(IFV [×10−7 m/s])	1.82 ± 0.89
IFV Skewness	0.79
IFV Kurtosis	3.84

Figure 4 exhibits the scatter plot of pre-TX total tumor volume versus mean CFM-estimated IFP. Significant correlation (Spearman) was found between pre-TX total tumor volume and the mean IFP (ρ = 0.5, P = .004). No significant correlation was found between mean tumor IFV and total tumor volume (ρ = −0.1, P = .5).

Figure 4.

Plot of total tumor volume, calculated from T2-weighted image, and mean tumor IFP as estimated by COMSOL simulation, for all pretreatment neck lymph node (LN) metastases. A significant correlation between tumor volume and the mean intratumor pressure is found (ρ = 0.5, P = .004).

Figure 5 shows the of IFP and IFV profiles from the sensitivity analysis as carried out for representative patients #1 and 2. The parameter with the greatest effect on amplitude of IFP values was the tumor vessel permeability (Lp). For the necrotic tumor in patient #1, changes to parameters affected the magnitude of IFP but the shape and slope of the IFP field did not vary greatly and hence had little influence on the internal IFV. The non-necrotic node in patient #2 showed greater IFP sensitivity to all changes in parameters, with more variability in the magnitude of predicted IFV, and maintaining IFV appreciably higher.

Figure 5.

Sensitivity analyses of IFP and IFV profiles for patients #1 (A, B) and #2 (C, D) as a result of varying model parameters; in each case, the high value is 2× the default parameter (as listed in Table 2), and the low value is 0.5× the default. Subscript t, n denote values in tumor and normal, respectively.

Discussion

This is the first feasibility study to apply CFM using DCE data to estimate IFP/IFV in HN LN metastases. A 3D mesh tumor ROI, derived from the contours on DCE imaging, and the continuity partial differential equation (to) describe the physiological processes (Starling law and Darcy law), were combined and solved with the finite element method (COMSOL Multiphysics computing environment) using the 3D values. The voxel values in the CFM were normalized by the mean values in each tumor ROI. CFM-generated IFP and IFV maps depict heterogeneity within the tumor, stemming from incorporation of to inform fluid influx rate variability from the spatial CA extravasation, similar to Pishko et al. (23).

The CFM-estimated mean tumor IFP showed a moderate but significant correlation with the total tumor volume at pre-TX. This finding is in good agreement with previous clinical studies in HN, cervical, and breast cancer that directly measured IFP by using an invasive WIN method (8, 48, 49). Across different tumor types, the mean IFP value may vary but the relationship between IFP and tumor volume still holds (8, 48, 49). IFP is elevated within the tumor and decays sharply at its boundary, and in accordance with the Darcy law, IFV attains a maximum value at the tumor boundary and movement is directed outward from the tumor (1, 50, 51). Most of the high-velocity fluid movement exists near the periphery of the tumor owing to a precipitous fall-off of pressure (52). Tumors with mild-to-severe necrotic fractions in their interior, as in patients #1 (severe) and 4 (mild), tend to quickly reach maximum IFP, starting off a gradual plateau and varying less inside the tumor compared with patients #2 and 3, which attained a central local maximum of pressure within the core (Figure 3). The physical impact of cystic necrosis on the IFV is shown in these contrasting profiles. The variability of IFP in the non-necrotic tumor in contrast to necrotic tumor may lead to a significant change in IFV of tumor core.

Sensitivity analyses illustrated the influence that tissue-descriptive parameters, such as hydraulic conductivity, KH, or the effectiveness of lymphatic vessels, have on estimates of IFP and IFV in the CFM. The analysis also assesses the stability of the results, in spite of estimated values from literature. In the core of the necrotic node, IFV was not changed from the baseline despite changes to the model parameters. In contrast, the IFV in the node with a large viable fraction showed greater tendency for variation depending on parameter values. Differing from the conditions of the main analysis, we considered lymphatic functionality for the tumor core. Tumor drainage introduces a distinct change in curvature of both patients' IFP profiles and is the only factor that introduces variability to IFV values in the core of the necrotic tumor. In the patient with necrotic tumor, increasing hydraulic conductivity ratio KH,t/KH,n has no noticeable effect from baseline value, but in the non-necrotic case, a lower-ratio KH,t/KH,n results in lower overall IFP with the greatest IFV compared with all other parameter changes. In contrast, for the non-necrotic tumor, a lower lymphatic clearance in normal tissue has no impact on IFP/IFV from baseline, but it effectively increases the IFP in the necrotic tumor. In general, any fluid delivered into a tumor will be mostly reliant on diffusive transport or be trapped owing to the approximately zero internal IFV. Intravenous drug delivery effectiveness may be affected by the presence of necrotic core and clearance mechanisms therein.

The effect of mesh resolution or partitioning size on resulting finite element method simulations was shown to be minimal. CFM studies in Pishko et al. (23) and Magdoom et al. (53) reanalyzed the same preclinical DCE data set using different quantization schemes for the mesh; their result shows the resulting pressure maps were comparable. Although the absolute IFP magnitude was affected by differences in spatial discretization due to changes in voxel size, the overall profile trends in IFP maps were not changed.

Known limitations of this study include literature-based parameter values in tumor and normal tissue (23), a shortcoming of all similarly constructed IFP simulations (1, 54–57). Physiologically accurate approaches in the future would aim to include additional imaging-based biomarkers to characterize the model parameters. A B1 nonuniformity acquisition is required for high-field DCE-MRI. The accuracy of ETM-derived values would be improved by inclusion of a B1 correction (58). Individual AIF measurements for large patient cohort studies can be a challenge, which can be addressed by using an average AIF (34). Lastly, the repeatability of the metric and comparison to a gold standard WIN IFP measurement needs to be assessed before the IFP derived from CFM method can be established as a prognostic/predictive biomarker.

In conclusion, our initial findings establish that DCE-MRI-based CFM can noninvasively provide estimates of tumor IFP and IFV in neck nodal metastases. After validation, IFP estimates from CFM may help treating physicians better understand and plan optimal radiation treatment dose for tumor control, and/or assess early treatment response of novel targeted therapies.