Novel Radiomic Measurements of Tumor-Associated Vasculature Morphology on Clinical Imaging as a Biomarker of Treatment Response in Multiple Cancers

Neoadjuvant chemotherapy, or chemotherapy administered prior to surgical intervention, often constitutes first-line intervention in a number of cancer domains (1–4). When successful, neoadjuvant chemotherapy can offer substantial benefits for patients by reducing tumor burden and increasing a patient's surgical options (5). However, many patients ultimately fail to respond and, accordingly, will endure financial burden and dangerous side effects without tangible benefit (6). Furthermore, in many cancers, including breast (BRCA) and non–small cell lung cancer (NSCLC), there is a current lack of validated predictive and prognostic biomarkers capable of definitively guiding first-line chemotherapeutic interventions (7–9).

Tumor angiogenesis has long been shown to be crucial in cancer progression. Through influence over the body's machinery for synthesizing vasculature, a tumor will initiate the rapid formation of new blood vessels from preexisting vessels in its surrounding peritumoral environment. This newly formed vessel network, known as the tumor-associated vasculature (TAV), enables tumor growth by perfusing it with abundant oxygen and nutrients, as well as providing an avenue for metastatic spread (10). Histologic and molecular evidence of elevated tumor angiogenesis, such as increased density of microvessels measured via immunostaining or elevated VEGF expression (11), is associated with poor prognosis and therapeutic response.

However, the TAV also possesses crucial architectural differences from healthy blood vessels that are undetectable through routine clinical assessment. Excessive upregulation of angiogenesis creates vessels that are twisted, leaky, and chaotically organized (12–15). Previous work has shown that abnormalities in the shape of tumor vessels are detectable on CT and MRI scans and can distinguish cancer from benign lesions (16–18). This aberrant vessel morphology has been implicated in potentiating treatment refractoriness by reducing drug transfer to the tumor bed, thus leading to a lack of durable response (19). Conversely, successful normalization of TAV architecture through antiangiogenic therapy can promote the efficacy of therapeutic intervention (20). It is likely that tumors that are resistant to treatment will differ in the twistedness and arrangement of their vasculature relative to responsive tumors (21), which in turn could potentially be captured quantitatively on radiologic imaging (22). Consequently, computerized analysis of TAV morphology and spatial organization might enable better guidance of chemotherapy-based treatment by stratifying patients according to likely therapeutic benefit.

In this paper, we present a new computational imaging biomarker based on quantitative tumor-associated vasculature (QuanTAV) measurements to characterize the morphology and architecture of the vessel network surrounding a tumor on radiology scans. We present and evaluate a number of computationally extracted measurements of the twistedness and organization of tumor vessels on pretreatment contrast-enhanced MRI of patients with breast cancer and CT of patients with lung cancer. We further demonstrate the predictive and prognostic utility of QuanTAV measurements in the context of response to chemotherapy-related treatments for four cases involving breast and lung cancers across these modalities. In total, the prognostic and predictive utility of QuanTAV was evaluated on 558 patients, including 242 breast cancer patients receiving anthracycline-based neoadjuvant chemotherapy [BRCA-ACT], 129 patients with breast cancer receiving neoadjuvant chemotherapy with HER2-targeted therapy [BRCA-TCHP], 97 patients with NSCLC receiving platinum-based chemotherapy without surgery [NSCLC-PLAT], and 90 patients with NSCLC receiving a trimodality regimen of neoadjuvant chemoradiation followed by surgical intervention [NSCLC-TRI].

From 3-dimensional volumes delineating a tumor and its corresponding vasculature, our approach mathematically characterizes the complexity of the TAV for use in machine learning models to predict outcomes. Vessel volumes are algorithmically reduced (23) to centerlines and split into discrete branches prior to analysis. A set of 91 QuanTAV measurements are then computed, belonging to one of two categories:

• QuanTAV morphology (17) (61 features): Features describing the 3D shape of tumor vessels. Metrics measuring the twistedness of vessels across different length scales are calculated: torsion (twistedness across a full vessel branch) and curvature (local twistedness among adjacent points along a branch). Additional metrics such as vessel volume and length, and the proportion of vessels entering a tumor, are also derived.

• QuanTAV Spatial Organization (24) (30 features): Features quantifying the degree of heterogeneity in the architecture of the tumor vasculature. 2D projections of the tumor vasculature are generated across each dimension of the imaging plane and in a spherical coordinate system relative to the tumor centroid within a fixed radius of the tumor. The set of QuanTAV Spatial Organization features are statistics describing vessel orientations across each projection image.

For each treatment group, we derive QuanTAV response and risk scores from these metrics, then evaluate their ability to predict response and time to recurrence or progression. Our experimental workflow is summarized in Fig. 1. Code for performing QuanTAV analysis and a workable demo are made available at: https://github.com/ccipd/QuanTAV.

This Health Insurance Portability and Accountability Act of 1996 regulations–compliant study was approved by the institutional review boards at the University Hospitals Cleveland Medical Center, Cleveland, Ohio and the Cleveland Clinic Foundation and the need for informed consent was waived.

A total of 470 patients who received breast neoadjuvant chemotherapy with pretreatment dynamic contrast-enhanced (DCE) MRI were identified for this study. Each breast MRI exam consisted of several T1-weighted acquisitions, including a pre-contrast scan and several scans acquired following the injection of gadolinium-based contrast agent. Thirty-one patients were excluded due to poor image quality resulting in flawed vascular segmentation (including low spatial resolution, insufficient temporal scans or poor temporal resolution, severe artifacts, or inadequate vessel enhancement). Sixty-eight patients were HER2-positive, but received treatment prior to the introduction of anti-HER2 agents, and were thus excluded from analysis. The total number of patients for analysis was 371. Patient response was defined as pathologic complete response (pCR) following chemotherapy, the most commonly used surrogate endpoint in the breast neoadjuvant chemotherapy setting (25, 26) and defined as a lack of remaining invasive cancer cells within the breast or axilla based on pathologic examination of excised surgical samples (ypT0/isN0), 115 achieved pCR, while 256 retained the presence of residual disease following chemotherapy (non-pCR). Patients received different chemotherapeutic regimens based on the expression of the HER2 receptor protein, and patients were split into corresponding treatment groups for analysis.

