Correlation of quantitative dynamic contrast‐enhanced MRI with microvascular density in necrotic, partial necrotic, and viable liver tumors in a rabbit model

The purpose of this study was to examine the correlation of quantitative dynamic contrast‐enhanced (DCE) magnetic resonance imaging (MRI) with microvessel density (MVD) in necrotic, partial necrotic, and viable tumors using a rabbit VX2 liver tumor model. Nine rabbits were used for this study. The complete necrotic area (CNA), partial necrotic area (PNA), and viable tumor area (VTA) of liver tumors were experimentally induced by radiofrequency ablation (RFA). DCE‐MRI data were processed based on the extended Kety model to estimate Ktrans,ve and vp parameters. The boundaries among CNA, PNA, and VTA were delineated based on H&E stain images, and MVD was assessed for each subregion of each VX2 tumor based. There were no correlations between ph‐parameters (Ktrans,ve, and vp) and MVD for CNA. For PNA, the Ktrans values were positively correlated with the MVD (r=0.8124,p<0.0001). For VTA, we found a positive correlation between Ktrans values and the MVD (r=0.5743,p<0.05). Measuring from both the PNA and the VTA, mean Ktrans values were positively correlated with mean MVD (r=0.8470,p<0.0001). In a rabbit VX2 liver tumor model, Ktrans values correlated well with MVD counts of PNA and VTA in liver tumors. PACS number(s): 87.19.If MRI


I. INTRODUCTION
Dynamic contrast-enhanced (DCE) magnetic resonance imaging (MRI) can provide functional information about tumors such as vascular endothelial proliferation, vascular density, and presence of angiogenesis. (1,2) DCE-MRI has recently been used to monitor therapeutic responses to new treatment methods and to evaluate the effects of new antiangiogenic drugs in preclinical and clinical studies. (3)(4)(5)(6)(7) With respect to liver cancer, previous studies have demonstrated the feasibility of using DCE-MRI as a tool to monitor therapeutic response to antiangiogenic drugs via quantitative DCE-MRI measurements (7) and semiquantitative DCE-MRI measurements (i.e., maximal enhancement and peak enhancement). (8) Another study showed the possibility of using semiquantitative DCE-MRI measurements to identify biomarkers for early prediction of hepatocellular carcinoma (HCC) response after radiotherapy (i.e., slope and peak enhancement). (9) These studies in liver tumors incorporated the basic assumption that the parameters computed from DCE-MRI data are noninvasive measurements of tumor microenvironments such as microvessel density (MVD). However, no previous studies have investigated correlations between the perfusion parameters and microvessel density in liver cancer.
Previous studies from other cancers have shown conflicting results. (10)(11)(12)(13) In rectal cancer, Atkin et al. (10) showed a negative correlation between K trans values and MVD counts. In contrast, Zhang et al. (11) showed a negative correlation between time-to-peak enhancement of DCE-MRI data and MVD counts. In prostate cancer, Schlemmer et al. (12) found a correlation between k ep (= K trans / v e ) values and MVD counts, whereas van Niekerk et al. (13) found no relationship between DCE-MRI perfusion parameters (K trans , v e , and k ep ) and MVD counts, but observed positive correlations of k ep values with MVD counts when tumor-normal ratio values were used for correcting intersubject variability.
To examine the relationship between perfusion parameters and microvascular density in liver cancer, we quantified DEC-MRI data based on pharmacokinetic modeling methods using individual arterial input function (AIF) rather than population-averaged or model AIF, to reflect intersubject variability among parameters, such as decline in cardiac output and increase in atherosclerotic vascular changes. And considering the heterogeneity of tumors in terms of enhancement patterns in DCE-MRI data, we experimentally induced complete necrotic area (CNA), partial necrotic area (PNA), and viable tumor area (VTA) by performing radiofrequency ablation (RFA) on a rabbit VX2 liver tumor model, and examined the correlation between pharmacokinetic parameters and MVD counts for each subregion of each VX2 liver tumor.

