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 wickinneedle (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 contrastenhanced (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 (^{18}FFDG) PET scan. MRI can provide both anatomical and functional results, which can be used to model tumor IFP estimates (22, 23).
DCEMRI uses lowmolecularweight gadoliniumbased 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 DCEMRI, 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 DCEMRI, 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 Ktran, 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 DCEMRI 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 HPVnegative (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.
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 contrastenhanced imaging. The images are contoured by a neuroradiologist on all slices to properly demarcate the 3dimensional tumor structure. The images are then processed through extended Tofts model (ETM) using MRIQAMPER software. A representative patient's individualized arterial input function (AIF) and tissue signals are plotted here. The ETM Ktran map is generated and incorporated into the simulation, along with the 3dimensional 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.
DCEMRI Data Acquisition
MRI studies were performed using a neurovascular phasedarray coil on a 3Tesla (T) scanner (Ingenia; Philips Healthcare; The Netherlands). The standard magnetic resonance (MR) acquisition parameters were as follows: multiplanar (axial, coronal, and sagittal) T2weighted (T2w), fatsuppressed, fast spinecho 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 T1weighted (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, T_{10}]), 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 singleexcitation; NA = 1, pixel bandwidth = 250 Hz/pixel, FOV = 20–24 cm^{2}, 256 × 128 matrix zerofilled 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 gadoliniumbased clinically approved contrast at 2 mL/s, followed by saline flush using an MRcompatible programmable power injector (Spectris; Medrad, Indianola, PA) after acquiring the first 5–6 precontrast images. The entire node was covered contiguously with 5mmthick slices, zero gap, yielding 8–10 slices and 40 or 60 phases were acquired with ≤8second temporal resolution. Total acquisition time was ≤8 minutes.
DCEMRI 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,
M_{0} is the equilibrium magnetization of the protons, θ is the flip angle, TR is repetition time and TE is echo time.
R_{1}(t) (
R_{1} = 1/
T_{1}) and
R_{2}*(
t) (
R_{2}* = 1/
T*
_{2}) are the time courses of longitudinal and transverse relaxation rates, respectively.
For the R_{1}(t) calculation, TE ≪ T_{2}*, and so approximating eTER2*(t)≈1, then Equation [1] becomes the following equation (33):
The solution for the relaxation rate R_{1} (t) follows from Equation [2] and is given by:
In the limit of fast water exchange, the change in water proton relaxation rate (ie, ΔR_{1} =(R1R10) [s^{−1}]) owing to CA relaxivity, r_{1} = 4.0 [(mM)^{−1}^{−1}]_{,} is linearly related to the tissue CA concentration, C_{t} (mM) as follows:
where
R_{10} is the precontrast longitudinal relaxation rate.
The extended Tofts model (ETM) is based on a 2compartment model (vascular space and extravascular extracelluar space [EES]). The ETM expression for modeling C_{t}(t) is given by the following equation (24):
where,
Ktran (min
^{−1}) is the volume transfer constant of CA,
C_{p}(t) is the timecourse of plasma CA concentration, which is called arterial input function (AIF). The rate constant describing CA backflux into the vascular space from the EES is defined as
k_{ep} ≡
Ktran/
v_{e,} where
v_{e} and
v_{p} 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 ITKSNAP (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.
T_{10} values were calculated from multiple flip angle precontrast T1w images (36). The time course of tissue concentration data, Ct(t), 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 Ktran, v_{e}, and v_{p}, and parameter estimation values were bounded in the fitting routine as follows: Ktran ∈ [0, 5] (min^{−1}), v_{e}, and v_{p} ∈ [0, 1]. Automation of AIF findings, quantitative pharmacokinetic analysis of DCEMRI data, and generation of parametric maps were performed using inhouse MRIQAMPER software (Quantitative Analysis MultiParametric 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 massbalance 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 u is the fluid velocity vector (m/s), pi is the (interstitial) fluid pressure (Pa), K is the permeability of the porous membrane (m^{2}, describing the packing of the extracellular matrix), ρ is the interstitial fluid density (kg/m^{3}), μ is interstitial fluid viscosity (Pa*s), additional forces acting on the fluid volume element may be represented as F (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 steadystate velocity. Neglecting both friction within fluid and exchange of momentum between fluid and solid phases, the case of an incompressible fluid (∇·u=0) in the Navier–Stokes equation (Equation [6]) simplifies to Darcy's law where interstitial fluid velocity (IFV) is related to the gradient in IFP (∇pi):
This Darcy velocity u describes bulk fluid movement within the interstitial space and expresses velocity as the gradient of instantaneous pressure, proportional to the local hydraulic conductivity K_{H}≡Kμ:
Fluid sourcesink assumptions are as follows: fluid can enter the interstitium via the vascular compartment (ϕv); fluid is drained from the system via lymphatic pathways (ϕL) (1).
Blood flow transport ϕv across capillary walls was regulated according to the Starling law (1, 31):
Where LP is the hydraulic conductivity of the capillary wall, or vessel permeability, S/V is microvascular surface area per unit volume,pV is the blood pressure in the microvessel, pi is IFP, πV is osmotic pressure in microvasculature, πi is osmotic pressure in interstitial space, and σT 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,pi,and pressure of the lymphatic system, pL and the lymphatic clearance rate, LpLSLV:
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 
L_{p}

m Pa^{−1} s^{−1}

Vessel permeability (1, 23) 
2 × 10^{−11} (tumor) 3 × 10^{−12} (normal) 
L_{pL}S_{L}/V

Pa^{−1} s^{−1}

Lymphatic filtration coefficient (23) 
1 × 10^{−7}

K

m^{2} 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 × 10^{4} (tumor) 7 × 10^{3} (normal) 
p_{V}

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 ETMderived Ktran is incorporated into the Starling equation (Equation [11]), and scaled by the mean tumor value, <Ktran> 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, pi:
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 p_{i} 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 mm^{3} 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).
ETMestimated Ktran 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 Ktran values were interpreted in COMSOL as a scalar field over the mesh domain. The Ktran field values were normalized by the mean Ktran 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 noflux 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, L_{p}), lymphatic clearance (L_{pL}S_{L}/V in normal tissue and tumor), and ratio of hydraulic conductivity (allowing K_{H,tumor} to vary relative to K_{H,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 ([L_{pL}S_{L}/V]_{t}/[L_{pL}S_{L}/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 Pvalue < .05 was considered statistically significant.
Results
All 22 patients with HN cancer underwent MRI pretreatment (preTX), 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 cm^{3}).
Figure 2.
T1weighted postcontrast image of patient #1 with right lateral neck nodal metastasis (yellow outline) (A); T1weighted (T1w) postcontrast image with overlaid ROI for analysis (green denotes normal tissue, red denotes viable tumor tissue, blue denotes necrotic core) (B); Ktran 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) (V_{1} = 14.26 cm^{3}, V_{2} = 25.20 cm^{3}, V_{3} = 31.27 cm^{3}, V_{4} = 28.80 cm^{3}); Row 2: Estimated IFP map of tumor and surrounding normal tissue (mean tumor pressure P¯1= 1.73 kPa, P¯2 = 2.05 kPa, P¯3 = 2.09 kPa,P¯4 = 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 u¯1 = 1.60 × 10^{−7}m/s, u¯2 = 1.80 × 10^{−7}m/s, u¯3 = 1.67 × 10^{−7}m/s, u¯4 = 1.51 × 10^{−7}m/s); Row 5: IFV profile along vertical bisector line, with tumor boundary (dash).
In Figure 3, mean IFP (p_{mean}) for the cystic node in patient #1 was 1.73 kPa with maximum IFP at tumor core p_{max} = 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 p_{mean} and p_{max} for the metastatic node in patient #2 was 2.05 kPa, and 2.64 kPa, respectively. Values of p_{mean} and p_{max} for patients #3 and 4 were 2.09, 1.65 kPa (p_{mean}) and 2.67, 2.04 kPa (p_{max}), respectively. In the IFP profiles of patients #2 and 3 with nonnecrotic 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 (cm^{3}), Ktran (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 cm^{3} ranging from 1.36 to 40.02 cm^{3}.
Table 3.
Summary of DCEMRI K^{trans} and CFMestimated IFP and IFV from 22 Patients with Head and Neck Squamous Cell Carcinoma
Parameter (Unit) 
Value (Mean ± SD) 
Total Tumor Volume (cm^{3}) 
15.83 ± 10.80 
K^{trans} (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 preTX total tumor volume versus mean CFMestimated IFP. Significant correlation (Spearman) was found between preTX 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 T2weighted 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 (L_{p}). 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 nonnecrotic 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 Ktran values. The Ktranvoxel values in the CFM were normalized by the mean Ktran values in each tumor ROI. CFMgenerated IFP and IFV maps depict heterogeneity within the tumor, stemming from incorporation of Ktran to inform fluid influx rate variability from the spatial CA extravasation, similar to Pishko et al. (23).
The CFMestimated mean tumor IFP showed a moderate but significant correlation with the total tumor volume at preTX. 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 highvelocity fluid movement exists near the periphery of the tumor owing to a precipitous falloff of pressure (52). Tumors with mildtosevere 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 nonnecrotic tumor in contrast to necrotic tumor may lead to a significant change in IFV of tumor core.
Sensitivity analyses illustrated the influence that tissuedescriptive parameters, such as hydraulic conductivity, K_{H}, 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 K_{H,t}/K_{H,n} has no noticeable effect from baseline value, but in the nonnecrotic case, a lowerratio K_{H,t}/K_{H,n} results in lower overall IFP with the greatest IFV compared with all other parameter changes. In contrast, for the nonnecrotic 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 literaturebased 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 imagingbased biomarkers to characterize the model parameters. A B1 nonuniformity acquisition is required for highfield DCEMRI. The accuracy of ETMderived Ktran 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 DCEMRIbased 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.
Acknowledgments
Equal contribution: E.L. and R.P. contributed equally to this work.
This work was supported by U01 CA211205 and in part through the NIH/NCI Cancer Center Support Grant P30 CA008748. We would like to thank the MRI technologists for their great efforts in helping to perform the MRI examinations, and Mr. Christian Czmielewski for his helpful contribution to patient enrollment and data management. This work was supported by the MSKCC internal IMRAS grant, U01 CA211205, and in part through the NIH/NCI Cancer Center Support Grant: P30 CA008748.
References

Baxter LT, Jain RK. Transport of fluid and macromolecules in tumors. I. Role of interstitial pressure and convection. Microvasc Res. 1989;37:77–104.

Leu AJ, Berk DA, Lymboussaki A, Alitalo K, Jain RK. Absence of functional lymphatics within a murine sarcoma: a molecular and functional evaluation. Cancer Res. 2000;60:4324–4327.

Jain RK. Delivery of molecular medicine to solid tumors. Science. 1996;271:1079–1080.

Jain RK. Barriers to drug delivery in solid tumors. Sci Am. 1994;271:58–65.

Heldin CH, Rubin K, Pietras K, Ostman A. High interstitial fluid pressure  an obstacle in cancer therapy. Nat Rev Cancer. 2004;4:806–813.

Lunt SJ, Fyles A, Hill RP, Milosevic M. Interstitial fluid pressure in tumors: therapeutic barrier and biomarker of angiogenesis. Future Oncol. 2008;4:793–802.

Boucher Y, Kirkwood JM, Opacic D, Desantis M, Jain RK. Interstitial hypertension in superficial metastatic melanomas in humans. Cancer Res. 1991;51:6691–6694.

Gutmann R, Leunig M, Feyh J, Goetz AE, Messmer K, Kastenbauer E, Jain RK. Interstitial hypertension in head and neck tumors in patients: correlation with tumor size. Cancer Res. 1992;52:1993–1995.

Boucher Y, Salehi H, Witwer B, Harsh GR, Jain RK. Interstitial fluid pressure in intracranial tumours in patients and in rodents. Br J Cancer. 1997;75:829–836.

Taghian AG, AbiRaad R, Assaad SI, Casty A, Ancukiewicz M, Yeh E, Molokhia P, Attia K, Sullivan T, Kuter I, Boucher Y, Powell SN. Paclitaxel decreases the interstitial fluid pressure and improves oxygenation in breast cancers in patients treated with neoadjuvant chemotherapy: clinical implications. J Clin Oncol. 2005;23:1951–1961.

Fadnes HO, Reed RK, Aukland K. Interstitial fluid pressure in rats measured with a modified wick technique. Microvasc Res. 1977;14:27–36.

Wiig H, Reed RK, Aukland K. Measurement of interstitial fluid pressure in dogs: evaluation of methods. Am J Physiol. 1987;253(2 Pt 2):H283–90.

Argiris A, Li S, Savvides P, Ohr JP, Gilbert J, Levine MA, Chakravarti A, Haigentz Jr M, Saba NF, Ikpeazu CV, Schneider CJ, Pinto HA, Forastiere AA, Burtness B. Phase III randomized trial of chemotherapy with or without bevacizumab in patients with recurrent or metastatic head and neck cancer. J Clin Oncol. 2019;37:3266–3274.

Beckham TH, Barney C, Healy E, Wolfe AR, Branstetter A, Yaney A, Riaz N, McBride SM, Tsai CJ, Kang J, Yu Y, Chen L, Sherman E, Dunn L, Pfister DG, Tan J, Rupert R, Bonomi M, Zhang Z, Lobaugh SM, Grecula JC, Mitchell DL, Wobb JL, Miller ED, Blakaj DM, Diavolitsis VM, Lee N, Bhatt AD. Platinumbased regimens versus cetuximab in definitive chemoradiation for human papillomavirusunrelated head and neck cancer. Int J Cancer. 2019.

Rosenthal DI, Istenmaa DA, Glatstein E. A review of neoadjuvant chemotherapy for head and neck cancer: partially shrunken tumors may be both leaner and meaner. Int J Radiat Oncol Biol Phys. 1994;28:315–320.

Han JE, Yi SK, Wang S, Erman A, Bearelly S, Sindhu S, Robbins JR, Bauman J, Hsu CC. Neoadjuvant chemotherapy improves survival compared with concurrent chemoradiation alone in nasopharyngeal carcinoma patients with N3 disease. Head Neck. 2019;41:4076–4087.

Ang KK, Harris J, Garden AS, Trotti A, Jones CU, Carrascosa L, Cheng JD, Spencer SS, Forastiere A, Weber RS. Concomitant boost radiation plus concurrent cisplatin for advanced head and neck carcinomas: radiation therapy oncology group phase II trial 9914. J Clin Oncol. 2005;23:3008–3015.

Dilling TJ, Bae K, Paulus R, WatkinsBruner D, Garden AS, Forastiere A, Kian Ang K, Movsas B. Impact of gender, partner status, and race on locoregional failure and overall survival in head and neck cancer patients in three radiation therapy oncology group trials. Int J Radiat Oncol Biol Phys. 2011;81:e101–9.

Setton J, Caria N, Romanyshyn J, Koutcher L, Wolden SL, Zelefsky MJ, Rowan N, Sherman EJ, Fury MG, Pfister DG, Wong RJ, Shah JP, Kraus DH, Shi W, Zhang Z, Schupak KD, Gelblum DY, Rao SD, Lee NY. Intensitymodulated radiotherapy in the treatment of oropharyngeal cancer: an update of the Memorial SloanKettering Cancer Center experience. Int J Radiat Oncol Biol Phys. 2012;82:291–298.

Eisenhauer EA, Therasse P, Bogaerts J, Schwartz LH, Sargent D, Ford R, Dancey J, Arbuck S, Gwyther S, Mooney M, Rubinstein L, Shankar L, Dodd L, Kaplan R, Lacombe D, Verweij J. New response evaluation criteria in solid tumours: revised RECIST guideline (version 1.1). Eur J Cancer. 2009;45:228–247.

Wahl RL, Jacene H, Kasamon Y, Lodge MA. From RECIST to PERCIST: evolving considerations for PET response criteria in solid tumors. J Nucl Med. 2009;50:122S–150S.

Bhandari A, Bansal A, Singh A, Sinha N. Perfusion kinetics in human brain tumor with DCEMRI derived model and CFD analysis. J Biomech. 2017;59:80–89.

Pishko GL, Astary GW, Mareci TH, Sarntinoranont M. Sensitivity analysis of an imagebased solid tumor computational model with heterogeneous vasculature and porosity. Ann Biomed Eng. 2011;39:2360–2373.

Tofts PS, Brix G, Buckley DL, Evelhoch JL, Henderson E, Knopp MV, Larsson HBW, Lee TY, Mayr NA, Parker GJM, Port RE, Taylor J, Weisskoff RM. Estimating kinetic parameters from dynamic contrastenhanced T(1)weighted MRI of a diffusable tracer: standardized quantities and symbols. J Magn Reson Imaging. 1999;10:223–232.

ShuklaDave A, Lee NY, Jansen JFA, Thaler HT, Stambuk HE, Fury MG, Patel SG, Moreira AL, Sherman E, Karimi S, Wang Y, Kraus D, Shah JP, Pfister DG, Koutcher JA. Dynamic contrastenhanced magnetic resonance imaging as a predictor of outcome in headandneck squamous cell carcinoma patients with nodal metastases. Int J Radiat Oncol Biol Phys. 2012;82:1837–1844.

Kim S, Loevner LA, Quon H, Kilger A, Sherman E, Weinstein G, Chalian A, Poptani H. Prediction of response to chemoradiation therapy in squamous cell carcinomas of the head and neck using dynamic contrastenhanced MR imaging. AJNR Am J Neuroradiol. 2010;31:262–268.

Lee N, Schoder H, Beattie B, Lanning R, Riaz N, McBride S, Katabi N, Li D, Yarusi B, Chan S, Mitrani L, Zhang Z, Pfister DG, Sherman E, Baxi S, Boyle J, Morris LGT, Ganly I, Wong R, Humm J. Strategy of using intratreatment hypoxia imaging to selectively and safely guide radiation dose deescalation concurrent with chemotherapy for locoregionally advanced human papillomavirusrelated oropharyngeal carcinoma. Int J Radiat Oncol Biol Phys. 2016;96:9–17.

Zhao J, Salmon H, Sarntinoranont M. Effect of heterogeneous vasculature on interstitial transport within a solid tumor. Microvasc Res. 2007;73:224–236.

Steuperaert M, Debbaut C, Carlier C, De Wever O, Descamps B, Vanhove C, Ceelen W, Segers P. A 3D CFD model of the interstitial fluid pressure and drug distribution in heterogeneous tumor nodules during intraperitoneal chemotherapy. Drug Deliv. 2019;26:404–415.

Bhandari A, Bansal A, Singh A, Gupta RK, Sinha N. Comparison of transport of chemotherapeutic drugs in voxelized heterogeneous model of human brain tumor. Microvasc Res. 2019;124:76–90.

Starling EH. On the absorption of fluids from the connective tissue spaces. J Physiol. 1896;19:312–326. sp000596.

BagherEbadian H, Jain R, NejadDavarani SP, Mikkelsen T, Lu M, Jiang Q, Scarpace L, Arbab AS, Narang J, SoltanianZadeh H, Paudyal R, Ewing JR. Model selection for DCET1 studies in glioblastoma. Magn Reson Med. 2012;68:241–251.

Aryal MP, Nagaraja TN, Brown SL, Lu M, BagherEbadian H, Ding G, Panda S, Keenan K, Cabral G, Mikkelsen T, Ewing JR. Intratumor distribution and testretest comparisons of physiological parameters quantified by dynamic contrastenhanced MRI in rat U251 glioma. NMR Biomed. 2014;27:1230–1238.

ShuklaDave A, Lee N, Stambuk H, Wang Y, Huang W, Thaler HT, Patel SG, Shah JP, Koutcher JA. Average arterial input function for quantitative dynamic contrast enhanced magnetic resonance imaging of neck nodal metastases. BMC Med Phys. 2009;9:4.

Yushkevich PA, Piven J, Hazlett HC, Smith RG, Ho S, Gee JC, Gerig G. Userguided 3D active contour segmentation of anatomical structures: significantly improved efficiency and reliability. Neuroimage. 2006;31:1116–1128.

Deoni SC, Josseau MJ, Rutt BK, Peters TM. Visualization of thalamic nuclei on high resolution, multiaveraged T1 and T2 maps acquired at 1.5 T. Hum Brain Mapp. 2005;25:353–359.

Paudyal R, Lu Y, Hatzoglou V, Moreira A, Stambuk HE, Oh JH, Cunanan KM, Nunez DA, Mazaheri Y, Gonen M, HO A, Fagin JA, Wong RJ, Shaha A, Tuttle RM, ShuklaDave A. Dynamic contrastenhanced MRI model selection for predicting tumor aggressiveness in papillary thyroid cancers. NMR Biomed. 2020;33(1):e4166.

Dullien FAL. Porous media: fluid transport and pore structure. NY: Academic Press; 1979.

Yao W, Li Y, Ding G. Interstitial fluid flow: the mechanical environment of cells and foundation of meridians. Evid Based Complement Alternat Med. 2012;2012:853516.

Baxter LT, Jain RK. Vascular permeability and interstitial diffusion in superfused tissues: a twodimensional model. Microvasc Res. 1988;36:108–115.

Padera TP, Kadambi A, di Tomaso E, Carreira CM, Brown EB, Boucher Y, Choi NC, Mathisen D, Wain J, Mark EJ, Munn LL, Jain RK. Lymphatic metastasis in the absence of functional intratumor lymphatics. Science. 2002;296:1883–1886.


Sakai O, Curtin HD, Romo LV, Som PM. Lymph node pathology. Benign proliferative, lymphoma, and metastatic disease. Radiol Clin North Am. 2000;38:979–998.

Som PM. Detection of metastasis in cervical lymph nodes: CT and MR criteria and differential diagnosis. AJR Am J Roentgenol. 1992;158:961–969.

Magdoom KN, Pishko GL, Kim JH, Sarntinoranont M. Evaluation of a voxelized model based on DCEMRI for tracer transport in tumor. J Biomech Eng. 2012;134:091004.

Stohrer M, Boucher Y, Stangassinger M, Jain RK. Oncotic pressure in solid tumors is elevated. Cancer Res. 2000;60:4251–4255.

Tan WH, Wang F, Lee T, Wang CH. Computer simulation of the delivery of etanidazole to brain tumor from PLGA wafers: comparison between linear and double burst release systems. Biotechnol Bioeng. 2003;82:278–288.

Haider MA, Milosevic M, Fyles A, Sitartchouk I, Yeung I, Henderson E, Lockwood G, Lee TY, Roberts TPL. Assessment of the tumor microenvironment in cervix cancer using dynamic contrast enhanced CT, interstitial fluid pressure and oxygen measurements. Int J Radiat Oncol Biol Phys. 2005;62:1100–1107.

Nathanson D, Nelson L. Interstitial fluid pressure in breast cancer, benign breast conditions, and breast parenchyma. Ann Surg Oncol. 1994;1:333–338.

Baxter LT, Jain RK. Transport of fluid and macromolecules in tumors. II. Role of heterogeneous perfusion and lymphatics. Microvasc Res. 1990;40:246–263.

Simonsen TG, Lund KV, Hompland T, Kristensen GB, Rofstad EK. DCEMRIderived measures of tumor hypoxia and interstitial fluid pressure predict outcomes in cervical carcinoma. Int J Radiat Oncol Biol Phys. 2018;102:1193–1201.

Ewing JR, Nagaraja TN, Aryal MP, Keenan KA, Elmghirbi R, BagherEbadian H, Panda S, Lu M, Mikkelsen T, Cabral G, Brown SL. Peritumoral tissue compression is predictive of exudate flux in a rat model of cerebral tumor: an MRI study in an embedded tumor. NMR Biomed. 2015;28:1557–1569.

Magdoom KN, Pishko GL, Rice L, Pampo C, Siemann DW, Sarntinoranont M. MRIbased computational model of heterogeneous tracer transport following local infusion into a mouse hind limb tumor. PLoS One. 2014;9 e89594.

Sarntinoranont M, Rooney F, Ferrari M. Interstitial stress and fluid pressure within a growing tumor. Ann Biomed Eng. 2003;31:327–335.

Pozrikidis C, Farrow DA. A model of fluid flow in solid tumors. Ann Biomed Eng. 2003;31:181–194.

Smith JH, Humphrey JA. Interstitial transport and transvascular fluid exchange during infusion into brain and tumor tissue. Microvasc Res. 2007;73:58–73.

Liu LJ, Brown SL, Ewing JR, Schlesinger M. Phenomenological model of interstitial fluid pressure in a solid tumor. Phys Rev E Stat Nonlin Soft Matter Phys. 2011;84:021919(2 Pt 1).

Sengupta A, Gupta RK, Singh A. Evaluation of B1 inhomogeneity effect on DCEMRI data analysis of brain tumor patients at 3T. J Transl Med. 2017;15:242.

Boucher Y, Jain RK. Microvascular pressure is the principal driving force for interstitial hypertension in solid tumors: implications for vascular collapse. Cancer Res. 1992;52:5110–5114.
Research Articles
Download PDF (2.72 MB)
TOMOGRAPHY, June 2020, Volume 6, Issue 2:129138
DOI: 10.18383/j.tom.2020.00005
Computational Modeling of Interstitial Fluid Pressure and Velocity in Head and Neck Cancer Based on Dynamic ContrastEnhanced Magnetic Resonance Imaging: Feasibility Analysis
Eve LoCastro^{1}, Ramesh Paudyal^{1}, Yousef Mazaheri^{1}, Vaios Hatzoglou^{2}, Jung Hun Oh^{1}, Yonggang Lu^{3}, Amaresha Shridhar Konar^{1}, Kira vom Eigen^{1}, Alan Ho^{4}, James R. Ewing^{5}, Nancy Lee^{7}, Joseph O. Deasy^{1}, Amita ShuklaDave^{1}
Abstract
We developed and tested the feasibility of computational fluid modeling (CFM) based on dynamic contrastenhanced magnetic resonance imaging (DCEMRI) for quantitative estimation of interstitial fluid pressure (IFP) and velocity (IFV) in patients with head and neck (HN) cancer with locoregional lymph node metastases. Twentytwo patients with HN cancer, with 38 lymph nodes, underwent pretreatment standard MRI, including DCEMRI, on a 3Tesla 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 Ktrans 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. CFMestimated 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 drugtherapy 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 wickinneedle (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 contrastenhanced (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 (^{18}FFDG) PET scan. MRI can provide both anatomical and functional results, which can be used to model tumor IFP estimates (22, 23).
DCEMRI uses lowmolecularweight gadoliniumbased 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 DCEMRI, 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 DCEMRI, 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 useK t r a n , 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 DCEMRI 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 HPVnegative (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
Figure 1.
Workflow for magnetic resonance imaging (MRI)based computational fluid modeling (CFM) simulations. The patient undergoes a magnetic resonance (MR) examination with dynamic contrastenhanced imaging. The images are contoured by a neuroradiologist on all slices to properly demarcate the 3dimensional tumor structure. The images are then processed through extended Tofts model (ETM) using MRIQAMPER software. A representative patient's individualized arterial input function (AIF) and tissue signals are plotted here. The ETMK t r a n map is generated and incorporated into the simulation, along with the 3dimensional 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.
DCEMRI Data Acquisition
MRI studies were performed using a neurovascular phasedarray coil on a 3Tesla (T) scanner (Ingenia; Philips Healthcare; The Netherlands). The standard magnetic resonance (MR) acquisition parameters were as follows: multiplanar (axial, coronal, and sagittal) T2weighted (T2w), fatsuppressed, fast spinecho 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 T1weighted (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, T_{10}]), 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 singleexcitation; NA = 1, pixel bandwidth = 250 Hz/pixel, FOV = 20–24 cm^{2}, 256 × 128 matrix zerofilled 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 gadoliniumbased clinically approved contrast at 2 mL/s, followed by saline flush using an MRcompatible programmable power injector (Spectris; Medrad, Indianola, PA) after acquiring the first 5–6 precontrast images. The entire node was covered contiguously with 5mmthick slices, zero gap, yielding 8–10 slices and 40 or 60 phases were acquired with ≤8second temporal resolution. Total acquisition time was ≤8 minutes.
DCEMRI Pharmacokinetic Analysis
The signal intensity in a voxel, measured from a SPGR T1w acquisition, is given by the following equation (32):
[1]
For the R_{1}(t) calculation, TE ≪ T_{2}*, and so approximatinge  T E R 2 * ( t ) ≈ 1 , then Equation [1] becomes the following equation (33):
[2]
The solution for the relaxation rate R_{1} (t) follows from Equation [2] and is given by:
[3]
In the limit of fast water exchange, the change in water proton relaxation rate (ie, ΔR_{1}= ( R 1  R 10 ) [s^{−1}]) owing to CA relaxivity, r_{1} = 4.0 [(mM)^{−1}^{−1}]_{,} is linearly related to the tissue CA concentration, C_{t} (mM) as follows:
[4]
The extended Tofts model (ETM) is based on a 2compartment model (vascular space and extravascular extracelluar space [EES]). The ETM expression for modeling C_{t}(t) is given by the following equation (24):
[5]
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 ITKSNAP (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.
T_{10} values were calculated from multiple flip angle precontrast T1w images (36). The time course of tissue concentration data,C t ( t ) , 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 K t r a n , v_{e}, and v_{p}, and parameter estimation values were bounded in the fitting routine as follows: K t r a n ∈ [0, 5] (min^{−1}), v_{e}, and v_{p} ∈ [0, 1]. Automation of AIF findings, quantitative pharmacokinetic analysis of DCEMRI data, and generation of parametric maps were performed using inhouse MRIQAMPER software (Quantitative Analysis MultiParametric 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 massbalance 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).
[6]
Whereu is the fluid velocity vector (m/s), p i is the (interstitial) fluid pressure (Pa), K is the permeability of the porous membrane (m^{2}, describing the packing of the extracellular matrix), ρ is the interstitial fluid density (kg/m^{3}), μ is interstitial fluid viscosity (Pa*s), additional forces acting on the fluid volume element may be represented as F (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 steadystate velocity. Neglecting both friction within fluid and exchange of momentum between fluid and solid phases, the case of an incompressible fluid (∇ · u = 0 ) in the Navier–Stokes equation (Equation [6]) simplifies to Darcy's law where interstitial fluid velocity (IFV) is related to the gradient in IFP (∇ p i ):
[7]
This Darcy velocityu describes bulk fluid movement within the interstitial space and expresses velocity as the gradient of instantaneous pressure, proportional to the local hydraulic conductivity K_{H}≡K μ :
[8]
Fluid sourcesink assumptions are as follows: fluid can enter the interstitium via the vascular compartment (ϕ v ); fluid is drained from the system via lymphatic pathways (ϕ L ) (1).
[9]
Blood flow transportϕ v across capillary walls was regulated according to the Starling law (1, 31):
[10]
WhereL P is the hydraulic conductivity of the capillary wall, or vessel permeability, S/V is microvascular surface area per unit volume,p V is the blood pressure in the microvessel, p i is IFP, π V is osmotic pressure in microvasculature, π i is osmotic pressure in interstitial space, and σ T 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,p i , and pressure of the lymphatic system, p L and the lymphatic clearance rate, L p L S L V :
[11]
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):
[12]
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
3 × 10^{−12} (normal)
3.8 × 10^{−13} (normal)
7 × 10^{3} (normal)
1,330 (normal)
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 ETMderivedK t r a n is incorporated into the Starling equation (Equation [11]), and scaled by the mean tumor value, <K t r a n > 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, p i :
[13]
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 p_{i} 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 mm^{3} 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).
ETMestimatedK t r a n 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 K t r a n values were interpreted in COMSOL as a scalar field over the mesh domain. The K t r a n field values were normalized by the mean K t r a n 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 noflux 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, L_{p}), lymphatic clearance (L_{pL}S_{L}/V in normal tissue and tumor), and ratio of hydraulic conductivity (allowing K_{H,tumor} to vary relative to K_{H,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 ([L_{pL}S_{L}/V]_{t}/[L_{pL}S_{L}/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 Pvalue < .05 was considered statistically significant.
Results
All 22 patients with HN cancer underwent MRI pretreatment (preTX), 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 cm^{3}).
Figure 2.
T1weighted postcontrast image of patient #1 with right lateral neck nodal metastasis (yellow outline) (A); T1weighted (T1w) postcontrast image with overlaid ROI for analysis (green denotes normal tissue, red denotes viable tumor tissue, blue denotes necrotic core) (B);K t r a n 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) (V_{1} = 14.26 cm^{3}, V_{2} = 25.20 cm^{3}, V_{3} = 31.27 cm^{3}, V_{4} = 28.80 cm^{3}); Row 2: Estimated IFP map of tumor and surrounding normal tissue (mean tumor pressureP ¯ 1 = 1.73 kPa, P ¯ 2 = 2.05 kPa, P ¯ 3 = 2.09 kPa,P ¯ 4 = 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 u ¯ 1 = 1.60 × 10^{−7}m/s, u ¯ 2 = 1.80 × 10^{−7}m/s, u ¯ 3 = 1.67 × 10^{−7}m/s, u ¯ 4 = 1.51 × 10^{−7}m/s); Row 5: IFV profile along vertical bisector line, with tumor boundary (dash).
In Figure 3, mean IFP (p_{mean}) for the cystic node in patient #1 was 1.73 kPa with maximum IFP at tumor core p_{max} = 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 p_{mean} and p_{max} for the metastatic node in patient #2 was 2.05 kPa, and 2.64 kPa, respectively. Values of p_{mean} and p_{max} for patients #3 and 4 were 2.09, 1.65 kPa (p_{mean}) and 2.67, 2.04 kPa (p_{max}), respectively. In the IFP profiles of patients #2 and 3 with nonnecrotic 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 (cm^{3}),K t r a n (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 cm^{3} ranging from 1.36 to 40.02 cm^{3}.
Table 3.
Summary of DCEMRI K^{trans} and CFMestimated IFP and IFV from 22 Patients with Head and Neck Squamous Cell Carcinoma
(IFP [kPa])
(IFV [×10^{−7} m/s])
Figure 4 exhibits the scatter plot of preTX total tumor volume versus mean CFMestimated IFP. Significant correlation (Spearman) was found between preTX 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 T2weighted 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 (L_{p}). 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 nonnecrotic 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 3DK t r a n values. The K t r a n voxel values in the CFM were normalized by the mean K t r a n values in each tumor ROI. CFMgenerated IFP and IFV maps depict heterogeneity within the tumor, stemming from incorporation of K t r a n to inform fluid influx rate variability from the spatial CA extravasation, similar to Pishko et al. (23).
The CFMestimated mean tumor IFP showed a moderate but significant correlation with the total tumor volume at preTX. 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 highvelocity fluid movement exists near the periphery of the tumor owing to a precipitous falloff of pressure (52). Tumors with mildtosevere 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 nonnecrotic tumor in contrast to necrotic tumor may lead to a significant change in IFV of tumor core.
Sensitivity analyses illustrated the influence that tissuedescriptive parameters, such as hydraulic conductivity, K_{H}, 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 K_{H,t}/K_{H,n} has no noticeable effect from baseline value, but in the nonnecrotic case, a lowerratio K_{H,t}/K_{H,n} results in lower overall IFP with the greatest IFV compared with all other parameter changes. In contrast, for the nonnecrotic 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 literaturebased 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 imagingbased biomarkers to characterize the model parameters. A B1 nonuniformity acquisition is required for highfield DCEMRI. The accuracy of ETMderivedK t r a n 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 DCEMRIbased 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.
Notes
[1] Abbreviations:
AIF
arterial input function
CA
contrast agent
CFM
computational fluid modeling
DCEMRI
dynamic contrastenhanced magnetic resonance imaging
EES
extravascular extracellular space
ETM
extended Tofts model
FOV
field of view
HN
head and neck
IFP
interstitial fluid pressure
IFV
interstitial fluid volume
LN
lymph node
PDE
partial differential equation
ROIs
regions of interest
SPGR
spoiled gradient recalled echo
SSE
sum of squared errors
STL
stereolithography
TE
echo time
TR
repetition time
NA
number of averages
WIN
wickinneedle
Acknowledgments
Equal contribution: E.L. and R.P. contributed equally to this work.
This work was supported by U01 CA211205 and in part through the NIH/NCI Cancer Center Support Grant P30 CA008748. We would like to thank the MRI technologists for their great efforts in helping to perform the MRI examinations, and Mr. Christian Czmielewski for his helpful contribution to patient enrollment and data management. This work was supported by the MSKCC internal IMRAS grant, U01 CA211205, and in part through the NIH/NCI Cancer Center Support Grant: P30 CA008748.
References
Journal Information
Journal ID (nlmta): tom
Journal ID (publisherid): TOMOG
Title: Tomography
Subtitle: A Journal for Imaging Research
Abbreviated Title: Tomog.
ISSN (print): 23791381
ISSN (electronic): 2379139X
Publisher: Grapho Publications, LLC (Ann Abor, Michigan)
Article Information
Self URI: media/vol6/issue2/images/GPTOMJ200009.pdf
Copyright statement: © 2020 The Authors. Published by Grapho Publications, LLC
Copyright: 2020, Grapho Publications, LLC
License (openaccess, http://creativecommons.org/licenses/byncnd/4.0/):
This is an open access article under the CC BYNCND license (http://creativecommons.org/licenses/byncnd/4.0/).
Publication date (print): June 2020
Volume: 6
Issue: 2
Pages: 129138
Publisher ID: TOMO.2020.00005
DOI: 10.18383/j.tom.2020.00005
PDF
Download the article PDF (2.72 MB)
Download the full issue PDF (12.51 MB)
Mobileready Flipbook
View the full issue as a flipbook (Desktop and Mobileready)