• BRCA-ACT: 242 patients were HER2-negative, and received an anthracycline-based regimen with or without a taxane. The cohort consisted of 85 patients from University Hospitals Cleveland Medical Center and 157 patients available publicly through the Cancer Imaging Archive (27–29). Following chemotherapy, 48 patients achieved pCR and 194 retained the presence of residual disease (non-pCR). This cohort included patients from the ISPY1 (n = 109) and Breast-NAC Pilot (n = 48) studies that also had recurrence-free survival (RFS) information available. We considered RFS from the initiation of neoadjuvant chemotherapy (RFS for the Breast-NAC MRI Pilot study was recorded following completion of chemotherapy, but was adjusted based on the duration of treatment according to the study protocol; ref. 30).

• BRCA-TCHP: A multi-institutional cohort of 129 HER2-positive patients who received targeted neoadjuvant therapy at University Hospitals Cleveland Medical Center (n = 28) or Cleveland Clinic (n = 101) was also assessed. The majority of patients received neoadjuvant chemotherapy supplemented with trastumuzab and pertuzumab (n = 125), while 5 patients from University Hospitals Cleveland Medical Center received only trastuzumab. Sixty-seven patients achieved pCR and 62 non-pCR. No BRCA-TCHP patients had survival information available.

A total of 187 standard dose, non-contrast lung CT volumes collected prior to treatment were included for analysis. Patients were treated and imaged at University Hospitals Cleveland Medical Center, and were divided into two groups depending on the type of therapeutic regimen that they received (i.e., trimodality or pemetrexed chemotherapy).

• NSCLC-PLAT: A total of 97 patients who received platinum-based chemotherapy without surgical intervention at Cleveland Clinic with available pretreatment CT scans were identified. In the absence of posttreatment surgical samples, response was determined from imaging by RECIST criteria based on size changes between pre- and posttreatment CT. Forty-patients patients were identified as responders, indicated by response or stable disease following platinum-based chemotherapy, while 49 had progression on imaging and were deemed nonresponders. Ninety-two patients had progression-free survival (PFS) information available, which was defined as the time from initiation of treatment to the detection of progressive disease or death, whichever occurred earlier, and was censored at the date of last follow-up for those alive without progression.

• NSCLC-TRI: 90 patients received trimodality therapy, consisting of neoadjuvant chemoradiation followed by surgical intervention. The endpoint for response was major pathologic response (MPR), defined as 10% or less residual viable tumor after neoadjuvant chemoradiation and the recommended surrogate endpoint in resectable NSCLC (31). Thirty-six patients achieved MPR. RFS was measured from the date of surgery to the date of recurrence or the date of death, whichever occurred earlier, and censored at the date of last follow-up for those alive without disease recurrence.

For each treatment group, patients were divided into training and testing sets. Models were developed and optimized on the training set, then applied to the testing set. Three of the treatment groups (BRCA-TCHP, NSCLC-PLAT, NSCLC-TRI) had response rates of approximately 50%, and were accordingly divided randomly in half for training and testing, when possible using the same splits from prior studies (32, 33). Relative to these treatment strategies, the rate of response to BRCA-ACT is substantially lower (34). Given the potential of training data imbalanced between categories to negatively impact classifier performance and robustness (35), a BRCA-ACT training cohort was randomly chosen containing 50% of responders and enough nonresponders to enforce a 3:1 class balance [previously shown to limit the negative effects of class imbalance for an linear discriminant analysis (LDA) classifier; ref. 36, 37]. The composition of each training and testing set, along with availability of response and survival endpoints, is summarized in Supplementary Table S1.

All lung CT volumes were resized to an isotropic resolution of 1 mm³. Tumor boundaries were manually annotated in 3D by an experienced reader. Automatic segmentation of the tumor vasculature was then performed, as depicted in Supplementary Fig. S1. Next, the TAV was extracted in several steps with a protocol previously shown to effectively segment pulmonary vasculature on non-contrast CT (38, 39). Each CT image was masked to the lungs by thresholding at a value of –550 HU followed by morphologic processing (ref. 40; Supplementary Fig. S1B). Next, an open-source (41), multi-scale 3D vessel enhancement filter (42) was applied to emphasize tubular vessel-like structures (see Supplementary Methods - implementation details for parameters), as illustrated in Supplementary Fig. S1C. Thresholding was applied to the vessel enhancement image via Otsu's method (43) to isolate pixels belonging to the vasculature, and then morphologic operations were applied to remove noise and non-vessel artifacts (Supplementary Fig. S1D). A box containing the tumor and an additional 5 cm in each direction was extracted for further analysis (Supplementary Fig. D1E). An open-source fast marching algorithm (44) was applied to the segmented vasculature to identify the center lines of vessels (45) and divide the vessel network into discrete constituent branches (Fig. 1).

The first post-contrast scan was spatially aligned to the pre-contrast scan via affine registration (46) and the difference in image intensities before and after contrast enhancement was then computed, yielding a subtraction image (Supplementary Fig. S2A). Volumes were resized to an isotropic resolution of 1 mm³. 3D tumor boundaries were obtained with a combination of manual annotation and automated segmentation techniques. First, partial tumor annotations on several adjacent axial slices were manually delineated by an experienced reader or derived from segmentations provided for publicly available data (28, 29). A 3-dimensional active contour segmentation algorithm (ref. 47; the ‘chenvese’ function in MATLAB; ref. 47) was applied to expand the annotated 2D slices to a full volumetric segmentation of the tumor in 3D. Vessel segmentation is depicted in Supplementary Fig. S2. The heart and posterior torso were automatically detected and removed (Supplementary Fig. S2B), and a vessel enhancement filter (41, 42) was again applied (Supplementary Fig. S2C) to detect vessel-like objects within the volume (see supplementary implementation details for parameters). Given the lack of true quantitative values in MRI as compared with CT, the vessel enhancement volume was segmented at multiple thresholds derived by Otsu's method (43), which were each refined by morphologic operations (Supplementary Fig. S2D). The resulting segmentations (Supplementary Fig. S2E) were assessed for alignment with vessel enhancement on maximum intensity projections and in 3D. The threshold that best captured the enhancing vasculature within each scan was selected manually by a single reader blinded to clinical data and therapeutic outcome for further analysis. Volumes were cropped 5 cm from the tumor in each dimension (Supplementary Fig. S2F), and center line coordinates and branches of the final vessel network were computed by fast marching (Fig. 1; ref. 23, 44).