II. MATERIALS AND METHODS
A. Rabbit VX2 liver tumor model All animal work was performed under a license issued by our Institutional Animal Care and Use Committee. A total of nine female New Zealand white rabbits, each weighing between 2.5 and 3.6 kg (mean, 2.7 kg), were used in this study. Before all procedures, rabbits were sedated with intramuscular injections of 1.0 mg/kg acepromazine maleate (Fermenta Animal Health, Kansas City, MO) and 50 mg/kg of ketamine hydrochloride (Ketaject; Phoenix Scientific, St. Joseph, MO). For the rabbit VX2 liver tumor model, VX2 carcinomas were prepared as described previously. (14) Briefly, a midline laparotomy was first performed in which 1-2 VX2 carcinoma fragments approximately 1 mm 3 in volume were injected into each lobe of the exposed liver using a 16-gauge needle. Electrocoagulation was then performed at the implantation surface to prevent tumor spillage. A total of 34 VX2 carcinomas were implanted into the livers of nine rabbits. Fifteen of 34 VX2 tumors were used, that were greater than 5 mm in size. We performed percutaneous RFA 14 days after VX2 liver tumor surgery.

B. Radiofrequency ablation
To prepare for RFA experiments, the epigastrium and back of each rabbit were shaved and sterilized. A wire-mesh ground pad (10 × 15 cm) and conductive gel were then placed on the back of the rabbit. Next, a 0.5 cm incision was made for insertion of the RFA needle. The Cool-tip RFA System (Valleylab; COVIDIEN, Mansfield, MA) with a 1 cm active tip was placed in the presumed target area of the liver under ultrasonography guidance with a 4V1 Acuson probe (Siemens Medical Solutions USA, Malvern, PA). RFA was performed over a period of 3-10 s using a 500-kHz RF generator (series CC-3; Radionics, Burlington, MA). To create partial necrotic areas within single tumors, incomplete ablation was performed for the partial volume of certain tumors visualized by ultrasonography. The power output was set to 50W and the applied current and power output of the generator was continuously monitored and recorded during RFA.

C. MR imaging
MR experiments were performed four days after RFA. All MR images were acquired using a 3.0 T MRI scanner (Intera Achieva 3T, Philips Medical Systems, Best, The Netherlands) with an extremity coil. Before MR experiments, rabbits were sedated and intravenous access was acquired via a marginal ear vein for injection of contrast agent. High spatial resolution coronal T2-weighted MR images were acquired using a turbo spin echo (TSE) sequence (TE/TR = 80/3300 ms, turbo factor = 10, FOV = 120 mm × 120 mm, matrix = 512 × 512, slice thickness = 2 mm, and acquisition time = 5-6 min). DCE-MR images including whole liver regions were acquired in the coronal plane with a T1-weighted turbo field echo (TFE) sequence (TE / TR = 4.5 / 2.3 ms, turbo factor = 30, FOV = 120 mm × 120 mm, matrix = 160 × 160, slice thickness = 2 mm, flip angle = 12°, and sampling interval = 2.4 s). Baseline images were acquired for 5 s, followed by automatic injection of 0.2 ml/kg Gd-DOTA (Dotarem; Guerbet, France) at 2 mL/s and then a 5 mL normal saline flush with additional acquisition over a total time of 4.8 min (total 120 dynamic images: 2 baseline images, and 118 postcontrast images). For T1 mapping, four precontrast images were acquired with the same imaging parameters using different flip angles (2°, 5°, 10°, and 12°).

D. DCE-MRI processing
MR signals did not exhibit a linear relationship with concentration of the contrast agent, and thus MR signal versus time curves were converted into concentration of contrast agent versus time curves. Specifically, the concentration of contrast agent in a given tissue voxel was estimated by determining the difference in longitudinal relaxation rate as follows: (1) where T 1 (t) and T 10 are the post-and precontrast T1 values, respectively, and r 1 denotes the longitudinal relaxivity (4.39 s -1 mM -1 for blood). (15) The T 10 value of each voxel was estimated using a variable flip angle method (2°, 5°, 10°, and 12°). (16) After estimating S 0 and T 1 values for precontrast images, the postcontrast T 1 value was obtained for the postcontrast image with 12° flip angles based on the fact that T 1 (t) can be estimated as a function of time from the SI(t).
For the quantification of DCE-MR data, we computed pharmacokinetic parameters by nonlinear fitting of concentration time curves into the extended Kety two-compartment model for each voxel: where C t is the concentration of contrast agent at the observed tissue, C p is the concentration in blood plasma, t 0 is the time delay between C p and C t , K trans is the volume transfer constant, v p is the fractional blood plasma volume per unit volume of tissue, and v e is the fractional extravascular extracellular space (EES) per unit volume of tissue. The AIF was manually measured by visual inspection of concentration time curves near the heart in blood vessels C b through the hematocrit for each rabbit: where the H ct is the hematocrit, estimated in rabbits to be 0.45. (17) For estimation of K trans , v p , v e , and t 0 , nonlinear least squares fitting was implemented using the lsqnonlin function from the MATLAB Optimization Toolbox (MathWorks Inc., Natick, MA).
Region of interest (ROI) was manually defined on high-resolution T2-weighted images by one radiologist, who also took into account the pathologic findings of overall tumor extent. Each manually defined ROI was resampled into a low-resolution image (DCE-MR image space) and used for quantitative analysis of DCE-MR images. For each ROI, we computed K trans , v p , and v e parameters from DCE-MRI data in voxel-by-voxel manner. The ph-parameters were averaged across 2nd neighborhood voxels (9 voxels) for each CNA, PNA, and VTA regions based on the K trans parametric map.

E. Pathologic analysis
After MR experiments, rabbits were euthanized by intravenous bolus injection of potassium chloride. Livers were subsequently removed and sectioned in the coronal plane similar to that used for MR imaging. Tumor specimens were divided into anterior and posterior sections. One specimen of each pair was fixed in 10% neutral buffered formalin, embedded in paraffin and then sliced into 4 μm sections for hematoxylin and eosin (H&E) staining and CD31 immunohistochemistry analysis.
For immunohistochemical staining with anti-CD31 antibody, paraffin-embedded sections were deparaffinized in xylene, rehydrated in graded alcohol, and transferred to 0.01 M phosphatebuffered saline (pH 7.4). After heat induced epitope retrieval with citrate buffer (pH 6.0; Dako, Carpinteria, CA) for 3 min at 121°C to reveal hidden antigen epitopes, endogenous peroxidase was blocked with 3% hydrogen peroxide in PBS for 10 min at room temperature. After washing in PBS buffer, sections were treated with serum-free protein blocking solution (Dako) for 20 min at room temperature to block nonspecific binding. Subsequently, sections were incubated with anti-CD31 rabbit polyclonal antibody (1/200; Novus Biologicals, Littleton, CO) for 60 min at room temperature. After washing in PBS, the sections were incubated for 30 min at room temperature with HRP-labeled polymer conjugated secondary antibodies against rabbit IgG (Dako). The color reaction was developed using the ready-to-use DAB (3,3 ′ -diaminobenzidine) substrate-chromogen solution (Dako) for 5 min and then washed with distilled water. Finally, sections were lightly counterstained with Mayer's hematoxylin for 1 min. All sections were also stained with H&E using standard techniques for histological analysis.

F. Microvessel density analysis
In order to examine the relationships between pharmacokinetic parameters estimated from the DCE-MR data and MVD estimated from the histopathological data, an experienced pathologist delineated boundaries among the CNA, the PNA, and the VTA on H&E stain images. We defined the CNA as the regions composed of only necrotic tissues, the PNA as the transition zone from necrotic tissues to viable tumors, and the VTA as regions composed only of viable tumors. Segmented histopathologic images were used as reference images to assess microvessel density.
For MVD analysis, each side was imaged at low magnification (×12.5) to identify CNA, PNA and VTA. For each subregion, five areas were randomly selected. These regions were then imaged at high magnification (×200) and individual microvessels were counted. Any clearly brown stained endothelial cells were counted as individual vessels. The average value of the five areas was used as the representative value of MVD for each CNA, PNA, and VTA.