From 3D vessel skeletons, 61 quantitative vessel tortuosity features, expanded from a set of 35 introduced previously (17), were computed. The full set of QuanTAV Morphology features is described in Supplementary Table S2. At each point within a branch, curvature was computed as the inverse of the radius of the circle containing that point and the two adjacent points within the branch. Distribution of curvature was summarized along the entire vasculature and each branch through first order statistics (mean, standard deviation, max, skewness, and kurtosis), and branch-level statistics were summarized at the patient level with the same statistics. For each branch, torsion was computed as one minus the ratio of the Euclidean distance between the first and last points of a branch to the branch's length and summarized at the patient level via first order statistics. The distributions of curvature and torsion across the full vasculature were further summarized via 10-bin histograms. Additional vessel metrics—including vessel volume, length, number of vessels entering the tumor, and percentage of vessels in the vessel network feeding the tumor—were also computed.

A set of 30 features describing the organization of the TAV, previously introduced (24) and listed in Supplementary Table S3, were computed. The steps for extracting QuanTAV Spatial Organization features are depicted in Supplementary Fig. S3. From vessel centerlines, a set of 2D vessel projection images are generated, across which statistics summarizing the local orientation of vessels are computed. Along a projection image, the five most prominent vessel orientations are computed within a local window of fixed size via the Hough transform, a mathematical operation for the detection of lines within an image. The window of analysis is moved incrementally along the image to obtain a distribution of vessel orientations across the entire vessel image. The overall distribution of vessel orientations is then summarized by five first order statistics (mean, median, standard deviation, skewness, and kurtosis), which constitute the set of QuanTAV Spatial Organization features.

This process is applied to six distinct projection images. A set of three Cartesian projections is obtained by flattening the vasculature along one of the three spatial dimensions, in the axial, sagittal, or coronal planes. In addition to analyzing the TAV in the original coordinate system, each point within the 3D vasculature is also converted to a spherical coordinate system in order to capture vessel position relative to the tumor. Rather than (X, Y, Z) position, spherical coordinates correspond to elevation from the tumor center, rotation about the tumor center, and distance from the tumor surface. As with Cartesian views, the tumor vasculature is projected along each of these dimensions to obtain three 2D projection images: elevation with respect to rotation, rotation with respect to distance, and elevation with respect to distance. QuanTAV organization features are computed with two tunable parameters: maximum vessel distance from the tumor to include and size of the sliding window used to compute vessel orientations. These parameters were optimized within each imaging modality/cancer domain, and the process and results are described in greater detail within the expanded implementation details located in the supplementary methods.

The set of top features that best predicted therapeutic response for each use case was identified in two rounds of Wilcoxon feature selection in 3-fold cross-validation within the training set. The size of this feature set was determined per cohort based on performance within the training set in cross-validation (see supplementary implementation details). For each cohort, top vessel features were incorporated into a LDA classifier and trained across the full training cohort to predict response in the testing set. The output of this classifier was a score between 0 and 1, in turn corresponding to the level of confidence that a patient would achieve a response following the conclusion of therapy.

For each cohort with survival information available, a survival model was derived in the training set to generate QuanTAV risk scores using a strategy inspired by Bhargava and colleagues (48). All observation times were censored at a maximum of 10 years. Features that were highly correlated and likely redundant were pruned from the feature set, retaining the feature with the highest absolute coefficient value in a multivariable proportional hazards model. A Cox regression model was trained using the remaining vessel features via 10-fold elastic net regularization using the Glmnet for MATLAB package (49). The coefficient values for the model were then applied to training and testing sets to derive patient risk scores. A risk score threshold to optimally stratify patients into high and low risk groups based on maximizing the HR was derived in the training set for each cohort. Further implementation and optimization details are described in the Supplementary Methods.

The primary metric used to evaluate response prediction models were OR and area under the receiver operating-characteristic curve (AUC). Significance level and 95% confidence intervals (CI) of the AUC were computed via permutation testing with Monte Carlo sampling (50, 51) across 50,000 iterations, described in detail in the Supplementary Material of a previous manuscript (52). The univariable and multivariable association of QuanTAV response score and clinical variables with response were assessed on the basis of OR in a logistic regression model. Clinical variables with univariable significance in the training set were incorporated into a clinical feature only logistic regression model, as well as a logistic regression model combining clinical variables with QuanTAV response score. The univariable and multivariable association of QuanTAV response score and clinical variables with response were assessed on the basis of OR in a logistic regression model containing all features.

For prognostic models, both the QuanTAV risk score and categorical QuanTAV risk groups were assessed in univariable and multivariable settings along with baseline clinical variables, as well as pathologic and treatment response information available at the completion of chemotherapeutic regimen. Cox proportional hazards models and risk groups were derived from only baseline clinical variables for comparison against QuanTAV risk groups. The primary metrics used to evaluate association with survival were HR and concordance index (C-index). Univariable and multivariable testing for significant association with prognosis was assessed on the basis of the coefficients of a Cox model within the cohort of interest.