G. Statistical analysis
For statistical analysis, we used Wilcoxon signed-rank test to examine differences in pharmacokinetic parameters, and differences of MVD between CNA, PNA, and VTA. Pearson correlation analysis was used to examine the relationships between pharmacokinetic parameters and MVD counts in each CNA, PNA, and VTA. P-values of 0.05 or less were considered statistically significant.

III. RESULTS
A. VX2 tumors after RFA Two rabbits died before MR experiments. We excluded VX2 tumors with nonenhanced patterns in DCE-MRI data or that were smaller than 5 mm in size. Fifteen of the total experimentally induced 34 VX2 tumors were included, and the mean size of 15 VX2 tumors was 14 mm, with a range of 8 mm to 22 mm. Figure 1 shows a representative VX2 tumor on T2-weighted MR ( Fig. 1(a)), the corresponding histopathologic specimen ( Fig. 1(b)), and an H&E stain image ( Fig. 1(c)). We did not detect heterogeneity of VX2 tumors on T2-weighted MR images, but the H&E stain images clearly showed heterogeneous VX2 tumors after RFA, composed of completely necrotic tissues and viable tumors. Based on the H&E stain images, boundaries were delineated among the CNA, PNA, and VTA for microvessel density analysis ( Fig. 1(d)).

C. Pharmacokinetic parameters from DCE-MRI data
For the quantification of DCE-MRI data, we used the extended Kety two-compartment model to estimate K trans , v e and v p parameters. The tumor boundaries were manually delineated on T2-weighted MR images as shown in Fig 2(a). The manually defined ROIs on T2-weighted MR images were resampled into the DCE-MRI data for DCE-MRI processing, as shown in Fig. 2(b). Figure 2(c) shows an example of the manually defined AIF near the heart, and Fig. 2(d) shows the AIFs for each rabbit. VX2 tumors exhibited heterogeneous properties in concentration profiles. Figure 2 Figure 3 shows pharmacokinetic parametric maps for K trans , v e , and v p . Statistical analysis revealed that there were significant differences in K trans value between the CNA and the PNA (p < 0.001), between the PNA and the VTA (p < 0.001), and between the CNA and the VTA (p < 0.001). We observed significant differences in the v p parameter between the CNA and the PNA (p < 0.05), and between the CNA and the VTA (p < 0.05), but not between the PNA and the VTA. We did not observe significant differences of v e between the CNA, the PNA, and the VTA. Table 1 shows a summary of the K trans , v e , and v p parameters for each tumor.

D. Microvessel density of subregions for VX2 tumors
To examine the relationship between pharmacokinetic parameters and microvessel density, we manually defined boundaries among the CNA, the PNA, and the VTA based on H&E stain images for each VX2 tumor. MVD was counted for each subregion of a VX2 tumor (Fig. 4). We observed significant differences in MVD count between the CNA and the PNA (p < 0.001), between the CNA and the VTA (p < 0.001), and between the PNA and the VTA (p < 0.001). Table 1 summarizes microvessel density for CNA, PNA, and VTA for each VX2 tumor.

E. Correlations between ph-parameters and MVD for subregions
As assessed using CNA, the ph-parameters were not significantly correlated with MVD count for K trans (r = 0.0075, p = 0.9788), v e (r = -0.1153, p = 0.6824), or v p (r = 0.0245, p = 0.9309). For PNA, the K trans values were positively correlated with the MVD count (r = 0.8124, p < 0.001, Fig. 5(a)), but not with v e (r = 0.4228, p = 0.1164) or v p (r = -0.2750, p = 0.3212). For VTA, we also observed a positive correlation between the K trans values and the MVD count (r = 0.5743, p < 0.05, Fig. 5B), but not for v e (r = 0.1743, p = 0.5345) or v p (r = -0.1340, p = 0.6340). Measuring from both the PNA and the VTA, the mean K trans values were positively correlated with the mean MVD count (r = 0.8470, p < 0.0001, Fig. 5(c)), but not v e (r = 0.3043, p = 0.2701) or v p (r = -0.1943, p = 0.4877). In correlation analyses, we found positive correlations between the K trans values and the MVD count in the PNA and the VTA, but not in the  CNA. When measuring the mean K trans values and the mean MVD count from the PNA and the VTA, the correlation coefficient was higher than when measuring from the PNA or the VTA.