A portion of the data used in this study is publicly available through the Cancer Imaging Archive (TCIA; ref. 27) and include data from the ISPY1-TRIAL (https://wiki.cancerimagingarchive.net/display/Public/ISPY1) and the Breast MRI NACT Pilot (https://wiki.cancerimagingarchive.net/display/Public/Breast-MRI-NACT-Pilot). Datasets from the University Hospitals Cleveland Medical Center (https://www.uhhospitals.org/uh-research/for-researchers) were used with permission for this study and any requests for this data must be directed to the institution.

For the majority of patients with breast cancer who receive neoadjuvant treatment, a chemotherapy-only regimen followed by surgery is standard of care (53). A multi-institutional cohort of 242 patients who received anthracycline-cyclophosphamide alone or followed by a taxane (BRCA-ACT) was assembled and divided into subsets for training (D_tr¹) and independent testing (D_te¹). 19.8% achieved pathologic response on surgical samples following chemotherapy. One hundred fifty-seven BRCA-ACT recipients additionally had 10-year RFS information available. Clinical details are summarized in Table 1.

Via cross-validation, a set of features discriminative of pathologic response were selected and used to train a classification model to yield a QuanTAV response score (Supplementary Table S4) that maximized performance in D_tr¹ (Supplementary Table S5). An increase in average torsion across vessels was identified as the feature most strongly associated with failure to achieve complete response (Fig. 2,A and B). Torsion is defined as the complement of the ratio of the Euclidean distance between a vessel's start and end points to its total length, and is elevated in vessels with internal looping or “U”-shaped vessels that terminate near their origin (13, 54). The presence of such patterns in the vessels surrounding nonresponsive tumors (Fig. 2B) could impede the delivery of systemic therapy to the tumor and subsequently contribute to poor therapeutic response.

When applied to D_te¹, QuanTAV response score identified pathologic response with AUC = 0.65 (95% CI, 0.54–0.76; P = 0.009) and was independently significant in a multivariable comparison with clinico-pathologic variables (Supplementary Table S6). Of three available clinical parameters, only hormone receptor positivity had univariable significance in D_tr¹. A model combining this variable and QuanTAV response score yielded an AUC = 0.78 (95% CI, 0.63–0.87; P = 2e-5) in D_te¹, an increase over hormone receptor status only performance (AUC = 0.69; 95% CI, 0.58–0.80; P = 0.0316). ROC curves for all models in D_te¹ are depicted in Fig. 3A.

A regularized Cox proportional hazards model of QuanTAV features (Supplementary Table S7) was trained to derive a QuanTAV risk score via cross-validation in D_tr¹ (n = 63). A risk score threshold for optimally stratifying patients into low- and high-risk groups was also derived in D_tr¹. Performance of QuanTAV risk score and risk groups in D_tr¹ and D_te¹ are listed in Supplementary Table S8. In D_te¹ (n = 94), the model was significantly prognostic as both a continuous score (P = 0.0022; HR = 1.25; 95% CI, 1.08–1.44; C-index = 0.66) and categorical low- and high-risk groups (P = 0.0096; HR = 4.25; 95% CI, 1.29–14.07; C-index = 0.62). Despite only using measurements from pretreatment imaging, QuanTAV risk groups (Fig. 3E) achieved similar prognostic performance to pathologic response on surgical sample after chemotherapy (Fig. 3F).

We assessed the multivariable significance of QuanTAV risk predictions when compared with baseline clinical variables (age, size, and hormone receptor positivity) and functional tumor volume (FTV; ref. 55), the volume of tumor that is actively vascularized on the basis of contrast agent kinetics, which had previously been assessed for association with survival by Hylton and colleagues in data comprising a portion of the BRCA-ACT cohort. The QuanTAV model remained independently prognostic as both a continuous risk score (HR = 1.20; 95% CI, 1.04–1.40; P = 0.014) and categorical risk groups (HR = 5.51; 95% CI, 1.41–21.49; P = 0.014), along with the majority of clinical variables and FTV (Supplementary Table S9). We also assessed the correlation of each individual feature of the QuanTAV model with FTV. While several features were found to be significantly associated with FTV (Supplementary Table S10), such as features characterizing the quantity of vessels feeding the tumor, the large majority (10 of 14) were independent (P > 0.05). This finding is consistent with the mutual independence observed between risk score and FTV (Supplementary Table S9), and suggests that QuanTAV provides prognostic information beyond clinical measures of perfusion.

Breast cancers with overexpression of the HER2 surface protein are highly aggressive, but can often be effectively combated through a targeted therapeutic strategy supplementing chemotherapy with monoclonal antibodies targeting the HER2 receptor. A second QuanTAV model was trained to predict response to a neoadjuvant regimen combining chemo- plus targeted therapy among patients with HER2-positive breast cancer. The cohort (Table 1) consisted of 129 patients who were HER2-positive and received treatment with docetaxel, carboplatin, trastuzumab, and/or pertuzumab (TCHP), denoted as the BRCA-TCHP treatment group and divided into training (D_tr²) and testing sets (D_te²). Rate of pathologic response was 51.9%.

A QuanTAV response score model (Supplementary Table S11) was trained within D_tr² (Supplementary Table S5) to predict pathologic response to BRCA-TCHP. As was observed in BRCA-ACT, poor response to BRCA-TCHP was associated with elevated vessel torsion, as well as increased skewness of vessel orientations within the XY plane. Within D_te² (Fig. 3B), the vessel model significantly predicted pathologic response (AUC = 0.63; 95% CI, 0.47–0.76; P = 0.042). As a covariate in logistic regression models (Supplementary Table S12), however, QuanTAV response score was not significant in D_te² alone (OR = 0.17; 95% CI, 0.01–2.38; P = 0.188). Hormone receptor status was significant in all subsets, with a testing set AUC of 0.64 (95% CI, 0.52–0.75; P = 0.017). Combining QuanTAV response score with hormone receptor status increased testing AUC to 0.70 (95% CI, 0.53–0.82; P = 0.0036).

In advanced NSCLC, platinum-based chemotherapy is standard-of-care first-line treatment for patients lacking actionable mutations. The NSCLC-PLAT cohort consisted of 97 patients with NSCLC who received a pemetrexed-based platinum doublet regimen and CT imaging before and after treatment at a single institution. In the absence of surgical samples, response was assessed on posttreatment CT based on change from baseline in longest lesion diameter according to RECIST criteria (56). Forty-eight percent had responsive or stable disease on posttreatment imaging, and were categorized as responders, while the remaining patients experienced progression. Fifty-three patients were used for training (D_tr³) and 44 for testing (D_te³).

The NSCLC-PLAT response score derived in D_tr³ (Supplementary Table S5) consisted entirely of QuanTAV spatial organization features (Supplementary Table S13). According to this signature, progression was distinguishable by QuanTAV spatial organization features that corresponded to heterogeneous distribution of vessel orientations (particularly in the region immediately surrounding the tumor). While many NSCLC tumors shared high vascular density regardless of therapeutic outcome (Fig. 4,A–D), QuanTAV spatial organization features reveal crucial architectural differences between responders (Fig. 4B) and progressors (Fig. 4D) at the tumor-vasculature interface. Vessel positions were converted to a spherical coordinate system (Fig. 4,E and F), which were used to derive projection images of vessel organization relative to the tumor (Fig. 4,G and H). For instance, elevated standard deviation of vessel orientations on projection images reflecting rotation and elevation with respect to distance from the tumor were strongly associated with progression (Fig. 4H). Conversely, vessels surrounding responsive tumors maintained a consistent orientation towards the tumor's surface (Fig. 4G).

When applied to D_te³ (Fig. 3C), QuanTAV response scores significantly predicted response on posttreatment imaging with AUC = 0.70 (95% CI, 0.54–0.85; P = 0.024). Only age (P = 0.048) and QuanTAV response score (P = 0.010) significantly differed between responders and nonresponders in D_tr³. However, age was not predictive in D_te³ (P = 0.232) and did not improve the performance of the QuanTAV response score (Fig. 3C). In contrast, QuanTAV response score was the only variable found to be independently significant (P < 0.045) in a multivariable comparison with six clinical variables in D_te³ (Supplementary Table S14).

A QuanTAV risk score model (Supplementary Table S15) and corresponding low/high risk groups were derived in D_tr³ (Supplementary Table S8) to predict PFS: the time from initiation of chemotherapy until progression on imaging, metastasis, or death. Within D_te³ (n = 39), risk group (P = 0.034; HR = 2.29; 95% CI, 1.07–4.94; C-index = 0.62), but not risk score (P = 0.141; HR = 1.12; 95% CI, 0.96–1.31; C-index = 0.61), was significantly associated with PFS. When assessed for independence in a multivariable cox proportional hazards model (Supplementary Table S16) with clinical variables, QuanTAV risk group was the only variable found to be significant (P = 0.028). KM plots for QuanTAV risk group and posttreatment RECIST response are depicted in Fig. 3G and 3H.

For patients with stage III resectable NSCLC, survival can be significantly improved by supplementing platinum-based chemotherapy with radiotherapy and surgical intervention (57), known as tri-modality therapy and denoted here as NSCLC-TRI. Ninety patients received pretreatment CT, followed by neoadjuvant chemoradiation and surgery (Table 2). 41.1% of trimodality recipients achieved pathologic response, and longitudinal outcome data was available for all patients. Patients were divided randomly into training (D_tr⁴) and held-out testing cohorts (D_te⁴).

A NSCLC-TRI QuanTAV response score (Supplementary Table S17) was derived within D_tr⁴ (Supplementary Table S5). In D_te⁴ (Fig. 3D), QuanTAV response score distinguished pathologic response with AUC = 0.71 (95% CI, 0.51–0.84; P = 0.0093). Of eight clinical and treatment-related variables, only Histology (adenocarcinoma vs. squamous cell carcinoma/other) was individually significant (P = 0.0075) in D_tr⁴ and predicted pathologic response with AUC = 0.73 (95% CI, 0.59–0.86; P = 0.0002) in D_te³. The combination of QuanTAV response score and histology outperformed either alone (AUC = 0.85; 95% CI, 0.69–0.94; P = 2E-5). QuanTAV response score remained significantly associated with pathologic response in a multivariable comparison with all available clinical variables in D_te⁴ (OR = 0.0004; 95% CI, 0.00–0.18; P = 0.012), as did histology (Supplementary Table S18).

Next, we assessed the capability of QuanTAV measures to predict RFS from date of surgery in recipients of trimodality therapy. A QuanTAV risk model (Supplementary Table S19) and corresponding risk groups were derived in D_tr⁴ to stratify patients by RFS (Supplementary Table S8). Increases in the standard deviation of curvature across the length of the vessel was associated with elevated risk of recurrence (Fig. 2,C), whereas tumors achieving durable response possessed fewer local variations in curvature due to bends and twists. Similar to the risk score derived for patients with NSCLC receiving chemotherapy alone, QuanTAV Spatial Organization features measuring standard deviation of vessel orientation relative to the tumor centroid were also associated with recurrence or metastasis following surgery.

When applied to D_te⁴, QuanTAV risk score (P = 0.039; HR = 1.28; 95% CI, 1.01–1.62; C-index = 0.66) and categorical risk groups (P = 0.036; HR = 3.77; 95% CI, 1.09–13.00; C-index = 0.64) were significantly prognostic. Kaplan–Meier curves illustrate the ability of pretreatment QuanTAV risk groups (Fig. 3I) and posttreatment pathologic response (Fig. 3J) to stratify patients by RFS.

We assessed QuanTAV risk score and groups for independent prognostic value in D_te⁴ in a multivariable comparison (Supplementary Table S20) including baseline clinical variables (age, sex, histology, clinical stage, largest lesion diameter, ECOG performance status, chemotherapy regimen, radiotherapy induction dose, and surgical procedure type), as well as features of pathology (vascular invasion and lymphatic invasion). Of these, only the QuanTAV model (P = 0.037) and induction dose (P = 0.035) were significant in a comparison with continuous risk score. Categorical QuanTAV risk group was similarly significant in a multivariable setting (P = 0.013), along with several variables including induction dose, surgical procedure, lesion diameter, and presence of vascular invasion.

We sought to understand the impact of vessel segmentation errors on resulting QuanTAV models. To assess the robustness of our approach, we evaluated the performance of QuanTAV response scores in two testing sets (NSCLC-TRI, n = 45 and BRCA-ACT, n = 144) following various levels of perturbation to vessel masks (see Supplementary Methods - additional experiments). The vessel segmentations for each patient in the testing set were degraded through multiple iterations of randomized morphologic operations at each branchpoint and endpoint in the vessel skeletons (Supplementary Fig. S4). Constituent QuanTAV features and corresponding response scores were recomputed from vessel segmentations after 5, 10, 15, and 20 iterations of perturbation. When QuanTAV response score generated from perturbed vessel segmentations to the original response score values via Delong test of paired ROC curves (58), no significant difference in AUC was found at any perturbation level in either the NSCLC-TRI (P = 0.12–0.65) or BRCA-ACT (P = 0.11–0.30) cohorts.

In addition, we sought to assess the generalizability of QuanTAV analysis across institutions. Of the cohorts used in this study, only BRCA-ACT had sufficient data from multiple institutions to assess external generalizability. We conducted a post hoc experiment after the finalization of our primary results within this cohort, where we repeated the training and validation of a pathologic response prediction model, but instead split our dataset according to its source: public (ISPY1-TRIAL and UCSF PILOT; ref. 28, 29; n = 158) or private institution (University Hospitals, n = 84). When QuanTAV response score was trained on public data and tested on private data, response prediction improved to AUC = 0.71 (95% CI, 0.56–0.84; P = 0.01). Similarly, when trained on private data and tested on public data, performance was slightly reduced to AUC = 0.63 (95% CI, 0.51–0.72; P = 0.006). These findings are consistent with our primary findings without separation by institution (AUC = 0.65; 95% CI, 0.54–0.76; P = 0.009; n = 144). Accordingly, we believe that QuanTAV model performs robustly across institutions, with only slight variations in performance that seem to correspond roughly with the number of patients used in training (training n = 85, testing AUC = 0.63; training n = 98, testing AUC = 0.65; training n = 144, testing AUC = 0.71).

In this study, we presented a novel radiomic biomarker that was associated with prognosis and treatment response for two different cancer and three different therapy types. This new category of computational imaging biomarkers leverages morphologic measurements of the twistedness and architecture of the TAV. These vessel-based measurements were found to predict response and survival following intervention across two cancers, two imaging modalities, four chemotherapy-related treatment strategies, and a total of 558 patients. The construction of the tumor's vascular network through neo-angiogenesis plays a crucial role in the determination of patient outcomes by fostering a tumor microenvironment that promotes tumor progression and therapeutic resistance (14, 59). The structural abnormality of the resultant vasculature directly opposes successful therapeutic intervention, possibly owing to poorer delivery of therapeutic agents to the tumor bed (19), while also encouraging the formation of hypoxic regions (60) that reduce efficacy of therapeutic agents and accelerate the development of drug-resistant subclones (61). Consistent with the known deleterious role of abnormal tumor vascularization (10, 11, 59), we found that the expression of features reflecting erratic vascular shape and arrangement were predictive of poor response and elevated risk following chemotherapeutic intervention. Our findings suggest the critical role played by the tumor-associated vessel network across cancer domains in promoting therapeutic response and outcome.

In breast cancer, QuanTAV measurements on pretreatment DCE-MRI predicted patient outcomes following neoadjuvant treatment with two standard-of-care therapeutic strategies in need of validated predictive markers (7, 53). QuanTAV response scores were developed to predict pathologic response following a chemotherapy-only regimen of BRCA-ACT and a BRCA-TCHP regimen of chemotherapy and targeted therapy for patients with a targetable HER2 receptor status. In recipients of BRCA-ACT, QuanTAV-derived models were shown to strongly predict response and RFS independent of clinical variables including hormone receptor status: one of the few predictive markers available across this large and heterogeneous patient population. Additional biomarkers of response and survival for BRCA-ACT is of high clinical interest, because only roughly a quarter of recipients will achieve a complete pathologic response (34). Our findings are consistent with prior research demonstrating that the formation of hook-like vessels feeding the tumor, known as adjacent vessel sign, is associated unfavorable prognosis and tumor phenotype (62). We also investigated the ability of QuanTAV to predict response to a targeted BRCA-TCHP regimen. We observed that the BRCA-TCHP QuanTAV response score achieved a statistically significant ROC AUC and improved performance when incorporated into a clinical model; however, it was not found to be independent as a logistic regression coefficient. In contrast to the other therapies explored in this work, HER2-targeted therapy is mechanistically antiangiogenic (63, 64) and helps normalize the TAV, thus potentially reducing the prognostic value of vascular shape and architecture within this treatment group.

These findings were mirrored in advanced NSCLC, where QuanTAV measures extracted from pretreatment CT volumes were associated with both response and survival following two intervention strategies. First, for patients with advanced NSCLC without actionable mutation, a platinum-based chemotherapy regimen is standard first-line intervention (NSCLC-PLAT). However, only 24% to 31% of patients will achieve response and there are no clinically validated biomarkers for the guidance of platinum-based chemotherapy by benefit (65). QuanTAV measures were predictive of response on posttreatment imaging according to RECIST criteria (56), as well as PFS. Second, for patients with stage III resectable NSCLC, survival can be significantly improved by a trimodality regimen supplementing platinum-based chemotherapy with radiotherapy and surgical intervention. However, trimodality therapy lacks predictive pretreatment markers of benefit and bears a high rate of mortality between 5% and 15% (57). Elevated disorganization and twistedness of the TAV on imaging was associated with a failure to achieve pathologic response and poorer 10-year RFS. Our findings are in agreement with the crucial role of the TAV in NSCLC outcomes, evidenced by the importance of lymphovascular invasion (66) as a prognostic marker and the benefits of TAV-normalization via antiangiogenic therapy for many patients with NSCLC (67). The discriminability of QuanTAV in a regimen including radiotherapy is also consistent with the known role of abnormal vessel geometry in creating a low blood flow, poorly oxygenated tumor microenvironment that that facilitates radioresistance (68).

Critically, we observed that our measurements offered prognostic value independent of measures of functional volume on DCE-MRI (Supplementary Tables S9 and S10), suggesting that discriminative attributes of tumor vascular network architecture may not be captured by contrast agent-based perfusion imaging. Morphologic aberrations of the TAV on radiology have previously been shown to be elevated in the case of breast (16) and NSCLC (17) malignancy, as compared with benign lesions. Reduction in the tortuosity of the TAV on high-resolution brain MR angiography throughout treatment has been shown to be associated with favorable treatment outcomes in metastatic breast cancer (69, 70). Conversely, vessel tortuosity that has not normalized following treatment provide an earlier indication of treatment failure than monitoring tumor growth (71). To our knowledge, this is the first study to date demonstrating the potential of 3D vascular morphology for predicting therapeutic outcomes prior to treatment, as well as the most comprehensive investigation of its role as potential predictive and prognostic biomarkers across cancers, imaging modalities, and treatment types.

QuanTAV analysis represents a new addition to an expanding body of work suggesting the potential of quantitative imaging features mined from radiology to provide predictive biomarkers: (72) an approach known as radiomics. One of the most frequently deployed families of radiomic features is image texture, which quantifies the heterogeneity or spatial arrangement of image signals. Across numerous cancers and imaging modalities, texture-based features of the tumor and its environment have allowed for stratification of tumors into clinically significant biology- and outcome-associated groups (73–75). In breast cancer, textural patterns of the tumor (76, 77), peri-tumoral surroundings (52, 78), bulk parenchyma (79, 80), and lymph nodes (81) on imaging has shown associations with risk and responsiveness to neoadjuvant therapy. Likewise in NSCLC, textural analysis of the tumor and peri-tumoral lung parenchyma has shown promise in predicting benefit of a number of therapeutic approaches, including chemoradiation with and without surgery (32, 33, 82), targeted therapy (83), and immunotherapy (84, 85). Consistently across both cancers, evidence suggests that elevated textural heterogeneity portends poor prognosis and increased risk of nonresponse (52, 73, 78). Tortuous tumor vasculature plays an established role in fostering a heterogeneous, treatment-resistant tumor microenvironment, and, in turn, that heterogeneity fuels further chaotic tumor angiogenesis (19). Thus, a disorganized TAV may be to some degree intertwined with the development of a texturally complex tumor microenvironment on imaging that forms the basis of such prognostic radiomic signatures. To investigate a potential explanatory relationship between the TAV and prognostic texture signatures, we performed a comparison (Supplementary Table S21) of QuanTAV features and risk score with a previously published intra- and peri-tumoral texture-based risk score (33) that was derived within the same NSCLC-TRI cohort. QuanTAV and texture-based risk scores were found to be significantly correlated (r = 0.23, P = 0.030). Of the five most prognostic individual QuanTAV features, reduced variability of curvature along vessels was inversely correlated with texture-derived risk score (r = -0.41, P = 0.0001). This result warrants additional study of the role of angiogenesis as a potential basis for image texture-based biomarkers. Despite potential interactions of vascularity and image texture, we crucially also found that QuanTAV was independent and complementary to texture-based analysis. When these two signatures were combined, they better stratified patients by RFS (C-index = 0.70) than risk scores from texture (C-index = 0.61) or vessel (C-index = 0.66) features alone: suggesting the potential of computational vessel features to complement and improve traditional radiomic analysis. We also repeated a multivariable comparison including QuanTAV and texture risk scores along with the clinical variables previously examined in Supplementary Table S20. We found QuanTAV (P = 0.02) risk score but not texture (P = 0.08) to be independently prognostic in this setting.

Our findings suggest that patients with convoluted vasculature at time of treatment are less likely to derive benefit from systemic therapeutic intervention. Resultantly, patients flagged as nonresponders based upon analysis of the TAV may benefit from antiangiogenic therapy. In NSCLC, bevacizumab, an anti-antiangiogenic targeted therapy, in combination with chemotherapy provides therapeutic benefit by blocking the VEGF receptor, downregulating tumor angiogenesis (86–88), and facilitating delivery of other systemic therapeutics (20). However, bevacizumab is currently prescribed conservatively in NSCLC due to its toxicity and current lack of validated predictive markers of therapeutic benefit (88). QuanTAV measurements could potentially identify patients with NSCLC who would benefit from vascular normalization through the addition of anti-VEGF therapy to their therapeutic regimen. The role of bevacizumab in breast cancer remains controversial, having been previously approved and subsequently revoked for treatment of metastatic breast cancer by the FDA due to safety concerns (89). However, its use in the neoadjuvant setting in combination with chemotherapy has been shown to improve rate of pathologic response (90, 91) and overall survival (92) in breast cancer subsets with specific receptor status and genotype. These results illustrate the important role of patient selection in success of vascular normalization in breast cancer and raise the question of whether antiangiogenic therapies could still be an effective therapeutic option for these patients given more effective tools for targeting their application. Future work should explore the potential association of QuanTAV phenotype and benefit of antiangiogenic therapy in NSCLC and breast cancer.

Our study did have its limitations. First, our segmentation protocol was formulated to achieve a balance of accuracy and efficiency in order to enable analysis within such a large cohort. To assess the robustness of our approach, we evaluated the performance of QuanTAV-based response scores in a breast MRI and lung CT dataset following various levels of disruption to vessel masks and found QuanTAV signatures to be robust to noise in vessel segmentations (Supplementary Fig. S4). Beyond this experiment, it is encouraging that even accounting for segmentation errors, our approach was found to be predictive and prognostic in a wide number of use cases. We did not explore more sophisticated methods of isolating the tumor vasculature in this work, such as specialized deep learning segmentation strategies (93, 94). However, we have shown in a subsequent, preliminary study (95) that fully automated deep learning-based vessel segmentation also enables prognostic QuanTAV analysis in an additional disease and treatment domain—liver metastases treated with CDK 4/6 inhibitors—and validated this signature across institutions. Future work should compare various strategies of vessel segmentation in the context of QuanTAV performance. Second, breast MRI datasets were assembled across institutions and trials, and, consequently, imaging data was highly heterogeneous in acquisition protocol—a confounder we attempted to minimize through preprocessing strategies. Encouragingly, unsupervised clusterings of top QuanTAV features did not reveal site-based batch effects in either breast cancer cohort (Supplementary Fig. S5). Third, across the datasets analyzed, it was necessary to use different clinical endpoints for response (pCR for BRCA-ACT and BRCA-TCHP, MPR for NSCLC-TRI, and RECIST response for NSCLC-PLAT) and survival (RFS for BRCA-ACT and NSCLC-TRI and PFS for NSCLC-PLAT), due to differing accepted and feasible clinical endpoints in the various clinical contexts (See Experimental Design for further definition). For instance, pathologic response could not be assessed in NSCLC-PLAT because patients did not receive surgery. Finally, further validation of our approach is required in a prospective setting prior to clinical adoption. QuanTAV-based measurements should next be evaluated for their ability to predict well-defined clinical endpoints such as pathologic response among patients enrolled in clinical trials including chemotherapy.

N. Braman reports grants from NIH NCI and NIH National Institute of Biomedical Imaging and Bioengineering during the conduct of the study; personal fees and other support from Tempus Labs and Picture Health and personal fees from IBM Research outside the submitted work; in addition, N. Braman has a patent 10,861,152 issued, licensed, and with royalties paid from Picture Health and a patent 10,064,594 issued, licensed, and with royalties paid from Picture Health. P. Prasanna reports a patent for US 10861152 B2 issued. P. Leo reports personal fees from Genentech and other support from Roche outside the submitted work. N.A. Pennell reports personal fees from Merck, AstraZeneca, Genentech, Eli Lilly, Sanofi Genzyme, Mirati, Amgen, Boehringer Ingelheim, Xencor, Janssen, G1 Therapeutics, Pfizer and BMS outside the submitted work. V. Velcheti reports personal fees from Picture Health during the conduct of the study; personal fees from Bristol Myers-Squibb, Merck, AstraZeneca, Bayer, and lFoundation Medicine outside the submitted work; in addition, V. Velcheti has a patent for Radiomics issued. A. Madabhushi reports personal fees and other support from Picture Health and SimBioSys; other support from Elucid Bioimaging; grants from AstraZeneca, Bristol Myers-Squibb, Boehringer Ingelheim, and Eli Lilly; and personal fees from Biohme, Castle Biosciences and Aiforia during the conduct of the study; in addition, A. Madabhushi has a patent 10861152 issued and licensed to Picture Health. No disclosures were reported by the other authors.

The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH, the U.S. Department of Veterans Affairs, the Department of Defense, or the United States Government.

N. Braman: Conceptualization, resources, data curation, software, formal analysis, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. P. Prasanna: Conceptualization, software, methodology, writing–review and editing. K. Bera: Resources, data curation, writing–review and editing. M. Alilou: Software, methodology. M. Khorrami: Data curation, software. P. Leo: Software, methodology. M. Etesami: Resources, data curation, writing–review and editing. M. Vulchi: Resources, data curation. P. Turk: Resources, data curation. A. Gupta: Resources, data curation. P. Jain: Resources, data curation. P. Fu: Validation, writing–review and editing. N. Pennell: Resources. V. Velcheti: Resources. J. Abraham: Resources. D. Plecha: Resources, data curation, writing–original draft. A. Madabhushi: Conceptualization, resources, supervision, funding acquisition, writing–original draft, writing–review and editing.

Research reported in this publication performed by Nathaniel Braman was supported by:

• the NCI under award number F31CA221383-01A1

• the National Institute for Biomedical Imaging and Bioengineering under award numbers T32EB007509

• the Hartwell Foundation.

Research reported in this publication performed by Prateek Prasanna was supported by the NCI under award number 1R21CA258493-01A1.

Research reported in this publication performed by Anant Madabhushi was supported by:

• the NCI under award numbers 1U24CA199374-01, R01CA202752-01A1, R01CA208236-01A1, R01CA216579-01A1, R01CA220581-01A1, 1U01CA239055-01, 1U01CA248226-01, 1U54CA254566-01

• the National Heart, Lung and Blood Institute under award number 1R01HL15127701A1

• the National Institute for Biomedical Imaging and Bioengineering under award numbers 1R43EB028736-01

• the National Center for Research Resources under award number 1 C06 RR12463-01

• the National Institute of Diabetes and Digestive and Kidney Diseases through the Kidney Precision Medicine Project (KPMP) Glue Grant

• the United States Department of Veterans Affairs Biomedical Laboratory Research and Development Service under VA Merit Review Award IBX004121A

• the Office of the Assistant Secretary of Defense for Health Affairs, through

   ◦ the Breast Cancer Research Program (W81XWH-19-1-0668)

   ◦ the Prostate Cancer Research Program (W81XWH-15-1-0558, W81XWH-20-1-0851)

   ◦ the Lung Cancer Research Program (W81XWH-18-1-0440, W81XWH-20-1-0595)

   ◦ the Peer Reviewed Cancer Research Program (W81XWH-18-1-0404)

• the Ohio Third Frontier Technology Validation Fund

• the Clinical and Translational Science Collaborative of Cleveland (UL1TR0002548) from the National Center for Advancing Translational Sciences (NCATS) component of the NIH and NIH roadmap for Medical Research

• The Wallace H. Coulter Foundation Program in the Department of Biomedical Engineering at Case Western Reserve University

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.