IV. DISCUSSION
In this study, we examine the relationship between pharmacokinetic parameters and MVD in a rabbit VX2 liver tumor model, and demonstrate significant positive correlations of K trans values and MVD counts in the PNA and the VTA. K trans is the volume transfer constant between the intravascular (v p ) and extravascular extracellular compartment (v e ). This exchange rate between v p and v e spaces depends on vascular perfusion factors such as blood flow and permeability surface area product. Thus, K trans may be interpreted as tumor blood flow or vascular permeability, depending on whether the underlying tumor condition is flow-limited or permeability-limited. (18) In this study, the increased numbers of blood vessels counted from CD31 stain images indicate increased blood flow of contrast agent or increased permeability surface area in the PNA and the VTA of VX2 tumors, reflecting increased K trans value in DCE-MRI data, in turn showing positive correlations between MVD and K trans values.
We observed strong correlations when measuring both the PNA and the VTA rather than the PNA or the VTA alone. Interestingly, the correlation was stronger when measured from the PNA than the VTA. The ph-parameters of DCE-MRI data and the MVD counts of CD31 stained data are strongly affected by the definition of the tumor ROI. (19) Inclusion of voxels in the CNA artificially lowered the mean of the values, whereas inclusion of the voxels in VTA artificially increased the mean of the values while ignoring tumor heterogeneity. As shown in Table 1, we found that intersubject variability is higher in VTA than PNA for both ph-parameters and MVD counts. This large intersubject variability may explain the lower correlation in VTA compared with that in PNA. For robust correlations, our findings suggest that ph-parameters and the MVD counts should be measured from both the PNA and the VTA to reflect tumor heterogeneity, while excluding the CNA.
There may have been possible errors in the individual AIF measurement such as inflow effects and inaccurate T 1 measurement caused by B 1 field inhomogeneity. There may be also a possibility of the underestimation of the individual AIF due to the partial volume effects. In this study, to minimize the partial volume effects, we carefully defined the AIF voxels near heart by visual inspection of the shape of concentration time curves (i.e., high peak, narrow shape, and quick wash-in and wash-out).
DCE-MRI data obtained from the liver during free breathing leads to artifacts in pharmacokinetic parameters, particularly in areas close to liver boundaries (such as the diaphragm). Conventional registration methods are likely to fail with DCE-MRI data because of important local intensity changes across different time-points. New methods have been proposed using tracer kinetic model-driven registration, (20) progressive principal component registration, (21) and robust principal component registration, (22) but it is difficult to apply these methods into practical and clinical setting due to its complexity of algorithm. To minimize motion artifacts, in this study, we included only large VX2 tumors (> 5 mm), which are less affected by motion.
More sophisticated pharmacokinetic models are required for the quantification of DCE-MRI data in the liver. Considering the dual bloods flows (arterial and venous) into tissues of the liver, the extended Kety model (two compartments with arterial input function) could not perfectly explain the hemodynamic changes of contrast agent in the liver. In our study, we assumed that a single arterial input (extended Kety model) may suffice to approximate the vascular behavior in VX2 tumors, because hypervascular VX2 tumors are predominantly supplied by hepatic arterial vessels. Many previous studies have demonstrated the usefulness of the extended Kety model for studying liver tumors, in which the predominant blood supply comes from the hepatic arterial blood flow. (23)(24)(25) Possible correlations of pharmacokinetic parameters based on the two compartments with arterial and venous blood flows with histopathologic microvascular parameters should be the focus of future studies.

V. CONCLUSION
We demonstrated the strong correlations of K trans values with MVD counts in the PNA and the VTA in a rabbit VX2 liver tumor model. Thus, our results provide direct evidence that K trans values reflect microvascular information regarding liver tumors, suggesting the possibility of DCE-MRI tool to generate image biomarkers for prediction of response to therapeutic treatments such as antiangiogenic drugs, RFA, and radiotherapy in patients with HCC.